Skip to content

Commit

Permalink
refactored intersect-polygon
Browse files Browse the repository at this point in the history
  • Loading branch information
@DanielJDufour
DanielJDufour committed Apr 28, 2024
1 parent 5191d49 commit 4fcae79
Show file tree
Hide file tree
Showing 2 changed files with 108 additions and 35 deletions.
66 changes: 48 additions & 18 deletions src/intersect-polygon/intersect-polygon.module.js
View file
Original file line number Diff line number Diff line change
@@ -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";
Expand All @@ -11,23 +13,63 @@ import snap from "snap-bbox";
import wrapParse from "../wrap-parse";
// import writePng from "@danieljdufour/write-png";

const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level = 0, vrm = [1, 1] } = {}) => {
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") vrm = [vrm, vrm]
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 = [
Math.ceil(fastMax(geometry_bboxes_multipliers.map(([xmul, ymul]) => xmul))),
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;
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 (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;
Expand All @@ -46,17 +88,6 @@ const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level =

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);

Expand Down Expand Up @@ -118,8 +149,6 @@ 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,
Expand All @@ -129,7 +158,7 @@ const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level =
pixel_width: georaster.pixelWidth / xvrm,
geometry
};
console.log("intersect_params:", intersect_params);
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));
Expand Down Expand Up @@ -164,6 +193,7 @@ const intersectPolygon = (georaster, geometry, perPixelFunction, { debug_level =
});
});
});
return { vrm };
});
};

Expand Down
77 changes: 60 additions & 17 deletions src/intersect-polygon/intersect-polygon.test.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,12 +13,17 @@ import intersectPolygon from "../intersect-polygon";

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;
const result = await intersectPolygon(
georaster,
geom,
value => {
if (includeNoData || value !== georaster.noDataValue) {
numberOfIntersectingPixels++;
}
},
{ debug_level, vrm }
);
return { numberOfIntersectingPixels, ...result };
}

async function fetch_json(url) {
Expand Down Expand Up @@ -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);
});

Expand All @@ -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);
});

Expand Down Expand Up @@ -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");
Expand All @@ -193,25 +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", async ({ eq }) => {
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 = await countIntersectingPixels(georaster, geom, false, { debug_level: 10, vrm: 100 });
// same as rasterstats
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", async ({ eq }) => {
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 = await countIntersectingPixels(georaster, geom, false, { debug_level: 10, vrm: 100 });
// same as rasterstats
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);
})
});

0 comments on commit 4fcae79

Please sign in to comment.