Coverage for gws-app/gws/lib/crs/__init_

1"""Coordinate reference systems.

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.

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.

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

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``

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

23The package provides:

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.

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

32Example::

34 import gws.lib.crs

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

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

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

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

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

43"""

45from typing import Optional

47import math

48import re

49import warnings

51import pyproj.crs

52import pyproj.exceptions

53import pyproj.transformer

55import gws

58##

61class Object(gws.Crs):

62 """Coordinate reference system."""

64 def __init__(self, **kwargs):

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

67 Args:

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

69 """

70 vars(self).update(kwargs)

72 # crs objects with the same srid must be equal

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

75 def __hash__(self):

76 return self.srid

78 def __eq__(self, other):

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

81 def __repr__(self):

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

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)

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)

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

95 tr = self.transformer(crs_to)

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 }

197##

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 return _get_crs(crs_name)

321

322

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

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

325

326 Args:

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

328

329 Returns:

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

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

332 """

333 fmt, srid = _parse(crs_name)

334 if not fmt:

335 return gws.CrsFormat.none, None

336 return fmt, _get_crs(srid)

337

338

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

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

341

342 Args:

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

344

345 Returns:

346 The CRS object.

347

348 Raises:

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

350 """

351 crs = _get_crs(crs_name)

352 if not crs:

353 raise Error(f'invalid CRS {crs_name!r}')

354 return crs

357##

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

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

362

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

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

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

366 or the given CRS if the list is empty.

367

368 Args:

369 crs: Target CRS.

370 supported_crs: List of supported CRSs.

371

372 Returns:

373 A CRS object.

374 """

375

376 if crs in supported_crs:

377 return crs

378

379 bst = _best_match(crs, supported_crs)

380 if not bst:

381 bst = supported_crs[0] if supported_crs else crs

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

383 return bst

384

385

386def _best_match(crs, supported_crs):

387 # @TODO find a projection with less errors

388 # @TODO find a projection with same units

389

390 if crs.isProjected:

391 # for a projected crs, find webmercator

392 for sup in supported_crs:

393 if sup.srid == WEBMERCATOR.srid:

394 return sup

395

396 # not found, return the first projected crs

397 for sup in supported_crs:

398 if sup.isProjected:

399 return sup

400

401 if crs.isGeographic:

402 # for a geographic crs, try wgs first

403 for sup in supported_crs:

404 if sup.srid == WGS84.srid:

405 return sup

406

407 # not found, return the first geographic crs

408 for sup in supported_crs:

409 if sup.isGeographic:

410 return sup

413##

416def _get_crs(crs_name):

417 if crs_name in _obj_cache:

418 return _obj_cache[crs_name]

419

420 fmt, srid = _parse(crs_name)

421 if not fmt:

422 gws.log.warning(f'CRS: cannot parse {crs_name!r}')

423 _obj_cache[crs_name] = None

424 return None

425

426 if srid in _obj_cache:

427 _obj_cache[crs_name] = _obj_cache[srid]

428 return _obj_cache[srid]

429

430 obj = _get_new_crs(srid)

431 _obj_cache[crs_name] = _obj_cache[srid] = obj

432 return obj

433

434

435def _get_new_crs(srid):

436 pp = _pyproj_crs_object(srid)

437 if not pp:

438 _warnings[srid] = f'CRS: unknown srid {srid!r}'

439 return None

440

441 au = _axis_and_unit(pp)

442 if not au:

443 _warnings[srid] = f'CRS: unsupported srid {srid!r}'

444 return None

445

446 axis, uom = au

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

448 _warnings[srid] = f'CRS: unsupported unit {uom!r} for {srid!r}'

449 return None

450

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

452

453

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

455 if srid in _pyproj_cache:

456 return _pyproj_cache[srid]

457

458 try:

459 pp = pyproj.CRS.from_epsg(srid)

460 except pyproj.exceptions.CRSError:

461 return None

462

463 _pyproj_cache[srid] = pp

464 return _pyproj_cache[srid]

465

466

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

468 key = srid_from, srid_to

469

470 if key in _transformer_cache:

471 return _transformer_cache[key]

472

473 pa = _pyproj_crs_object(srid_from)

474 pb = _pyproj_crs_object(srid_to)

475

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

477 return _transformer_cache[key]

478

479

480def _transform_extent_check(ext, srid_from, srid_to):

481 ext_nor = _normalize_extent(ext)

482

483 if srid_from == WGS84.srid:

484 ext_wgs = ext_nor

485 else:

486 tr_to_wgs = _pyproj_transformer(srid_from, WGS84.srid)

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

488

489 if srid_to == WGS84.srid:

490 return _normalize_extent(ext_wgs)

491

492 pp = _pyproj_crs_object(srid_to)

493 if not pp:

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

495

496 tr_from_wgs = _pyproj_transformer(WGS84.srid, srid_to)

497 au = pp.area_of_use

498

499 if au:

500 ext_au = au.bounds

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

502 ext_to = _transform_extent_sampled(ext_wgs, tr_from_wgs, au)

503 if ext_to:

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

505 return _normalize_extent(ext_to)

506

507 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]:

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

509

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

511

512 return _normalize_extent(ext_to)

513

514

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

516 x0, y0, x1, y1 = ext_wgs

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

518 mid_lon = (ax0 + ax1) / 2

519

520 xys = []

521

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

523 for i in range(samples + 1):

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

525 for j in range(samples + 1):

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

527 try:

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

529 except Exception:

530 pass

531

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

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

534 for i in range(samples + 1):

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

536 try:

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

538 except Exception:

539 pass

540

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

542 # to capture the full X spread at various latitudes

543 for i in range(samples + 1):

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

545 for j in range(samples + 1):

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

547 try:

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

549 except Exception:

550 pass

551

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

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

554

555 if not xs:

556 return

557

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

559

560

561def _is_big_extent(ext_wgs):

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

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

564 return dx > 350 or dy > 160

565

566

567def _transform_extent_direct(ext, srid_from, srid_to):

568 tr = _pyproj_transformer(srid_from, srid_to)

569

570 ext_nor = _normalize_extent(ext)

571

572 res = tr.transform_bounds(

573 left=ext_nor[0],

574 bottom=ext_nor[1],

575 right=ext_nor[2],

576 top=ext_nor[3],

577 errcheck=True,

578 )

579 return _normalize_extent(res)

580

581

582def _normalize_extent(ext):

583 return (

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

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

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

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

588 )

589

590

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

592 crs = Object()

593

594 crs.srid = srid

595

596 with warnings.catch_warnings():

597 warnings.simplefilter('ignore')

598 try:

599 crs.proj4text = pp.to_proj4()

600 except pyproj.exceptions.CRSError:

601 gws.log.error(f'CRS: cannot convert {srid!r} to proj4')

602 return None

603

604 crs.wkt = pp.to_wkt()

605

606 crs.axis = axis

607 crs.uom = uom

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

609

610 crs.isGeographic = pp.is_geographic

611 crs.isProjected = pp.is_projected

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

613

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

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

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

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

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

619

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

621 d = pp.to_json_dict()

622

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

624

625 def _datum(x):

626 if 'datum_ensemble' in x:

627 return x['datum_ensemble']['name']

628 if 'datum' in x:

629 return x['datum']['name']

630 return ''

631

632 def _bbox(d):

633 b = d.get('bbox')

634 if b:

635 # pyproj 3.6

636 return b

637 if d.get('usages'):

638 # pyproj 3.7

639 for u in d['usages']:

640 b = u.get('bbox')

641 if b:

642 return b

643

644 b = d.get('base_crs')

645 if b:

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

647 crs.datum = _datum(b)

648 else:

649 crs.base = 0

650 crs.datum = _datum(d)

651

652 b = _bbox(d)

653 if not b:

654 _warnings[srid] = f'CRS: no bbox for {crs.srid!r}'

655 return

656

657 crs.wgsExtent = (

658 b['west_longitude'],

659 b['south_latitude'],

660 b['east_longitude'],

661 b['north_latitude'],

662 )

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

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

665

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

667

668 return crs

669

670

671def _wgs_max_extent(pp, wgs_extent):

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

673

674 geo = pp.geodetic_crs

675 au = geo.area_of_use if geo else None

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

677 return _normalize_extent(au.bounds)

678

679 x0, y0, x1, y1 = wgs_extent

680 w = x1 - x0

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

682

683

684_AXES_AND_UNITS = {

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

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

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

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

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

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

691}

692

693

694def _axis_and_unit(pp):

695 ax = []

696 for a in pp.axis_info:

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

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

699

700

701##

702

703"""

704Projections can be referenced by:

705

706 - int/numeric SRID: 4326

707 - EPSG Code: EPSG:4326

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

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

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

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

712

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

714"""

715

716_WRITE_FORMATS = {

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

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

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

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

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

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

723}

724

725_PARSE_FORMATS = {

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

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

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

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

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

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

732}

733

734# @TODO

735

736_aliases = {

737 'crs:84': 4326,

738 'crs84': 4326,

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

740 'wgs84': 4326,

741 'epsg:900913': 3857,

742 'epsg:102100': 3857,

743 'epsg:102113': 3857,

744}

745

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

747# EPSG:4326 longitude/latitude

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

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

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

751

752_AXIS_FOR_FORMAT = {

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

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

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

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

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

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

759}

760

761

762def _parse(crs_name):

763 if isinstance(crs_name, int):

764 return gws.CrsFormat.epsg, crs_name

765

766 if isinstance(crs_name, bytes):

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

768

769 if isinstance(crs_name, str):

770 crs_name = crs_name.lower()

771

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

773 return gws.CrsFormat.crs, 4326

774

775 if crs_name in _aliases:

776 crs_name = _aliases[crs_name]

777 if isinstance(crs_name, int):

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

779

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

781 m = re.match(r, crs_name)

782 if m:

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

784

785 return None, 0

786

787

788def _unparse(srid, fmt):

789 return _WRITE_FORMATS[fmt].format(srid)

792##

795_obj_cache: dict = {

796 WGS84.srid: WGS84,

797 WEBMERCATOR.url: WEBMERCATOR,

798}

799

800_pyproj_cache: dict = {}

801

802_transformer_cache: dict = {}

803

804_warnings: dict = {}

Coverage for gws-app/gws/lib/crs/init.py: 92%

359 statements