stac-utils / pystac-client

Python client for searching STAC APIs
https://pystac-client.readthedocs.io
Other
154 stars 47 forks source link

Add special handling for "Feature" dictionaries #694

Closed johntruckenbrodt closed 3 months ago

johntruckenbrodt commented 3 months ago

Hi there, I am trying to use the __geo_interface__ concept for defining a search geometry.
From how I interpret Sean Gillies' Post, I can define it as a regular GeoJSON:

geometry = {'type': 'Feature',
            'geometry': {'type': 'Polygon',
                         'coordinates': [[[10.0, 50.0], [10.0, 55.0],
                                          [15.0, 55.0], [15.0, 50.0], [10.0, 50.0]]]
                         }
            }

...which can then for example be read into geopandas:

df = GeoDataFrame.from_features(features=geometry)

However, when I pass this geometry definition to Client.search:

catalog = Client.open(url='https://planetarycomputer.microsoft.com/api/stac/v1')
result = catalog.search(collections='sentinel-1-grd',
                        max_items=None,
                        intersects=geometry,
                        datetime=['2022-01-01', '2022-01-02'])
result = list(result.items())
print(len(result))

...I get this error:

...
body -> intersects -> type
  unexpected value; permitted: 'Polygon' (type=value_error.const; given=Feature; permitted=['Polygon'])
...

So, I need to define my geometry like this to get a working example:

geometry = {'type': 'Polygon',
            'coordinates': [[[10.0, 50.0], [10.0, 55.0],
                             [15.0, 55.0], [15.0, 50.0],
                             [10.0, 50.0]]]}

Is this expected behavior? Is the latter a correct geometry as intended for __geo_interface__?

TomAugspurger commented 3 months ago

It's a little unclear to me whether "type": "Feature" is allowed in the STAC API spec or not: https://github.com/radiantearth/stac-api-spec/tree/release/v1.0.0/item-search#query-parameter-table. But looking at the jsonschema: https://github.com/radiantearth/stac-api-spec/blob/604ade6158de15b8ab068320ca41e25e2bf0e116/core/commons.yaml#L730-L738, I don't see anything for "type": "Feature"

So if you want your objects to be automatically usable as an intersects parameter then you should probably avoid the Feature type.

gadomski commented 3 months ago

I think the "GeoJSON Geometry" type unambiguously resolves to https://datatracker.ietf.org/doc/html/rfc7946#section-3.1, aka not Feature.

johntruckenbrodt commented 3 months ago

Thanks a lot for your answers @TomAugspurger and @gadomski.
Sorry for digging deeper, but I think this doesn't quite solve it...

Sean Gillies defined the __geo_interface__ with

type (required) A string indicating the geospatial type. Possible values are "Feature" or a geometry type: "Point", "LineString", "Polygon", etc. ... properties (optional) A mapping of feature properties (labels, populations ... you name it. Dependent on the data). Valid for "Feature" types only.

So the "type": "Feature" can be expected for __geo_interface__ returns.
Furthermore, the GeoJSON spec chapter 1.4 reads

  • Inside this document, the term "geometry type" refers to seven case-sensitive strings: "Point", "MultiPoint", "LineString", "MultiLineString", "Polygon", "MultiPolygon", and "GeometryCollection".
  • As another shorthand notation, the term "GeoJSON types" refers to nine case-sensitive strings: "Feature", "FeatureCollection", and the geometry types listed above.

__geo_interface__ thus supports every allowed GeoJSON type except FeatureCollection. So basically I think if pystac-client does not support "type": "Feature", it cannot entirely support the __geo_interface__ protocol.

Since a Feature is basically just a single simple geometry with a properties attribute, couldn't pystac-client just extract the geometry attribute from the object and use that for search?

gadomski commented 3 months ago

couldn't pystac-client just extract the geometry attribute from the object and use that for search?

It does, but only if the the value is not a dictionary or a string: https://github.com/stac-utils/pystac-client/blob/918d97eb5b0268894a9ae9c835924e5015e4fe4d/pystac_client/item_search.py#L645-L653

In your example, the value is a dictionary, so it is passed straight through. It would be a (probably good) enhancement to add handing for "Feature" ~and "FeatureCollection"~ to _format_intersects ... I'll update this issue's title to reflect that request.

johntruckenbrodt commented 3 months ago

Great, this is exactly what I had in mind. Thanks.