Coverage for gws-app/gws/lib/crs/__init__.py: 93%

1"""Coordinate reference systems. 

2 

3This package provides ``gws.Crs`` objects, which describe coordinate reference systems 

4and transform extents, points and resolutions between them. CRS objects are created 

5from EPSG definitions with ``pyproj``, and are cached, so that there is one object per SRID. 

6Only CRSs with meter or degree units are supported. CRS objects with the same SRID compare equal. 

7 

8``transform_resolution`` samples nine points of the extent (corners, edge midpoints, centre), 

9transforms each together with a one-pixel step in x and in y, and returns the smallest 

10finite positive step length in the target CRS. 

11 

12A CRS can be referenced by a ``gws.CrsName`` in one of these formats (see ``gws.CrsFormat``): 

13 

14- numeric SRID: ``4326`` 

15- EPSG code: ``EPSG:4326`` 

16- OGC HTTP URL: ``http://www.opengis.net/gml/srs/epsg.xml#4326`` 

17- OGC experimental URN: ``urn:x-ogc:def:crs:EPSG:4326`` 

18- OGC URN: ``urn:ogc:def:crs:EPSG::4326`` 

19- OGC HTTP URI: ``http://www.opengis.net/def/crs/EPSG/0/4326`` 

20 

21Names are case-insensitive. Some aliases, like ``CRS:84`` or ``EPSG:900913``, are also recognized. 

22 

23The package provides: 

24 

25- predefined CRS objects ``WGS84`` and ``WEBMERCATOR`` and related constants, 

26- ``get``, ``require`` and ``parse`` to look up a CRS by name, 

27- ``best_match`` to pick a CRS from a list of supported CRSs, 

28- ``qgis_extent_width`` to compute the width of a geographic extent the way QGIS does. 

29 

30The ``gws.Crs`` interface itself is defined in ``types.pyinc``. 

31 

32Example:: 

33 

34 import gws.lib.crs 

35 

36 crs = gws.lib.crs.require('EPSG:25832') 

37 crs.to_string(gws.CrsFormat.urn) # 'urn:ogc:def:crs:EPSG::25832' 

38 

39 ext = gws.lib.crs.WGS84.transform_extent((5.0, 47.0, 15.0, 55.0), crs) 

40 

41 fmt, crs = gws.lib.crs.parse('urn:ogc:def:crs:EPSG::4326') 

42 crs.axis_for_format(fmt) # gws.Axis.yx 

43""" 

44 

45from typing import Optional 

46 

47import math 

48import re 

49import warnings 

50 

51import pyproj.crs 

52import pyproj.exceptions 

53import pyproj.transformer 

54 

55import gws 

56 

57 

58## 

59 

60 

61class Object(gws.Crs): 

62 """Coordinate reference system.""" 

63 

64 def __init__(self, **kwargs): 

65 """Create a CRS object with the given attributes. 

66 

67 Args: 

68 **kwargs: Attribute values, see ``gws.Crs``. 

69 """ 

70 vars(self).update(kwargs) 

71 

72 # crs objects with the same srid must be equal 

73 # (despite caching, they can be different due to pickling) 

74 

75 def __hash__(self): 

76 return self.srid 

77 

78 def __eq__(self, other): 

79 return isinstance(other, Object) and other.srid == self.srid 

80 

81 def __repr__(self): 

82 return f'<crs:{self.srid}>' 

83 

84 def axis_for_format(self, fmt): 

85 if not self.isYX: 

86 return self.axis 

87 return _AXIS_FOR_FORMAT.get(fmt, self.axis) 

88 

89 def transform_extent(self, ext, crs_to): 

90 if crs_to == self: 

91 return ext 

92 return _transform_extent_check(ext, self.srid, crs_to.srid) 

93 

94 def transform_resolution(self, extent, res, crs_to): 

95 tr = self.transformer(crs_to) 

96 

97 x0, y0, x1, y1 = extent 

98 xm = (x0 + x1) / 2 

99 ym = (y0 + y1) / 2 

100 

101 points = [ 

102 (x0, y0), 

103 (xm, y0), 

104 (x1, y0), 

105 (x0, ym), 

106 (xm, ym), 

107 (x1, ym), 

108 (x0, y1), 

109 (xm, y1), 

110 (x1, y1), 

111 ] 

112 

113 ds = [] 

114 for x, y in points: 

115 ax, ay = tr(x, y) 

116 bx, by = tr(x + res, y) 

117 cx, cy = tr(x, y + res) 

118 ds.append(math.hypot(bx - ax, by - ay)) 

119 ds.append(math.hypot(cx - ax, cy - ay)) 

120 

121 ds = [d for d in ds if math.isfinite(d) and d > 0] 

122 return min(ds) if ds else 0.0 

123 

124 def clip_wgs_extent(self, wgs_extent): 

125 a = wgs_extent 

126 b = self.wgsMaxExtent 

127 x0, y0, x1, y1 = max(a[0], b[0]), max(a[1], b[1]), min(a[2], b[2]), min(a[3], b[3]) 

128 if x0 >= x1 or y0 >= y1: 

129 return None 

130 return x0, y0, x1, y1 

131 

132 def transformer(self, crs_to): 

133 tr = _pyproj_transformer(self.srid, crs_to.srid) 

134 return tr.transform 

135 

136 def extent_size_in_meters(self, extent): 

137 x0, y0, x1, y1 = extent 

138 

139 if self.isProjected: 

140 if self.uom != gws.Uom.m: 

141 # @TODO support non-meter crs 

142 raise Error(f'unsupported unit: {self.uom}') 

143 return abs(x1 - x0), abs(y1 - y0) 

144 

145 geod = pyproj.Geod(ellps='WGS84') 

146 

147 mid_lat = (y0 + y1) / 2 

148 _, _, w = geod.inv(x0, mid_lat, x1, mid_lat) 

149 mid_lon = (x0 + x1) / 2 

150 _, _, h = geod.inv(mid_lon, y0, mid_lon, y1) 

151 

152 return w, h 

153 

154 def point_offset_in_meters(self, xy, dist, az): 

155 x, y = xy 

156 

157 if self.isProjected: 

158 if self.uom != gws.Uom.m: 

159 # @TODO support non-meter crs 

160 raise Error(f'unsupported unit: {self.uom}') 

161 

162 if az == 0: 

163 return x, y + dist 

164 if az == 90: 

165 return x + dist, y 

166 if az == 180: 

167 return x, y - dist 

168 if az == 270: 

169 return x - dist, y 

170 

171 az_rad = math.radians(90 - az) 

172 return ( 

173 x + dist * math.cos(az_rad), 

174 y + dist * math.sin(az_rad), 

175 ) 

176 

177 geod = pyproj.Geod(ellps='WGS84') 

178 x, y, _ = geod.fwd(x, y, dist=dist, az=az) 

179 return x, y 

180 

181 def to_string(self, fmt=None): 

182 fmt = fmt or gws.CrsFormat.epsg 

183 if fmt == gws.CrsFormat.srid: 

184 return str(self.srid) 

185 return getattr(self, str(fmt).lower()) 

186 

187 def to_geojson(self): 

188 # https://geojson.org/geojson-spec#named-crs 

189 return { 

190 'type': 'name', 

191 'properties': { 

192 'name': self.urn, 

193 }, 

194 } 

195 

196 

197## 

198 

199 

200def qgis_extent_width(extent: gws.Extent) -> float: 

201 """Compute the width of a geographic extent in meters, the way QGIS does. 

202 

203 This is a port of ``QgsScaleCalculator::calculateGeographicDistance`` from QGIS. 

204 The distance is measured along the middle latitude of the extent. 

205 

206 Args: 

207 extent: Extent in degrees. 

208 

209 Returns: 

210 The width in meters. 

211 """ 

212 # straight port from QGIS/src/core/qgsscalecalculator.cpp QgsScaleCalculator::calculateGeographicDistance 

213 x0, y0, x1, y1 = extent 

214 

215 lat = (y0 + y1) * 0.5 

216 RADS = (4.0 * math.atan(1.0)) / 180.0 

217 a = math.pow(math.cos(lat * RADS), 2) 

218 c = 2.0 * math.atan2(math.sqrt(a), math.sqrt(1.0 - a)) 

219 RA = 6378000 

220 E = 0.0810820288 

221 radius = RA * (1.0 - E * E) / math.pow(1.0 - E * E * math.sin(lat * RADS) * math.sin(lat * RADS), 1.5) 

222 return (x1 - x0) / 180.0 * radius * c 

223 

224 

225## 

226 

227# enough precision to represent 1cm 

228COORDINATE_PRECISION_DEG = 7 

229COORDINATE_PRECISION_M = 2 

230 

231WGS84: gws.Crs = Object( 

232 srid=4326, 

233 proj4text='+proj=longlat +datum=WGS84 +no_defs +type=crs', 

234 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]]', 

235 axis=gws.Axis.yx, 

236 uom=gws.Uom.deg, 

237 isGeographic=True, 

238 isProjected=False, 

239 isYX=True, 

240 epsg='EPSG:4326', 

241 urn='urn:ogc:def:crs:EPSG::4326', 

242 urnx='urn:x-ogc:def:crs:EPSG:4326', 

243 url='http://www.opengis.net/gml/srs/epsg.xml#4326', 

244 uri='http://www.opengis.net/def/crs/epsg/0/4326', 

245 name='WGS 84', 

246 base=0, 

247 datum='World Geodetic System 1984 ensemble', 

248 wgsExtent=(-180, -90, 180, 90), 

249 extent=(-180, -90, 180, 90), 

250 wgsMaxExtent=(-180, -90, 180, 90), 

251 coordinatePrecision=COORDINATE_PRECISION_DEG, 

252) 

253"""WGS 84 geographic CRS (EPSG:4326).""" 

254 

255WGS84.bounds = gws.Bounds(crs=WGS84, extent=WGS84.extent) 

256 

257WEBMERCATOR: gws.Crs = Object( 

258 srid=3857, 

259 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', 

260 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]]', 

261 axis=gws.Axis.xy, 

262 uom=gws.Uom.m, 

263 isGeographic=False, 

264 isProjected=True, 

265 isYX=False, 

266 epsg='EPSG:3857', 

267 urn='urn:ogc:def:crs:EPSG::3857', 

268 urnx='urn:x-ogc:def:crs:EPSG:3857', 

269 url='http://www.opengis.net/gml/srs/epsg.xml#3857', 

270 uri='http://www.opengis.net/def/crs/epsg/0/3857', 

271 name='WGS 84 / Pseudo-Mercator', 

272 base=4326, 

273 datum='World Geodetic System 1984 ensemble', 

274 wgsExtent=(-180, -85.06, 180, 85.06), 

275 extent=( 

276 -20037508.342789244, 

277 -20048966.104014598, 

278 20037508.342789244, 

279 20048966.104014598, 

280 ), 

281 wgsMaxExtent=(-180, -85.06, 180, 85.06), 

282 coordinatePrecision=COORDINATE_PRECISION_M, 

283) 

284"""WGS 84 / Pseudo-Mercator CRS (EPSG:3857).""" 

285 

286WEBMERCATOR.bounds = gws.Bounds(crs=WEBMERCATOR, extent=WEBMERCATOR.extent) 

287 

288WEBMERCATOR_RADIUS = 6378137 

289"""Radius of the web mercator sphere (the WGS84 semi-major axis), metres.""" 

290 

291METERS_PER_DEGREE = 2 * math.pi * WEBMERCATOR_RADIUS / 360 

292"""Metres per degree at the equator; the OGC convention for scale denominators in geographic CRS (WMTS 1.0, 6.1).""" 

293 

294WEBMERCATOR_SQUARE = ( 

295 -math.pi * WEBMERCATOR_RADIUS, 

296 -math.pi * WEBMERCATOR_RADIUS, 

297 +math.pi * WEBMERCATOR_RADIUS, 

298 +math.pi * WEBMERCATOR_RADIUS, 

299) 

300"""Square web mercator extent that covers the whole world width, in meters.""" 

301 

302 

303class Error(gws.Error): 

304 """CRS error.""" 

305 

306 pass 

307 

308 

309def get(crs_name: Optional[gws.CrsName]) -> Optional[gws.Crs]: 

310 """Get the CRS for a given CRS name or SRID. 

311 

312 Args: 

313 crs_name: CRS name in any supported format, or an SRID. 

314 

315 Returns: 

316 The CRS object, or ``None`` if the name is empty, cannot be parsed or refers to an unsupported CRS. 

317 """ 

318 if not crs_name: 

319 return None 

320 crs, err = _get_crs(crs_name) 

321 if err: 

322 gws.log.warning(err) 

323 return crs 

324 

325 

326def is_valid(crs_name: Optional[gws.CrsName]) -> bool: 

327 """Check if a CRS name or SRID refers to a supported CRS, without logging. 

328 

329 Args: 

330 crs_name: CRS name in any supported format, or an SRID. 

331 

332 Returns: 

333 ``True`` if ``get`` would return a CRS. 

334 """ 

335 if not crs_name: 

336 return False 

337 crs, _ = _get_crs(crs_name) 

338 return crs is not None 

339 

340 

341def parse(crs_name: gws.CrsName) -> tuple[gws.CrsFormat, Optional[gws.Crs]]: 

342 """Parse a CRS name into its format and the CRS itself. 

343 

344 Args: 

345 crs_name: CRS name in any supported format, or an SRID. 

346 

347 Returns: 

348 A tuple of the name format and the CRS object. If the name cannot be parsed, 

349 ``(CrsFormat.none, None)``. If the CRS is unknown or unsupported, the CRS is ``None``. 

350 """ 

351 fmt, srid = _parse(crs_name) 

352 if not fmt: 

353 return gws.CrsFormat.none, None 

354 crs, err = _get_crs(srid) 

355 if err: 

356 gws.log.warning(err) 

357 return fmt, crs 

358 

359 

360def require(crs_name: gws.CrsName) -> gws.Crs: 

361 """Get the CRS for a given CRS name or SRID, and fail if there is none. 

362 

363 Args: 

364 crs_name: CRS name in any supported format, or an SRID. 

365 

366 Returns: 

367 The CRS object. 

368 

369 Raises: 

370 ``Error``: If the name cannot be parsed or refers to an unknown or unsupported CRS. 

371 """ 

372 crs, err = _get_crs(crs_name) 

373 if not crs: 

374 raise Error(err) 

375 return crs 

376 

377 

378## 

379 

380 

381def best_match(crs: gws.Crs, supported_crs: list[gws.Crs]) -> gws.Crs: 

382 """Return a CRS from the list that most closely matches the given CRS. 

383 

384 If the CRS is in the list, it is returned. Otherwise, for a projected CRS, web mercator 

385 or the first projected CRS from the list is preferred, and for a geographic CRS, 

386 WGS84 or the first geographic CRS. Failing that, the first CRS from the list is returned, 

387 or the given CRS if the list is empty. 

388 

389 Args: 

390 crs: Target CRS. 

391 supported_crs: List of supported CRSs. 

392 

393 Returns: 

394 A CRS object. 

395 """ 

396 

397 if crs in supported_crs: 

398 return crs 

399 

400 bst = _best_match(crs, supported_crs) 

401 if not bst: 

402 bst = supported_crs[0] if supported_crs else crs 

403 gws.log.debug(f'CRS: best_crs: using {bst.srid!r} for {crs.srid!r}') 

404 return bst 

405 

406 

407def _best_match(crs, supported_crs): 

408 # @TODO find a projection with less errors 

409 # @TODO find a projection with same units 

410 

411 if crs.isProjected: 

412 # for a projected crs, find webmercator 

413 for sup in supported_crs: 

414 if sup.srid == WEBMERCATOR.srid: 

415 return sup 

416 

417 # not found, return the first projected crs 

418 for sup in supported_crs: 

419 if sup.isProjected: 

420 return sup 

421 

422 if crs.isGeographic: 

423 # for a geographic crs, try wgs first 

424 for sup in supported_crs: 

425 if sup.srid == WGS84.srid: 

426 return sup 

427 

428 # not found, return the first geographic crs 

429 for sup in supported_crs: 

430 if sup.isGeographic: 

431 return sup 

432 

433 

434## 

435 

436 

437def _get_crs(crs_name) -> tuple[Optional[gws.Crs], str]: 

438 if crs_name in _obj_cache: 

439 return _obj_cache[crs_name] 

440 

441 fmt, srid = _parse(crs_name) 

442 if not fmt: 

443 res = None, f'CRS: cannot parse {crs_name!r}' 

444 _obj_cache[crs_name] = res 

445 return res 

446 

447 if srid in _obj_cache: 

448 _obj_cache[crs_name] = _obj_cache[srid] 

449 return _obj_cache[srid] 

450 

451 res = _load_crs(srid) 

452 _obj_cache[crs_name] = _obj_cache[srid] = res 

453 return res 

454 

455 

456def _load_crs(srid) -> tuple[Optional[gws.Crs], str]: 

457 pp = _pyproj_crs_object(srid) 

458 if not pp: 

459 return None, f'CRS: unknown srid {srid!r}' 

460 

461 au = _axis_and_unit(pp) 

462 if not au: 

463 return None, f'CRS: unsupported srid {srid!r}' 

464 

465 axis, uom = au 

466 if uom not in (gws.Uom.m, gws.Uom.deg): 

467 return None, f'CRS: unsupported unit {uom!r} for {srid!r}' 

468 

469 return _make_crs(srid, pp, axis, uom) 

470 

471 

472def _pyproj_crs_object(srid) -> Optional[pyproj.CRS]: 

473 if srid in _pyproj_cache: 

474 return _pyproj_cache[srid] 

475 

476 try: 

477 pp = pyproj.CRS.from_epsg(srid) 

478 except pyproj.exceptions.CRSError: 

479 return None 

480 

481 _pyproj_cache[srid] = pp 

482 return _pyproj_cache[srid] 

483 

484 

485def _pyproj_transformer(srid_from, srid_to) -> pyproj.transformer.Transformer: 

486 key = srid_from, srid_to 

487 

488 if key in _transformer_cache: 

489 return _transformer_cache[key] 

490 

491 pa = _pyproj_crs_object(srid_from) 

492 pb = _pyproj_crs_object(srid_to) 

493 

494 _transformer_cache[key] = pyproj.transformer.Transformer.from_crs(pa, pb, always_xy=True) 

495 return _transformer_cache[key] 

496 

497 

498def _transform_extent_check(ext, srid_from, srid_to): 

499 ext_nor = _normalize_extent(ext) 

500 

501 if srid_from == WGS84.srid: 

502 ext_wgs = ext_nor 

503 else: 

504 tr_to_wgs = _pyproj_transformer(srid_from, WGS84.srid) 

505 ext_wgs = tr_to_wgs.transform_bounds(ext_nor[0], ext_nor[1], ext_nor[2], ext_nor[3]) 

506 

507 if srid_to == WGS84.srid: 

508 return _normalize_extent(ext_wgs) 

509 

510 pp = _pyproj_crs_object(srid_to) 

511 if not pp: 

512 raise Error(f'_transform_extent: unknown {srid_to=}') 

513 

514 tr_from_wgs = _pyproj_transformer(WGS84.srid, srid_to) 

515 au = pp.area_of_use 

516 

517 if au: 

518 ext_au = au.bounds 

519 if _is_big_extent(ext_wgs) and not _is_big_extent(ext_au): 

520 ext_to = _transform_extent_sampled(ext_wgs, tr_from_wgs, au) 

521 if ext_to: 

522 gws.log.debug(f'transform_extent: {ext=} {srid_from!r}->{srid_to!r}: big extent {ext_to=} ') 

523 return _normalize_extent(ext_to) 

524 

525 if ext_wgs[2] < ext_au[0] or ext_wgs[0] > ext_au[2] or ext_wgs[3] < ext_au[1] or ext_wgs[1] > ext_au[3]: 

526 gws.log.debug(f'transform_extent: {ext=} {srid_from!r}->{srid_to!r}: outside of AoU ') 

527 

528 ext_to = tr_from_wgs.transform_bounds(ext_wgs[0], ext_wgs[1], ext_wgs[2], ext_wgs[3]) 

529 

530 return _normalize_extent(ext_to) 

531 

532 

533def _transform_extent_sampled(ext_wgs, tr, au, samples=50): 

534 x0, y0, x1, y1 = ext_wgs 

535 ax0, ay0, ax1, ay1 = au.bounds 

536 mid_lon = (ax0 + ax1) / 2 

537 

538 xys = [] 

539 

540 # 1) Dense grid sampling within the area of use (accurate core extent) 

541 for i in range(samples + 1): 

542 lon = ax0 + (ax1 - ax0) * i / samples 

543 for j in range(samples + 1): 

544 lat = ay0 + (ay1 - ay0) * j / samples 

545 try: 

546 xys.append(tr.transform(lon, lat, errcheck=True)) 

547 except Exception: 

548 pass 

549 

550 # 2) Sample the full latitude range along the central meridian of the AoU 

551 # This captures the full Y extent for "the world" in the target projection 

552 for i in range(samples + 1): 

553 lat = y0 + (y1 - y0) * i / samples 

554 try: 

555 xys.append(tr.transform(mid_lon, lat, errcheck=True)) 

556 except Exception: 

557 pass 

558 

559 # 3) Sample several meridians across the full longitude range 

560 # to capture the full X spread at various latitudes 

561 for i in range(samples + 1): 

562 lon = x0 + (x1 - x0) * i / samples 

563 for j in range(samples + 1): 

564 lat = y0 + (y1 - y0) * j / samples 

565 try: 

566 xys.append(tr.transform(lon, lat, errcheck=True)) 

567 except Exception: 

568 pass 

569 

570 xs = [x for x, y in xys if math.isfinite(x) and math.isfinite(y)] 

571 ys = [y for x, y in xys if math.isfinite(x) and math.isfinite(y)] 

572 

573 if not xs: 

574 return 

575 

576 return (min(xs), min(ys), max(xs), max(ys)) 

577 

578 

579def _is_big_extent(ext_wgs): 

580 dx = abs(ext_wgs[2] - ext_wgs[0]) 

581 dy = abs(ext_wgs[3] - ext_wgs[1]) 

582 return dx > 350 or dy > 160 

583 

584 

585def _transform_extent_direct(ext, srid_from, srid_to): 

586 tr = _pyproj_transformer(srid_from, srid_to) 

587 

588 ext_nor = _normalize_extent(ext) 

589 

590 res = tr.transform_bounds( 

591 left=ext_nor[0], 

592 bottom=ext_nor[1], 

593 right=ext_nor[2], 

594 top=ext_nor[3], 

595 errcheck=True, 

596 ) 

597 return _normalize_extent(res) 

598 

599 

600def _normalize_extent(ext): 

601 return ( 

602 min(ext[0], ext[2]), 

603 min(ext[1], ext[3]), 

604 max(ext[0], ext[2]), 

605 max(ext[1], ext[3]), 

606 ) 

607 

608 

609def _make_crs(srid, pp, axis, uom): 

610 crs = Object() 

611 

612 crs.srid = srid 

613 

614 with warnings.catch_warnings(): 

615 warnings.simplefilter('ignore') 

616 try: 

617 crs.proj4text = pp.to_proj4() 

618 except pyproj.exceptions.CRSError: 

619 return None, f'CRS: cannot convert {srid!r} to proj4' 

620 

621 crs.wkt = pp.to_wkt() 

622 

623 crs.axis = axis 

624 crs.uom = uom 

625 crs.coordinatePrecision = COORDINATE_PRECISION_M if uom == gws.Uom.m else COORDINATE_PRECISION_DEG 

626 

627 crs.isGeographic = pp.is_geographic 

628 crs.isProjected = pp.is_projected 

629 crs.isYX = crs.axis == gws.Axis.yx 

630 

631 crs.epsg = _unparse(crs.srid, gws.CrsFormat.epsg) 

632 crs.urn = _unparse(crs.srid, gws.CrsFormat.urn) 

633 crs.urnx = _unparse(crs.srid, gws.CrsFormat.urnx) 

634 crs.url = _unparse(crs.srid, gws.CrsFormat.url) 

635 crs.uri = _unparse(crs.srid, gws.CrsFormat.uri) 

636 

637 # see https://proj.org/schemas/v0.5/projjson.schema.json 

638 d = pp.to_json_dict() 

639 

640 crs.name = d.get('name') or str(crs.srid) 

641 

642 def _datum(x): 

643 if 'datum_ensemble' in x: 

644 return x['datum_ensemble']['name'] 

645 if 'datum' in x: 

646 return x['datum']['name'] 

647 return '' 

648 

649 def _bbox(d): 

650 b = d.get('bbox') 

651 if b: 

652 # pyproj 3.6 

653 return b 

654 if d.get('usages'): 

655 # pyproj 3.7 

656 for u in d['usages']: 

657 b = u.get('bbox') 

658 if b: 

659 return b 

660 

661 b = d.get('base_crs') 

662 if b: 

663 crs.base = b['id']['code'] 

664 crs.datum = _datum(b) 

665 else: 

666 crs.base = 0 

667 crs.datum = _datum(d) 

668 

669 b = _bbox(d) 

670 if not b: 

671 return None, f'CRS: no bbox for {crs.srid!r}' 

672 

673 crs.wgsExtent = ( 

674 b['west_longitude'], 

675 b['south_latitude'], 

676 b['east_longitude'], 

677 b['north_latitude'], 

678 ) 

679 crs.extent = _transform_extent_check(crs.wgsExtent, WGS84.srid, srid) 

680 crs.bounds = gws.Bounds(extent=crs.extent, crs=crs) 

681 

682 crs.wgsMaxExtent = _wgs_max_extent(pp, crs.wgsExtent) 

683 

684 return crs, '' 

685 

686 

687def _wgs_max_extent(pp, wgs_extent): 

688 """The datum's area of use, or, for global datums, the own extent widened by its width on each side.""" 

689 

690 geo = pp.geodetic_crs 

691 au = geo.area_of_use if geo else None 

692 if au and not _is_big_extent(au.bounds): 

693 return _normalize_extent(au.bounds) 

694 

695 x0, y0, x1, y1 = wgs_extent 

696 w = x1 - x0 

697 return max(x0 - w, -180), y0, min(x1 + w, 180), y1 

698 

699 

700_AXES_AND_UNITS = { 

701 'Easting/metre,Northing/metre': (gws.Axis.xy, gws.Uom.m), 

702 'Northing/metre,Easting/metre': (gws.Axis.yx, gws.Uom.m), 

703 'Geodetic latitude/degree,Geodetic longitude/degree': (gws.Axis.yx, gws.Uom.deg), 

704 'Geodetic longitude/degree,Geodetic latitude/degree': (gws.Axis.xy, gws.Uom.deg), 

705 'Easting/US survey foot,Northing/US survey foot': (gws.Axis.xy, gws.Uom.us_ft), 

706 'Easting/foot,Northing/foot': (gws.Axis.xy, gws.Uom.ft), 

707} 

708 

709 

710def _axis_and_unit(pp): 

711 ax = [] 

712 for a in pp.axis_info: 

713 ax.append(a.name + '/' + a.unit_name) 

714 return _AXES_AND_UNITS.get(','.join(ax)) 

715 

716 

717## 

718 

719""" 

720Projections can be referenced by: 

721 

722 - int/numeric SRID: 4326 

723 - EPSG Code: EPSG:4326 

724 - OGC HTTP URL: http://www.opengis.net/gml/srs/epsg.xml#4326 

725 - OGC Experimental URN: urn:x-ogc:def:crs:EPSG:4326 

726 - OGC URN: urn:ogc:def:crs:EPSG::4326 

727 - OGC HTTP URI: http://www.opengis.net/def/crs/EPSG/0/4326 

728 

729# https://docs.geoserver.org/stable/en/user/services/wfs/webadmin.html#gml 

730""" 

731 

732_WRITE_FORMATS = { 

733 gws.CrsFormat.srid: '{:d}', 

734 gws.CrsFormat.epsg: 'EPSG:{:d}', 

735 gws.CrsFormat.url: 'http://www.opengis.net/gml/srs/epsg.xml#{:d}', 

736 gws.CrsFormat.uri: 'http://www.opengis.net/def/crs/epsg/0/{:d}', 

737 gws.CrsFormat.urnx: 'urn:x-ogc:def:crs:EPSG:{:d}', 

738 gws.CrsFormat.urn: 'urn:ogc:def:crs:EPSG::{:d}', 

739} 

740 

741_PARSE_FORMATS = { 

742 gws.CrsFormat.srid: r'^(\d+)$', 

743 gws.CrsFormat.epsg: r'^epsg:(\d+)$', 

744 gws.CrsFormat.url: r'^http://www.opengis.net/gml/srs/epsg.xml#(\d+)$', 

745 gws.CrsFormat.uri: r'http://www.opengis.net/def/crs/epsg/0/(\d+)$', 

746 gws.CrsFormat.urnx: r'^urn:x-ogc:def:crs:epsg:(\d+)$', 

747 gws.CrsFormat.urn: r'^urn:ogc:def:crs:epsg:[0-9.]*:(\d+)$', 

748} 

749 

750# @TODO 

751 

752_aliases = { 

753 'crs:84': 4326, 

754 'crs84': 4326, 

755 'urn:ogc:def:crs:ogc:1.3:crs84': 'urn:ogc:def:crs:epsg::4326', 

756 'wgs84': 4326, 

757 'epsg:900913': 3857, 

758 'epsg:102100': 3857, 

759 'epsg:102113': 3857, 

760} 

761 

762# https://docs.geoserver.org/latest/en/user/services/wfs/axis_order.html 

763# EPSG:4326 longitude/latitude 

764# http://www.opengis.net/gml/srs/epsg.xml#xxxx longitude/latitude 

765# urn:x-ogc:def:crs:EPSG:xxxx latitude/longitude 

766# urn:ogc:def:crs:EPSG::4326 latitude/longitude 

767 

768_AXIS_FOR_FORMAT = { 

769 gws.CrsFormat.srid: gws.Axis.xy, 

770 gws.CrsFormat.epsg: gws.Axis.xy, 

771 gws.CrsFormat.url: gws.Axis.xy, 

772 gws.CrsFormat.uri: gws.Axis.xy, 

773 gws.CrsFormat.urnx: gws.Axis.yx, 

774 gws.CrsFormat.urn: gws.Axis.yx, 

775} 

776 

777 

778def _parse(crs_name): 

779 if isinstance(crs_name, int): 

780 return gws.CrsFormat.epsg, crs_name 

781 

782 if isinstance(crs_name, bytes): 

783 crs_name = crs_name.decode('ascii').lower() 

784 

785 if isinstance(crs_name, str): 

786 crs_name = crs_name.lower() 

787 

788 if crs_name in {'crs84', 'crs:84'}: 

789 return gws.CrsFormat.crs, 4326 

790 

791 if crs_name in _aliases: 

792 crs_name = _aliases[crs_name] 

793 if isinstance(crs_name, int): 

794 return gws.CrsFormat.epsg, int(crs_name) 

795 

796 for fmt, r in _PARSE_FORMATS.items(): 

797 m = re.match(r, crs_name) 

798 if m: 

799 return fmt, int(m.group(1)) 

800 

801 return None, 0 

802 

803 

804def _unparse(srid, fmt): 

805 return _WRITE_FORMATS[fmt].format(srid) 

806 

807 

808## 

809 

810