diff --git a/README.md b/README.md index 24276ab..6706a15 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,5 @@ -[![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*** diff --git a/data/virtual-resampling/virtual-resampling-intersect.geojson b/data/virtual-resampling/virtual-resampling-intersect.geojson new file mode 100644 index 0000000..3662829 --- /dev/null +++ b/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 ] ] ] } } +] +} diff --git a/data/virtual-resampling/virtual-resampling-one.geojson b/data/virtual-resampling/virtual-resampling-one.geojson new file mode 100644 index 0000000..c1dd305 --- /dev/null +++ b/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 ] ] ] } } +] +} diff --git a/package.json b/package.json index 80dc6b0..82b59ac 100644 --- a/package.json +++ b/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", diff --git a/src/intersect-polygon/intersect-polygon.module.js b/src/intersect-polygon/intersect-polygon.module.js index 9b8aef6..e702130 100644 --- a/src/intersect-polygon/intersect-polygon.module.js +++ b/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 }; }); }; diff --git a/src/intersect-polygon/intersect-polygon.test.js b/src/intersect-polygon/intersect-polygon.test.js index 247e295..0c3fa7a 100644 --- a/src/intersect-polygon/intersect-polygon.test.js +++ b/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); +}); diff --git a/src/max/index.js b/src/max/index.js index 4ef0bf5..b127dec 100644 --- a/src/max/index.js +++ b/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" }); } diff --git a/src/max/max.test.js b/src/max/max.test.js index 826fa23..c76aaac 100644 --- a/src/max/max.test.js +++ b/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]); +}); diff --git a/src/mean/index.js b/src/mean/index.js index ae2912e..f751bd4 100644 --- a/src/mean/index.js +++ b/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" }); } diff --git a/src/mean/mean.test.js b/src/mean/mean.test.js index 7c9907d..ad3659f 100644 --- a/src/mean/mean.test.js +++ b/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,32 @@ 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, [38]); +}); + +test("virtual resampling, intersecting 4 pixels (sync)", 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 georaster = await load(url); + const result = await mean(georaster, geojson); + eq(result, [38]); +}); + +test("virtual resampling with bbox", async ({ eq }) => { + const url = "http://localhost:3000/data/geotiff-test-data/nz_habitat_anticross_4326_1deg.tif"; + const geom = [166.06811846012403, -46.42907237702467, 166.23189176587618, -46.40887433963862]; + const result = await mean(url, geom); + eq(result, [38]); +}); diff --git a/src/median/index.js b/src/median/index.js index 7fbd4a6..5323c1f 100644 --- a/src/median/index.js +++ b/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: "minimal" }); } diff --git a/src/median/median.test.js b/src/median/median.test.js index d830ceb..e415387 100644 --- a/src/median/median.test.js +++ b/src/median/median.test.js @@ -4,7 +4,7 @@ import { serve } from "srvd"; import load from "../load"; import median from "."; -serve({ debug: true, max: 1, port: 3000 }); +serve({ debug: true, max: 6, port: 3000 }); const url = "http://localhost:3000/data/test.tiff"; @@ -78,3 +78,17 @@ test("Get Median from Whole Raster from url", async ({ eq }) => { const value = result[0]; eq(value, 0); }); + +test("Get Median with 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]); +}); diff --git a/src/min/index.js b/src/min/index.js index 8449dea..7fac1d3 100644 --- a/src/min/index.js +++ b/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: "minimal" }); } diff --git a/src/min/min.test.js b/src/min/min.test.js index e87f53a..99e0825 100644 --- a/src/min/min.test.js +++ b/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, wait: 120 }); 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, [38]); +}); diff --git a/src/mode/index.js b/src/mode/index.js index 88b5aed..ada8d6b 100644 --- a/src/mode/index.js +++ b/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: "minimal" }); } diff --git a/src/mode/mode.test.js b/src/mode/mode.test.js index 9c568fd..76dfee6 100644 --- a/src/mode/mode.test.js +++ b/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]); +}); diff --git a/src/modes/index.js b/src/modes/index.js index 448de33..b2e52b3 100644 --- a/src/modes/index.js +++ b/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: "minimal" }); } diff --git a/src/modes/modes.test.js b/src/modes/modes.test.js index d355d05..77c0d9d 100644 --- a/src/modes/modes.test.js +++ b/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]]); +}); diff --git a/src/range/index.js b/src/range/index.js index 56a4894..99c2a28 100644 --- a/src/range/index.js +++ b/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: "minimal" }); } diff --git a/src/range/range.test.js b/src/range/range.test.js index 608aeb7..86456ba 100644 --- a/src/range/range.test.js +++ b/src/range/range.test.js @@ -2,7 +2,9 @@ import test from "flug"; import { serve } from "srvd"; import range from "."; -serve({ debug: true, max: 20, port: 3000 }); +const sleep = ms => new Promise(resolve => setTimeout(() => resolve(), ms)); + +serve({ debug: true, max: 22, port: 3000, wait: 160 }); const url = "http://localhost:3000/data/test.tiff"; const bbox = [80.63, 7.42, 84.21, 10.1]; @@ -40,4 +42,22 @@ test("(Modern) Get Range from Polygon", async ({ eq }) => { test("(Modern) Get Range from whole Raster", async ({ eq }) => { const results = await range(url); eq(results, [8131.2]); + + // weird hack to give time for main thread and/or srvd to clear up tasks before next request is sent + await sleep(1000); +}); + +test("virtual resampling, contained", async ({ eq }) => { + const geojson_res = await fetch("http://localhost:3000/data/virtual-resampling/virtual-resampling-one.geojson"); + const geojson = await geojson_res.json(); + const url = "http://localhost:3000/data/geotiff-test-data/nz_habitat_anticross_4326_1deg.tif"; + 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, [0]); // range is zero because by default doing minimal (i.e. bare minimum) resampling }); diff --git a/src/stat/stat.module.js b/src/stat/stat.module.js index 1b2db34..4c87243 100644 --- a/src/stat/stat.module.js +++ b/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])); } diff --git a/src/stats/stats.core.js b/src/stats/stats.core.js index 6103d64..f5343de 100644 --- a/src/stats/stats.core.js +++ b/src/stats/stats.core.js @@ -1,12 +1,16 @@ import calcStats from "calc-stats"; import QuickPromise from "quick-promise"; +import polygon from "bbox-fns/polygon"; +import validate from "bbox-fns/validate"; import get from "../get"; import utils from "../utils"; 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 VRM_NO_RESAMPLING = [1, 1]; + +const stats = (georaster, geometry, calcStatsOptions, test, { debug_level = 0, include_meta = false, rescale = false, vrm = VRM_NO_RESAMPLING } = {}) => { try { // shallow clone calcStatsOptions = { ...calcStatsOptions }; @@ -18,6 +22,13 @@ const stats = (georaster, geometry, calcStatsOptions, test, { debug_level = 0 } calcStatsOptions.noData = noDataValue; } + 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]; + } + if (test) { if (calcStatsOptions && calcStatsOptions.filter) { const original_filter = calcStatsOptions.filter; @@ -30,11 +41,25 @@ const stats = (georaster, geometry, calcStatsOptions, test, { debug_level = 0 } const flat = true; const getStatsByBand = values => values.map(band => calcStats(band, calcStatsOptions)); - if (geometry === null || geometry === undefined) { + const resample = vrm === "minimal" || (Array.isArray(vrm) && (vrm[0] !== 1 || vrm[1] !== 1)); + const geometry_is_nullish = geometry === null || geometry === undefined; + + if (resample === true) { + if (geometry_is_nullish) { + geometry = polygon([georaster.xmin, georaster.ymin, georaster.xmax, georaster.ymax]); + } else if (validate(geometry)) { + geometry = polygon(geometry); + } else if (utils.isBboxObj(geometry)) { + // convert { xmin: 20, xmax: 32, ymin: -3, ymax: 0 } to geojson polygon + geometry = polygon([geometry.xmin, geometry.ymin, geometry.xmax, geometry.ymax]); + } + } + + if (geometry_is_nullish && resample === false) { if (debug_level >= 2) console.log("[geoblaze] geometry is nullish"); const values = get(georaster, undefined, flat); return QuickPromise.resolve(values).then(getStatsByBand); - } else if (utils.isBbox(geometry)) { + } else if (utils.isBbox(geometry) && resample === false) { if (debug_level >= 2) console.log("[geoblaze] geometry is a rectangle"); geometry = convertBbox(geometry); const values = get(georaster, geometry, flat); @@ -59,13 +84,88 @@ const stats = (georaster, geometry, calcStatsOptions, test, { debug_level = 0 } values[bandIndex] = [value]; } }, - { debug_level } + { debug_level, vrm } ); - return QuickPromise.resolve(done).then(() => { + return QuickPromise.resolve(done).then(({ vrm }) => { + // check if the user doesn't want the number of valid pixels returned + const want_valid = calcStatsOptions.stats ? calcStatsOptions.stats.includes("valid") : calcStatsOptions.calcValid !== false; + + const use_virtual_resampling = vrm[0] !== 1 && vrm[1] !== 1; + const bands = values.filter(band => band.length !== 0); - if (bands.length > 0) return bands.map(band => calcStats(band, calcStatsOptions)); - else throw "No Values were found in the given geometry"; + + if (bands.length > 0) { + // if we need to know the number of valid pixels for intermediate calculations, + // but the user didn't ask for it, calculate it anyway + // and then remove the valid count from the returned results + if (use_virtual_resampling) { + if (calcStatsOptions.stats) { + if (calcStatsOptions.stats.includes("product") && calcStatsOptions.stats.includes("valid") === false) { + calcStatsOptions.stats = [...calcStatsOptions.stats, "valid"]; + if (debug_level >= 2) console.log('[geoblaze] added "valid" to stats'); + } + } else { + if (calcStatsOptions.calcProduct === true && calcStatsOptions.calcValid === false) { + calcStatsOptions.calcValid = true; + if (debug_level >= 2) console.log('[geoblaze] set "calcValid" to true'); + } + } + } + + const stats = bands.map(band => calcStats(band, calcStatsOptions)); + if (debug_level >= 2) console.log("[geoblaze] stats (before rescaling):", stats); + + // only rescaling results if virtual resampling is on + if (use_virtual_resampling && rescale) { + if (debug_level >= 2) console.log("[geoblaze] rescaling results based on relative size of virtual pixels"); + + // if vrm is [2, 4] then area_multiplier will be 8, + // meaning there are 8 virtual pixels for every actual pixel + const area_multiplier = vrm[0] * vrm[1]; + if (debug_level >= 2) console.log("[geoblaze] area_multiplier:", area_multiplier); + + stats.forEach(band => { + const { valid } = band; + if (typeof band.count === "number") band.count /= area_multiplier; + if (typeof band.invalid === "number") band.invalid /= area_multiplier; + if (typeof band.sum === "number") band.sum /= area_multiplier; + if (typeof band.valid === "number") band.valid /= area_multiplier; + + if (band.histogram) { + for (const key in band.histogram) { + // this will lead to fractions of pixels + // for example, a pixel could appear 3.75 times if vrm is [2, 2] + band.histogram[key].ct /= area_multiplier; + } + } + + if (typeof band.product === "number") { + band.product /= Math.pow(area_multiplier, valid); + } + }); + } + + // if the user asked not to have the valid stat, + // exclude if from the results + if (want_valid === false) { + stats.forEach(band => delete band.valid); + } + + if (include_meta) { + stats.forEach(band => { + band._meta = { + vph: georaster.pixelHeight / vrm[1], + vpw: georaster.pixelWidth / vrm[0], + vrm: [vrm[0], vrm[1]] + }; + }); + } + + return stats; + } else { + throw "No Values were found in the given geometry"; + } }); } else { throw "Geometry Type is Not Supported"; diff --git a/src/stats/stats.test.js b/src/stats/stats.test.js index 3025307..585e651 100644 --- a/src/stats/stats.test.js +++ b/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: 30, port: 3000, wait: 240 }); +serve({ debug: true, max: 35, port: 3000, wait: 240 }); const url = "http://localhost:3000/data/test.tiff"; @@ -79,8 +79,9 @@ test("(Sync) Stats without Geometry", async ({ eq }) => { const georaster = await load(url); const results = stats(georaster, undefined); results.forEach(band => { - delete band.uniques; + delete band.frequency; delete band.histogram; + delete band.uniques; }); eq(results, EXPECTED_RASTER_STATS); }); @@ -88,6 +89,7 @@ test("(Sync) Stats without Geometry", async ({ eq }) => { test("(Async) Stats without Geometry", async ({ eq }) => { const results = await stats(url, undefined); results.forEach(band => { + delete band.frequency; delete band.histogram; delete band.uniques; }); @@ -98,6 +100,7 @@ test("(Sync) Stats with Bounding Box", async ({ eq }) => { const georaster = await load(url); const results = stats(georaster, bbox); results.forEach(band => { + delete band.frequency; delete band.histogram; delete band.uniques; }); @@ -107,6 +110,7 @@ test("(Sync) Stats with Bounding Box", async ({ eq }) => { test("(Async) Stats with Bounding Box", async ({ eq }) => { const results = await stats(url, bbox); results.forEach(band => { + delete band.frequency; delete band.histogram; delete band.uniques; }); @@ -117,6 +121,7 @@ test("(Sync) Stats with Polygon", async ({ eq }) => { const georaster = await load(url); const results = stats(georaster, polygon); results.forEach(band => { + delete band.frequency; delete band.histogram; delete band.uniques; }); @@ -126,6 +131,7 @@ test("(Sync) Stats with Polygon", async ({ eq }) => { test("(Async) Stats with Polygon", async ({ eq }) => { const results = await stats(url, polygon); results.forEach(band => { + delete band.frequency; delete band.histogram; delete band.uniques; }); @@ -234,13 +240,93 @@ test("multipolygon vs 2 polygons", async ({ eq }) => { 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, { stats: _stats }, undefined, { debug_level: 0 }), 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 }); + const results1 = await stats(georaster, poly1, { debug_level: 0, 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 }); + const results2 = await stats(georaster, poly2, { debug_level: 0, stats: _stats }); eq(results2, [{ count: 576, valid: 4, invalid: 572, min: 0, max: 0, sum: 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 stats(url, geojson, undefined, undefined, { debug_level: 0, rescale: true, vrm: "minimal" }); + eq(result, [ + { + count: 0.007936507936507936, + frequency: { 38: { n: 38, freq: 1 } }, + invalid: 0, + median: 38, + min: 38, + max: 38, + product: 0.022738725119677502, + sum: 0.30158730158730157, + range: 0, + mean: 38, + valid: 0.007936507936507936, + variance: 0, + std: 0, + histogram: { 38: { n: 38, ct: 0.007936507936507936 } }, + modes: [38], + mode: 38, + uniques: [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 = JSON.parse(readFileSync("./data/virtual-resampling/virtual-resampling-intersect.geojson", "utf-8")); + const result = await stats(url, geojson, null, null, { include_meta: false, rescale: true, vrm: [10, 10] }); + eq(result, [ + { + count: 0.18, + invalid: 0, + median: 38, + min: 1, + max: 38, + product: 1.5122637183654097e-10, + sum: 6.17, + range: 37, + mean: 34.27777777777778, + valid: 0.18, + variance: 112.20061728395062, + std: 10.592479279373201, + frequency: { + 1: { + freq: 0.05555555555555555, + n: 1 + }, + 8: { + freq: 0.05555555555555555, + n: 8 + }, + 38: { + freq: 0.8888888888888888, + n: 38 + } + }, + histogram: { + 1: { + n: 1, + ct: 0.01 + }, + 8: { + n: 8, + ct: 0.01 + }, + 38: { + n: 38, + ct: 0.16 + } + }, + modes: [38], + mode: 38, + uniques: [1, 8, 38] + } + ]); +}); diff --git a/src/sum/index.js b/src/sum/index.js index 38f3405..1751194 100644 --- a/src/sum/index.js +++ b/src/sum/index.js @@ -18,5 +18,5 @@ import stat from "../stat"; * [217461, 21375, 57312, 457125] */ export default function sum(georaster, geometry, test) { - return stat(georaster, geometry, "sum", test); + return stat(georaster, geometry, "sum", test, { rescale: true, vrm: "minimal" }); } diff --git a/src/sum/sum.test.js b/src/sum/sum.test.js index a9e85c2..7b17461 100644 --- a/src/sum/sum.test.js +++ b/src/sum/sum.test.js @@ -13,7 +13,7 @@ import sum from "."; const { fetchJson, fetchJsons } = utils; -const { port } = serve({ debug: true, max: 100, port: 8888, wait: 120 }); +const { port } = serve({ debug: true, max: 31, port: 8888, wait: 120 }); const urlRwanda = `http://localhost:${port}/data/RWA_MNH_ANC.tif`; const bboxRwanda = require("../../data/RwandaBufferedBoundingBox.json"); @@ -59,6 +59,7 @@ test("(Legacy) Get Sum from Veneto Geonode", async ({ eq }) => { ]; const [georaster, geojson] = values; const results = sum(georaster, geojson); + eq(Array.isArray(results), true); const actualValue = Number(results[0].toFixed(2)); const expectedValue = 24_943.31; // rasterstats says 24,963.465454101562 eq(actualValue, expectedValue); @@ -305,3 +306,13 @@ test("(Modern) Get Sum from Polygon Above X Value", async ({ eq }) => { const value = Number(sum(georaster, polygon, v => v > 3000)[0].toFixed(2)); eq(value, 1_454_066); }); + +test("Virtual Resampling", async ({ eq }) => { + const values = [ + await load(`http://localhost:${port}/data/geotiff-test-data/nz_habitat_anticross_4326_1deg.tif`), + await fetchJson(`http://localhost:${port}/data/virtual-resampling/virtual-resampling-one.geojson`) + ]; + const [georaster, geojson] = values; + const results = await sum(georaster, geojson, undefined); + eq(results, [0.30158730158730157]); +}); diff --git a/src/utils/utils.module.js b/src/utils/utils.module.js index 91255e6..5e5034b 100644 --- a/src/utils/utils.module.js +++ b/src/utils/utils.module.js @@ -122,6 +122,18 @@ const utils = { return Array.isArray(parsed) ? parsed[0] : parsed; }, + // checks if bbox in form { xmin: 20, xmax: 32, ymin: -3, ymax: 0 } + isBboxObj(geometry) { + return ( + typeof geometry === "object" && + Array.isArray(geometry) === false && + typeof geometry.xmin === "number" && + typeof geometry.xmax === "number" && + typeof geometry.ymin === "number" && + typeof geometry.ymax === "number" + ); + }, + isPolygonal(geometry) { const polys = mpoly.get(geometry);