Merge pull request #217 from GeoTIFF/antimeridian-2023-08-24

improved antimeridian support

data/antimeridian/across.geojson

@@ -0,0 +1,35 @@
+{
+  "type": "Feature",
+  "geometry": {
+    "type": "Polygon",
+    "coordinates": 
+      [
+        [
+          [
+            -182,
+            -18
+          ],
+          [
+            -178,
+            -18
+          ],
+          [
+            -178,
+            -20
+          ],
+          [
+            -182,
+            -20
+          ],
+          [
+            -182,
+            -18
+          ]
+        ]
+      ]
+  },
+  "properties": {
+    "description": "rectangular polygon that cross the anti-meridian (on the left side), technically same area as split.geojson"
+  },
+  "srs": 4326
+}

data/antimeridian/across.prj

@@ -0,0 +1 @@
+PROJCS["WGS_1984_EASE-Grid_2.0_Global",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Behrmann"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",30.0],UNIT["Meter",1.0]]

data/create_expected_truth_data.py

@@ -26,15 +26,19 @@
   # https://github.com/perrygeo/python-rasterstats/issues/26
   # ogr2ogr right-edge.shp right-edge.geojson -t_srs EPSG:6933
   ["./antimeridian/right-edge.shp", "gfwfiji_6933_COG.tiff"],
-  ["./antimeridian/split.shp", "gfwfiji_6933_COG_Binary.tif"]
+  ["./antimeridian/right-edge.shp", "gfwfiji_6933_COG.tiff", { "boundless": False, "nodata": -9999 }],
+  ["./antimeridian/split.shp", "gfwfiji_6933_COG_Binary.tif"],
+  ["./antimeridian/across.shp", "gfwfiji_6933_COG_Binary.tif"]
 ]
 
-for i, (geom, raster) in enumerate(test_cases):
+for i, (geom, raster, *opts) in enumerate(test_cases):
   print("\n\ncase: " + str(i + 1))
   print("  vector: " + geom)
   print("  raster: " + raster)
+  print("  opts: " + str(opts))
 
-  feature_stats = zonal_stats(geom, raster, stats=VALID_STATS, band=1)
+
+  feature_stats = zonal_stats(geom, raster, stats=VALID_STATS, band=1, **(opts[0] if opts else {}))
 
   try:
     # merge feature_stats (doesn't handle overlapping polygons)
@@ -67,7 +71,7 @@
 
 
 
-country_files = [[f.split(".")[0], f]for f in os.listdir("./gadm/geojsons") if f.endswith(".geojson")]
+country_files = sorted([[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

data/expected_data.txt

@@ -3,6 +3,7 @@
 case: 1
   vector: part-of-india.geojson
   raster: test.tiff
+  opts: []
   result:
         count: 1,672
         min: 0.0
@@ -13,6 +14,7 @@ case: 1
 case: 2
   vector: part-of-india-2.geojson
   raster: test.tiff
+  opts: []
   result:
         count: 109
         min: 0.0
@@ -23,6 +25,7 @@ case: 2
 case: 3
   vector: part-of-india-3.geojson
   raster: test.tiff
+  opts: []
   result:
         count: 67
         min: 0.0
@@ -33,6 +36,7 @@ case: 3
 case: 4
   vector: ./gadm/geojsons/Akrotiri and Dhekelia.geojson
   raster: test.tiff
+  opts: []
   result:
         count: 1
         min: 0.0
@@ -43,6 +47,7 @@ case: 4
 case: 5
   vector: ./gadm/geojsons/Afghanistan.geojson
   raster: ./mapspam/spam2005v3r2_harvested-area_wheat_total.tiff
+  opts: []
   result:
         count: 7,112
         min: 0.0
@@ -53,6 +58,7 @@ case: 5
 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
+  opts: []
   result:
         count: 229
         min: 0.0
@@ -63,6 +69,7 @@ case: 6
 case: 7
   vector: ./gadm/geojsons/Ukraine.geojson
   raster: ./mapspam/spam2005v3r2_harvested-area_wheat_total.tiff
+  opts: []
   result:
         count: 8,596
         min: 0.0
@@ -73,6 +80,7 @@ case: 7
 case: 8
   vector: ./veneto/veneto.geojson
   raster: ./veneto/geonode_atlanteil.tif
+  opts: []
   result:
         count: 59,409
         min: 0.0
@@ -83,6 +91,7 @@ case: 8
 case: 9
   vector: ./santa-maria/santa-maria-mpa.geojson
   raster: ./geotiff-test-data/gfw-azores.tif
+  opts: []
   result:
         count: 2
         min: 9.936111450195312
@@ -93,6 +102,7 @@ case: 9
 case: 10
   vector: ./antimeridian/right-edge.shp
   raster: gfwfiji_6933_COG.tiff
+  opts: []
   result:
         count: 314,930
         min: 0.20847222208976746
@@ -101,15 +111,38 @@ case: 10
 
 
 case: 11
+  vector: ./antimeridian/right-edge.shp
+  raster: gfwfiji_6933_COG.tiff
+  opts: [{'boundless': False, 'nodata': -9999}]
+  result:
+        count: 314,930
+        min: 0.20847222208976746
+        max: 492.3219299316406
+        sum: 12,783,873.0
+
+
+case: 12
   vector: ./antimeridian/split.shp
   raster: gfwfiji_6933_COG_Binary.tif
+  opts: []
   result:
         count: 327,972
         min: 1.0
         max: 1.0
         sum: 327,972.0
 
 
+case: 13
+  vector: ./antimeridian/across.shp
+  raster: gfwfiji_6933_COG_Binary.tif
+  opts: []
+  result:
+        count: 29,211,252
+        min: 1.0
+        max: 1.0
+        sum: 29,211,252.0
+
+
 sum
 test.tiff
 108343045.40000004
@@ -128,8 +161,13 @@ sum
 sum
 ./veneto/geonode_atlanteil.tif
 27953.119140625
+Afghanistan:  27,024,636.0
+Akrotiri and Dhekelia:  23,469.134765625
+Croatia:  4,685,367.0
+Cyprus:  790,242.0625
+Jamaica:  2,332,581.75
 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
+Ukraine:  12,697,956.0
+Uruguay:  3,303,090.0

package.json

@@ -59,20 +59,23 @@
   "dependencies": {
     "@turf/boolean-clockwise": "^6.5.0",
     "@turf/combine": "^6.5.0",
-    "bbox-fns": "^0.13.0",
+    "@turf/rewind": "^6.5.0",
+    "bbox-fns": "^0.17.0",
     "calc-stats": "^2.2.0",
     "cross-fetch": "^4.0.0",
     "dufour-peyton-intersection": "0.2.0",
     "fast-max": "^0.4.0",
     "fast-min": "^0.3.0",
     "faster-median": "^1.0.0",
+    "geojson-antimeridian-cut": "^0.1.0",
     "georaster": "^1.6.0",
     "get-depth": "^0.0.3",
     "mathjs": "^11.9.1",
     "mpoly": "^0.2.0",
     "preciso": "^0.12.0",
     "proj4": "^2.9.0",
     "proj4-fully-loaded": "^0.2.0",
+    "quick-promise": "^0.1.0",
     "quick-resolve": "^0.0.1",
     "reproject-bbox": "^0.12.0",
     "reproject-geojson": "^0.3.0",