Multipolygon producing unexpected sum - antimeridian crossing scenario

[twelch]

Describe the bug
Geoblaze.sum is returning an unexpected result.

I have a raster of sample data with data near Fiji, called fiji_anticross_meridian_test, EPSG:4326, where the data spans the 180 degree antimeridian. It is uncompressed.

And I have a multipolygon clip-test-clipped focused on Fiji that also crosses the 180 degree antimeridian and is already split (see screenshots of both pieces). Note, this multipolygon has a lot of holes from cutting out Fiji land.

To Reproduce
Zipped test input is attached - geoblaze-2.4-clip-test.zip

Steps to reproduce the behavior:

1. Run geoblaze.sum() with the test raster and multipolygon

> await gb.sum(rast, mp)
[ 111 ]

Expected behavior
I expected the result of geoblaze.sum() to be 3, but receive 111. From visual inspection, I believe the multipolygon overlaps with only a single cell of the raster, and that cell has a value of 3. However, a geoblaze.sum() produces a value of 111.

I'm not sure the antimeridian has anything to do with this. What I do notice is that 111/3 = 37. And there seems to be ~36 holes of the multipolygon overlapping with that one raster cell. I wonder if the hole overlap may be being counted the sum in addition to the outer ring. If that is expected behavior let me know.

Screenshots
Western portion of multipolygon

Easter portion of multipolygon (overlapping 1 cell with value of 3)

[DanielJDufour]

Hello. Appreciate you catching this. I can confirm the issue. The problem is here: https://github.com/GeoTIFF/geoblaze/blob/master/src/intersect-polygon/intersect-polygon.module.js#L54

Because geotiffs that cross the antimeridian generally use less than 1% of pixels in the overall extent, this triggers geoblaze to disaggregate the extent and sample around each sub-geometry. This is done to preserve memory.

There currently is no process for checking for overlap between sample areas because I have naively assumed the sub-geometries+samples would be far away from each other and not covering the same pixel.

Put simply, this wasn't a test case I had thought about and I'm very appreciative for you bringing this to my attention :-)

This can definitely be fixed by adding some union/de-confliction logic for sample areas, but may take a week or two.

Let me know if you need a fix on a quicker time table and I can push a patch that allows you to over-ride the fancy sub-sampling memory-safe logic.

[DanielJDufour]

technically got the union of some of the sub-geometry bboxes, but issue is that if the geometries are really small (not covering half a pixel), the size of the bbox will be 0-width or 0-height. I think we can get around the issue with this particular step by padding the snapping to the image space by 1 pixel on each side.

[twelch]

A week or two should be fine, thanks Daniel

[DanielJDufour]

Got tests fixed on my local. Solution was to add a few more logical steps:

• filtering out geometries that don't cross the bounds of the raster (e.g., the right-hand geometry in the test case technically is just outside the bounds of the raster, so is filtered out)
• union/combine sampling areas expressed in pixel coordinates (avoiding double counting pixels)
• recompute bboxes of sampling areas in spatial reference system coordinates after previous step

Still will probably take a week or so to clean up and republish a new version of geoblaze. Just wanted to keep you up to date :-)

[twelch]

Confirmed fixed in 2.5.0