TopologyException happens while loading data for a polygon region

[Bis-Bala]

Background

I'm using agdc_statistics to process data from all of time for polygonal regions around Australia.

Expected behaviour

The datacube.utils.geometry.unary_union() function is being used to Union the valid_data polygons from my source datasets (~1900), it should return a valid polygon for all of the sources.

The relevant line of code is:

sources_union = geometry.unary_union(source.extent.to_crs(geobox.crs) for source in sources)

Actual behaviour

A TopologyException is raised instead of a valid shape being returned.

ERROR 1: TopologyException: found non-noded intersection between LINESTRING (1.7e+06 -2.20922e+06, 1.7e+06 -2.23041e+06) and LINESTRING (1.7e+06 -2.3e+06, 1.7e+06 -2.2087e+06) at 1699999.9999999993 -2210399.5855833939

Steps to reproduce the behaviour

Perform a unary union on the valid data regions from these dataset ids: datasetids.txt

OR

Perform a unary union on these geometries: source_geoms_wkt.txt

The exception occurs at geometry.py#L644 which is a call to the underlying ogr library.

union = geom.UnionCascaded()

Environment information

• Which datacube --version are you using? 1.3.2
• What datacube deployment/enviornment are you running against? Running on the NCI, loaded from agdc-py3-prod

[sixy6e]

The reported line segment looks like it would be a self intersecting line, hence the bad topology error. You could try simplifying the geometry first.

I tried the cascaded_union function direct from shapely and it succeeded, so maybe switch to using that?

from shapely.ops import cascaded_union
from shapely import wkt
with open('source_geoms_wkt.txt', 'r') as src:
    geometry = wkt.loads(src.readlines()[0])

polygon = cascaded_union(geometry)
print(polygon.wkt)

[Bis-Bala]

Currently I am using ConvexHull of the geometry instead of UnionCascaded and seems to be working. I tried cascaded_union from shapely and didn't get any topolgy exception. But downstream ogr's Geometry_Intersection complains of a type error as it is looking for type OGRGeometryShadow while creating geom from ogr.

[omad]

There's a couple of usable workarounds here.

At some point we may improve the vector geometry tools. They do currently work, but with some unavoidable limitations due to the nature of vector geometries.