Merge pull request #221 from GeoTIFF/issue-212

Issue 212

data/gadm/download_gadm.sh

@@ -1,5 +1,6 @@
 if [ ! -f gadm28_levels.shp.zip ]; then
-    wget http://biogeo.ucdavis.edu/data/gadm2.8/gadm28_levels.shp.zip
+    # originally from http://biogeo.ucdavis.edu/data/gadm2.8/gadm28_levels.shp.zip
+    wget https://geoblaze.s3.amazonaws.com/gadm28_levels.shp.zip
 fi
 
 if [ ! -f gadm28_adm0.shp ]; then

data/geotiff-test-data/setup.sh

@@ -2,6 +2,7 @@
 # download from https://github.com/GeoTIFF/test-data/
 wget https://github.com/GeoTIFF/test-data/archive/refs/heads/main.zip -O geotiff-test-data.zip
 unzip -j -o geotiff-test-data.zip "test-data-*/files/*" -d .
-rm geotiff-test-data.zip
+unzip spam2005v3r2_harvested-area_wheat_total.tiff.zip
+rm geotiff-test-data.zip spam2005v3r2_harvested-area_wheat_total.tiff.zip
 
 gdalwarp -t_srs EPSG:4326 gadas.tif gadas-4326.tif

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

@@ -9,8 +9,6 @@ import wrapParse from "../wrap-parse";
 // import writePng from "@danieljdufour/write-png";
 
 const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level = 0 } = {}) => {
-  const { noDataValue } = georaster;
-
   const georaster_bbox = [georaster.xmin, georaster.ymin, georaster.xmax, georaster.ymax];
 
   const precisePixelHeight = georaster.pixelHeight.toString();
@@ -41,6 +39,7 @@ const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level =
     samples = [{ intersections, sample, offset: [0, 0] }];
   } else if (georaster.getValues) {
     const geometry_bboxes = calcAll(geometry);
+    if (debug_level >= 2) console.log("[geoblaze] geometry_bboxes:", geometry_bboxes);
 
     if (!geometry_bboxes.some(bbox => booleanIntersects(bbox, georaster_bbox))) return;
 
@@ -121,9 +120,7 @@ const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level =
                 const row = band[irow];
                 if (row) {
                   const value = row[icol];
-                  if (value != undefined && value !== noDataValue) {
-                    perPixelFunction(value, iband, yoff + irow, xoff + icol);
-                  }
+                  perPixelFunction(value, iband, yoff + irow, xoff + icol);
                 }
               });
             }

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

@@ -11,6 +11,16 @@ import parse from "../parse";
 import { convertMultiPolygon } from "../convert-geometry";
 import intersectPolygon from "../intersect-polygon";
 
+async function countIntersectingPixels(georaster, geom, includeNoData) {
+  let numberOfIntersectingPixels = 0;
+  await intersectPolygon(georaster, geom, value => {
+    if (includeNoData || value !== georaster.noDataValue) {
+      numberOfIntersectingPixels++;
+    }
+  });
+  return numberOfIntersectingPixels;
+}
+
 async function fetch_json(url) {
   let res;
   let text;
@@ -42,7 +52,7 @@ test("(Legacy) Testing intersection of box", async ({ eq }) => {
   values.forEach(band => {
     band.forEach(row => {
       row.forEach(value => {
-        if (value != georaster.noDataValue) {
+        if (value !== georaster.noDataValue) {
           expectedNumberOfIntersectingPixels++;
         }
       });
@@ -60,18 +70,16 @@ test("(Legacy) Testing intersection of box", async ({ eq }) => {
     georaster.ymax - 0.5 * pixelHeight
   ]);
   const coords = geom.geometry.coordinates;
-  let numberOfIntersectingPixels = 0;
-  intersectPolygon(georaster, coords, () => numberOfIntersectingPixels++);
+  const numberOfIntersectingPixels = await countIntersectingPixels(georaster, coords, false);
   eq(numberOfIntersectingPixels, expectedNumberOfIntersectingPixels);
 });
 
 test("(Legacy) Test intersection/sum calculations for Country with Multiple Rings", async ({ eq }) => {
   const georaster = await load(urlToData + "ghsl/tiles/GHS_POP_GPW42015_GLOBE_R2015A_54009_1k_v1_0_4326_30_40.tif");
   const response = await fetch(urlToGeojson);
   const country = await response.json();
-  let numberOfIntersectingPixels = 0;
   const geom = convertMultiPolygon(country);
-  intersectPolygon(georaster, geom, () => numberOfIntersectingPixels++);
+  const numberOfIntersectingPixels = await countIntersectingPixels(georaster, geom, true);
   eq(numberOfIntersectingPixels, EXPECTED_NUMBER_OF_INTERSECTING_PIXELS);
 });
 
@@ -85,9 +93,7 @@ test("(Modern) Testing intersection of box", async ({ eq }) => {
   values.forEach(band => {
     band.forEach(row => {
       row.forEach(value => {
-        if (value != georaster.noDataValue) {
-          expectedNumberOfIntersectingPixels++;
-        }
+        expectedNumberOfIntersectingPixels++;
       });
     });
   });
@@ -153,7 +159,11 @@ test("antimerdian #1", async ({ eq }) => {
   // reproject geometry to projection of raster
   geojson = reprojectGeoJSON(geojson, { from: 4326, to: georaster.projection });
   const geom = convertMultiPolygon(geojson);
-  await intersectPolygon(georaster, geom, () => numberOfIntersectingPixels++);
+  await intersectPolygon(georaster, geom, value => {
+    if (value !== georaster.noDataValue) {
+      numberOfIntersectingPixels++;
+    }
+  });
   // rasterstats says 314,930
   eq(numberOfIntersectingPixels, 314_930);
 });
@@ -182,9 +192,8 @@ test("parse no overlap", async ({ eq }) => {
 test("more testing", async ({ eq }) => {
   const georaster = await parse(urlToData + "test.tiff");
   const geojson = await fetch_json(urlToData + "part-of-india.geojson");
-  let numberOfIntersectingPixels = 0;
   const geom = convertMultiPolygon(geojson);
-  await intersectPolygon(georaster, geom, () => numberOfIntersectingPixels++);
+  const numberOfIntersectingPixels = await countIntersectingPixels(georaster, geom, false);
   // same as rasterstats
   eq(numberOfIntersectingPixels, 1_672);
 });

src/stats/stats.core.js

@@ -49,13 +49,18 @@ const stats = (georaster, geometry, calcStatsOptions, test, { debug_level = 0 }
       // the third argument of this function is a function which
       // runs for every pixel in the polygon. Here we add them to
       // an array, so we can later on calculate stats for each band separately
-      const done = intersectPolygon(georaster, geometry, (value, bandIndex) => {
-        if (values[bandIndex]) {
-          values[bandIndex].push(value);
-        } else {
-          values[bandIndex] = [value];
-        }
-      });
+      const done = intersectPolygon(
+        georaster,
+        geometry,
+        (value, bandIndex) => {
+          if (values[bandIndex]) {
+            values[bandIndex].push(value);
+          } else {
+            values[bandIndex] = [value];
+          }
+        },
+        { debug_level }
+      );
 
       return QuickPromise.resolve(done).then(() => {
         const bands = values.filter(band => band.length !== 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: 26, port: 3000, wait: 240 });
+serve({ debug: true, max: 30, port: 3000, wait: 240 });
 
 const url = "http://localhost:3000/data/test.tiff";
 
@@ -58,8 +58,8 @@ const EXPECTED_BBOX_STATS = [
 
 const EXPECTED_POLYGON_STATS = [
   {
-    count: 1672, // rasterstats says 1672
-    invalid: 0, // expected because geoblaze ignores no data values (but maybe this isn't ideal?)
+    count: 1722,
+    invalid: 50,
     max: 7807.4,
     mean: 1853.7104066985623,
     median: 1637.35,
@@ -70,7 +70,7 @@ const EXPECTED_POLYGON_STATS = [
     range: 7807.4,
     std: 1507.3255322956486,
     sum: 3_099_403.799999996, // rasterstats says 3,099,403.8
-    valid: 1672,
+    valid: 1672, // same as rasterstats
     variance: 2_272_030.2603103607
   }
 ];
@@ -187,11 +187,13 @@ test("antimerdian #2 (split at antimeridian)", async ({ eq }) => {
 test("antimeridian New Zealand EEZ and Habitat", async ({ eq }) => {
   const georaster = await parse("http://localhost:3000/data/geotiff-test-data/nz_habitat_anticross_4326_1deg.tif");
   const geojson = JSON.parse(readFileSync("./data/geojson-test-data/eez_land_union/EEZ_Land_v3_202030_New_Zealand.geojson", "utf-8"));
-  const results = await stats(georaster, geojson, { stats: ["count", "min", "max", "sum"] });
-  eq(results, [{ count: 454, min: 1, max: 71, sum: 4512 }]);
+  const results = await stats(georaster, geojson, { stats: ["count", "invalid", "min", "max", "sum", "valid"] });
+  eq(results, [{ count: 487, invalid: 33, min: 1, max: 71, sum: 4512, valid: 454 }]);
 
-  const results2 = await Promise.all(geojson.features[0].geometry.coordinates.map(geom => stats(georaster, geom, { stats: ["count", "min", "max", "sum"] })));
-  eq(results2, [[{ count: 322, min: 1, max: 71, sum: 3931 }], [{ count: 132, min: 1, max: 22, sum: 581 }]]);
+  const results2 = await Promise.all(
+    geojson.features[0].geometry.coordinates.map(geom => stats(georaster, geom, { stats: ["count", "invalid", "min", "max", "sum", "valid"] }))
+  );
+  eq(results2, [[{ count: 354, valid: 322, invalid: 32, min: 1, max: 71, sum: 3931 }], [{ count: 133, valid: 132, invalid: 1, min: 1, max: 22, sum: 581 }]]);
 });
 
 test("antimerdian #3 (across antimeridian on left-side)", async ({ eq }) => {
@@ -216,3 +218,20 @@ test("edge", async ({ eq }) => {
   const results = await stats(georaster, geojson, { stats: ["count", "min", "max", "sum"] });
   eq(results, [{ count: 2, min: 9.936111450195312, max: 19.24805450439453, sum: 29.184165954589844 }]);
 });
+
+test("multipolygon vs 2 polygons", async ({ eq }) => {
+  const geojson = JSON.parse(readFileSync("./data/antimeridian/split.geojson", "utf-8"));
+  const georaster = await parse("http://localhost:3000/data/geotiff-test-data/spam2005v3r2_harvested-area_wheat_total.tiff");
+  const _stats = ["count", "invalid", "min", "max", "sum", "valid"];
+
+  const expected = [{ count: 1152, valid: 4, invalid: 1148, min: 0, max: 0, sum: 0 }];
+  eq(await stats(georaster, geojson, { stats: _stats }, undefined, { debug_level: 5 }), expected);
+  eq(await stats(georaster, geojson.features[0], { stats: _stats }), expected);
+
+  const [poly1, poly2] = geojson.features[0].geometry.coordinates;
+  const results1 = await stats(georaster, poly1, { debug_level: 5, stats: _stats });
+  eq(results1, [{ count: 576, valid: 0, invalid: 576, sum: 0, min: undefined, max: undefined }]);
+
+  const results2 = await stats(georaster, poly2, { debug_level: 5, stats: _stats });
+  eq(results2, [{ count: 576, valid: 4, invalid: 572, min: 0, max: 0, sum: 0 }]);
+});