Implemented support for the filter-crs query parameter

docs/source/cql.rst

@@ -40,6 +40,15 @@ Using Elasticsearch the following type of queries are supported right now:
 * Logical ``and`` query with ``between`` and ``eq`` expression
 * Spatial query with ``bbox``
 
+Note that when using a spatial operator in your filter expression, geometries are by default interpreted as being
+in the OGC:CRS84 Coordinate Reference System. If you wish to provide geometries in other CRS, use the ``filter-crs``
+query parameter with a suitable value.
+
+Alternatively, a geometry's CRS may also be included using Extended Well-Known Text, in which case it will override
+the value of ``filter-crs`` (if any) - this can be useful if your filtering expression is complex enough to
+need multiple geometries expressed in different CRSs. The standard way of providing ``filter-crs`` as an additional
+query parameter is preferable for most cases.
+
 Examples
 ^^^^^^^^
 
@@ -93,4 +102,20 @@ A ``CROSSES`` example via an HTTP GET request.  The CQL text is passed via the `
 
   curl "http://localhost:5000/collections/hot_osm_waterways/items?f=json&filter=CROSSES(foo_geom,%20LINESTRING(28%20-2,%2030%20-4))"
 
+A ``DWITHIN`` example via HTTP GET and using a custom CRS for the filter geometry:
+
+.. code-block:: bash
+
+  curl "http://localhost:5000/collections/beni/items?filter=DWITHIN(geometry,POINT(1392921%205145517),100,meters)&filter-crs=http://www.opengis.net/def/crs/EPSG/0/3857"
+
+
+The same example, but this time providing a geometry in EWKT format:
+
+.. code-block:: bash
+
+  curl "http://localhost:5000/collections/beni/items?filter=DWITHIN(geometry,SRID=3857;POINT(1392921%205145517),100,meters)"
+
+
+
+
 Note that the CQL text has been URL encoded. This is required in curl commands but when entering in a browser, plain text can be used e.g. ``CROSSES(foo_geom, LINESTRING(28 -2, 30 -4))``.

pygeoapi/api.py

@@ -1376,7 +1376,7 @@ def get_collection_items(
         reserved_fieldnames = ['bbox', 'bbox-crs', 'crs', 'f', 'lang', 'limit',
                                'offset', 'resulttype', 'datetime', 'sortby',
                                'properties', 'skipGeometry', 'q',
-                               'filter', 'filter-lang']
+                               'filter', 'filter-lang', 'filter-crs']
 
         collections = filter_dict_by_key_value(self.config['resources'],
                                                'type', 'collection')
@@ -1589,15 +1589,18 @@ def get_collection_items(
         else:
             skip_geometry = False
 
+        LOGGER.debug('Processing filter-crs parameter')
+        filter_crs_uri = request.params.get('filter-crs', DEFAULT_CRS)
         LOGGER.debug('processing filter parameter')
         cql_text = request.params.get('filter')
         if cql_text is not None:
             try:
                 filter_ = parse_ecql_text(cql_text)
                 filter_ = modify_pygeofilter(
                     filter_,
+                    filter_crs_uri=filter_crs_uri,
                     storage_crs_uri=provider_def.get('storage_crs'),
-                    geometry_column_name=provider_def.get('geom_field')
+                    geometry_column_name=provider_def.get('geom_field'),
                 )
             except Exception as err:
                 LOGGER.error(err)
@@ -1616,7 +1619,6 @@ def get_collection_items(
             return self.get_exception(
                 HTTPStatus.BAD_REQUEST, headers, request.format,
                 'InvalidParameterValue', msg)
-
         # Get provider locale (if any)
         prv_locale = l10n.get_plugin_locale(provider_def, request.raw_locale)
 
@@ -1637,6 +1639,7 @@ def get_collection_items(
         LOGGER.debug(f'cql_text: {cql_text}')
         LOGGER.debug(f'filter_: {filter_}')
         LOGGER.debug(f'filter-lang: {filter_lang}')
+        LOGGER.debug(f'filter-crs: {filter_crs_uri}')
 
         try:
             content = p.query(offset=offset, limit=limit,
@@ -1806,7 +1809,7 @@ def post_collection_items(
         reserved_fieldnames = ['bbox', 'f', 'limit', 'offset',
                                'resulttype', 'datetime', 'sortby',
                                'properties', 'skipGeometry', 'q',
-                               'filter-lang']
+                               'filter-lang', 'filter-crs']
 
         collections = filter_dict_by_key_value(self.config['resources'],
                                                'type', 'collection')
@@ -1969,6 +1972,8 @@ def post_collection_items(
         else:
             skip_geometry = False
 
+        LOGGER.debug('Processing filter-crs parameter')
+        filter_crs = request.params.get('filter-crs', DEFAULT_CRS)
         LOGGER.debug('Processing filter-lang parameter')
         filter_lang = request.params.get('filter-lang')
         if filter_lang != 'cql-json':  # @TODO add check from the configuration
@@ -1988,6 +1993,7 @@ def post_collection_items(
         LOGGER.debug(f'skipGeometry: {skip_geometry}')
         LOGGER.debug(f'q: {q}')
         LOGGER.debug(f'filter-lang: {filter_lang}')
+        LOGGER.debug(f'filter-crs: {filter_crs}')
 
         LOGGER.debug('Processing headers')
 
@@ -2027,6 +2033,7 @@ def post_collection_items(
                 filter_ = parse_cql_json(data)
                 filter_ = modify_pygeofilter(
                     filter_,
+                    filter_crs_uri=filter_crs,
                     storage_crs_uri=provider_def.get('storage_crs'),
                     geometry_column_name=provider_def.get('geom_field')
                 )

pygeoapi/util.py

@@ -884,6 +884,7 @@ def bbox2geojsongeometry(bbox: list) -> dict:
 def modify_pygeofilter(
         ast_tree: pygeofilter.ast.Node,
         *,
+        filter_crs_uri: str,
         storage_crs_uri: Optional[str] = None,
         geometry_column_name: Optional[str] = None
 ) -> pygeofilter.ast.Node:
@@ -892,6 +893,8 @@ def modify_pygeofilter(
 
     :param ast_tree: `pygeofilter.ast.Node` representing the
     already parsed pygeofilter expression
+    :param filter_crs_uri: URI of the CRS being used in the filtering
+    expression
     :param storage_crs_uri: An optional string containing the URI of
     the provider's storage CRS
     :param geometry_column_name: An optional string containing the
@@ -915,14 +918,16 @@ def modify_pygeofilter(
     new_tree = deepcopy(ast_tree)
     if storage_crs_uri:
         storage_crs = get_crs_from_uri(storage_crs_uri)
-        _inplace_transform_filter_geometries(new_tree, storage_crs)
+        filter_crs = get_crs_from_uri(filter_crs_uri)
+        _inplace_transform_filter_geometries(new_tree, filter_crs, storage_crs)
     if geometry_column_name:
         _inplace_replace_geometry_filter_name(new_tree, geometry_column_name)
     return new_tree
 
 
 def _inplace_transform_filter_geometries(
         node: pygeofilter.ast.Node,
+        filter_crs: pyproj.CRS,
         storage_crs: pyproj.CRS
 ):
     """
@@ -941,33 +946,44 @@ def _inplace_transform_filter_geometries(
             is_geometry_node = isinstance(
                 sub_node, pygeofilter.values.Geometry)
             if is_geometry_node:
-                # NOTE: We specify a default CRS using a URI of type URN
+                # NOTE1: To be flexible, and since pygeofilter
+                # already supports it, in addition to supporting
+                # the `filter-crs` parameter, we also support having a
+                # geometry defined in EWKT, meaning the CRS is provided
+                # inline, like this `SRID=<CRS_CODE>;<WKT>` - If provided,
+                # this overrides the value of `filter-crs`. This enables
+                # supporting, for example, an exotic filter expression with
+                # multiple geometries specified in different CRSs
+
+                # NOTE2: We specify a default CRS using a URI of type URN
                 # because this is what pygeofilter uses internally too
-                crs = get_crs_from_uri(
-                    sub_node.geometry.get(
-                        'crs', {}
-                    ).get(
-                        'properties', {}
-                    ).get('name', 'urn:ogc:def:crs:OGC:1.3:CRS84')
-                )
+
+                crs_urn_provided_in_ewkt = sub_node.geometry.get(
+                    'crs', {}).get('properties', {}).get('name')
+                if crs_urn_provided_in_ewkt is not None:
+                    crs = get_crs_from_uri(crs_urn_provided_in_ewkt)
+                else:
+                    crs = filter_crs
+                print(f"{crs=}")
+                print(f"{storage_crs=}")
                 if crs != storage_crs:
                     # convert geometry coordinates to storage crs
                     geom = geojson_to_geom(sub_node.geometry)
                     coord_transformer = pyproj.Transformer.from_crs(
                         crs_from=crs, crs_to=storage_crs).transform
                     transformed_geom = ops.transform(coord_transformer, geom)
-                    authority, code = storage_crs.to_authority()
-                    sub_node.geometry = {
-                        **geom_to_geojson(transformed_geom),
-                        'crs': {
-                            'properties': {
-                                'name': f'urn:ogc:def:crs:{authority}::{code}'
-                            }
-                        }
+                    sub_node.geometry = geom_to_geojson(transformed_geom)
+                # ensure the crs is encoded in the sub-node, otherwise
+                # pygeofilter will assign it its own default CRS
+                authority, code = storage_crs.to_authority()
+                sub_node.geometry['crs'] = {
+                    'properties': {
+                        'name': f'urn:ogc:def:crs:{authority}::{code}'
                     }
+                }
             else:
                 _inplace_transform_filter_geometries(
-                    sub_node, storage_crs)
+                    sub_node, filter_crs, storage_crs)
 
 
 def _inplace_replace_geometry_filter_name(

tests/test_util.py

@@ -331,9 +331,10 @@ def test_prefetcher():
     assert headers.get('content-type') == 'image/png'
 
 
-@pytest.mark.parametrize('original_filter, storage_crs, geometry_colum_name, expected', [  # noqa
+@pytest.mark.parametrize('original_filter, filter_crs, storage_crs, geometry_colum_name, expected', [  # noqa
     pytest.param(
-        "INTERSECTS(geometry, POINT(1 1))",
+        'INTERSECTS(geometry, POINT(1 1))',
+        'http://www.opengis.net/def/crs/OGC/1.3/CRS84',
         None,
         None,
         pygeofilter.ast.GeometryIntersects(
@@ -344,6 +345,7 @@ def test_prefetcher():
     ),
     pytest.param(
         "INTERSECTS(geometry, POINT(1 1))",
+        'http://www.opengis.net/def/crs/OGC/1.3/CRS84',
         None,
         'custom_geom_name',
         pygeofilter.ast.GeometryIntersects(
@@ -354,6 +356,7 @@ def test_prefetcher():
     ),
     pytest.param(
         "some_attribute = 10 AND INTERSECTS(geometry, POINT(1 1))",
+        'http://www.opengis.net/def/crs/OGC/1.3/CRS84',
         None,
         'custom_geom_name',
         pygeofilter.ast.And(
@@ -369,6 +372,7 @@ def test_prefetcher():
     pytest.param(
         "(some_attribute = 10 AND INTERSECTS(geometry, POINT(1 1))) OR "
         "DWITHIN(geometry, POINT(2 2), 10, meters)",
+        'http://www.opengis.net/def/crs/OGC/1.3/CRS84',
         None,
         'custom_geom_name',
         pygeofilter.ast.Or(
@@ -391,6 +395,7 @@ def test_prefetcher():
     ),
     pytest.param(
         "INTERSECTS(geometry, POINT(12.512829 41.896698))",
+        'http://www.opengis.net/def/crs/OGC/1.3/CRS84',
         'http://www.opengis.net/def/crs/EPSG/0/3004',
         None,
         pygeofilter.ast.GeometryIntersects(
@@ -401,6 +406,7 @@ def test_prefetcher():
     ),
     pytest.param(
         "some_attribute = 10 AND INTERSECTS(geometry, POINT(12.512829 41.896698))",  # noqa
+        'http://www.opengis.net/def/crs/OGC/1.3/CRS84',
         'http://www.opengis.net/def/crs/EPSG/0/3004',
         None,
         pygeofilter.ast.And(
@@ -416,6 +422,7 @@ def test_prefetcher():
     pytest.param(
         "(some_attribute = 10 AND INTERSECTS(geometry, POINT(12.512829 41.896698))) OR "  # noqa
         "DWITHIN(geometry, POINT(12 41), 10, meters)",
+        'http://www.opengis.net/def/crs/OGC/1.3/CRS84',
         'http://www.opengis.net/def/crs/EPSG/0/3004',
         None,
         pygeofilter.ast.Or(
@@ -438,16 +445,40 @@ def test_prefetcher():
     ),
     pytest.param(
         "INTERSECTS(geometry, SRID=3857;POINT(1392921 5145517))",
+        'http://www.opengis.net/def/crs/OGC/1.3/CRS84',
+        'http://www.opengis.net/def/crs/EPSG/0/3004',
+        None,
+        pygeofilter.ast.GeometryIntersects(
+            pygeofilter.ast.Attribute(name='geometry'),
+            Geometry({'type': 'Point', 'coordinates': (2313681.8086284213, 4641307.939955416)})  # noqa
+        ),
+        id='unnested-geometry-transformed-coords-explicit-input-crs-ewkt'
+    ),
+    pytest.param(
+        "INTERSECTS(geometry, POINT(1392921 5145517))",
+        'http://www.opengis.net/def/crs/EPSG/0/3857',
+        'http://www.opengis.net/def/crs/EPSG/0/3004',
+        None,
+        pygeofilter.ast.GeometryIntersects(
+            pygeofilter.ast.Attribute(name='geometry'),
+            Geometry({'type': 'Point', 'coordinates': (2313681.8086284213, 4641307.939955416)})  # noqa
+        ),
+        id='unnested-geometry-transformed-coords-explicit-input-crs-filter-crs'
+    ),
+    pytest.param(
+        "INTERSECTS(geometry, SRID=3857;POINT(1392921 5145517))",
+        'http://www.opengis.net/def/crs/OGC/1.3/CRS84',
         'http://www.opengis.net/def/crs/EPSG/0/3004',
         None,
         pygeofilter.ast.GeometryIntersects(
             pygeofilter.ast.Attribute(name='geometry'),
             Geometry({'type': 'Point', 'coordinates': (2313681.8086284213, 4641307.939955416)})  # noqa
         ),
-        id='unnested-geometry-transformed-coords-explicit-input-crs'
+        id='unnested-geometry-transformed-coords-ewkt-crs-overrides-filter-crs'
     ),
     pytest.param(
         "INTERSECTS(geometry, POINT(12.512829 41.896698))",
+        'http://www.opengis.net/def/crs/OGC/1.3/CRS84',
         'http://www.opengis.net/def/crs/EPSG/0/3004',
         'custom_geom_name',
         pygeofilter.ast.GeometryIntersects(
@@ -458,10 +489,16 @@ def test_prefetcher():
     ),
 ])
 def test_modify_pygeofilter(
-        original_filter, storage_crs, geometry_colum_name, expected):
+        original_filter,
+        filter_crs,
+        storage_crs,
+        geometry_colum_name,
+        expected
+):
     parsed_filter = parse(original_filter)
     result = util.modify_pygeofilter(
         parsed_filter,
+        filter_crs_uri=filter_crs,
         storage_crs_uri=storage_crs,
         geometry_column_name=geometry_colum_name
     )