Merge pull request #206 from GeoTIFF/update-2023-06-10

Update 2023 06 10

.eslintignore

@@ -0,0 +1 @@
+*.test.js

.github/workflows/test-main.yml

@@ -10,6 +10,9 @@ jobs:
       - uses: actions/setup-node@v3
         with:
           node-version: 18
+      - uses: actions/setup-python@v4
+        with:
+          python-version: '3.10'
       - run: npm install
       - run: npm run lint
       - run: npm run setup

.gitignore

@@ -106,3 +106,4 @@ pnpm-lock.yaml
 out
 data/geotiff-test-data/*.cog
 data/geotiff-test-data/*.tif*
+Pipfile*

data/create_expected_truth_data.py

@@ -0,0 +1,69 @@
+import csv
+import os
+
+import rasterio
+from rasterio.windows import from_bounds
+
+from rasterstats import zonal_stats
+from rasterstats.utils import VALID_STATS
+
+# known limitations
+# - just for first band
+# - doesn't handle overlapping polygons
+
+test_cases = [
+  ["part-of-india.geojson", "test.tiff"],
+  ["part-of-india-2.geojson", "test.tiff"],
+  ["part-of-india-3.geojson", "test.tiff"],
+  ["./gadm/geojsons/Akrotiri and Dhekelia.geojson", "test.tiff"],
+  ["./gadm/geojsons/Afghanistan.geojson", "./mapspam/spam2005v3r2_harvested-area_wheat_total.tiff"],
+  ["./gadm/geojsons/Akrotiri and Dhekelia.geojson", "./ghsl/tiles/GHS_POP_GPW42015_GLOBE_R2015A_54009_1k_v1_0_4326_30_40.tif"],
+  ["./gadm/geojsons/Ukraine.geojson", "./mapspam/spam2005v3r2_harvested-area_wheat_total.tiff"],
+  ["./veneto/veneto.geojson", "./veneto/geonode_atlanteil.tif"],
+]
+
+for i, (geom, raster) in enumerate(test_cases):
+  print("\n\ncase: " + str(i + 1))
+  feature_stats = zonal_stats(geom, raster, stats=VALID_STATS, band=1)
+
+  # merge feature_stats (doesn't handle overlapping polygons)
+  results = {
+    "count": sum(f['count'] for f in feature_stats),
+    "min": min(f['min'] for f in feature_stats),
+    "max": min(f['max'] for f in feature_stats),
+    "sum": sum(f['sum'] for f in feature_stats),
+  }
+  print("  vector: " + geom)
+  print("  raster: " + raster)
+  print("  result:")
+  for key, value in results.items():
+    print(f"        {key}: {value:,}")
+
+
+sum_test_cases = [
+  "test.tiff",
+  "./mapspam/spam2005v3r2_harvested-area_wheat_total.tiff",
+  "./ghsl/tiles/GHS_POP_GPW42015_GLOBE_R2015A_54009_1k_v1_0_4326_30_40.tif",
+  "./veneto/geonode_atlanteil.tif"
+]
+
+for filepath in sum_test_cases:
+  with rasterio.open(filepath) as src:
+    arr = src.read(1)
+    result = arr[arr != src.meta['nodata']].sum()
+    print(f"\n\nsum\n{filepath}\n{result}")
+
+
+
+country_files = [[f.split(".")[0], f]for f in os.listdir("./gadm/geojsons") if f.endswith(".geojson")]
+ghsl_tiles = [f for f in os.listdir("./ghsl/tiles") if f.endswith(".tif")]
+
+# calculate populations
+with open("country_populations.csv", "w") as csv_file:
+  writer = csv.DictWriter(csv_file, fieldnames=["country", "population"])
+  writer.writeheader()
+  for country, file in country_files:
+    pop = 0
+    for tile in ghsl_tiles:
+      pop += zonal_stats(f"./gadm/geojsons/{file}", f"./ghsl/tiles/{tile}", stats="sum", band=1)[0]['sum']
+    print(f"{country}:  {pop:,}")

data/expected_data.txt

@@ -0,0 +1,110 @@
+
+
+case: 1
+  vector: part-of-india.geojson
+  raster: test.tiff
+  result:
+        count: 1,672
+        min: 0.0
+        max: 7,807.4
+        sum: 3,099,403.8
+
+
+case: 2
+  vector: part-of-india-2.geojson
+  raster: test.tiff
+  result:
+        count: 109
+        min: 0.0
+        max: 6,624.4
+        sum: 291,490.5
+
+
+case: 3
+  vector: part-of-india-3.geojson
+  raster: test.tiff
+  result:
+        count: 67
+        min: 0.0
+        max: 4,124.6
+        sum: 129,465.70000000001
+
+
+case: 4
+  vector: ./gadm/geojsons/Akrotiri and Dhekelia.geojson
+  raster: test.tiff
+  result:
+        count: 1
+        min: 0.0
+        max: 0.0
+        sum: 0.0
+
+
+case: 5
+  vector: ./gadm/geojsons/Afghanistan.geojson
+  raster: ./mapspam/spam2005v3r2_harvested-area_wheat_total.tiff
+  result:
+        count: 7,112
+        min: 0.0
+        max: 5,012.60009765625
+        sum: 2,227,296.0
+
+
+case: 6
+  vector: ./gadm/geojsons/Akrotiri and Dhekelia.geojson
+  raster: ./ghsl/tiles/GHS_POP_GPW42015_GLOBE_R2015A_54009_1k_v1_0_4326_30_40.tif
+  result:
+        count: 229
+        min: 0.0
+        max: 3,256.270263671875
+        sum: 23,469.134765625
+
+
+case: 7
+  vector: ./gadm/geojsons/Ukraine.geojson
+  raster: ./mapspam/spam2005v3r2_harvested-area_wheat_total.tiff
+  result:
+        count: 8,596
+        min: 0.0
+        max: 3,476.800048828125
+        sum: 5,877,552.0
+
+
+case: 8
+  vector: ./veneto/veneto.geojson
+  raster: ./veneto/geonode_atlanteil.tif
+  result:
+        count: 59,409
+        min: 0.0
+        max: 1.3135790824890137
+        sum: 24,963.465454101562
+
+
+sum
+test.tiff
+108343045.40000004
+
+
+sum
+./mapspam/spam2005v3r2_harvested-area_wheat_total.tiff
+215356768.0
+
+
+sum
+./ghsl/tiles/GHS_POP_GPW42015_GLOBE_R2015A_54009_1k_v1_0_4326_30_40.tif
+120669560.0
+
+
+sum
+./veneto/geonode_atlanteil.tif
+27953.119140625
+Lebanon:  5,554,060.0
+Uruguay:  3,303,090.0
+Macedonia:  2,239,499.25
+Ukraine:  12,697,956.0
+Nicaragua:  5,066,313.5
+Afghanistan:  27,024,636.0
+Akrotiri and Dhekelia:  23,469.134765625
+Cyprus:  790,242.0625
+Croatia:  4,685,367.0
+Jamaica:  2,332,581.75

data/holes/hole-2.geojson

@@ -0,0 +1,24 @@
+{
+  "type": "Feature",
+  "properties": {},
+  "geometry": {
+    "type": "Polygon",
+    "coordinates": [
+      [
+        [0, 0],
+        [20, 0],
+        [20, 20],
+        [0, 20],
+        [0, 0]
+      ],
+      [
+        [4, 4],
+        [3.9, 15.9],
+        [15.8, 15.8],
+        [15.9, 3.9],
+        [4, 4]
+      ]
+    ]
+  },
+  "bbox": [0, 0, 20, 20]
+}

data/part-of-india-2.geojson

@@ -0,0 +1,43 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+  {
+    "type": "Feature",
+    "properties": {},
+    "geometry": {
+    "type": "Polygon",
+    "coordinates": [
+      [
+      [
+        86.220703125,
+        22.156883186860703
+      ],
+      [
+        86.341552734375,
+        21.43261686447735
+      ],
+      [
+        86.912841796875,
+        21.70847301324597
+      ],
+      [
+        87.000732421875,
+        22.39071391683855
+      ],
+      [
+        85.968017578125,
+        22.49225722008518
+      ],
+      [
+        85.726318359375,
+        21.912470952680266
+      ],
+      [
+        86.220703125,
+        22.156883186860703
+      ]
+      ]
+    ]
+    }
+  }]
+}

data/part-of-india-3.geojson

@@ -0,0 +1,43 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+  {
+    "type": "Feature",
+    "properties": {},
+    "geometry": {
+    "type": "Polygon",
+    "coordinates": [
+      [
+      [
+        85.869140625,
+        20.64306554672647
+      ],
+      [
+        86.12182617187499,
+        21.022982546427425
+      ],
+      [
+        85.330810546875,
+        21.361013117950915
+      ],
+      [
+        84.44091796875,
+        21.3303150734318
+      ],
+      [
+        85.594482421875,
+        21.074248926792812
+      ],
+      [
+        85.067138671875,
+        20.715015145512087
+      ],
+      [
+        85.869140625,
+        20.64306554672647
+      ]
+      ]
+    ]
+    }
+  }]
+}

data/part-of-india.geojson

@@ -0,0 +1,23 @@
+{
+  "type": "Feature",
+  "properties": {},
+  "geometry": {
+    "type": "Polygon",
+    "coordinates":[
+      [
+        [83.12255859375, 22.49225722008518],
+        [82.96875, 21.57571893245848],
+        [81.58447265624999, 21.207458730482642],
+        [83.07861328125, 20.34462694382967],
+        [83.8037109375, 19.497664168139053],
+        [84.814453125, 19.766703551716976],
+        [85.078125, 21.166483858206583],
+        [86.044921875, 20.838277806058933],
+        [86.98974609375, 22.49225722008518],
+        [85.58349609375, 24.54712317973075],
+        [84.6826171875, 23.36242859340884],
+        [83.12255859375, 22.49225722008518]
+      ]
+    ]
+  }
+}