fixed issues whereby calculating stats in a rectangle would pull all …

…pixels that touch the rectangle and not just those whose centroid is inside

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

@@ -2,11 +2,11 @@
 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
@@ -20,24 +20,33 @@
   ["./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"],
+
+  # 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"]
 ]
 
 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:,}")
+
+  feature_stats = zonal_stats(geom, raster, stats=VALID_STATS, band=1)
+
+  try:
+    # 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": max(f['max'] for f in feature_stats),
+      "sum": sum(f['sum'] for f in feature_stats),
+    }
+    print("  result:")
+    for key, value in results.items():
+      print(f"        {key}: {value:,}")
+  except Exception as e:
+    print(feature_stats)
+    raise e
 
 
 sum_test_cases = [

data/expected_data.txt

@@ -76,10 +76,20 @@ case: 8
   result:
         count: 59,409
         min: 0.0
-        max: 1.3135790824890137
+        max: 5.398769378662109
         sum: 24,963.465454101562
 
 
+case: 9
+  vector: ./antimeridian/right-edge.shp
+  raster: gfwfiji_6933_COG.tiff
+  result:
+        count: 314,930
+        min: 0.20847222208976746
+        max: 492.3219299316406
+        sum: 12,783,873.0
+
+
 sum
 test.tiff
 108343045.40000004

package.json

@@ -22,6 +22,7 @@
     "build:prod:lite": "npm run build:prod:node:lite && npm run build:prod:web:lite",
     "build:prod:node:lite": "GEOBLAZE_WEBPACK_ENTRY='./lite/index.js' GEOBLAZE_WEBPACK_OUTPUT_FILENAME='geoblaze.lite.node.min.js' webpack --mode production --target node",
     "build:prod:web:lite": "GEOBLAZE_WEBPACK_ENTRY='./lite/index.js' GEOBLAZE_WEBPACK_OUTPUT_FILENAME='geoblaze.lite.web.min.js' webpack --mode production --target web",
+    "create-expected-data": "cd data && pipenv run python3 create_expected_truth_data.py > expected_data.txt",
     "dev": "webpack --mode development --target node --watch",
     "document": "npx documentation build src/** --config documentation.yml -f html -o docs && echo 'docs.geoblaze.io' > docs/CNAME",
     "fix": "eslint lite --fix && eslint src --fix",

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

@@ -3,14 +3,15 @@ import { serve } from "srvd";
 import fetch from "cross-fetch";
 import parseGeoraster from "georaster";
 import bboxPolygon from "@turf/bbox-polygon";
+import reprojectGeoJSON from "reproject-geojson";
 
 import get from "../get";
 import load from "../load";
 import parse from "../parse";
 import { convertMultiPolygon } from "../convert-geometry";
 import intersectPolygon from "../intersect-polygon";
 
-serve({ debug: true, max: 5, port: 3000 });
+serve({ debug: true, max: 10, port: 3000 });
 
 const urlToGeojson = "http://localhost:3000/data/gadm/geojsons/Akrotiri and Dhekelia.geojson";
 
@@ -131,3 +132,16 @@ test("Hole Test", async ({ eq }) => {
   await intersectPolygon(georaster, geom, () => numberOfIntersectingPixels++);
   eq(numberOfIntersectingPixels, 12);
 });
+
+test("antimerdian #1", async ({ eq }) => {
+  const georaster = await parse(urlToData + "gfwfiji_6933_COG.tiff");
+  const geojson = await fetch(urlToData + "antimeridian/right-edge.geojson").then(res => res.json());
+  console.dir({ georaster, geojson });
+  let numberOfIntersectingPixels = 0;
+  let geom = convertMultiPolygon(geojson);
+  // reproject geometry to projection of raster
+  geom = reprojectGeoJSON(geom, { from: 4326, to: georaster.projection });
+  await intersectPolygon(georaster, geom, () => numberOfIntersectingPixels++);
+  // same as rasterstats
+  eq(numberOfIntersectingPixels, 314_930);
+});

src/stats/stats.test.js

@@ -1,10 +1,13 @@
-import test from "flug";
 import { readFileSync } from "fs";
+
+import test from "flug";
+import fetch from "cross-fetch";
 import { serve } from "srvd";
 import parseGeoraster from "georaster";
-import fetch from "cross-fetch";
+import reprojectGeoJSON from "reproject-geojson";
 
 import load from "../load";
+import parse from "../parse";
 import stats from "./stats.module";
 
 serve({ debug: true, max: 8, port: 3000, wait: 240 });
@@ -160,3 +163,11 @@ test("Hole Test", async ({ eq }) => {
   eq(results[0].count, 12);
   eq(results[0].sum, 12);
 });
+
+test("antimerdian #1", async ({ eq }) => {
+  const georaster = await parse("http://localhost:3000/data/gfwfiji_6933_COG.tiff");
+  let geom = JSON.parse(readFileSync("./data/antimeridian/right-edge.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: 314_930, min: 0.20847222208976746, max: 492.3219299316406, sum: 12_783_872.545041203 }]);
+});

src/sum/sum.test.js

@@ -1,8 +1,10 @@
 import { readFileSync } from "fs";
+import fetch from "cross-fetch";
 import test from "flug";
 import { serve } from "srvd";
 import turfBbox from "@turf/bbox";
 import { polygon as turfPolygon } from "@turf/helpers";
+import reprojectGeoJSON from "reproject-geojson";
 
 import utils from "../utils";
 import load from "../load";

src/utils/utils.module.js

@@ -16,6 +16,11 @@ function fetchJsons(urls) {
     throw error;
   }
 }
+
+function roundDown(n) {
+  return -1 * Math.round(-1 * n);
+}
+
 const utils = {
   convertToGeojsonIfNecessary(input, debug = false) {
     if (this.isEsriJson(input, debug)) {
@@ -45,10 +50,10 @@ const utils = {
     system is inverted along the y axis relative to the lat/long (geographic)
     coordinate system */
     return {
-      xmin: Math.floor((crsXMin - georaster.xmin) / georaster.pixelWidth),
-      ymin: Math.floor((georaster.ymax - crsYMax) / georaster.pixelHeight),
-      xmax: Math.ceil((crsXMax - georaster.xmin) / georaster.pixelWidth),
-      ymax: Math.ceil((georaster.ymax - crsYMin) / georaster.pixelHeight)
+      xmin: roundDown((crsXMin - georaster.xmin) / georaster.pixelWidth),
+      ymin: roundDown((georaster.ymax - crsYMax) / georaster.pixelHeight),
+      xmax: Math.round((crsXMax - georaster.xmin) / georaster.pixelWidth),
+      ymax: Math.round((georaster.ymax - crsYMin) / georaster.pixelHeight)
     };
   },
 
@@ -187,7 +192,9 @@ const utils = {
       results.push(i);
     }
     return results;
-  }
+  },
+
+  roundDown
 };
 
 export default utils;