Convert coordinates to storage crs when filtering via cql (#1489)

* Allow either URL or URN CRS URIs

* Implemented transformation of geometries in CQL filter

* Fixed flake8 issue

* Removed commented out code

* Implemented support for the filter-crs query parameter

* Remove unneeded print() calls

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

@@ -85,7 +85,8 @@
                            json_serial, render_j2_template, str2bool,
                            TEMPLATES, to_json, get_api_rules, get_base_url,
                            get_crs_from_uri, get_supported_crs_list,
-                           CrsTransformSpec, transform_bbox)
+                           modify_pygeofilter, CrsTransformSpec,
+                           transform_bbox)
 
 LOGGER = logging.getLogger(__name__)
 
@@ -1501,7 +1502,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')
@@ -1714,11 +1715,19 @@ 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'),
+                )
             except Exception as err:
                 LOGGER.error(err)
                 msg = f'Bad CQL string : {cql_text}'
@@ -1736,7 +1745,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)
 
@@ -1755,7 +1763,9 @@ def get_collection_items(
         LOGGER.debug(f'language: {prv_locale}')
         LOGGER.debug(f'q: {q}')
         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,
@@ -1925,7 +1935,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')
@@ -2012,19 +2022,21 @@ def post_collection_items(
         LOGGER.debug('Loading provider')
 
         try:
-            p = load_plugin('provider', get_provider_by_type(
-                collections[dataset]['providers'], 'feature'))
+            provider_def = get_provider_by_type(
+                collections[dataset]['providers'], 'feature')
         except ProviderTypeError:
             try:
-                p = load_plugin('provider', get_provider_by_type(
-                    collections[dataset]['providers'], 'record'))
+                provider_def = get_provider_by_type(
+                    collections[dataset]['providers'], 'record')
             except ProviderTypeError:
                 msg = 'Invalid provider type'
                 return self.get_exception(
                     HTTPStatus.BAD_REQUEST, headers, request.format,
                     'NoApplicableCode', msg)
+
+        try:
+            p = load_plugin('provider', provider_def)
         except ProviderGenericError as err:
-            LOGGER.error(err)
             return self.get_exception(
                 err.http_status_code, headers, request.format,
                 err.ogc_exception_code, err.message)
@@ -2086,6 +2098,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
@@ -2105,6 +2119,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')
 
@@ -2142,6 +2157,12 @@ def post_collection_items(
             LOGGER.debug('processing PostgreSQL CQL_JSON data')
             try:
                 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')
+                )
             except Exception as err:
                 LOGGER.error(err)
                 msg = f'Bad CQL string : {data}'

pygeoapi/provider/postgresql.py

@@ -55,7 +55,6 @@
 from geoalchemy2.functions import ST_MakeEnvelope
 from geoalchemy2.shape import to_shape
 from pygeofilter.backends.sqlalchemy.evaluate import to_filter
-import pygeofilter.ast
 import pyproj
 import shapely
 from sqlalchemy import create_engine, MetaData, PrimaryKeyConstraint, asc, desc
@@ -139,8 +138,7 @@ def query(self, offset=0, limit=10, resulttype='results',
 
         LOGGER.debug('Preparing filters')
         property_filters = self._get_property_filters(properties)
-        modified_filterq = self._modify_pygeofilter(filterq)
-        cql_filters = self._get_cql_filters(modified_filterq)
+        cql_filters = self._get_cql_filters(filterq)
         bbox_filter = self._get_bbox_filter(bbox)
         order_by_clauses = self._get_order_by_clauses(sortby, self.table_model)
         selected_properties = self._select_properties_clause(select_properties,
@@ -497,40 +495,3 @@ def _get_crs_transform(self, crs_transform_spec=None):
         else:
             crs_transform = None
         return crs_transform
-
-    def _modify_pygeofilter(
-            self,
-            ast_tree: pygeofilter.ast.Node,
-    ) -> pygeofilter.ast.Node:
-        """
-        Prepare the input pygeofilter for querying the database.
-
-        Returns a new ``pygeofilter.ast.Node`` object that can be used for
-        querying the database.
-        """
-        new_tree = deepcopy(ast_tree)
-        _inplace_replace_geometry_filter_name(new_tree, self.geom)
-        return new_tree
-
-
-def _inplace_replace_geometry_filter_name(
-        node: pygeofilter.ast.Node,
-        geometry_column_name: str
-):
-    """Recursively traverse node tree and rename nodes of type ``Attribute``.
-
-    Nodes of type ``Attribute`` named ``geometry`` are renamed to the value of
-    the ``geometry_column_name`` parameter.
-    """
-    try:
-        sub_nodes = node.get_sub_nodes()
-    except AttributeError:
-        pass
-    else:
-        for sub_node in sub_nodes:
-            is_attribute_node = isinstance(sub_node, pygeofilter.ast.Attribute)
-            if is_attribute_node and sub_node.name == "geometry":
-                sub_node.name = geometry_column_name
-            else:
-                _inplace_replace_geometry_filter_name(
-                    sub_node, geometry_column_name)