Skip to content

Commit

Permalink
sample twice when crossing antimeridian
Browse files Browse the repository at this point in the history
  • Loading branch information
@DanielJDufour
DanielJDufour committed Sep 17, 2023
1 parent 344d4c4 commit 19b1ff7
Show file tree
Hide file tree
Showing 6 changed files with 170 additions and 113 deletions.
5 changes: 4 additions & 1 deletion package.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,20 +59,23 @@
"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",
"mpoly": "^0.2.0",
"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",
Expand Down
153 changes: 91 additions & 62 deletions src/intersect-polygon/intersect-polygon.module.js
View file
Original file line number Diff line number Diff line change
@@ -1,36 +1,34 @@
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];

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,
Expand All @@ -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);
}
}
});
}
});
}
});
});
});
});
Expand Down
29 changes: 20 additions & 9 deletions src/intersect-polygon/intersect-polygon.test.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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";

Expand Down Expand Up @@ -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 }) => {
Expand All @@ -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);
});
58 changes: 30 additions & 28 deletions src/stats/stats.core.js
View file
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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";
Expand Down
Loading

0 comments on commit 19b1ff7

Please sign in to comment.