moradology/build_equal_area_test_tif.py Secret

## build_equal_area_test_tif.py
#!/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)

## calculate_expected_value.py
#!/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}")

## test_minmax.py
#!/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}")

## test_zone.geojson

      
Display the source blob

    
Display the rendered blob

    
    Raw
  

              test_zone.geojson
            
          
        Loading

      Sorry, something went wrong. Reload?
      Sorry, we cannot display this file.
      Sorry, this file is invalid so it cannot be displayed.
      
          Viewer requires iframe.
      
    
## windows_for_mosaic.py
#!/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)
	#!/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)
	#!/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}")
	#!/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)