Convert coordinates to storage crs when filtering via cql

[ricardogsilva]

Overview

This PR implements conversion of geometries supplied in CQL filters to the storage CRS of a feature collection.

The main part of the implementation follows this logic:

1. Iterate through the nodes of the previously-parsed pygeofilter filter
2. If a node is of type pygeofilter.values.Geometry, retrieve its CRS information or default to the OGC:CRS84 CRS if none is provided - CRS information is gotten from the filter-crs query parameter, if present - it can also be provided in the geometry definition, using EWKT notation
3. Depending on the CRS of the collection, as designated in the storage_crs configuration option of the relevant provider, determine whether it is necessary to transform the geometry to match the same CRS and the one used by the provider
4. If needed, convert the filter geometry's coordinates using a combination of pyproj and shapely, in similar fashion as the approach already being used elsewhere in the pygeoapi codebase
5. Update the pygeofilter node tree with the converted geometry, which will subsequently be used by the provider to query its datasource

This PR includes the following related modifications:

• pygeoapi.util.get_crs_from_uri() has been refactored in order to support both URN and URL types of URI. This was done in order to accomodate the fact that pygeofilter is using the URN form internally

• the logic to iterate through the pygeofilter node tree was previously only available inside the postgres provider. It has been moved out to the pygeoapi.util.modify_pygeofilter() function. This function now performs two checks:

  • Translates from the generic geometry name to whatever name the provider specifies in its geom_field configuration parameter - This is the same functionality introduced recently in Translate between generic geometry queryable name and the actual geom field name in postgres db #1453
  • Converts coordinates to the appropriate storage CRS, i.e. the main purpose of this PR

This PR enables the use cases mentioned in #1488, which is basically to use whatever CRS when defining a geometry in a CQL filter, and default to assuming OGC:CRS84.

Example requests:

# no CRS specified - pygeoapi assumes OGC:CRS84 by default
curl http://localhost:5000/collections/beni_confiscati/items?\
    filter=DWITHIN(geometry,POINT(12.626000006069875 41.8108000011414),100,meters)

# CRS specified by means of the filter-crs query parameter
curl http://localhost:5000/collections/beni_confiscati/items?\
    filter=DWITHIN(geometry,POINT(2313683.65%204641312.29),100,meters)&\
    filter-crs=http://www.opengis.net/def/crs/EPSG/0/3004

# CRS specified by defining the filter geometry as EWKT
curl http://localhost:5000/collections/beni_confiscati/items?\
    filter=DWITHIN(geometry,SRID=3004;POINT(2313683.65%204641312.29),100,meters)

Note that regardless of the CRS being defined in the request, this PR makes pygeoapi ensure that it will always use the CRS of the underlying collection (as defined in its storage_crs configuration parameter), converting the filter geometry to the correct CRS if needed.

Related issue / discussion

• fixes Enable CQL filtering to accept filter geometries in any CRS #1488

Dependency policy (RFC2)

• I have ensured that this PR meets RFC2 requirements

Updates to public demo

• I have ensured that breaking changes to the pygeoapi master demo server have been addressed
  • https://github.com/geopython/demo.pygeoapi.io/blob/master/services/pygeoapi_master/local.config.yml

Contributions and licensing

(as per https://github.com/geopython/pygeoapi/blob/master/CONTRIBUTING.md#contributions-and-licensing)

• I'd like to contribute [feature X|bugfix Y|docs|something else] to pygeoapi. I confirm that my contributions to pygeoapi will be compatible with the pygeoapi license guidelines at the time of contribution
• I have already previously agreed to the pygeoapi Contributions and Licensing Guidelines

[ricardogsilva]

The failing CI is unrelated to the changes contained in this PR, as mentioned in #1490

[walkermatt]

@ricardogsilva Hi 👋. I think the standards based approach to specifying the CRS of geometries in a CQL filter is via the filter-crs parameter as per https://portal.ogc.org/files/96288#filter-filter-crs.

[ricardogsilva]

hi @walkermatt thanks for the pointer 👍

I'll take a look at the spec and implement accordingly

[walkermatt]

@ricardogsilva fantastic. Thanks for implementing support for EWKT in geopython/pygeofilter, it's a great help.

[ricardogsilva]

@walkermatt I've now updated this PR to use a filter-crs query parameter, as mentioned in the OAPIF - Part 3 spec you mentioned

[walkermatt]

@walkermatt I've now updated this PR to use a filter-crs query parameter, as mentioned in the OAPIF - Part 3 spec you mentioned

I've just updated my local copy of pygeoapi to use this PR branch and it's working well for me. I have a table in Postgres/PostGIS in EPSG:27700 with storage_crs set to http://www.opengis.net/def/crs/EPSG/0/27700. I'm now able to use the WITHIN CQL function like so:

http://localhost:5000/collections/licensed_area/items?org=example-com&f=json&filter=WITHIN(POINT(558760.2888089775%20139992.53573017253),geometry)&filter-crs=http://www.opengis.net/def/crs/EPSG/0/27700&crs=http://www.opengis.net/def/crs/EPSG/0/27700&limit=1

[francbartoli]

@ricardogsilva found a small bug if the SRID is explicitly sent in the EWKT:

# no CRS specified - pygeoapi assumes OGC:CRS84 by default
curl http://localhost:5000/collections/beni_confiscati/items?\
    filter=DWITHIN(geometry,SRID=4326;POINT(12.626000006069875 41.8108000011414),100,meters)

returns 0 results, even if the recommended query parameter &filter-crs=https://www.opengis.net/def/crs/OGC/1.3/CRS84 is added to the request