-
-
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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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)