`lonlat_bounds` sometimes causes error "linearring requires at least 4 coordinates"

Ariana-B commented 4 months ago

Calling lonlat_bounds on a valid Polygon sometimes correctly outputs a BoundingBox, but other times causes shapely to error: ValueError: A linearring requires at least 4 coordinates. Running with python 3.11, odc-geo 0.4.5

>>> import odc.geo.geom as geom
>>> p = geom.polygon([[-1632000.0, -960000.0], [-1632000.0, -1056000.0], [-1536000.0, -1056000.0], [-1536000.0, -960000.0], [-1632000.0, -960000.0]], crs='epsg:3577')
>>> p
Geometry(POLYGON ((-1632000 -960000, -1632000 -1056000, -1536000 -1056000, -1536000 -960000, -1632000 -960000)), EPSG:3577)
>>> p.area
9216000000.0
>>> geom.lonlat_bounds(p)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/ubuntu/mambaforge/envs/cubeenv/lib/python3.11/site-packages/odc/geo/geom.py", line 1415, in lonlat_bounds
    bbox = geom.to_crs("EPSG:4326", check_and_fix=True).boundingbox
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/ubuntu/mambaforge/envs/cubeenv/lib/python3.11/site-packages/odc/geo/geom.py", line 742, in to_crs
    return maybe_fix(geom._to_crs(crs))
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/ubuntu/mambaforge/envs/cubeenv/lib/python3.11/site-packages/odc/geo/geom.py", line 729, in maybe_fix
    g = g.dropna()
        ^^^^^^^^^^
  File "/home/ubuntu/mambaforge/envs/cubeenv/lib/python3.11/site-packages/odc/geo/geom.py", line 1017, in dropna
    return self.filter(lambda x, y: math.isfinite(x) and math.isfinite(y))
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/ubuntu/mambaforge/envs/cubeenv/lib/python3.11/site-packages/odc/geo/geom.py", line 986, in filter
    return polygon(pts, self.crs, *inners)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/ubuntu/mambaforge/envs/cubeenv/lib/python3.11/site-packages/odc/geo/geom.py", line 1184, in polygon
    return Geometry({"type": "Polygon", "coordinates": (outer,) + inners}, crs=crs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/ubuntu/mambaforge/envs/cubeenv/lib/python3.11/site-packages/odc/geo/geom.py", line 503, in __init__
    self.geom = _geojson_to_shapely(geom)
                ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/ubuntu/mambaforge/envs/cubeenv/lib/python3.11/site-packages/odc/geo/geom.py", line 436, in _geojson_to_shapely
    return to_geom(xx)
           ^^^^^^^^^^^
  File "/home/ubuntu/mambaforge/envs/cubeenv/lib/python3.11/site-packages/odc/geo/geom.py", line 420, in to_geom
    return geometry.shape(force_2d(x))
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/ubuntu/mambaforge/envs/cubeenv/lib/python3.11/site-packages/shapely/geometry/geo.py", line 102, in shape
    return Polygon(ob["coordinates"][0], ob["coordinates"][1:])
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/ubuntu/mambaforge/envs/cubeenv/lib/python3.11/site-packages/shapely/geometry/polygon.py", line 230, in __new__
    shell = LinearRing(shell)
            ^^^^^^^^^^^^^^^^^
  File "/home/ubuntu/mambaforge/envs/cubeenv/lib/python3.11/site-packages/shapely/geometry/polygon.py", line 104, in __new__
    geom = shapely.linearrings(coordinates)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/ubuntu/mambaforge/envs/cubeenv/lib/python3.11/site-packages/shapely/decorators.py", line 77, in wrapped
    return func(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^
  File "/home/ubuntu/mambaforge/envs/cubeenv/lib/python3.11/site-packages/shapely/creation.py", line 173, in linearrings
    return lib.linearrings(coords, out=out, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: A linearring requires at least 4 coordinates.

But on second and subsequent tries (no change to p), it works:

>>> geom.lonlat_bounds(p)
BoundingBox(left=117.57871740949655, bottom=-9.311556020945595, right=118.50974150539429, top=-8.321682036216568, crs=CRS('EPSG:4326'))

Kirill888 commented 4 months ago

@Ariana-B that sounds like broken python environment to me, pyproj shipped libs clashing with rasterio maybe? This creates all sorts of nasty import order dependent issues. Could also be missing proj database or proj config pointing to wrong proj data dir.

Those points are nowhere near any discontinuities. So there is an error in code: we should check that enough points survive transform from one projection to another before attempting to reconstruct original shape class in .dropna() method, BUT equally, this should not be happening for this data and suggests an issue with pyproj in your environment.

Real issue is: .dropna()/.filter(..) should return empty geometry when not enough points are finite for geometry type being processed.

Ariana-B commented 3 months ago

You're right, it seems to have been a pyproj issue; upgrading to 3.6.1 appears to have resolved the error. Perhaps it would be worth updating the requirements, as they list pyproj>=3.0.0.

Kirill888 commented 3 months ago

We are not using any new features of pyproj, most likely it's pyproj vs PROJ lib version that got fixed in the update.

opendatacube / odc-geo

`lonlat_bounds` sometimes causes error "linearring requires at least 4 coordinates" #160