Coverage for pygeoapi/crs.py: 48%
157 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-10-07 08:15 +0000
1# =================================================================
2#
3# Authors: Tom Kralidis <tomkralidis@gmail.com>
4# Just van den Broecke <justb4@gmail.com>
5#
6# Copyright (c) 2025 Tom Kralidis
7# Copyright (c) 2025 Just van den Broecke
8#
9# Permission is hereby granted, free of charge, to any person
10# obtaining a copy of this software and associated documentation
11# files (the "Software"), to deal in the Software without
12# restriction, including without limitation the rights to use,
13# copy, modify, merge, publish, distribute, sublicense, and/or sell
14# copies of the Software, and to permit persons to whom the
15# Software is furnished to do so, subject to the following
16# conditions:
17#
18# The above copyright notice and this permission notice shall be
19# included in all copies or substantial portions of the Software.
20#
21# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
22# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
23# OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
24# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
25# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
26# WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
28# OTHER DEALINGS IN THE SOFTWARE.
29#
30# =================================================================
32"""Generic CRS functions used in the code"""
34from copy import deepcopy
35import functools
36from functools import partial
37from dataclasses import dataclass
38import logging
39from typing import Union, Optional, Callable
41import pyproj
42import pygeofilter.ast
43import pygeofilter.values
44from pyproj.exceptions import CRSError
45from shapely import ops, Geometry
46from shapely.geometry import (
47 shape as geojson_to_geom,
48 mapping as geom_to_geojson
49)
51LOGGER = logging.getLogger(__name__)
53DEFAULT_CRS_LIST = [
54 'http://www.opengis.net/def/crs/OGC/1.3/CRS84',
55 'http://www.opengis.net/def/crs/OGC/1.3/CRS84h',
56]
58DEFAULT_CRS = 'http://www.opengis.net/def/crs/OGC/1.3/CRS84'
59DEFAULT_STORAGE_CRS = DEFAULT_CRS
62@dataclass
63class CrsTransformSpec:
64 source_crs_uri: str
65 source_crs_wkt: str
66 target_crs_uri: str
67 target_crs_wkt: str
68 always_xy: bool = False
71def get_srid(crs: Union[str, pyproj.CRS]) -> Union[int, None]:
72 """
73 Helper function to attempt to extract an EPSG SRID from
74 a `pyproj.CRS` object.
76 :param crs: `pyproj.CRS` object
78 :returns: int of EPSG SRID, if found
79 """
80 if isinstance(crs, str):
81 crs = get_crs(crs)
83 if crs.to_epsg():
84 return crs.to_epsg()
86 try:
87 return pyproj.CRS(crs.to_proj4()).to_epsg()
88 except KeyError:
89 LOGGER.debug('Unable to extract SRID from proj4 string')
92def get_supported_crs_list(
93 provider_def: dict, default_crs_list: list = DEFAULT_CRS_LIST
94) -> list:
95 """
96 Helper function to get a complete list of supported CRSs
97 from a (Provider) config dict. Result should always include
98 a default CRS according to OAPIF Part 2 OGC Standard.
99 This will be the default when no CRS list in config or
100 added when (partially) missing in config.
102 Author: @justb4
104 :param provider_def: dictionary with or without a list of CRSs
105 :param default_crs_list: default CRS alternatives, first is default
107 :returns: list of supported CRSs
108 """
109 supported_crs_list = provider_def.get('crs', list())
110 contains_default = False
111 for uri in supported_crs_list: 111 ↛ 112line 111 didn't jump to line 112 because the loop on line 111 never started
112 if uri in default_crs_list:
113 contains_default = True
114 break
116 # A default CRS is missing: add the first which is the default
117 if not contains_default: 117 ↛ 120line 117 didn't jump to line 120 because the condition on line 117 was always true
118 supported_crs_list.append(default_crs_list[0])
120 return supported_crs_list
123def get_crs(crs: Union[str, pyproj.CRS]) -> pyproj.CRS:
124 """
125 Get a `pyproj.CRS` instance from a CRS.
126 Author: @MTachon
128 :param crs: Uniform resource identifier of the coordinate
129 reference system. In accordance with
130 https://docs.ogc.org/pol/09-048r5.html#_naming_rule
131 URIs can take either the form of a URL or a URN
132 or `pyproj.CRS` object
133 :raises `CRSError`: Error raised if no CRS could be identified from the
134 URI.
136 :returns: `pyproj.CRS` instance matching the input URI.
137 """
139 if isinstance(crs, pyproj.CRS): 139 ↛ 140line 139 didn't jump to line 140 because the condition on line 139 was never true
140 return crs
142 # normalize the input `uri` to a URL first
143 uri = str(crs)
144 url = uri.replace(
145 'urn:ogc:def:crs', 'http://www.opengis.net/def/crs'
146 ).replace(':', '/')
147 try:
148 authority, code = url.rsplit('/', maxsplit=3)[1::2]
149 crs = pyproj.CRS.from_authority(authority, code)
150 except ValueError:
151 msg = (
152 f'CRS could not be identified from URI {uri!r}. CRS URIs must '
153 'follow one of two formats: '
154 '"http://www.opengis.net/def/crs/{authority}/{version}/{code}" or '
155 '"urn:ogc:def:crs:{authority}:{version}:{code}" '
156 '(see https://docs.opengeospatial.org/is/18-058r1/18-058r1.html#crs-overview).' # noqa
157 )
158 LOGGER.error(msg)
159 raise CRSError(msg)
160 except CRSError:
161 msg = f"CRS could not be identified from URI {uri!r}"
162 LOGGER.error(msg)
163 raise CRSError(msg)
165 return crs
168def get_transform_from_spec(
169 crs_transform_spec: CrsTransformSpec
170) -> Callable[[Geometry], Geometry]:
171 """ Get transformation function from a `CrsTransformSpec` instance.
173 :param crs_transform_spec: `CrsTransformSpec`
175 :returns: `callable` Function to transform the coordinates of a `Geometry`.
176 """
177 if crs_transform_spec is None:
178 LOGGER.warning('No transform spec found')
179 return None
181 return get_transform_from_crs(
182 pyproj.CRS.from_wkt(crs_transform_spec.source_crs_wkt),
183 pyproj.CRS.from_wkt(crs_transform_spec.target_crs_wkt),
184 crs_transform_spec.always_xy
185 )
188def get_transform_from_crs(
189 crs_in: pyproj.CRS, crs_out: pyproj.CRS, always_xy: bool = False
190) -> Callable[[Geometry], Geometry]:
191 """ Get transformation function from two `pyproj.CRS` instances.
193 Get function to transform the coordinates of a Shapely geometrical object
194 from one coordinate reference system to another.
196 :param crs_in: `pyproj.CRS` Input Coordinate Reference System
197 :param crs_out: `pyproj.CRS` Output Coordinate Reference System
198 :param always_xy: 'bool' should axis order be forced to x,y (lon, lat)
199 even if CRSdeclares y,x (lat,lon)
201 :returns: `callable` Function to transform the coordinates of a `Geometry`.
202 """
203 crs_transform = pyproj.Transformer.from_crs(
204 crs_in, crs_out, always_xy=always_xy,
205 ).transform
206 return partial(ops.transform, crs_transform)
209def crs_transform(func):
210 """Decorator that transforms the geometry's/geometries' coordinates of a
211 Feature/FeatureCollection.
213 This function can be used to decorate another function which returns either
214 a Feature or a FeatureCollection (GeoJSON-like `dict`). For a
215 FeatureCollection, the Features are stored in a ´list´ available at the
216 'features' key of the returned `dict`. For each Feature, the geometry is
217 available at the 'geometry' key. The decorated function may take a
218 'crs_transform_spec' parameter, which accepts a `CrsTransformSpec` instance
219 as value. If the `CrsTransformSpec` instance represents a coordinates
220 transformation between two different CRSs, the coordinates of the
221 Feature's/FeatureCollection's geometry/geometries will be transformed
222 before returning the Feature/FeatureCollection. If the 'crs_transform_spec'
223 parameter is not given, passed `None` or passed a `CrsTransformSpec`
224 instance which does not represent a coordinates transformation, the
225 Feature/FeatureCollection is returned unchanged. This decorator can for
226 example be use to help supporting coordinates transformation of
227 Feature/FeatureCollection `dict` objects returned by the `get` and `query`
228 methods of (new or with no native support for transformations) providers of
229 type 'feature'.
231 :param func: Function to decorate.
233 :returns: Decorated function.
234 """
235 @functools.wraps(func)
236 def get_geojsonf(*args, **kwargs):
237 crs_transform_spec = kwargs.get('crs_transform_spec')
238 result = func(*args, **kwargs)
239 if crs_transform_spec is None: 239 ↛ 246line 239 didn't jump to line 246 because the condition on line 239 was always true
240 # No coordinates transformation for feature(s) returned by the
241 # decorated function.
242 LOGGER.debug('crs_transform: NOT applying coordinate transforms')
243 return result
244 # Create transformation function and transform the output feature(s)'
245 # coordinates before returning them.
246 transform_func = get_transform_from_spec(crs_transform_spec)
247 LOGGER.debug(f'crs_transform: transforming features CRS '
248 f'from {crs_transform_spec.source_crs_uri} '
249 f'to {crs_transform_spec.target_crs_uri}')
251 features = result.get('features')
252 # Decorated function returns a single Feature
253 if features is None:
254 # Transform the feature's coordinates
255 crs_transform_feature(result, transform_func)
256 # Decorated function returns a FeatureCollection
257 else:
258 # Transform all features' coordinates
259 for feature in features:
260 crs_transform_feature(feature, transform_func)
261 return result
262 return get_geojsonf
265def crs_transform_feature(feature: dict, transform_func: Callable):
266 """Transform the coordinates of a Feature.
268 :param feature: Feature (GeoJSON-like `dict`) to transform.
269 :param transform_func: Function that transforms the coordinates of a
270 `Geometry` instance.
272 :returns: None
273 """
274 json_geometry = feature.get('geometry')
275 if json_geometry is not None:
276 feature['geometry'] = geom_to_geojson(
277 transform_func(geojson_to_geom(json_geometry))
278 )
281def transform_bbox(bbox: list, from_crs: Union[str, pyproj.CRS],
282 to_crs: Union[str, pyproj.CRS]) -> list:
283 """
284 helper function to transform a bounding box (bbox) from
285 a source to a target CRS. CRSs in URI str format.
286 Uses pyproj Transformer.
288 :param bbox: list of coordinates in 'from_crs' projection
289 :param from_crs: CRS to transform from
290 :param to_crs: CRS to transform to
291 :raises `CRSError`: Error raised if no CRS could be identified from an
292 URI.
294 :returns: list of 4 or 6 coordinates
295 """
297 from_crs_obj = get_crs(from_crs)
298 to_crs_obj = get_crs(to_crs)
299 transform_func = pyproj.Transformer.from_crs(
300 from_crs_obj, to_crs_obj).transform
302 n_dims = len(bbox) // 2
303 return list(transform_func(*bbox[:n_dims]) + transform_func(
304 *bbox[n_dims:]))
307def modify_pygeofilter(
308 ast_tree: pygeofilter.ast.Node,
309 *,
310 filter_crs_uri: str,
311 storage_crs_uri: Optional[str] = None,
312 geometry_column_name: Optional[str] = None
313) -> pygeofilter.ast.Node:
314 """
315 Modifies the input pygeofilter with information from the provider.
317 :param ast_tree: `pygeofilter.ast.Node` representing the
318 already parsed pygeofilter expression
319 :param filter_crs_uri: URI of the CRS being used in the filtering
320 expression
321 :param storage_crs_uri: An optional string containing the URI of
322 the provider's storage CRS
323 :param geometry_column_name: An optional string containing the
324 actual name of the provider's geometry field
325 :returns: A new pygeofilter.ast.Node, with the modified filter
326 expression
328 This function modifies the parsed pygeofilter that contains the raw
329 filter expression provided by an external client. It performs the
330 following modifications:
332 - if the filter includes any spatial coordinates and they are being
333 provided in a different CRS from the provider's storage CRS, the
334 corresponding geometries are transformed into the storage CRS
336 - if the filter includes the generic 'geometry' name as a reference to
337 the actual geometry of features, it is replaced by the actual name
338 of the geometry field, as specified by the provider
340 """
341 new_tree = deepcopy(ast_tree)
342 if storage_crs_uri: 342 ↛ 343line 342 didn't jump to line 343 because the condition on line 342 was never true
343 _inplace_transform_filter_geometries(
344 new_tree, get_crs(filter_crs_uri), get_crs(storage_crs_uri)
345 )
346 if geometry_column_name: 346 ↛ 347line 346 didn't jump to line 347 because the condition on line 346 was never true
347 _inplace_replace_geometry_filter_name(new_tree, geometry_column_name)
348 return new_tree
351def _inplace_transform_filter_geometries(
352 node: pygeofilter.ast.Node,
353 filter_crs: pyproj.CRS,
354 storage_crs: pyproj.CRS
355) -> None:
356 """
357 Recursively traverse node tree and convert coordinates to the storage CRS.
359 This function modifies nodes in the already-parsed filter in order to find
360 any geometry literals that may be used in the filter and, if necessary,
361 proceeds to convert spatial coordinates to the CRS used by the provider.
362 """
363 try:
364 sub_nodes = node.get_sub_nodes()
365 except AttributeError:
366 pass
367 else:
368 for sub_node in sub_nodes:
369 is_geometry_node = isinstance(
370 sub_node, pygeofilter.values.Geometry)
371 if is_geometry_node:
372 # NOTE1: To be flexible, and since pygeofilter
373 # already supports it, in addition to supporting
374 # the `filter-crs` parameter, we also support having a
375 # geometry defined in EWKT, meaning the CRS is provided
376 # inline, like this `SRID=<CRS_CODE>;<WKT>` - If provided,
377 # this overrides the value of `filter-crs`. This enables
378 # supporting, for example, an exotic filter expression with
379 # multiple geometries specified in different CRSs
381 # NOTE2: We specify a default CRS using a URI of type URN
382 # because this is what pygeofilter uses internally too
384 crs_urn_provided_in_ewkt = sub_node.geometry.get(
385 'crs', {}).get('properties', {}).get('name')
386 if crs_urn_provided_in_ewkt is not None:
387 crs = get_crs(crs_urn_provided_in_ewkt)
388 else:
389 crs = filter_crs
390 if crs != storage_crs:
391 # convert geometry coordinates to storage crs
392 geom = geojson_to_geom(sub_node.geometry)
393 coord_transformer = pyproj.Transformer.from_crs(
394 crs_from=crs, crs_to=storage_crs).transform
395 transformed_geom = ops.transform(coord_transformer, geom)
396 sub_node.geometry = geom_to_geojson(transformed_geom)
397 # ensure the crs is encoded in the sub-node, otherwise
398 # pygeofilter will assign it its own default CRS
399 authority, code = storage_crs.to_authority()
400 sub_node.geometry['crs'] = {
401 'properties': {
402 'name': f'urn:ogc:def:crs:{authority}::{code}'
403 }
404 }
405 else:
406 _inplace_transform_filter_geometries(
407 sub_node, filter_crs, storage_crs)
410def _inplace_replace_geometry_filter_name(
411 node: pygeofilter.ast.Node,
412 geometry_column_name: str
413) -> None:
414 """Recursively traverse node tree and rename nodes of type ``Attribute``.
416 Nodes of type ``Attribute`` named ``geometry`` are renamed to the value of
417 the ``geometry_column_name`` parameter.
418 """
419 try:
420 sub_nodes = node.get_sub_nodes()
421 except AttributeError:
422 pass
423 else:
424 for sub_node in sub_nodes:
425 is_attribute_node = isinstance(sub_node, pygeofilter.ast.Attribute)
426 if is_attribute_node and sub_node.name == "geometry":
427 sub_node.name = geometry_column_name
428 else:
429 _inplace_replace_geometry_filter_name(
430 sub_node, geometry_column_name)
433def create_crs_transform_spec(
434 provider_def: dict, query_crs_uri: Optional[str] = None
435) -> Union[None, CrsTransformSpec]:
436 """
437 Create a `CrsTransformSpec` instance based on provider config and
438 *crs* query parameter.
440 :param provider_def: Provider config dictionary.
441 :param query_crs_uri: Uniform resource identifier of the coordinate
442 reference system (CRS) specified in query parameter (if specified).
444 :raises ValueError: Error raised if the CRS specified in the query
445 parameter is not in the list of supported CRSs of the provider.
446 :raises `CRSError`: Error raised if no CRS could be identified from the
447 query *crs* parameter (URI).
449 :returns: `CrsTransformSpec` instance if the CRS specified in query
450 parameter differs from the storage CRS, else `None`.
451 """
453 # Get storage/default CRS for Collection.
454 always_xy = provider_def.get('always_xy', False)
455 storage_crs_uri = provider_def.get('storage_crs', DEFAULT_STORAGE_CRS)
456 storage_crs = get_crs(storage_crs_uri)
458 if not query_crs_uri:
459 if storage_crs_uri in DEFAULT_CRS_LIST: 459 ↛ 464line 459 didn't jump to line 464 because the condition on line 459 was always true
460 # Could be that storage_crs is
461 # http://www.opengis.net/def/crs/OGC/1.3/CRS84h
462 query_crs_uri = storage_crs_uri
463 else:
464 query_crs_uri = DEFAULT_CRS
465 LOGGER.debug(f'no crs parameter, using default: {query_crs_uri}')
467 supported_crs_list = get_supported_crs_list(provider_def)
468 # Check that the crs specified by the query parameter is supported.
469 if query_crs_uri not in supported_crs_list: 469 ↛ 470line 469 didn't jump to line 470 because the condition on line 469 was never true
470 raise ValueError(
471 f'CRS {query_crs_uri!r} not supported for this '
472 'collection. List of supported CRSs: '
473 f'{", ".join(supported_crs_list)}.'
474 )
475 crs_out = get_crs(query_crs_uri)
477 # Check if the crs specified in query parameter differs from the
478 # storage crs.
479 if storage_crs == crs_out: 479 ↛ 483line 479 didn't jump to line 483 because the condition on line 479 was always true
480 LOGGER.debug('No CRS transformation')
481 return None
483 LOGGER.debug(
484 f'CRS transformation: {storage_crs} -> {crs_out}'
485 )
486 return CrsTransformSpec(
487 source_crs_uri=storage_crs_uri,
488 source_crs_wkt=storage_crs.to_wkt(),
489 target_crs_uri=query_crs_uri,
490 target_crs_wkt=crs_out.to_wkt(),
491 always_xy=always_xy
492 )
495def set_content_crs_header(
496 headers: dict, config: dict, query_crs_uri: Optional[str] = None,
497) -> None:
498 """Set the *Content-Crs* header in responses from providers of Feature
499 type.
501 :param headers: Response headers dictionary.
502 :param config: Provider config dictionary.
503 :param query_crs_uri: Uniform resource identifier of the coordinate
504 reference system specified in query parameter (if specified).
506 :returns: None
507 """
509 if query_crs_uri:
510 content_crs_uri = query_crs_uri
511 else:
512 # If empty use default CRS
513 storage_crs_uri = config.get('storage_crs', DEFAULT_STORAGE_CRS)
514 if storage_crs_uri in DEFAULT_CRS_LIST: 514 ↛ 517line 514 didn't jump to line 517 because the condition on line 514 was always true
515 content_crs_uri = storage_crs_uri
516 else:
517 content_crs_uri = DEFAULT_CRS
519 headers['Content-Crs'] = f'<{content_crs_uri}>'