Coverage for gws-app / gws / lib / crs / __init__.py: 81%
286 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-03 10:12 +0100
1from typing import Optional
3import math
4import re
5import warnings
7import pyproj.crs
8import pyproj.exceptions
9import pyproj.transformer
11import gws
14##
17class Object(gws.Crs):
18 def __init__(self, **kwargs):
19 vars(self).update(kwargs)
21 # crs objects with the same srid must be equal
22 # (despite caching, they can be different due to pickling)
24 def __hash__(self):
25 return self.srid
27 def __eq__(self, other):
28 return isinstance(other, Object) and other.srid == self.srid
30 def __repr__(self):
31 return f'<crs:{self.srid}>'
33 def axis_for_format(self, fmt):
34 if not self.isYX:
35 return self.axis
36 return _AXIS_FOR_FORMAT.get(fmt, self.axis)
38 def transform_extent(self, ext, crs_to):
39 if crs_to == self:
40 return ext
41 return _transform_extent_constrained(ext, self.srid, crs_to.srid)
43 def transformer(self, crs_to):
44 tr = _pyproj_transformer(self.srid, crs_to.srid)
45 return tr.transform
47 def extent_size_in_meters(self, extent):
48 x0, y0, x1, y1 = extent
50 if self.isProjected:
51 if self.uom != gws.Uom.m:
52 # @TODO support non-meter crs
53 raise Error(f'unsupported unit: {self.uom}')
54 return abs(x1 - x0), abs(y1 - y0)
56 geod = pyproj.Geod(ellps='WGS84')
58 mid_lat = (y0 + y1) / 2
59 _, _, w = geod.inv(x0, mid_lat, x1, mid_lat)
60 mid_lon = (x0 + x1) / 2
61 _, _, h = geod.inv(mid_lon, y0, mid_lon, y1)
63 return w, h
65 def point_offset_in_meters(self, xy, dist, az):
66 x, y = xy
68 if self.isProjected:
69 if self.uom != gws.Uom.m:
70 # @TODO support non-meter crs
71 raise Error(f'unsupported unit: {self.uom}')
73 if az == 0:
74 return x, y + dist
75 if az == 90:
76 return x + dist, y
77 if az == 180:
78 return x, y - dist
79 if az == 270:
80 return x - dist, y
82 az_rad = math.radians(90 - az)
83 return (
84 x + dist * math.cos(az_rad),
85 y + dist * math.sin(az_rad),
86 )
88 geod = pyproj.Geod(ellps='WGS84')
89 x, y, _ = geod.fwd(x, y, dist=dist, az=az)
90 return x, y
92 def to_string(self, fmt=None):
93 fmt = fmt or gws.CrsFormat.epsg
94 if fmt == gws.CrsFormat.srid:
95 return str(self.srid)
96 return getattr(self, str(fmt).lower())
98 def to_geojson(self):
99 # https://geojson.org/geojson-spec#named-crs
100 return {
101 'type': 'name',
102 'properties': {
103 'name': self.urn,
104 },
105 }
108#
111def qgis_extent_width(extent: gws.Extent) -> float:
112 # straight port from QGIS/src/core/qgsscalecalculator.cpp QgsScaleCalculator::calculateGeographicDistance
113 x0, y0, x1, y1 = extent
115 lat = (y0 + y1) * 0.5
116 RADS = (4.0 * math.atan(1.0)) / 180.0
117 a = math.pow(math.cos(lat * RADS), 2)
118 c = 2.0 * math.atan2(math.sqrt(a), math.sqrt(1.0 - a))
119 RA = 6378000
120 E = 0.0810820288
121 radius = RA * (1.0 - E * E) / math.pow(1.0 - E * E * math.sin(lat * RADS) * math.sin(lat * RADS), 1.5)
122 return (x1 - x0) / 180.0 * radius * c
125#
127# enough precision to represent 1cm
128COORDINATE_PRECISION_DEG = 7
129COORDINATE_PRECISION_M = 2
131WGS84: gws.Crs = Object(
132 srid=4326,
133 proj4text='+proj=longlat +datum=WGS84 +no_defs +type=crs',
134 wkt='GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]',
135 axis=gws.Axis.yx,
136 uom=gws.Uom.deg,
137 isGeographic=True,
138 isProjected=False,
139 isYX=True,
140 epsg='EPSG:4326',
141 urn='urn:ogc:def:crs:EPSG::4326',
142 urnx='urn:x-ogc:def:crs:EPSG:4326',
143 url='http://www.opengis.net/gml/srs/epsg.xml#4326',
144 uri='http://www.opengis.net/def/crs/epsg/0/4326',
145 name='WGS 84',
146 base=0,
147 datum='World Geodetic System 1984 ensemble',
148 wgsExtent=(-180, -90, 180, 90),
149 extent=(-180, -90, 180, 90),
150 coordinatePrecision=COORDINATE_PRECISION_DEG,
151)
153WGS84.bounds = gws.Bounds(crs=WGS84, extent=WGS84.extent)
155WEBMERCATOR: gws.Crs = Object(
156 srid=3857,
157 proj4text='+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs',
158 wkt='PROJCRS["WGS 84 / Pseudo-Mercator",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["Popular Visualisation Pseudo-Mercator",METHOD["Popular Visualisation Pseudo Mercator",ID["EPSG",1024]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Web mapping and visualisation."],AREA["World between 85.06°S and 85.06°N."],BBOX[-85.06,-180,85.06,180]],ID["EPSG",3857]]',
159 axis=gws.Axis.xy,
160 uom=gws.Uom.m,
161 isGeographic=False,
162 isProjected=True,
163 isYX=False,
164 epsg='EPSG:3857',
165 urn='urn:ogc:def:crs:EPSG::3857',
166 urnx='urn:x-ogc:def:crs:EPSG:3857',
167 url='http://www.opengis.net/gml/srs/epsg.xml#3857',
168 uri='http://www.opengis.net/def/crs/epsg/0/3857',
169 name='WGS 84 / Pseudo-Mercator',
170 base=4326,
171 datum='World Geodetic System 1984 ensemble',
172 wgsExtent=(-180, -85.06, 180, 85.06),
173 extent=(
174 -20037508.342789244,
175 -20048966.104014598,
176 20037508.342789244,
177 20048966.104014598,
178 ),
179 coordinatePrecision=COORDINATE_PRECISION_M,
180)
182WEBMERCATOR.bounds = gws.Bounds(crs=WEBMERCATOR, extent=WEBMERCATOR.extent)
184WEBMERCATOR_RADIUS = 6378137
185WEBMERCATOR_SQUARE = (
186 -math.pi * WEBMERCATOR_RADIUS,
187 -math.pi * WEBMERCATOR_RADIUS,
188 +math.pi * WEBMERCATOR_RADIUS,
189 +math.pi * WEBMERCATOR_RADIUS,
190)
193class Error(gws.Error):
194 pass
197def get(crs_name: Optional[gws.CrsName]) -> Optional[gws.Crs]:
198 """Returns the CRS for a given CRS-code or SRID."""
199 if not crs_name:
200 return None
201 return _get_crs(crs_name)
204def parse(crs_name: gws.CrsName) -> tuple[gws.CrsFormat, Optional[gws.Crs]]:
205 """Parses a CRS to a tuple of CRS-format and the CRS itself."""
206 fmt, srid = _parse(crs_name)
207 if not fmt:
208 return gws.CrsFormat.none, None
209 return fmt, _get_crs(srid)
212def require(crs_name: gws.CrsName) -> gws.Crs:
213 """Raises an error if no correct CRS is given."""
214 crs = _get_crs(crs_name)
215 if not crs:
216 raise Error(f'invalid CRS {crs_name!r}')
217 return crs
220##
223def best_match(crs: gws.Crs, supported_crs: list[gws.Crs]) -> gws.Crs:
224 """Return a crs from the list that most closely matches the given crs.
226 Args:
227 crs: target CRS
228 supported_crs: CRS list
230 Returns:
231 A CRS object
232 """
234 if crs in supported_crs:
235 return crs
237 bst = _best_match(crs, supported_crs)
238 if not bst:
239 bst = supported_crs[0] if supported_crs else crs
240 gws.log.debug(f'CRS: best_crs: using {bst.srid!r} for {crs.srid!r}')
241 return bst
244def _best_match(crs, supported_crs):
245 # @TODO find a projection with less errors
246 # @TODO find a projection with same units
248 if crs.isProjected:
249 # for a projected crs, find webmercator
250 for sup in supported_crs:
251 if sup.srid == WEBMERCATOR.srid:
252 return sup
254 # not found, return the first projected crs
255 for sup in supported_crs:
256 if sup.isProjected:
257 return sup
259 if crs.isGeographic:
260 # for a geographic crs, try wgs first
261 for sup in supported_crs:
262 if sup.srid == WGS84.srid:
263 return sup
265 # not found, return the first geographic crs
266 for sup in supported_crs:
267 if sup.isGeographic:
268 return sup
271##
274def _get_crs(crs_name):
275 if crs_name in _obj_cache:
276 return _obj_cache[crs_name]
278 fmt, srid = _parse(crs_name)
279 if not fmt:
280 gws.log.warning(f'CRS: cannot parse {crs_name!r}')
281 _obj_cache[crs_name] = None
282 return None
284 if srid in _obj_cache:
285 _obj_cache[crs_name] = _obj_cache[srid]
286 return _obj_cache[srid]
288 obj = _get_new_crs(srid)
289 _obj_cache[crs_name] = _obj_cache[srid] = obj
290 return obj
293def _get_new_crs(srid):
294 pp = _pyproj_crs_object(srid)
295 if not pp:
296 gws.log.warning(f'CRS: unknown srid {srid!r}')
297 return None
299 au = _axis_and_unit(pp)
300 if not au:
301 gws.log.warning(f'CRS: unsupported srid {srid!r}')
302 return None
304 axis, uom = au
305 if uom not in (gws.Uom.m, gws.Uom.deg):
306 gws.log.warning(f'CRS: unsupported unit {uom!r} for {srid!r}')
307 return None
309 return _make_crs(srid, pp, axis, uom)
312def _pyproj_crs_object(srid) -> Optional[pyproj.CRS]:
313 if srid in _pyproj_cache:
314 return _pyproj_cache[srid]
316 try:
317 pp = pyproj.CRS.from_epsg(srid)
318 except pyproj.exceptions.CRSError:
319 return None
321 _pyproj_cache[srid] = pp
322 return _pyproj_cache[srid]
325def _pyproj_transformer(srid_from, srid_to) -> pyproj.transformer.Transformer:
326 key = srid_from, srid_to
328 if key in _transformer_cache:
329 return _transformer_cache[key]
331 pa = _pyproj_crs_object(srid_from)
332 pb = _pyproj_crs_object(srid_to)
334 _transformer_cache[key] = pyproj.transformer.Transformer.from_crs(pa, pb, always_xy=True)
335 return _transformer_cache[key]
338def _transform_extent_constrained(ext, srid_from, srid_to):
339 ext_nor = _normalize_extent(ext)
341 if srid_from == WGS84.srid:
342 ext_wgs = ext_nor
343 else:
344 tr_to_wgs = _pyproj_transformer(srid_from, WGS84.srid)
345 ext_wgs = tr_to_wgs.transform_bounds(
346 left=ext_nor[0],
347 bottom=ext_nor[1],
348 right=ext_nor[2],
349 top=ext_nor[3],
350 errcheck=True,
351 )
353 if srid_to == WGS84.srid:
354 return _normalize_extent(ext_wgs)
356 pp = _pyproj_crs_object(srid_to)
357 if not pp:
358 raise Error(f'_transform_extent: unknown {srid_to=}')
360 au = pp.area_of_use
361 if not au:
362 ext_constrained = ext_wgs
363 else:
364 ext_use = au.bounds
365 x0 = max(ext_wgs[0], ext_use[0])
366 y0 = max(ext_wgs[1], ext_use[1])
367 x1 = min(ext_wgs[2], ext_use[2])
368 y1 = min(ext_wgs[3], ext_use[3])
369 if x0 >= x1 or y0 >= y1:
370 raise Error(f'transform_extent: no overlap {ext=} {srid_from!r}->{srid_to!r}')
371 ext_constrained = (x0, y0, x1, y1)
373 tr_from_wgs = _pyproj_transformer(WGS84.srid, srid_to)
374 ext_to = tr_from_wgs.transform_bounds(
375 left=ext_constrained[0],
376 bottom=ext_constrained[1],
377 right=ext_constrained[2],
378 top=ext_constrained[3],
379 errcheck=True,
380 )
381 return _normalize_extent(ext_to)
384def _transform_extent_direct(ext, srid_from, srid_to):
385 tr = _pyproj_transformer(srid_from, srid_to)
387 ext_nor = _normalize_extent(ext)
389 res = tr.transform_bounds(
390 left=ext_nor[0],
391 bottom=ext_nor[1],
392 right=ext_nor[2],
393 top=ext_nor[3],
394 errcheck=True,
395 )
396 return _normalize_extent(res)
399def _normalize_extent(ext):
400 return (
401 min(ext[0], ext[2]),
402 min(ext[1], ext[3]),
403 max(ext[0], ext[2]),
404 max(ext[1], ext[3]),
405 )
408def _make_crs(srid, pp, axis, uom):
409 crs = Object()
411 crs.srid = srid
413 with warnings.catch_warnings():
414 warnings.simplefilter('ignore')
415 try:
416 crs.proj4text = pp.to_proj4()
417 except pyproj.exceptions.CRSError:
418 gws.log.error(f'CRS: cannot convert {srid!r} to proj4')
419 return None
421 crs.wkt = pp.to_wkt()
423 crs.axis = axis
424 crs.uom = uom
425 crs.coordinatePrecision = COORDINATE_PRECISION_M if uom == gws.Uom.m else COORDINATE_PRECISION_DEG
427 crs.isGeographic = pp.is_geographic
428 crs.isProjected = pp.is_projected
429 crs.isYX = crs.axis == gws.Axis.yx
431 crs.epsg = _unparse(crs.srid, gws.CrsFormat.epsg)
432 crs.urn = _unparse(crs.srid, gws.CrsFormat.urn)
433 crs.urnx = _unparse(crs.srid, gws.CrsFormat.urnx)
434 crs.url = _unparse(crs.srid, gws.CrsFormat.url)
435 crs.uri = _unparse(crs.srid, gws.CrsFormat.uri)
437 # see https://proj.org/schemas/v0.5/projjson.schema.json
438 d = pp.to_json_dict()
440 crs.name = d.get('name') or str(crs.srid)
442 def _datum(x):
443 if 'datum_ensemble' in x:
444 return x['datum_ensemble']['name']
445 if 'datum' in x:
446 return x['datum']['name']
447 return ''
449 def _bbox(d):
450 b = d.get('bbox')
451 if b:
452 # pyproj 3.6
453 return b
454 if d.get('usages'):
455 # pyproj 3.7
456 for u in d['usages']:
457 b = u.get('bbox')
458 if b:
459 return b
461 b = d.get('base_crs')
462 if b:
463 crs.base = b['id']['code']
464 crs.datum = _datum(b)
465 else:
466 crs.base = 0
467 crs.datum = _datum(d)
469 b = _bbox(d)
470 if not b:
471 gws.log.error(f'CRS: no bbox for {crs.srid!r}')
472 return
474 crs.wgsExtent = (
475 b['west_longitude'],
476 b['south_latitude'],
477 b['east_longitude'],
478 b['north_latitude'],
479 )
480 crs.extent = _transform_extent_constrained(crs.wgsExtent, WGS84.srid, srid)
481 crs.bounds = gws.Bounds(extent=crs.extent, crs=crs)
483 return crs
486_AXES_AND_UNITS = {
487 'Easting/metre,Northing/metre': (gws.Axis.xy, gws.Uom.m),
488 'Northing/metre,Easting/metre': (gws.Axis.yx, gws.Uom.m),
489 'Geodetic latitude/degree,Geodetic longitude/degree': (gws.Axis.yx, gws.Uom.deg),
490 'Geodetic longitude/degree,Geodetic latitude/degree': (gws.Axis.xy, gws.Uom.deg),
491 'Easting/US survey foot,Northing/US survey foot': (gws.Axis.xy, gws.Uom.us_ft),
492 'Easting/foot,Northing/foot': (gws.Axis.xy, gws.Uom.ft),
493}
496def _axis_and_unit(pp):
497 ax = []
498 for a in pp.axis_info:
499 ax.append(a.name + '/' + a.unit_name)
500 return _AXES_AND_UNITS.get(','.join(ax))
503##
505"""
506Projections can be referenced by:
508 - int/numeric SRID: 4326
509 - EPSG Code: EPSG:4326
510 - OGC HTTP URL: http://www.opengis.net/gml/srs/epsg.xml#4326
511 - OGC Experimental URN: urn:x-ogc:def:crs:EPSG:4326
512 - OGC URN: urn:ogc:def:crs:EPSG::4326
513 - OGC HTTP URI: http://www.opengis.net/def/crs/EPSG/0/4326
515# https://docs.geoserver.org/stable/en/user/services/wfs/webadmin.html#gml
516"""
518_WRITE_FORMATS = {
519 gws.CrsFormat.srid: '{:d}',
520 gws.CrsFormat.epsg: 'EPSG:{:d}',
521 gws.CrsFormat.url: 'http://www.opengis.net/gml/srs/epsg.xml#{:d}',
522 gws.CrsFormat.uri: 'http://www.opengis.net/def/crs/epsg/0/{:d}',
523 gws.CrsFormat.urnx: 'urn:x-ogc:def:crs:EPSG:{:d}',
524 gws.CrsFormat.urn: 'urn:ogc:def:crs:EPSG::{:d}',
525}
527_PARSE_FORMATS = {
528 gws.CrsFormat.srid: r'^(\d+)$',
529 gws.CrsFormat.epsg: r'^epsg:(\d+)$',
530 gws.CrsFormat.url: r'^http://www.opengis.net/gml/srs/epsg.xml#(\d+)$',
531 gws.CrsFormat.uri: r'http://www.opengis.net/def/crs/epsg/0/(\d+)$',
532 gws.CrsFormat.urnx: r'^urn:x-ogc:def:crs:epsg:(\d+)$',
533 gws.CrsFormat.urn: r'^urn:ogc:def:crs:epsg:[0-9.]*:(\d+)$',
534}
536# @TODO
538_aliases = {
539 'crs:84': 4326,
540 'crs84': 4326,
541 'urn:ogc:def:crs:ogc:1.3:crs84': 'urn:ogc:def:crs:epsg::4326',
542 'wgs84': 4326,
543 'epsg:900913': 3857,
544 'epsg:102100': 3857,
545 'epsg:102113': 3857,
546}
548# https://docs.geoserver.org/latest/en/user/services/wfs/axis_order.html
549# EPSG:4326 longitude/latitude
550# http://www.opengis.net/gml/srs/epsg.xml#xxxx longitude/latitude
551# urn:x-ogc:def:crs:EPSG:xxxx latitude/longitude
552# urn:ogc:def:crs:EPSG::4326 latitude/longitude
554_AXIS_FOR_FORMAT = {
555 gws.CrsFormat.srid: gws.Axis.xy,
556 gws.CrsFormat.epsg: gws.Axis.xy,
557 gws.CrsFormat.url: gws.Axis.xy,
558 gws.CrsFormat.uri: gws.Axis.xy,
559 gws.CrsFormat.urnx: gws.Axis.yx,
560 gws.CrsFormat.urn: gws.Axis.yx,
561}
564def _parse(crs_name):
565 if isinstance(crs_name, int):
566 return gws.CrsFormat.epsg, crs_name
568 if isinstance(crs_name, bytes):
569 crs_name = crs_name.decode('ascii').lower()
571 if isinstance(crs_name, str):
572 crs_name = crs_name.lower()
574 if crs_name in {'crs84', 'crs:84'}:
575 return gws.CrsFormat.crs, 4326
577 if crs_name in _aliases:
578 crs_name = _aliases[crs_name]
579 if isinstance(crs_name, int):
580 return gws.CrsFormat.epsg, int(crs_name)
582 for fmt, r in _PARSE_FORMATS.items():
583 m = re.match(r, crs_name)
584 if m:
585 return fmt, int(m.group(1))
587 return None, 0
590def _unparse(srid, fmt):
591 return _WRITE_FORMATS[fmt].format(srid)
594##
597_obj_cache: dict = {
598 WGS84.srid: WGS84,
599 WEBMERCATOR.url: WEBMERCATOR,
600}
602_pyproj_cache: dict = {}
604_transformer_cache: dict = {}