Merge pull request #214 from GeoTIFF/rect-2023-07-22

fixed rectangular intersecting

.github/workflows/test-main.yml

@@ -13,9 +13,14 @@ jobs:
       - uses: actions/setup-python@v4
         with:
           python-version: '3.10'
+      - run: sudo apt-add-repository ppa:ubuntugis/ubuntugis-unstable
+      - run: sudo apt-get update
+      - run: sudo apt-get install gdal-bin libgdal-dev
       - run: npm install
       - run: npm run lint
       - run: npm run setup
+      # list all files in the data folder for debugging purposes
+      - run: find ./data
       - run: npm run test
       # see https://github.com/webpack/webpack/issues/14532
       - run: NODE_OPTIONS=--openssl-legacy-provider npm run build

data/antimeridian/right-edge.geojson

@@ -0,0 +1,15 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "properties": {},
+      "geometry": {
+        "coordinates": [
+         [ [175, -18],[175, -20],[179, -20],[179, -18],[175, -18]]
+        ],
+        "type": "Polygon"
+      }
+    }
+  ]
+}

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",

setup.sh

@@ -1,3 +1,5 @@
+set -e
+
 # start in data directory
 cd data
 

src/histogram/histogram.test.js

@@ -29,13 +29,13 @@ const ratioEIGeorasterResults = {
 };
 
 const ratioQuantileBboxResults = {
-  "0 - 93.2": 31,
-  ">93.2 - 272": 31,
-  ">272 - 724": 31,
-  ">724 - 1141.9": 31,
-  ">1141.9 - 1700.6": 31,
-  ">1700.6 - 2732.4": 31,
-  ">2732.4 - 5166.7": 28
+  "0 - 97.1": 27,
+  ">97.1 - 272": 27,
+  ">272 - 733.5": 27,
+  ">733.5 - 1176.9": 27,
+  ">1176.9 - 1592.1": 27,
+  ">1592.1 - 2672.4": 27,
+  ">2672.4 - 5166.7": 27
 };
 
 const ratioQuantilePolygonResults = {

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

@@ -1,16 +1,29 @@
 import test from "flug";
 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 });
+async function fetch_json(url) {
+  let res;
+  let text;
+  try {
+    res = await fetch(url);
+    text = await res.text();
+    return JSON.parse(text);
+  } catch (error) {
+    console.log("text:", text);
+    throw error;
+  }
+}
+
+serve({ debug: true, max: 11, port: 3000, wait: 60 });
 
 const urlToGeojson = "http://localhost:3000/data/gadm/geojsons/Akrotiri and Dhekelia.geojson";
 
@@ -131,3 +144,15 @@ 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_json(urlToData + "antimeridian/right-edge.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/mean/mean.test.js

@@ -12,7 +12,7 @@ const { round } = utils;
 const url = "http://localhost:3000/data/test.tiff";
 
 const bbox = [80.63, 7.42, 84.21, 10.1];
-const expectedBboxValue = 1232.47;
+const expectedBboxValue = 1257.64;
 
 const polygon = JSON.parse(readFileSync("./data/part-of-india.geojson", "utf-8"));
 const expectedPolygonValue = 1853.71;

src/median/median.test.js

@@ -9,7 +9,7 @@ serve({ debug: true, max: 1, port: 3000 });
 const url = "http://localhost:3000/data/test.tiff";
 
 const bbox = [80.63, 7.42, 84.21, 10.1];
-const expectedBboxValue = 906.7;
+const expectedBboxValue = 915.85;
 
 const bboxGeojson = `{
   "type": "FeatureCollection",
@@ -26,7 +26,7 @@ const bboxGeojson = `{
     }
   }]
 }`;
-const expectedBboxGeojsonValue = 1849.4;
+const expectedBboxGeojsonValue = 1831.9;
 
 const polygon = JSON.parse(readFileSync("./data/part-of-india.geojson", "utf-8"));
 const expectedPolygonValue = 1637.35;

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 });
@@ -35,19 +38,19 @@ const EXPECTED_RASTER_STATS = [
 
 const EXPECTED_BBOX_STATS = [
   {
-    count: 213,
+    count: 188,
     invalid: 0,
     max: 5166.7,
-    mean: 1232.4718309859154,
-    median: 906.7,
+    mean: 1257.6351063829786,
+    median: 915.85,
     min: 0,
     mode: 0,
     modes: [0],
     range: 5166.7,
-    std: 1195.3529916721104,
-    sum: 262516.5,
-    valid: 213,
-    variance: 1428868.7746994647
+    std: 1216.4677587709607,
+    sum: 236435.4,
+    valid: 188,
+    variance: 1479793.808129244
   }
 ];
 
@@ -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";
@@ -29,7 +31,7 @@ const bbox = [80.63, 7.42, 84.21, 10.1];
 //   arr = src.read(1, window=from_bounds(80.63, 7.42, 84.21, 10.1, src.transform))
 //   arr[arr > 0].sum()
 // 236_435.40000000002
-const expectedBboxValue = 262_516.5;
+const expectedBboxValue = 236_435.4;
 
 const polygon = JSON.parse(readFileSync("./data/part-of-india.geojson", "utf-8"));
 

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;