Merge pull request #228 from GeoTIFF/issue-56-virtual-resampling

Virtual Resampling

README.md

@@ -1,8 +1,5 @@
 <img src="http://geoblaze.io/assets/img/logo.svg">
 
-[![Maintainability](https://api.codeclimate.com/v1/badges/a99a88d28ad37a79dbf6/maintainability)](https://codeclimate.com/github/codeclimate/codeclimate/maintainability)
-[![Test Coverage](https://api.codeclimate.com/v1/badges/a99a88d28ad37a79dbf6/test_coverage)](https://codeclimate.com/github/codeclimate/codeclimate/test_coverage)
-
 # GeoBlaze
 
 ***A blazing fast javascript raster processing engine***

data/virtual-resampling/virtual-resampling-intersect.geojson

@@ -0,0 +1,8 @@
+{
+"type": "FeatureCollection",
+"name": "virtual-resampling-intersect",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
+"features": [
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 165.939001402474787, -46.486236933721152 ], [ 166.303394074387001, -46.121247183868078 ], [ 166.581806452926685, -46.392534642381364 ], [ 165.996321598056511, -46.660679777789575 ], [ 165.939001402474787, -46.486236933721152 ] ] ] } }
+]
+}

data/virtual-resampling/virtual-resampling-one.geojson

@@ -0,0 +1,8 @@
+{
+"type": "FeatureCollection",
+"name": "virtual-resampling-one",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
+"features": [
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 166.128374676391331, -46.408874339638622 ], [ 166.231891765876185, -46.50319516961985 ], [ 166.051123117074297, -46.506566757293193 ], [ 166.068118460124026, -46.429072377024667 ], [ 166.128374676391331, -46.408874339638622 ] ] ] } }
+]
+}

package.json

@@ -29,7 +29,17 @@
     "format": "npx prettier --arrow-parens=avoid --print-width=160 --trailing-comma=none --write lite/*.js src/*.js src/*/*.js ",
     "lint": "eslint lite && eslint src",
     "serve": "npx srvd --debug --port=3000 --wait=120",
-    "test": "set -e; for f in src/*/*test*.js; do echo \"\nrunning $f\" && sleep 5 && node -r esm $f; done",
+    "test": "set -e; for f in src/*/*test*.js; do echo \"\nrunning $f\" && sleep 5 && TEST_GAP_TIME=1 node -r esm $f; done",
+    "test:intersect-polygon": "node -r esm ./src/intersect-polygon/intersect-polygon.test.js",
+    "test:max": "node -r esm ./src/max/max.test.js",
+    "test:mean": "node -r esm ./src/mean/mean.test.js",
+    "test:median": "node -r esm ./src/median/median.test.js",
+    "test:min": "node -r esm ./src/min/min.test.js",
+    "test:mode": "node -r esm ./src/mode/mode.test.js",
+    "test:modes": "node -r esm ./src/modes/modes.test.js",
+    "test:range": "node -r esm ./src/range/range.test.js",
+    "test:stats": "node -r esm ./src/stats/stats.test.js",
+    "test:sum": "node -r esm ./src/sum/sum.test.js",
     "test-loading-builds": "node test-loading-builds.js",
     "setup": "bash setup.sh"
   },
@@ -114,7 +124,7 @@
     "eslint-plugin-standard": "^3.1.0",
     "esm": "3.2.25",
     "find-and-read": "^1.2.0",
-    "flug": "^2.6.0",
+    "flug": "^2.8.2",
     "jsdoc": "^3.6.7",
     "memory-fs": "^0.4.1",
     "null-loader": "^4.0.1",

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

@@ -1,7 +1,9 @@
 import bboxArea from "bbox-fns/bbox-area.js";
+import bboxSize from "bbox-fns/bbox-size.js";
 import booleanIntersects from "bbox-fns/boolean-intersects.js";
 import calcAll from "bbox-fns/calc-all.js";
 import dufour_peyton_intersection from "dufour-peyton-intersection";
+import fastMax from "fast-max";
 import merge from "bbox-fns/merge.js";
 import union from "bbox-fns/union.js";
 import reproject from "bbox-fns/precise/reproject.js";
@@ -11,18 +13,64 @@ import snap from "snap-bbox";
 import wrapParse from "../wrap-parse";
 // import writePng from "@danieljdufour/write-png";
 
-const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level = 0 } = {}) => {
+const VRM_NO_RESAMPLING = [1, 1];
+
+const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level = 0, vrm = VRM_NO_RESAMPLING } = {}) => {
   const georaster_bbox = [georaster.xmin, georaster.ymin, georaster.xmax, georaster.ymax];
 
   const precisePixelHeight = georaster.pixelHeight.toString();
   const precisePixelWidth = georaster.pixelWidth.toString();
 
+  if (typeof vrm === "number") {
+    if (vrm <= 0 || vrm !== Math.round(vrm)) {
+      throw new Error("[geoblaze] vrm can only be defined as a positive integer");
+    }
+    vrm = [vrm, vrm];
+  }
+
+  let geometry_bboxes = union(calcAll(geometry));
+  if (debug_level >= 2) console.log("[geoblaze] geometry_bboxes:", geometry_bboxes);
+
+  // filter out geometry bboxes that don't intersect raster
+  // assuming anti-meridian geometries have been normalized beforehand
+  geometry_bboxes = geometry_bboxes.filter(bbox => booleanIntersects(bbox, georaster_bbox));
+  if (debug_level >= 2) console.log("[geoblaze] intersecting geometry bboxes:", geometry_bboxes);
+
+  // no intersecting pixels
+  if (geometry_bboxes.length === 0) return;
+
+  // calculate size of bboxes in pixels
+  const geometry_bboxes_sizes = geometry_bboxes.map(bbox => bboxSize(bbox));
+  if (debug_level >= 2) console.log("[geoblaze] size of geometry bboxes:", geometry_bboxes_sizes);
+
+  const geometry_bbox_size_ratios = geometry_bboxes_sizes.map(([width, height]) => [width / georaster.pixelWidth, height / georaster.pixelHeight]);
+  if (debug_level >= 2) console.log("[geoblaze] geometry_bbox_size_ratios:", geometry_bbox_size_ratios);
+
+  // resample just enough to ensure at least one resampled pixel intersects
+  if (vrm === "minimal") {
+    // check if any geometry is smaller than half a pixel
+    if (geometry_bbox_size_ratios.some(([xratio, yratio]) => xratio <= 1 || yratio <= 1)) {
+      const geometry_bboxes_multipliers = geometry_bbox_size_ratios.map(([xratio, yratio]) => [2 / xratio, 2 / yratio]);
+      vrm = [
+        // don't drop more than 10,000 sample lines per pixel
+        Math.min(10000, Math.ceil(fastMax(geometry_bboxes_multipliers.map(([xmul, ymul]) => xmul)))),
+        Math.min(10000, Math.ceil(fastMax(geometry_bboxes_multipliers.map(([xmul, ymul]) => ymul))))
+      ];
+    } else {
+      vrm = VRM_NO_RESAMPLING;
+    }
+  }
+
+  if (debug_level >= 2) console.log("[geoblaze] vrm:", vrm);
+  const [xvrm, yvrm] = vrm;
+
   // run intersect for each sample
   // each sample is a multi-dimensional array of numbers
   // in the following xdim format [b][r][c]
   let samples;
 
   if (georaster.values) {
+    if (debug_level >= 2) console.log("[geoblaze] georaster already has all values loaded into memory");
     // if we have already loaded all the values into memory,
     // just pass those along and avoid using up more memory
     const sample = georaster.values;
@@ -32,26 +80,15 @@ const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level =
     const intersections = dufour_peyton_intersection.calculate({
       debug: false,
       raster_bbox: georaster_bbox,
-      raster_height: georaster.height,
-      raster_width: georaster.width,
-      pixel_height: georaster.pixelHeight,
-      pixel_width: georaster.pixelWidth,
+      raster_height: georaster.height * yvrm,
+      raster_width: georaster.width * xvrm,
+      pixel_height: georaster.pixelHeight / yvrm,
+      pixel_width: georaster.pixelWidth / xvrm,
       geometry
     });
 
     samples = [{ intersections, sample, offset: [0, 0] }];
   } else if (georaster.getValues) {
-    let geometry_bboxes = union(calcAll(geometry));
-    if (debug_level >= 2) console.log("[geoblaze] geometry_bboxes:", geometry_bboxes);
-
-    // filter out geometry bboxes that don't intersect raster
-    // assuming anti-meridian geometries have been normalized beforehand
-    geometry_bboxes = geometry_bboxes.filter(bbox => booleanIntersects(bbox, georaster_bbox));
-    if (debug_level >= 2) console.log("[geoblaze] intersecting geometry bboxes:", geometry_bboxes);
-
-    // no intersecting pixels
-    if (geometry_bboxes.length === 0) return;
-
     const geotransfrom = PreciseGeotransform([georaster.xmin.toString(), precisePixelWidth, "0", georaster.ymax.toString(), "0", "-" + precisePixelHeight]);
     if (debug_level >= 2) console.log("[geoblaze] geotransfrom:", geotransfrom);
 
@@ -116,12 +153,13 @@ const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level =
       const intersect_params = {
         debug: true,
         raster_bbox: sample_bbox,
-        raster_height: sample_height,
-        raster_width: sample_width,
-        pixel_height: georaster.pixelHeight,
-        pixel_width: georaster.pixelWidth,
+        raster_height: sample_height * yvrm,
+        raster_width: sample_width * xvrm,
+        pixel_height: georaster.pixelHeight / yvrm,
+        pixel_width: georaster.pixelWidth / xvrm,
         geometry
       };
+      if (debug_level >= 3) console.log("[geoblaze] intersect_params:", intersect_params);
 
       const intersections = dufour_peyton_intersection.calculate(intersect_params);
       if (debug_level >= 3) console.log("[geoblaze] intersections:", JSON.stringify(intersections, undefined, 2));
@@ -145,9 +183,9 @@ const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level =
           row.forEach(([start, end], irange) => {
             for (let icol = start; icol <= end; icol++) {
               imageBands.forEach((band, iband) => {
-                const row = band[irow];
+                const row = band[Math.floor(irow / yvrm)];
                 if (row) {
-                  const value = row[icol];
+                  const value = row[Math.floor(icol / xvrm)];
                   perPixelFunction(value, iband, yoff + irow, xoff + icol);
                 }
               });
@@ -156,6 +194,7 @@ const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level =
         });
       });
     });
+    return { vrm };
   });
 };
 

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

@@ -11,14 +11,19 @@ import parse from "../parse";
 import { convertMultiPolygon } from "../convert-geometry";
 import intersectPolygon from "../intersect-polygon";
 
-async function countIntersectingPixels(georaster, geom, includeNoData) {
+async function countIntersectingPixels(georaster, geom, includeNoData, { debug_level, vrm } = {}) {
   let numberOfIntersectingPixels = 0;
-  await intersectPolygon(georaster, geom, value => {
-    if (includeNoData || value !== georaster.noDataValue) {
-      numberOfIntersectingPixels++;
-    }
-  });
-  return numberOfIntersectingPixels;
+  const result = await intersectPolygon(
+    georaster,
+    geom,
+    value => {
+      if (includeNoData || value !== georaster.noDataValue) {
+        numberOfIntersectingPixels++;
+      }
+    },
+    { debug_level, vrm }
+  );
+  return { numberOfIntersectingPixels, ...result };
 }
 
 async function fetch_json(url) {
@@ -34,7 +39,7 @@ async function fetch_json(url) {
   }
 }
 
-serve({ debug: true, max: 25, port: 3000, wait: 60 });
+serve({ debug: true, max: 30, port: 3000, wait: 60 });
 
 const urlToGeojson = "http://localhost:3000/data/gadm/geojsons/Akrotiri and Dhekelia.geojson";
 
@@ -70,7 +75,7 @@ test("(Legacy) Testing intersection of box", async ({ eq }) => {
     georaster.ymax - 0.5 * pixelHeight
   ]);
   const coords = geom.geometry.coordinates;
-  const numberOfIntersectingPixels = await countIntersectingPixels(georaster, coords, false);
+  const { numberOfIntersectingPixels } = await countIntersectingPixels(georaster, coords, false);
   eq(numberOfIntersectingPixels, expectedNumberOfIntersectingPixels);
 });
 
@@ -79,7 +84,7 @@ test("(Legacy) Test intersection/sum calculations for Country with Multiple Ring
   const response = await fetch(urlToGeojson);
   const country = await response.json();
   const geom = convertMultiPolygon(country);
-  const numberOfIntersectingPixels = await countIntersectingPixels(georaster, geom, true);
+  const { numberOfIntersectingPixels } = await countIntersectingPixels(georaster, geom, true);
   eq(numberOfIntersectingPixels, EXPECTED_NUMBER_OF_INTERSECTING_PIXELS);
 });
 
@@ -168,6 +173,26 @@ test("antimerdian #1", async ({ eq }) => {
   eq(numberOfIntersectingPixels, 314_930);
 });
 
+test("antimerdian #1, vrm=10", async ({ eq }) => {
+  const georaster = await parse(urlToData + "gfwfiji_6933_COG.tiff");
+  let geojson = await fetch_json(urlToData + "antimeridian/right-edge.geojson");
+  let numberOfIntersectingPixels = 0;
+  // reproject geometry to projection of raster
+  geojson = reprojectGeoJSON(geojson, { from: 4326, to: georaster.projection });
+  const geom = convertMultiPolygon(geojson);
+  await intersectPolygon(
+    georaster,
+    geom,
+    value => {
+      if (value !== georaster.noDataValue) {
+        numberOfIntersectingPixels++;
+      }
+    },
+    { vrm: 10 }
+  );
+  eq(numberOfIntersectingPixels, 31_501_375);
+});
+
 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");
@@ -193,7 +218,43 @@ test("more testing", async ({ eq }) => {
   const georaster = await parse(urlToData + "test.tiff");
   const geojson = await fetch_json(urlToData + "part-of-india.geojson");
   const geom = convertMultiPolygon(geojson);
-  const numberOfIntersectingPixels = await countIntersectingPixels(georaster, geom, false);
+  const { numberOfIntersectingPixels } = await countIntersectingPixels(georaster, geom, false);
   // same as rasterstats
   eq(numberOfIntersectingPixels, 1_672);
 });
+
+test("virtual resampling with vrm = minimal", async ({ eq }) => {
+  const georaster = await parse(urlToData + "/geotiff-test-data/nz_habitat_anticross_4326_1deg.tif");
+  const geojson = await fetch_json(urlToData + "/virtual-resampling/virtual-resampling-one.geojson");
+  const geom = convertMultiPolygon(geojson);
+  const { numberOfIntersectingPixels, vrm } = await countIntersectingPixels(georaster, geom, false, { debug_level: 0, vrm: "minimal" });
+  eq(vrm, [12, 21]);
+  eq(numberOfIntersectingPixels, 2);
+});
+
+test("virtual resampling with vrm = 100", async ({ eq }) => {
+  const georaster = await parse(urlToData + "/geotiff-test-data/nz_habitat_anticross_4326_1deg.tif");
+  const geojson = await fetch_json(urlToData + "/virtual-resampling/virtual-resampling-one.geojson");
+  const geom = convertMultiPolygon(geojson);
+  const { numberOfIntersectingPixels, vrm } = await countIntersectingPixels(georaster, geom, false, { debug_level: 0, vrm: 100 });
+  eq(vrm, [100, 100]);
+  eq(numberOfIntersectingPixels, 104);
+});
+
+test("virtual resampling intersect with vrm = minimal", async ({ eq }) => {
+  const georaster = await parse(urlToData + "/geotiff-test-data/nz_habitat_anticross_4326_1deg.tif");
+  const geojson = await fetch_json(urlToData + "/virtual-resampling/virtual-resampling-intersect.geojson");
+  const geom = convertMultiPolygon(geojson);
+  const { numberOfIntersectingPixels, vrm } = await countIntersectingPixels(georaster, geom, false, { debug_level: 0, vrm: "minimal" });
+  eq(vrm, [4, 4]);
+  eq(numberOfIntersectingPixels, 2);
+});
+
+test("virtual resampling intersect with vrm = 100", async ({ eq }) => {
+  const georaster = await parse(urlToData + "/geotiff-test-data/nz_habitat_anticross_4326_1deg.tif");
+  const geojson = await fetch_json(urlToData + "/virtual-resampling/virtual-resampling-intersect.geojson");
+  const geom = convertMultiPolygon(geojson);
+  const { numberOfIntersectingPixels, vrm } = await countIntersectingPixels(georaster, geom, false, { debug_level: 0, vrm: 100 });
+  eq(vrm, [100, 100]);
+  eq(numberOfIntersectingPixels, 1577);
+});

src/max/index.js

@@ -19,5 +19,5 @@ import stat from "../stat";
  */
 
 export default function max(georaster, geometry, test) {
-  return stat(georaster, geometry, "max", test);
+  return stat(georaster, geometry, "max", test, { vrm: "minimal" });
 }

src/max/max.test.js

@@ -5,7 +5,7 @@ import { serve } from "srvd";
 import load from "../load";
 import max from ".";
 
-serve({ debug: true, max: 1, port: 3000 });
+serve({ debug: true, max: 7, port: 3000, wait: 240 });
 
 const url = "http://localhost:3000/data/test.tiff";
 
@@ -61,3 +61,17 @@ test("Max with Web Mercator Bounding Box and GeoRaster URL", async ({ eq }) => {
   const value = Number(result[0].toFixed(2));
   eq(value, expectedBboxValue);
 });
+
+test("virtual resampling, contained", async ({ eq }) => {
+  const url = "http://localhost:3000/data/geotiff-test-data/nz_habitat_anticross_4326_1deg.tif";
+  const geojson = await fetch("http://localhost:3000/data/virtual-resampling/virtual-resampling-one.geojson").then(res => res.json());
+  const result = await max(url, geojson);
+  eq(result, [38]);
+});
+
+test("virtual resampling, intersecting 4 pixels", async ({ eq }) => {
+  const url = "http://localhost:3000/data/geotiff-test-data/nz_habitat_anticross_4326_1deg.tif";
+  const geojson = await fetch("http://localhost:3000/data/virtual-resampling/virtual-resampling-intersect.geojson").then(res => res.json());
+  const result = await max(url, geojson);
+  eq(result, [38]);
+});

src/mean/index.js

@@ -19,5 +19,5 @@ import stat from "../stat";
  */
 
 export default function mean(georaster, geometry, test) {
-  return stat(georaster, geometry, "mean", test);
+  return stat(georaster, geometry, "mean", test, { vrm: "minimal" });
 }