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

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(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 return getattr(self, str(fmt or gws.CrsFormat.epsg).lower())

95 def to_geojson(self):

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

97 return {

98 'type': 'name',

99 'properties': {

100 'name': self.urn,

101 },

102 }

105#

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

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

110 x0, y0, x1, y1 = extent

111

112 lat = (y0 + y1) * 0.5

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

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

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

116 RA = 6378000

117 E = 0.0810820288

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

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

120

121

122#

123

124# enough precision to represent 1cm

125COORDINATE_PRECISION_DEG = 7

126COORDINATE_PRECISION_M = 2

127

128WGS84 = Object(

129 srid=4326,

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

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

132 axis=gws.Axis.yx,

133 uom=gws.Uom.deg,

134 isGeographic=True,

135 isProjected=False,

136 isYX=True,

137 epsg='EPSG:4326',

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

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

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

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

142 name='WGS 84',

143 base=0,

144 datum='World Geodetic System 1984 ensemble',

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

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

147 coordinatePrecision=COORDINATE_PRECISION_DEG,

148)

149

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

151

152WEBMERCATOR = Object(

153 srid=3857,

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

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

156 axis=gws.Axis.xy,

157 uom=gws.Uom.m,

158 isGeographic=False,

159 isProjected=True,

160 isYX=False,

161 epsg='EPSG:3857',

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

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

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

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

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

167 base=4326,

168 datum='World Geodetic System 1984 ensemble',

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

170 extent=(

171 -20037508.342789244,

172 -20048966.104014598,

173 20037508.342789244,

174 20048966.104014598,

175 ),

176 coordinatePrecision=COORDINATE_PRECISION_M,

177)

178

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

180

181WEBMERCATOR_RADIUS = 6378137

182WEBMERCATOR_SQUARE = (

183 -math.pi * WEBMERCATOR_RADIUS,

184 -math.pi * WEBMERCATOR_RADIUS,

185 +math.pi * WEBMERCATOR_RADIUS,

186 +math.pi * WEBMERCATOR_RADIUS,

187)

188

189

190class Error(gws.Error):

191 pass

192

193

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

195 """Returns the CRS for a given CRS-code or SRID."""

196 if not crs_name:

197 return None

198 return _get_crs(crs_name)

199

200

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

202 """Parses a CRS to a tuple of CRS-format and the CRS itself."""

203 fmt, srid = _parse(crs_name)

204 if not fmt:

205 return gws.CrsFormat.none, None

206 return fmt, _get_crs(srid)

207

208

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

210 """Raises an error if no correct CRS is given."""

211 crs = _get_crs(crs_name)

212 if not crs:

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

214 return crs

217##

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

221 """Return a crs from the list that most closely matches the given crs.

222

223 Args:

224 crs: target CRS

225 supported_crs: CRS list

226

227 Returns:

228 A CRS object

229 """

230

231 if crs in supported_crs:

232 return crs

233

234 bst = _best_match(crs, supported_crs)

235 if not bst:

236 bst = supported_crs[0] if supported_crs else crs

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

238 return bst

239

240

241def _best_match(crs, supported_crs):

242 # @TODO find a projection with less errors

243 # @TODO find a projection with same units

244

245 if crs.isProjected:

246 # for a projected crs, find webmercator

247 for sup in supported_crs:

248 if sup.srid == WEBMERCATOR.srid:

249 return sup

250

251 # not found, return the first projected crs

252 for sup in supported_crs:

253 if sup.isProjected:

254 return sup

255

256 if crs.isGeographic:

257 # for a geographic crs, try wgs first

258 for sup in supported_crs:

259 if sup.srid == WGS84.srid:

260 return sup

261

262 # not found, return the first geographic crs

263 for sup in supported_crs:

264 if sup.isGeographic:

265 return sup

268##

271def _get_crs(crs_name):

272 if crs_name in _obj_cache:

273 return _obj_cache[crs_name]

274

275 fmt, srid = _parse(crs_name)

276 if not fmt:

277 _obj_cache[crs_name] = None

278 return None

279

280 if srid in _obj_cache:

281 _obj_cache[crs_name] = _obj_cache[srid]

282 return _obj_cache[srid]

283

284 obj = _get_new_crs(srid)

285 _obj_cache[crs_name] = _obj_cache[srid] = obj

286 return obj

287

288

289def _get_new_crs(srid):

290 pp = _pyproj_crs_object(srid)

291 if not pp:

292 gws.log.error(f'CRS: unknown srid {srid!r}')

293 return None

294

295 au = _axis_and_unit(pp)

296 if not au:

297 gws.log.error(f'CRS: unsupported srid {srid!r}')

298 return None

299

300 axis, uom = au

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

302 gws.log.error(f'CRS: unsupported unit {uom!r} for {srid!r}')

303 return None

304

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

306

307

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

309 if srid in _pyproj_cache:

310 return _pyproj_cache[srid]

311

312 try:

313 pp = pyproj.CRS.from_epsg(srid)

314 except pyproj.exceptions.CRSError:

315 return None

316

317 _pyproj_cache[srid] = pp

318 return _pyproj_cache[srid]

319

320

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

322 key = srid_from, srid_to

323

324 if key in _transformer_cache:

325 return _transformer_cache[key]

326

327 pa = _pyproj_crs_object(srid_from)

328 pb = _pyproj_crs_object(srid_to)

329

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

331 return _transformer_cache[key]

332

333

334def _transform_extent(ext, srid_from, srid_to):

335 tr = _pyproj_transformer(srid_from, srid_to)

336

337 e = tr.transform_bounds(

338 left=min(ext[0], ext[2]),

339 bottom=min(ext[1], ext[3]),

340 right=max(ext[0], ext[2]),

341 top=max(ext[1], ext[3]),

342 )

343

344 return (

345 min(e[0], e[2]),

346 min(e[1], e[3]),

347 max(e[0], e[2]),

348 max(e[1], e[3]),

349 )

350

351

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

353 crs = Object()

354

355 crs.srid = srid

356

357 with warnings.catch_warnings():

358 warnings.simplefilter('ignore')

359 try:

360 crs.proj4text = pp.to_proj4()

361 except pyproj.exceptions.CRSError:

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

363 return None

364

365 crs.wkt = pp.to_wkt()

366

367 crs.axis = axis

368 crs.uom = uom

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

370

371 crs.isGeographic = pp.is_geographic

372 crs.isProjected = pp.is_projected

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

374

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

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

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

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

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

380

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

382 d = pp.to_json_dict()

383

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

385

386 def _datum(x):

387 if 'datum_ensemble' in x:

388 return x['datum_ensemble']['name']

389 if 'datum' in x:

390 return x['datum']['name']

391 return ''

392

393 def _bbox(d):

394 b = d.get('bbox')

395 if b:

396 # pyproj 3.6

397 return b

398 if d.get('usages'):

399 # pyproj 3.7

400 for u in d['usages']:

401 b = u.get('bbox')

402 if b:

403 return b

404

405 b = d.get('base_crs')

406 if b:

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

408 crs.datum = _datum(b)

409 else:

410 crs.base = 0

411 crs.datum = _datum(d)

412

413 b = _bbox(d)

414 if not b:

415 gws.log.error(f'CRS: no bbox for {crs.srid!r}')

416 return

417

418 crs.wgsExtent = (

419 b['west_longitude'],

420 b['south_latitude'],

421 b['east_longitude'],

422 b['north_latitude'],

423 )

424 crs.extent = _transform_extent(crs.wgsExtent, WGS84.srid, srid)

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

426

427 return crs

428

429

430_AXES_AND_UNITS = {

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

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

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

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

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

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

437}

438

439

440def _axis_and_unit(pp):

441 ax = []

442 for a in pp.axis_info:

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

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

445

446

447##

448

449"""

450Projections can be referenced by:

451

452 - int/numeric SRID: 4326

453 - EPSG Code: EPSG:4326

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

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

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

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

458

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

460"""

461

462_WRITE_FORMATS = {

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

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

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

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

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

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

469}

470

471_PARSE_FORMATS = {

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

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

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

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

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

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

478}

479

480# @TODO

481

482_aliases = {

483 'crs:84': 4326,

484 'crs84': 4326,

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

486 'wgs84': 4326,

487 'epsg:900913': 3857,

488 'epsg:102100': 3857,

489 'epsg:102113': 3857,

490}

491

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

493# EPSG:4326 longitude/latitude

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

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

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

497

498_AXIS_FOR_FORMAT = {

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

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

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

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

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

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

505}

506

507

508def _parse(crs_name):

509 if isinstance(crs_name, int):

510 return gws.CrsFormat.epsg, crs_name

511

512 if isinstance(crs_name, bytes):

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

514

515 if isinstance(crs_name, str):

516 crs_name = crs_name.lower()

517

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

519 return gws.CrsFormat.crs, 4326

520

521 if crs_name in _aliases:

522 crs_name = _aliases[crs_name]

523 if isinstance(crs_name, int):

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

525

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

527 m = re.match(r, crs_name)

528 if m:

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

530

531 return None, 0

532

533

534def _unparse(srid, fmt):

535 return _WRITE_FORMATS[fmt].format(srid)

538##

541_obj_cache: dict = {

542 WGS84.srid: WGS84,

543 WEBMERCATOR.url: WEBMERCATOR,

544}

545

546_pyproj_cache: dict = {}

547

548_transformer_cache: dict = {}

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

254 statements