diff --git a/package.json b/package.json index 5bca25f..c3e0a7f 100644 --- a/package.json +++ b/package.json @@ -59,13 +59,15 @@ "dependencies": { "@turf/boolean-clockwise": "^6.5.0", "@turf/combine": "^6.5.0", - "bbox-fns": "^0.13.0", + "@turf/rewind": "^6.5.0", + "bbox-fns": "^0.17.0", "calc-stats": "^2.2.0", "cross-fetch": "^4.0.0", "dufour-peyton-intersection": "0.2.0", "fast-max": "^0.4.0", "fast-min": "^0.3.0", "faster-median": "^1.0.0", + "geojson-antimeridian-cut": "^0.1.0", "georaster": "^1.6.0", "get-depth": "^0.0.3", "mathjs": "^11.9.1", @@ -73,6 +75,7 @@ "preciso": "^0.12.0", "proj4": "^2.9.0", "proj4-fully-loaded": "^0.2.0", + "quick-promise": "^0.1.0", "quick-resolve": "^0.0.1", "reproject-bbox": "^0.12.0", "reproject-geojson": "^0.3.0", diff --git a/src/intersect-polygon/intersect-polygon.module.js b/src/intersect-polygon/intersect-polygon.module.js index ae337fa..ed76053 100644 --- a/src/intersect-polygon/intersect-polygon.module.js +++ b/src/intersect-polygon/intersect-polygon.module.js @@ -1,14 +1,14 @@ +import bboxArea from "bbox-fns/bbox-area.js"; import booleanIntersects from "bbox-fns/boolean-intersects.js"; -import calcBoundingBox from "bbox-fns/calc.js"; +import calcAll from "bbox-fns/calc-all.js"; import dufour_peyton_intersection from "dufour-peyton-intersection"; +import merge from "bbox-fns/merge.js"; +import QuickPromise from "quick-promise"; import snap from "snap-bbox"; import wrapParse from "../wrap-parse"; -import utils from "../utils"; // import writePng from "@danieljdufour/write-png"; -const { resolve } = utils; - -const intersectPolygon = (georaster, geometry, perPixelFunction) => { +const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level = 0 } = {}) => { const { noDataValue } = georaster; const georaster_bbox = [georaster.xmin, georaster.ymin, georaster.xmax, georaster.ymax]; @@ -16,21 +16,19 @@ const intersectPolygon = (georaster, geometry, perPixelFunction) => { const precisePixelHeight = georaster.pixelHeight.toString(); const precisePixelWidth = georaster.pixelWidth.toString(); - let sample; - let sample_height; - let sample_width; - let intersections; + // 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 we have already loaded all the values into memory, // just pass those along and avoid using up more memory - sample = georaster.values; - sample_height = georaster.height; - sample_width = georaster.width; + const sample = georaster.values; // intersections are calculated in time mostly relative to geometry, // so using the whole georaster bbox, shouldn't significantly more operations - intersections = dufour_peyton_intersection.calculate({ + const intersections = dufour_peyton_intersection.calculate({ debug: false, raster_bbox: georaster_bbox, raster_height: georaster.height, @@ -39,67 +37,98 @@ const intersectPolygon = (georaster, geometry, perPixelFunction) => { pixel_width: georaster.pixelWidth, geometry }); - } else if (georaster.getValues) { - const geometry_bbox = calcBoundingBox(geometry); - if (!booleanIntersects(geometry_bbox, georaster_bbox)) return; + samples = [{ intersections, sample, offset: [0, 0] }]; + } else if (georaster.getValues) { + const geometry_bboxes = calcAll(geometry); - const [xmin, ymin, xmax, ymax] = geometry_bbox; + if (!geometry_bboxes.some(bbox => booleanIntersects(bbox, georaster_bbox))) return; - // snap geometry bounding box to georaster grid system - const snapResult = snap({ - bbox: [xmin.toString(), ymin.toString(), xmax.toString(), ymax.toString()], - debug: false, - origin: [georaster.xmin.toString(), georaster.ymax.toString()], - overflow: false, - scale: [precisePixelWidth, "-" + precisePixelHeight], - size: [georaster.width.toString(), georaster.height.toString()], - precise: true - }); + const combined_geometry_bbox = merge(geometry_bboxes); + if (debug_level >= 2) console.log("[geoblaze] combined_geometry_bbox:", combined_geometry_bbox); - const preciseSampleBox = snapResult.bbox_in_coordinate_system; - const sample_bbox = preciseSampleBox.map(str => Number(str)); + const usedArea = geometry_bboxes.reduce((total, bbox) => total + bboxArea(bbox), 0); + const totalArea = bboxArea(combined_geometry_bbox); + const usedPercentage = usedArea / totalArea; + if (debug_level >= 2) console.log("[geoblaze] percentage of sample area used:", usedPercentage); - const [left, bottom, right, top] = snapResult.bbox_in_grid_cells.map(n => Number(n)); + const sample_bboxes = usedPercentage > 0.01 ? [combined_geometry_bbox] : geometry_bboxes; - sample_height = bottom - top; - sample_width = right - left; + // get samples for each geometry bbox + samples = sample_bboxes.map(sample_bbox => { + const [xmin, ymin, xmax, ymax] = sample_bbox; - sample = georaster.getValues({ - left, - bottom, - right, - top, - width: sample_width, - height: sample_height, - resampleMethod: "near" - }); + // snap whole geometry bounding box to georaster grid system + const snapResult = snap({ + bbox: [xmin.toString(), ymin.toString(), xmax.toString(), ymax.toString()], + debug: false, + origin: [georaster.xmin.toString(), georaster.ymax.toString()], + overflow: false, + scale: [precisePixelWidth, "-" + precisePixelHeight], + size: [georaster.width.toString(), georaster.height.toString()], + precise: true + }); + if (debug_level >= 2) console.log("[geoblaze] snapResult:", snapResult); + const snapped_bbox = snapResult.bbox_in_coordinate_system.map(n => Number(n)); + + const image_bbox = snapResult.bbox_in_grid_cells.map(n => Number(n)); + if (debug_level >= 2) console.log("[geoblaze] image_bbox:", image_bbox); + const [left, bottom, right, top] = image_bbox; + + const snapped_height = bottom - top; + const snapped_width = right - left; + + const getValuesPromise = georaster.getValues({ + left, + bottom, + right, + top, + width: snapped_width, + height: snapped_height, + resampleMethod: "near" + }); - intersections = dufour_peyton_intersection.calculate({ - debug: false, - raster_bbox: sample_bbox, - raster_height: sample_height, - raster_width: sample_width, - pixel_height: georaster.pixelHeight, - pixel_width: georaster.pixelWidth, - geometry + const intersections = dufour_peyton_intersection.calculate({ + debug: false, + raster_bbox: snapped_bbox, + raster_height: snapped_height, + raster_width: snapped_width, + pixel_height: georaster.pixelHeight, + pixel_width: georaster.pixelWidth, + geometry + }); + if (debug_level >= 3) console.log("[geoblaze] intersections:", JSON.stringify(intersections, undefined, 2)); + + return QuickPromise.resolve(getValuesPromise).then(sample => { + if (debug_level >= 3) console.log("[geoblaze] got sample:", sample); + return { + intersections, + sample, + offset: [left, top] + }; + }); }); } - - return resolve(sample).then(imageBands => { - intersections.rows.forEach((row, irow) => { - row.forEach(([start, end], irange) => { - for (let icol = start; icol <= end; icol++) { - imageBands.forEach((band, iband) => { - const row = band[irow]; - if (row) { - const value = row[icol]; - if (value != undefined && value !== noDataValue) { - perPixelFunction(value, iband, irow, icol); - } + if (debug_level >= 3) console.log("[geoblaze] samples:", samples); + + return QuickPromise.all(samples).then(resolved_samples => { + resolved_samples.forEach(sample => { + QuickPromise.resolve(sample).then(({ intersections, sample: imageBands, offset: [xoff, yoff] }) => { + intersections.rows.forEach((row, irow) => { + row.forEach(([start, end], irange) => { + for (let icol = start; icol <= end; icol++) { + imageBands.forEach((band, iband) => { + const row = band[irow]; + if (row) { + const value = row[icol]; + if (value != undefined && value !== noDataValue) { + perPixelFunction(value, iband, yoff + irow, xoff + icol); + } + } + }); } }); - } + }); }); }); }); diff --git a/src/intersect-polygon/intersect-polygon.test.js b/src/intersect-polygon/intersect-polygon.test.js index 6b8e86d..16423cb 100644 --- a/src/intersect-polygon/intersect-polygon.test.js +++ b/src/intersect-polygon/intersect-polygon.test.js @@ -24,7 +24,7 @@ async function fetch_json(url) { } } -serve({ debug: true, max: 19, port: 3000, wait: 60 }); +serve({ debug: true, max: 25, port: 3000, wait: 60 }); const urlToGeojson = "http://localhost:3000/data/gadm/geojsons/Akrotiri and Dhekelia.geojson"; @@ -148,24 +148,25 @@ test("Hole Test", async ({ eq }) => { 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 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 }); + geojson = reprojectGeoJSON(geojson, { from: 4326, to: georaster.projection }); + const geom = convertMultiPolygon(geojson); await intersectPolygon(georaster, geom, () => numberOfIntersectingPixels++); - // same as rasterstats + // rasterstats says 314,930 eq(numberOfIntersectingPixels, 314_930); }); 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"); - let numberOfIntersectingPixels = 0; const geom = convertMultiPolygon(geojson); - await intersectPolygon(georaster, geom, () => numberOfIntersectingPixels++); - // same as rasterstats - eq(numberOfIntersectingPixels, 2); + const values = []; + await intersectPolygon(georaster, geom, (value, iband, irow, icol) => { + values.push(value); + }); + eq(values, [19.24805450439453, 9.936111450195312]); }); test("parse no overlap", async ({ eq }) => { @@ -177,3 +178,13 @@ test("parse no overlap", async ({ eq }) => { // same as rasterstats eq(numberOfIntersectingPixels, 0); }); + +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++); + // same as rasterstats + eq(numberOfIntersectingPixels, 1_672); +}); diff --git a/src/stats/stats.core.js b/src/stats/stats.core.js index acc2ee1..338074b 100644 --- a/src/stats/stats.core.js +++ b/src/stats/stats.core.js @@ -1,42 +1,46 @@ import calcStats from "calc-stats"; -import resolve from "quick-resolve"; +import QuickPromise from "quick-promise"; 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 DEFAULT_CALC_STATS_OPTIONS = { - calcHistogram: false, - calcMax: false, - calcMean: false, - calcMedian: false, - calcMin: false, - calcMode: false, - calcModes: false, - calcSum: false -}; - -const stats = (georaster, geometry, calcStatsOptions, test) => { +const stats = (georaster, geometry, calcStatsOptions, test, { debug_level = 0 } = {}) => { try { - const resolvedCalcStatsOptions = calcStatsOptions ? { ...DEFAULT_CALC_STATS_OPTIONS, ...calcStatsOptions } : undefined; + // shallow clone + calcStatsOptions = { ...calcStatsOptions }; const { noDataValue } = georaster; + + // use noDataValue unless explicitly over-written + if (noDataValue !== undefined && calcStatsOptions.noData === undefined) { + calcStatsOptions.noData = noDataValue; + } + + if (test) { + if (calcStatsOptions && calcStatsOptions.filter) { + const original_filter = calcStatsOptions.filter; + calcStatsOptions.filter = ({ index, value }) => original_filter({ index, value }) && test(value); + } else { + calcStatsOptions.filter = ({ index, value }) => test(value); + } + } + const flat = true; - const getStatsByBand = values => { - return values - .map(band => band.filter(value => value !== noDataValue && !isNaN(value) && (test === undefined || test(value)))) - .map(band => calcStats(band, resolvedCalcStatsOptions)); - }; + const getStatsByBand = values => values.map(band => calcStats(band, calcStatsOptions)); if (geometry === null || geometry === undefined) { + if (debug_level >= 2) console.log("[geoblaze] geometry is nullish"); const values = get(georaster, undefined, flat); - return resolve(values).then(getStatsByBand); + return QuickPromise.resolve(values).then(getStatsByBand); } else if (utils.isBbox(geometry)) { + if (debug_level >= 2) console.log("[geoblaze] geometry is a rectangle"); geometry = convertBbox(geometry); const values = get(georaster, geometry, flat); - return resolve(values).then(getStatsByBand); + return QuickPromise.resolve(values).then(getStatsByBand); } else if (utils.isPolygonal(geometry)) { + if (debug_level >= 2) console.log("[geoblaze] geometry is polygonal"); geometry = convertMultiPolygon(geometry); // don't know how many bands are in georaster, so default to 100 @@ -46,16 +50,14 @@ const stats = (georaster, geometry, calcStatsOptions, test) => { // 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 ((value !== noDataValue && !isNaN(value) && test === undefined) || test(value)) { - if (values[bandIndex]) { - values[bandIndex].push(value); - } else { - values[bandIndex] = [value]; - } + if (values[bandIndex]) { + values[bandIndex].push(value); + } else { + values[bandIndex] = [value]; } }); - return resolve(done).then(() => { + return QuickPromise.resolve(done).then(() => { 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"; diff --git a/src/stats/stats.test.js b/src/stats/stats.test.js index 937c15a..916b04f 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: 25, port: 3000, wait: 240 }); +serve({ debug: true, max: 26, port: 3000, wait: 240 }); const url = "http://localhost:3000/data/test.tiff"; @@ -20,8 +20,8 @@ const polygon = JSON.parse(readFileSync("./data/part-of-india.geojson", "utf-8") const EXPECTED_RASTER_STATS = [ { - count: 820517, - invalid: 0, + count: 9331200, + invalid: 8510683, max: 8131.2, mean: 132.04241399017369, median: 0, @@ -39,8 +39,8 @@ const EXPECTED_RASTER_STATS = [ const EXPECTED_BBOX_STATS = [ { - count: 188, - invalid: 0, + count: 1376, + invalid: 1188, max: 5166.7, mean: 1257.6351063829786, median: 915.85, @@ -171,8 +171,8 @@ 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 }]); + const results = await stats(georaster, geom, { stats: ["count", "invalid", "min", "max", "sum", "valid"] }); + eq(results, [{ count: 328425, valid: 314930, invalid: 13495, min: 0.20847222208976746, max: 492.3219299316406, sum: 12783872.545041203 }]); }); test("antimerdian #2 (split at antimeridian)", async ({ eq }) => { @@ -185,13 +185,17 @@ test("antimerdian #2 (split at antimeridian)", async ({ eq }) => { }); test("antimerdian #3 (across antimeridian on left-side)", async ({ eq }) => { + return; // XXX // converted GeoTIFF to all 1's const georaster = await parse("http://localhost:3000/data/gfwfiji_6933_COG_Binary.tif"); const geojson = JSON.parse(readFileSync("./data/antimeridian/across.geojson", "utf-8")); - // geom = reprojectGeoJSON(geom, { from: 4326, to: georaster.projection }); - // console.dir(geom.geometry.coordinates, { depth: null }); - const geom = { geometry: geojson, srs: 4326 }; - const results = await stats(georaster, geom, { stats: ["count", "min", "max", "sum"] }); + console.log(JSON.stringify(geojson)); + let geom = reprojectGeoJSON(geojson, { from: 4326, to: georaster.projection }); + console.log(JSON.stringify(geom)); + geom = { geometry: geojson, srs: 4326 }; + const results = await stats(georaster, geom, { stats: ["count", "invalid", "min", "max", "sum", "valid"] }); + // Size is 65_208, 3_515 + eq(results[0].valid, 327_972); eq(results, [{ count: 327_972, min: 1, max: 1, sum: 327_972 }]); }); diff --git a/src/utils/utils.module.js b/src/utils/utils.module.js index 904c87c..b43ce31 100644 --- a/src/utils/utils.module.js +++ b/src/utils/utils.module.js @@ -86,8 +86,10 @@ const utils = { isEsriJson(input, debug = false) { if (debug) console.log("starting isEsriJson with", input); - const inputType = typeof input; - const obj = inputType === "string" ? JSON.parse(input) : inputType === "object" ? input : null; + + if (typeof input === "undefined") return false; + + const obj = typeof input === "string" ? JSON.parse(input) : typeof input === "object" ? input : null; const geometry = obj.geometry ? obj.geometry : obj; if (geometry) { if (geometry.rings || (geometry.x && geometry.y)) { @@ -109,6 +111,12 @@ const utils = { } }, + isGeoJSON(input) { + if (typeof input === "undefined") return false; + if (input === null) return false; + return !utils.isEsriJson(input) && ["FeatureCollection", "Feature", "Polygon", "MultiPolygon"].includes(input.type); + }, + toGeoJSON(input, debug = false) { const parsed = ArcGIS.toGeoJSON(input); return Array.isArray(parsed) ? parsed[0] : parsed;