-
-
Save moradology/531a4b044350bebf5c0fe8877cb74f3a to your computer and use it in GitHub Desktop.
testing equal area statistics
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/env python | |
import numpy as np | |
import rasterio | |
from rasterio.transform import Affine | |
from rasterio.crs import CRS | |
width, height = 10000, 10000 | |
encoded_values = np.zeros((height, width), dtype=np.int32) | |
for row in range(height): | |
encoded_values[row, :] = row + 1 | |
# (10000 * 200)meters == ~1242.7 miles | |
transform = Affine.translation(0, 0) * Affine.scale(200, -200) | |
# Cylindrical Equal Area projection | |
crs = CRS.from_epsg(6933) | |
with rasterio.open( | |
'row_encoded_raster_cea.tif', 'w', driver='GTiff', | |
height=height, width=width, | |
count=1, dtype=str(encoded_values.dtype), | |
crs=crs, | |
transform=transform, | |
) as dst: | |
dst.write(encoded_values, 1) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/env python | |
# OK, so we have this very simple case to work with. Areas are equal, row/column are the same. | |
# A diagonal will slice one cell in half per row, so we can get rough expected values from scratch | |
# Here, I chose the bottom of the diagonal that runs from bottom left to top right. | |
N = 10000 | |
row_values = [x+1 for x in range(N)] | |
total_weighted_sum = 0 | |
for i in range(N): | |
contribution = (i + 0.5) * row_values[i] | |
total_weighted_sum += contribution | |
# 333358332500.0 | |
print(f"Total weighted sum: {total_weighted_sum}") | |
cells_counted = (10000*10000)/2 | |
# 6667.16665 | |
print(f"Total weighted mean: {total_weighted_sum/cells_counted}") |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/env python | |
# OK, so we have this very simple case to work with. Areas are equal, row/column are the same. | |
# A diagonal will slice one cell in half per row, so we can get rough expected values from scratch | |
N = 10000 | |
row_values = [x+1 for x in range(N)] | |
total_weighted_sum = 0 | |
for i in range(N): | |
contribution = (i + 0.5) * row_values[i] | |
total_weighted_sum += contribution | |
# 333358332500.0 | |
print(f"Total weighted sum: {total_weighted_sum}") | |
cells_counted = (10000*10000)/2 | |
# 6667.16665 | |
print(f"Total weighted mean: {total_weighted_sum/cells_counted}") |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/env python | |
import rasterio | |
from rasterio.windows import Window | |
def create_subtiffs(tif_path, sub_width, sub_height): | |
with rasterio.open(tif_path) as src: | |
for i in range(0, src.width, sub_width): | |
for j in range(0, src.height, sub_height): | |
window = Window(i, j, sub_width, sub_height) | |
transform = src.window_transform(window) | |
sub_tif_data = src.read(window=window) | |
sub_tif_path = f'subs/sub_{i}_{j}.tif' | |
with rasterio.open( | |
sub_tif_path, 'w', | |
driver='COG', | |
height=sub_height, width=sub_width, | |
count=src.count, | |
dtype=sub_tif_data.dtype, | |
crs=src.crs, | |
transform=transform, | |
) as dst: | |
dst.write(sub_tif_data) | |
create_subtiffs('row_encoded_raster_wgs84.tif', 5000, 5000) |
Here's that equal area tif in 4326 (gdalwarp -t_srs EPSG:4326 -r bilinear -co COMPRESS=LZW row_encoded_raster_cea.tif row_encoded_raster_wgs84.tif
).
Note that the size has very much changed (as expected, moving away from an equal area projection)
Files: row_encoded_raster_wgs84.tif
Size is 11228, 8599
Coordinate System is:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (0.000000000000000,0.000000000000000)
Pixel Size = (0.001846156257629,-0.001846156257629)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 0.0000000, 0.0000000) ( 0d 0' 0.01"E, 0d 0' 0.01"N)
Lower Left ( 0.0000000, -15.8750977) ( 0d 0' 0.01"E, 15d52'30.35"S)
Upper Right ( 20.7286425, 0.0000000) ( 20d43'43.11"E, 0d 0' 0.01"N)
Lower Right ( 20.7286425, -15.8750977) ( 20d43'43.11"E, 15d52'30.35"S)
Center ( 10.3643212, -7.9375488) ( 10d21'51.56"E, 7d56'15.18"S)
Band 1 Block=11228x1 Type=Int32, ColorInterp=Gray
![image](https://private-user-images.githubusercontent.com/1977405/300927385-f9fddd28-bc3d-4716-bf03-5f5e4f3fa3b0.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MjA0Nzc0NTksIm5iZiI6MTcyMDQ3NzE1OSwicGF0aCI6Ii8xOTc3NDA1LzMwMDkyNzM4NS1mOWZkZGQyOC1iYzNkLTQ3MTYtYmYwMy01ZjVlNGYzZmEzYjAucG5nP1gtQW16LUFsZ29yaXRobT1BV1M0LUhNQUMtU0hBMjU2JlgtQW16LUNyZWRlbnRpYWw9QUtJQVZDT0RZTFNBNTNQUUs0WkElMkYyMDI0MDcwOCUyRnVzLWVhc3QtMSUyRnMzJTJGYXdzNF9yZXF1ZXN0JlgtQW16LURhdGU9MjAyNDA3MDhUMjIxOTE5WiZYLUFtei1FeHBpcmVzPTMwMCZYLUFtei1TaWduYXR1cmU9NGEzY2JkOTk0YzUxZWRmNWI2NmM1NmNlNjk1NjA3ZGMzYzE3M2YwNjgwYzViYzdmYjM0Y2Q5ZWM3YzQ3ZGVkYSZYLUFtei1TaWduZWRIZWFkZXJzPWhvc3QmYWN0b3JfaWQ9MCZrZXlfaWQ9MCZyZXBvX2lkPTAifQ.VJNCl6G9yOL4PHhMr_oMd1Sq33svp-9EMvMs7maIh-E)
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Here's the test tif. Looks big enough imo. Min value: 1, max value: 10000