WORK IN PROGRESS

package.json

@@ -30,6 +30,14 @@
     "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: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:range": "node -r esm ./src/range/range.test.js",
+    "test:sum": "node -r esm ./src/sum/sum.test.js",
     "test-loading-builds": "node test-loading-builds.js",
     "setup": "bash setup.sh"
   },

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

@@ -11,18 +11,23 @@ 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 intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level = 0, vrm = [1, 1] } = {}) => {
   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") vrm = [vrm, vrm]
+  const [xvrm, yvrm] = vrm;
+  console.log({xvrm, yvrm});
+
   // 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) {
+    console.log("georaster.values:", georaster.values);
     // if we have already loaded all the values into memory,
     // just pass those along and avoid using up more memory
     const sample = georaster.values;
@@ -113,15 +118,18 @@ const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level =
 
       const sample_bbox = precise_sample_bbox.map(str => Number(str));
 
+      console.log({vrm});
+      console.log(JSON.stringify(geometry));
       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
       };
+      console.log("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 +153,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);
                 }
               });

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

@@ -11,13 +11,13 @@ 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++;
     }
-  });
+  }, { debug_level, vrm });
   return numberOfIntersectingPixels;
 }
 
@@ -197,3 +197,21 @@ test("more testing", async ({ eq }) => {
   // same as rasterstats
   eq(numberOfIntersectingPixels, 1_672);
 });
+
+test("virtual resampling", 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 = await countIntersectingPixels(georaster, geom, false, { debug_level: 10, vrm: 100 });
+  // same as rasterstats
+  eq(numberOfIntersectingPixels, 104);
+})
+
+test("virtual resampling intersect", 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 = await countIntersectingPixels(georaster, geom, false, { debug_level: 10, vrm: 100 });
+  // same as rasterstats
+  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: 100 });
 }

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: 5, port: 3000 });
 
 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: 100 });
 }

src/mean/mean.test.js

@@ -5,7 +5,7 @@ import load from "../load";
 import mean from ".";
 import utils from "../utils";
 
-serve({ debug: true, max: 2, port: 3000 });
+serve({ debug: true, max: 4, port: 3000 });
 
 const { round } = utils;
 
@@ -91,3 +91,17 @@ test("(Modern) Mean from GeoJSON", async ({ eq }) => {
   const value = round(results[0]);
   eq(value, expectedPolygonGeojsonValue);
 });
+
+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 mean(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 mean(url, geojson)
+  eq(result, [33.61065313887127]);  
+})

src/median/index.js

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

src/median/median.test.js

@@ -78,3 +78,17 @@ test("Get Median from Whole Raster from url", async ({ eq }) => {
   const value = result[0];
   eq(value, 0);
 });
+
+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 median(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 median(url, geojson)
+  eq(result, [38]);  
+})

src/min/index.js

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

src/min/min.test.js

@@ -3,7 +3,7 @@ import { serve } from "srvd";
 import load from "../load";
 import min from ".";
 
-serve({ debug: true, max: 1, port: 3000 });
+serve({ debug: true, max: 5, port: 3000 });
 
 const url = "http://localhost:3000/data/test.tiff";
 const bbox = [80.63, 7.42, 84.21, 10.1];
@@ -61,3 +61,17 @@ test("(Modern) Get Min from whole Raster", async ({ eq }) => {
   const results = await min(url);
   eq(results, [expectedPolygonValue]);
 });
+
+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 min(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 min(url, geojson)
+  eq(result, [1]);  
+})

src/mode/index.js

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

src/mode/mode.test.js

@@ -4,7 +4,7 @@ import { serve } from "srvd";
 import load from "../load";
 import mode from ".";
 
-serve({ debug: true, max: 1, port: 3000 });
+serve({ debug: true, max: 5, port: 3000 });
 
 const url = "http://localhost:3000/data/test.tiff";
 
@@ -44,3 +44,17 @@ test("(Modern) Mode Polygon", async ({ eq }) => {
   const results = await mode(url, polygon);
   eq(results, [0]);
 });
+
+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 mode(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 mode(url, geojson)
+  eq(result, [38]);  
+})

src/modes/index.js

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

src/modes/modes.test.js

@@ -44,3 +44,17 @@ test("(Modern) Modes Polygon", async ({ eq }) => {
   const results = await modes(url, polygon);
   eq(results, [[0]]);
 });
+
+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 modes(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 modes(url, geojson)
+  eq(result, [38]);
+})

src/range/index.js

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

src/range/range.test.js

@@ -2,7 +2,7 @@ import test from "flug";
 import { serve } from "srvd";
 import range from ".";
 
-serve({ debug: true, max: 20, port: 3000 });
+serve({ debug: true, max: 25, port: 3000 });
 
 const url = "http://localhost:3000/data/test.tiff";
 const bbox = [80.63, 7.42, 84.21, 10.1];
@@ -41,3 +41,17 @@ test("(Modern) Get Range from whole Raster", async ({ eq }) => {
   const results = await range(url);
   eq(results, [8131.2]);
 });
+
+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 range(url, geojson)
+  eq(result, [0]);  
+})
+
+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 range(url, geojson)
+  eq(result, [37]);  
+})

src/stat/stat.module.js

@@ -1,6 +1,6 @@
 import QuickPromise from "quick-promise";
 import stats from "../stats";
 
-export default function stat(georaster, geometry, stat, test) {
-  return QuickPromise.resolve(stats(georaster, geometry, { stats: [stat] }, test)).then(stats => stats.map(it => it[stat]));
+export default function stat(georaster, geometry, stat, test, options) {
+  return QuickPromise.resolve(stats(georaster, geometry, { stats: [stat] }, test, options)).then(stats => stats.map(it => it[stat]));
 }

src/stats/stats.core.js

@@ -6,7 +6,7 @@ import wrap from "../wrap-parse";
 import { convertBbox, convertMultiPolygon } from "../convert-geometry";
 import intersectPolygon from "../intersect-polygon";
 
-const stats = (georaster, geometry, calcStatsOptions, test, { debug_level = 0 } = {}) => {
+const stats = (georaster, geometry, calcStatsOptions, test, { debug_level = 0, vrm } = {}) => {
   try {
     // shallow clone
     calcStatsOptions = { ...calcStatsOptions };
@@ -18,6 +18,20 @@ const stats = (georaster, geometry, calcStatsOptions, test, { debug_level = 0 }
       calcStatsOptions.noData = noDataValue;
     }
 
+    let xvrm;
+    let yvrm;
+    if (typeof vrm === "number") {
+      if (vrm !== Math.round(vrm)) {
+        throw new Error("[geoblaze] divisor must be an integer");
+      }
+      xvrm = vrm;
+      yvrm = vrm;
+    } else if (Array.isArray(vrm) && vrm.length === 2 && typeof vrm[0] === "number") {
+      [xvrm, yvrm] = vrm;
+    }
+
+    console.log("vrm:", [xvrm, yvrm]);
+
     if (test) {
       if (calcStatsOptions && calcStatsOptions.filter) {
         const original_filter = calcStatsOptions.filter;
@@ -37,6 +51,7 @@ const stats = (georaster, geometry, calcStatsOptions, test, { debug_level = 0 }
     } else if (utils.isBbox(geometry)) {
       if (debug_level >= 2) console.log("[geoblaze] geometry is a rectangle");
       geometry = convertBbox(geometry);
+      // if using multiplier, might need to pad get results or at least not round
       const values = get(georaster, geometry, flat);
       return QuickPromise.resolve(values).then(getStatsByBand);
     } else if (utils.isPolygonal(geometry)) {
@@ -59,7 +74,7 @@ const stats = (georaster, geometry, calcStatsOptions, test, { debug_level = 0 }
             values[bandIndex] = [value];
           }
         },
-        { debug_level }
+        { debug_level, vrm }
       );
 
       return QuickPromise.resolve(done).then(() => {

src/sum/index.js

@@ -17,6 +17,6 @@ import stat from "../stat";
  * // results is [red, green, blue, nir]
  * [217461, 21375, 57312, 457125]
  */
-export default function sum(georaster, geometry, test) {
-  return stat(georaster, geometry, "sum", test);
+export default function sum(georaster, geometry, test, { vrm } = {}) {
+  return stat(georaster, geometry, "sum", test, { vrm });
 }