Updates 2023 08 11

.gitignore

@@ -64,7 +64,7 @@ typings/
 *.swp
 
 # temporary backups
-*.bak
+*bak*
 
 # downloaded test data
 data/*.xml

data/antimeridian/split.geojson

@@ -0,0 +1,63 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "MultiPolygon",
+        "coordinates": [
+          [
+            [
+              [
+                -180,
+                -18
+              ],
+              [
+                -178,
+                -18
+              ],
+              [
+                -178,
+                -20
+              ],
+              [
+                -180,
+                -20
+              ],
+              [
+                -180,
+                -18
+              ]
+            ]
+          ],
+          [
+            [
+              [
+                180,
+                -20
+              ],
+              [
+                178,
+                -20
+              ],
+              [
+                178,
+                -18
+              ],
+              [
+                180,
+                -18
+              ],
+              [
+                180,
+                -20
+              ]
+            ]
+          ]
+        ]
+      },
+      "properties": {},
+      "srs": 4326
+    }
+  ]
+}

data/antimeridian/split.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

@@ -21,9 +21,12 @@
   ["./gadm/geojsons/Ukraine.geojson", "./mapspam/spam2005v3r2_harvested-area_wheat_total.tiff"],
   ["./veneto/veneto.geojson", "./veneto/geonode_atlanteil.tif"],
 
+  ["./santa-maria/santa-maria-mpa.geojson", "./geotiff-test-data/gfw-azores.tif"],
+
   # 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/right-edge.shp", "gfwfiji_6933_COG.tiff"],
+  ["./antimeridian/split.shp", "gfwfiji_6933_COG_Binary.tif"]
 ]
 
 for i, (geom, raster) in enumerate(test_cases):

data/expected_data.txt

@@ -81,6 +81,16 @@ case: 8
 
 
 case: 9
+  vector: ./santa-maria/santa-maria-mpa.geojson
+  raster: ./geotiff-test-data/gfw-azores.tif
+  result:
+        count: 2
+        min: 9.936111450195312
+        max: 19.24805450439453
+        sum: 29.184165954589844
+
+
+case: 10
   vector: ./antimeridian/right-edge.shp
   raster: gfwfiji_6933_COG.tiff
   result:
@@ -90,6 +100,16 @@ case: 9
         sum: 12,783,873.0
 
 
+case: 11
+  vector: ./antimeridian/split.shp
+  raster: gfwfiji_6933_COG_Binary.tif
+  result:
+        count: 327,972
+        min: 1.0
+        max: 1.0
+        sum: 327,972.0
+
+
 sum
 test.tiff
 108343045.40000004
@@ -113,8 +133,3 @@ 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

package.json

@@ -60,23 +60,23 @@
     "@turf/boolean-clockwise": "^6.5.0",
     "@turf/combine": "^6.5.0",
     "bbox-fns": "^0.13.0",
-    "calc-stats": "^2.1.0",
-    "cross-fetch": "^3.1.6",
+    "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",
-    "georaster": "^1.5.6",
+    "georaster": "^1.6.0",
     "get-depth": "^0.0.3",
-    "mathjs": "^11.8.0",
+    "mathjs": "^11.9.1",
     "mpoly": "^0.2.0",
     "preciso": "^0.12.0",
     "proj4": "^2.9.0",
     "proj4-fully-loaded": "^0.2.0",
     "quick-resolve": "^0.0.1",
     "reproject-bbox": "^0.12.0",
     "reproject-geojson": "^0.3.0",
-    "snap-bbox": "^0.4.1",
+    "snap-bbox": "^0.5.0",
     "terraformer-arcgis-parser": "^1.1.0",
     "underscore": "^1.13.6",
     "xdim": "^1.10.1"

src/intersect-polygon/intersect-polygon.module.js

@@ -1,3 +1,4 @@
+import booleanIntersects from "bbox-fns/boolean-intersects.js";
 import calcBoundingBox from "bbox-fns/calc.js";
 import dufour_peyton_intersection from "dufour-peyton-intersection";
 import snap from "snap-bbox";
@@ -10,6 +11,8 @@ const { resolve } = utils;
 const intersectPolygon = (georaster, geometry, perPixelFunction) => {
   const { noDataValue } = georaster;
 
+  const georaster_bbox = [georaster.xmin, georaster.ymin, georaster.xmax, georaster.ymax];
+
   const precisePixelHeight = georaster.pixelHeight.toString();
   const precisePixelWidth = georaster.pixelWidth.toString();
 
@@ -29,22 +32,28 @@ const intersectPolygon = (georaster, geometry, perPixelFunction) => {
     // so using the whole georaster bbox, shouldn't significantly more operations
     intersections = dufour_peyton_intersection.calculate({
       debug: false,
-      raster_bbox: [georaster.xmin, georaster.ymin, georaster.xmax, georaster.ymax],
+      raster_bbox: georaster_bbox,
       raster_height: georaster.height,
       raster_width: georaster.width,
       pixel_height: georaster.pixelHeight,
       pixel_width: georaster.pixelWidth,
       geometry
     });
   } else if (georaster.getValues) {
-    const [xmin, ymin, xmax, ymax] = calcBoundingBox(geometry);
+    const geometry_bbox = calcBoundingBox(geometry);
+
+    if (!booleanIntersects(geometry_bbox, georaster_bbox)) return;
+
+    const [xmin, ymin, xmax, ymax] = geometry_bbox;
 
     // snap geometry bounding box to georaster grid system
     const snapResult = snap({
       bbox: [xmin.toString(), ymin.toString(), xmax.toString(), ymax.toString()],
-      container: [georaster.xmin.toString(), georaster.ymin.toString(), georaster.xmax.toString(), georaster.ymax.toString()],
+      debug: false,
       origin: [georaster.xmin.toString(), georaster.ymax.toString()],
+      overflow: false,
       scale: [precisePixelWidth, "-" + precisePixelHeight],
+      size: [georaster.width.toString(), georaster.height.toString()],
       precise: true
     });
 

src/intersect-polygon/intersect-polygon.test.js

@@ -1,4 +1,5 @@
 import test from "flug";
+import fetch from "cross-fetch";
 import { serve } from "srvd";
 import parseGeoraster from "georaster";
 import bboxPolygon from "@turf/bbox-polygon";
@@ -23,7 +24,7 @@ async function fetch_json(url) {
   }
 }
 
-serve({ debug: true, max: 11, port: 3000, wait: 60 });
+serve({ debug: true, max: 19, port: 3000, wait: 60 });
 
 const urlToGeojson = "http://localhost:3000/data/gadm/geojsons/Akrotiri and Dhekelia.geojson";
 
@@ -156,3 +157,23 @@ test("antimerdian #1", async ({ eq }) => {
   // same as rasterstats
   eq(numberOfIntersectingPixels, 314_930);
 });
+
+test("parse", async ({ eq }) => {
+  const georaster = await parse(urlToData + "geotiff-test-data/gfw-azores.tif");
+  const geojson = await fetch_json(urlToData + "santa-maria/santa-maria-mpa.geojson");
+  let numberOfIntersectingPixels = 0;
+  const geom = convertMultiPolygon(geojson);
+  await intersectPolygon(georaster, geom, () => numberOfIntersectingPixels++);
+  // same as rasterstats
+  eq(numberOfIntersectingPixels, 2);
+});
+
+test("parse no overlap", async ({ eq }) => {
+  const georaster = await parse(urlToData + "geotiff-test-data/gfw-azores.tif");
+  const geojson = await fetch_json(urlToData + "santa-maria/santa-maria-mpa-offset.geojson");
+  let numberOfIntersectingPixels = 0;
+  const geom = convertMultiPolygon(geojson);
+  await intersectPolygon(georaster, geom, () => numberOfIntersectingPixels++);
+  // same as rasterstats
+  eq(numberOfIntersectingPixels, 0);
+});

src/stats/stats.test.js

@@ -10,7 +10,7 @@ import load from "../load";
 import parse from "../parse";
 import stats from "./stats.module";
 
-serve({ debug: true, max: 8, port: 3000, wait: 240 });
+serve({ debug: true, max: 25, port: 3000, wait: 240 });
 
 const url = "http://localhost:3000/data/test.tiff";
 
@@ -28,6 +28,7 @@ const EXPECTED_RASTER_STATS = [
     min: 0,
     mode: 0,
     modes: [0],
+    product: 0,
     range: 8131.2,
     std: 562.8169687364914,
     sum: 108343045.39997534,
@@ -46,6 +47,7 @@ const EXPECTED_BBOX_STATS = [
     min: 0,
     mode: 0,
     modes: [0],
+    product: 0,
     range: 5166.7,
     std: 1216.4677587709607,
     sum: 236435.4,
@@ -64,6 +66,7 @@ const EXPECTED_POLYGON_STATS = [
     min: 0,
     mode: 0,
     modes: [0],
+    product: 0,
     range: 7807.4,
     std: 1507.3255322956486,
     sum: 3_099_403.799999996, // rasterstats says 3,099,403.8
@@ -171,3 +174,20 @@ test("antimerdian #1", async ({ eq }) => {
   const results = await stats(georaster, geom, { stats: ["count", "min", "max", "sum"] });
   eq(results, [{ count: 314_930, min: 0.20847222208976746, max: 492.3219299316406, sum: 12_783_872.545041203 }]);
 });
+
+test("antimerdian #2 (split at antimeridian)", async ({ eq }) => {
+  // converted GeoTIFF to all 1's
+  const georaster = await parse("http://localhost:3000/data/gfwfiji_6933_COG_Binary.tif");
+  let geom = JSON.parse(readFileSync("./data/antimeridian/split.geojson", "utf-8"));
+  geom = reprojectGeoJSON(geom, { from: 4326, to: georaster.projection });
+  const results = await stats(georaster, geom, { stats: ["count", "min", "max", "sum"] });
+  eq(results, [{ count: 327_972, min: 1, max: 1, sum: 327_972 }]);
+});
+
+test("edge", async ({ eq }) => {
+  // converted GeoTIFF to all 1's
+  const georaster = await parse("http://localhost:3000/data/geotiff-test-data/gfw-azores.tif");
+  const geojson = JSON.parse(readFileSync("./data/santa-maria/santa-maria-mpa.geojson", "utf-8"));
+  const results = await stats(georaster, geojson, { stats: ["count", "min", "max", "sum"] });
+  eq(results, [{ count: 2, min: 9.936111450195312, max: 19.24805450439453, sum: 29.184165954589844 }]);
+});