testing equal area statistics

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

{
  "type": "Feature",
  "properties": {},
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [0.0000028, -15.8750977],
        [20.7286425, 0.0000028],
        [20.7286425, -15.8750977],
        [0.0000028, -15.8750977]
      ]
    ]
  }
}

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)

Here's the test tif. Looks big enough imo. Min value: 1, max value: 10000

Driver: GTiff/GeoTIFF
Files: row_encoded_raster_cea.tif
Size is 10000, 10000
Coordinate System is:
PROJCRS["WGS 84 / NSIDC EASE-Grid 2.0 Global",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["US NSIDC EASE-Grid 2.0 Global",
        METHOD["Lambert Cylindrical Equal Area",
            ID["EPSG",9835]],
        PARAMETER["Latitude of 1st standard parallel",30,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Environmental science - used as basis for EASE grid."],
        AREA["World between 86°S and 86°N."],
        BBOX[-86,-180,86,180]],
    ID["EPSG",6933]]
Data axis to CRS axis mapping: 1,2
Origin = (0.000000000000000,0.000000000000000)
Pixel Size = (200.000000000000000,-200.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Lower Left  (       0.000,-2000000.000) (  0d 0' 0.01"E, 15d52'27.63"S)
Upper Right ( 2000000.000,       0.000) ( 20d43'42.01"E,  0d 0' 0.01"N)
Lower Right ( 2000000.000,-2000000.000) ( 20d43'42.01"E, 15d52'27.63"S)
Center      ( 1000000.000,-1000000.000) ( 10d21'51.00"E,  7d51'45.47"S)
Band 1 Block=10000x1 Type=Int32, ColorInterp=Gray

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