GEOSException when reprojecting

[pjhartzell]

Summary

When loading and then reprojecting data in a sinusoidal (MODIS and VIIRS products) tile whose boundary spans the edge of the projection (and, hence, the antimeridian), a GEOS Exception is raised, e.g, this:

ds = odc.stac.load(
    [item],
    chunks={"x": 2048, "y": 2048},
).odc.reproject("EPSG:32603")

produces an error:

GEOSException: TopologyException: side location conflict at -140.51384784042935 59.245084530288445. This can occur if the input geometry is invalid.

• The data does not cross the antimeridian, only the tile bounds

Workarounds (unsatisfactory to current needs)

• If I remove the chunks option, reprojection is successful. However, we are using Dask and need to specify chunk size.
• If I reproject in the load call by passing the crs (rather than the separate .odc.reproject), there is no error. However, we moved to a load and then reproject order of operations to work around a small data gap when loading adjacent sinusoidal tiles.

Minimal Example

A minimal STAC Item and Jupyter notebook recreating the problem is here.

Versions, hardware

Pinned to: odc-stac==0.3.5 odc-geo==0.4.3 odc-algo==0.2.3

Tested with Python 3.9

Apple M1 (arm64)

[Kirill888]

@pjhartzell that would be odc-geo issue then, reproject code is there

[Kirill888]

That's a common family of errors we have no good solution for currently. Essentially we need a more robust mechanism for checking overlap of two geometries defined in two different projections, keeping in mind that not all vertices of A can be mapped to B and same applies to B vertices going to A.

This is a common operation and needs to be efficient in the common case, yet allowing for a more costly but more robust check if we detect issues with geometries after projection while doing a cheap version of the check.

[Kirill888]

STAC for that tile didn't get the footprint right either, as it spans lon=180 line.

[pjhartzell]

@Kirill888 Thanks for the info.

Yep, that geometry is incorrect. Should be a MultiPolygon with a Polygon on each side of the antimeridian.