Coverage for gws-app/gws/gis/crs/__init__.py: 25%

1from typing import Optional 

2 

3import math 

4import re 

5import warnings 

6 

7import pyproj.crs 

8import pyproj.exceptions 

9import pyproj.transformer 

10 

11import gws 

12 

13 

14## 

15 

16class Object(gws.Crs): 

17 def __init__(self, **kwargs): 

18 vars(self).update(kwargs) 

19 

20 # crs objects with the same srid must be equal 

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

22 

23 def __hash__(self): 

24 return self.srid 

25 

26 def __eq__(self, other): 

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

28 

29 def __repr__(self): 

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

31 

32 def axis_for_format(self, fmt): 

33 if not self.isYX: 

34 return self.axis 

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

36 

37 def transform_extent(self, ext, crs_to): 

38 if crs_to == self: 

39 return ext 

40 return _transform_extent(ext, self.srid, crs_to.srid) 

41 

42 def transformer(self, crs_to): 

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

44 return tr.transform 

45 

46 def to_string(self, fmt=None): 

47 return getattr(self, str(fmt or gws.CrsFormat.epsg).lower()) 

48 

49 def to_geojson(self): 

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

51 return { 

52 'type': 'name', 

53 'properties': { 

54 'name': self.urn, 

55 } 

56 } 

57 

58 

59# 

60 

61WGS84 = Object( 

62 srid=4326, 

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

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

65 axis=gws.Axis.yx, 

66 uom=gws.Uom.deg, 

67 isGeographic=True, 

68 isProjected=False, 

69 isYX=True, 

70 epsg='EPSG:4326', 

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

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

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

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

75 name='WGS 84', 

76 base=0, 

77 datum='World Geodetic System 1984 ensemble', 

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

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

80) 

81 

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

83 

84WEBMERCATOR = Object( 

85 srid=3857, 

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

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

88 axis=gws.Axis.xy, 

89 uom=gws.Uom.m, 

90 isGeographic=False, 

91 isProjected=True, 

92 isYX=False, 

93 epsg='EPSG:3857', 

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

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

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

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

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

99 base=4326, 

100 datum='World Geodetic System 1984 ensemble', 

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

102 extent=( 

103 -20037508.342789244, 

104 -20048966.104014598, 

105 20037508.342789244, 

106 20048966.104014598, 

107 ) 

108) 

109 

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

111 

112WEBMERCATOR_RADIUS = 6378137 

113WEBMERCATOR_SQUARE = ( 

114 -math.pi * WEBMERCATOR_RADIUS, 

115 -math.pi * WEBMERCATOR_RADIUS, 

116 +math.pi * WEBMERCATOR_RADIUS, 

117 +math.pi * WEBMERCATOR_RADIUS, 

118) 

119 

120 

121class Error(gws.Error): 

122 pass 

123 

124 

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

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

127 if not crs_name: 

128 return None 

129 return _get_crs(crs_name) 

130 

131 

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

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

134 fmt, srid = _parse(crs_name) 

135 if not fmt: 

136 return gws.CrsFormat.none, None 

137 return fmt, _get_crs(srid) 

138 

139 

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

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

142 crs = _get_crs(crs_name) 

143 if not crs: 

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

145 return crs 

146 

147 

148## 

149 

150 

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

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

153 

154 Args: 

155 crs: target CRS 

156 supported_crs: CRS list 

157 

158 Returns: 

159 A CRS object 

160 """ 

161 

162 if crs in supported_crs: 

163 return crs 

164 

165 bst = _best_match(crs, supported_crs) 

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

167 return bst 

168 

169 

170def _best_match(crs, supported_crs): 

171 # @TODO find a projection with less errors 

172 # @TODO find a projection with same units 

173 

174 if crs.isProjected: 

175 # for a projected crs, find webmercator 

176 for sup in supported_crs: 

177 if sup.srid == WEBMERCATOR.srid: 

178 return sup 

179 

180 # not found, return the first projected crs 

181 for sup in supported_crs: 

182 if sup.isProjected: 

183 return sup 

184 

185 if crs.isGeographic: 

186 

187 # for a geographic crs, try wgs first 

188 for sup in supported_crs: 

189 if sup.srid == WGS84.srid: 

190 return sup 

191 

192 # not found, return the first geographic crs 

193 for sup in supported_crs: 

194 if sup.isGeographic: 

195 return sup 

196 

197 # should never be here, but... 

198 

199 for sup in supported_crs: 

200 return sup 

201 

202 

203def best_bounds(crs: gws.Crs, supported_bounds: list[gws.Bounds]) -> gws.Bounds: 

204 """Return the best one from the list of supported bounds. 

205 

206 Args: 

207 crs: target CRS 

208 supported_bounds: Bounds list 

209 

210 Returns: 

211 A Bounds object 

212 """ 

213 

214 bst = best_match(crs, [b.crs for b in supported_bounds]) 

215 for b in supported_bounds: 

216 if b.crs == bst: 

217 return b 

218 

219 

220def best_axis( 

221 crs: gws.Crs, 

222 protocol: gws.OwsProtocol = None, 

223 protocol_version: str = None, 

224 crs_format: gws.CrsFormat = None, 

225 inverted_crs: Optional[list[gws.Crs]] = None 

226) -> gws.Axis: 

227 """Return the 'best guess' axis under given circumstances. 

228 

229 Args: 

230 crs: target CRS 

231 protocol: OWS protocol (WMS, WFS...) 

232 protocol_version: protocol version 

233 crs_format: the format the target_crs was obtained from 

234 inverted_crs: user-provided list of CRSes known to have an inverted (YX) axis 

235 

236 Returns: 

237 An axis 

238 """ 

239 

240 # @TODO some logic to guess the axis, based on crs, service protocol and version 

241 # see https://docs.geoserver.org/latest/en/user/services/wfs/basics.html#wfs-basics-axis 

242 

243 if inverted_crs and crs in inverted_crs: 

244 return gws.Axis.yx 

245 

246 return gws.Axis.xy 

247 

248 

249## 

250 

251def _get_crs(crs_name): 

252 if crs_name in _obj_cache: 

253 return _obj_cache[crs_name] 

254 

255 fmt, srid = _parse(crs_name) 

256 if not fmt: 

257 _obj_cache[crs_name] = None 

258 return None 

259 

260 if srid in _obj_cache: 

261 _obj_cache[crs_name] = _obj_cache[srid] 

262 return _obj_cache[srid] 

263 

264 pp = _pyproj_crs_object(srid) 

265 if not pp: 

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

267 _obj_cache[crs_name] = _obj_cache[srid] = None 

268 return None 

269 

270 au = _axis_and_unit(pp) 

271 if not au: 

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

273 _obj_cache[crs_name] = _obj_cache[srid] = None 

274 return None 

275 

276 _obj_cache[crs_name] = _obj_cache[srid] = _make_crs(srid, pp, au) 

277 return _obj_cache[srid] 

278 

279 

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

281 if srid in _pyproj_cache: 

282 return _pyproj_cache[srid] 

283 

284 try: 

285 pp = pyproj.crs.CRS.from_epsg(srid) 

286 except pyproj.exceptions.CRSError as exc: 

287 return None 

288 

289 _pyproj_cache[srid] = pp 

290 return _pyproj_cache[srid] 

291 

292 

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

294 key = srid_from, srid_to 

295 

296 if key in _transformer_cache: 

297 return _transformer_cache[key] 

298 

299 pa = _pyproj_crs_object(srid_from) 

300 pb = _pyproj_crs_object(srid_to) 

301 

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

303 return _transformer_cache[key] 

304 

305 

306def _transform_extent(ext, srid_from, srid_to): 

307 tr = _pyproj_transformer(srid_from, srid_to) 

308 

309 e = tr.transform_bounds( 

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

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

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

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

314 ) 

315 

316 return ( 

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

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

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

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

321 ) 

322 

323 

324def _make_crs(srid, pp, au): 

325 crs = Object() 

326 

327 crs.srid = srid 

328 

329 with warnings.catch_warnings(): 

330 warnings.simplefilter('ignore') 

331 crs.proj4text = pp.to_proj4() 

332 

333 crs.wkt = pp.to_wkt() 

334 

335 crs.axis = au[0] 

336 crs.uom = au[1] 

337 

338 crs.isGeographic = pp.is_geographic 

339 crs.isProjected = pp.is_projected 

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

341 

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

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

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

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

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

347 

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

349 d = pp.to_json_dict() 

350 

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

352 

353 def _datum(x): 

354 if 'datum_ensemble' in x: 

355 return x['datum_ensemble']['name'] 

356 if 'datum' in x: 

357 return x['datum']['name'] 

358 return '' 

359 

360 b = d.get('base_crs') 

361 if b: 

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

363 crs.datum = _datum(b) 

364 else: 

365 crs.base = 0 

366 crs.datum = _datum(d) 

367 

368 crs.wgsExtent = crs.extent = None 

369 

370 b = d.get('bbox') 

371 if b: 

372 crs.wgsExtent = ( 

373 b['west_longitude'], 

374 b['south_latitude'], 

375 b['east_longitude'], 

376 b['north_latitude'], 

377 ) 

378 else: 

379 crs.wgsExtent = WGS84.wgsExtent 

380 

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

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

383 

384 return crs 

385 

386 

387_AXES_AND_UNITS = { 

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

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

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

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

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

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

394} 

395 

396 

397def _axis_and_unit(pp): 

398 ax = [] 

399 for a in pp.axis_info: 

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

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

402 

403 

404## 

405 

406""" 

407Projections can be referenced by: 

408 

409 - int/numeric SRID: 4326 

410 - EPSG Code: EPSG:4326 

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

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

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

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

415 

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

417""" 

418 

419_WRITE_FORMATS = { 

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

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

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

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

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

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

426} 

427 

428_PARSE_FORMATS = { 

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

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

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

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

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

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

435} 

436 

437# @TODO 

438 

439_aliases = { 

440 'crs:84': 4326, 

441 'crs84': 4326, 

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

443 'wgs84': 4326, 

444 'epsg:900913': 3857, 

445 'epsg:102100': 3857, 

446 'epsg:102113': 3857, 

447} 

448 

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

450# EPSG:4326 longitude/latitude 

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

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

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

454 

455_AXIS_FOR_FORMAT = { 

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

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

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

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

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

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

462} 

463 

464 

465def _parse(crs_name): 

466 if isinstance(crs_name, int): 

467 return gws.CrsFormat.epsg, crs_name 

468 

469 if isinstance(crs_name, bytes): 

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

471 

472 if isinstance(crs_name, str): 

473 crs_name = crs_name.lower() 

474 

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

476 return gws.CrsFormat.crs, 4326 

477 

478 if crs_name in _aliases: 

479 crs_name = _aliases[crs_name] 

480 if isinstance(crs_name, int): 

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

482 

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

484 m = re.match(r, crs_name) 

485 if m: 

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

487 

488 return None, 0 

489 

490 

491def _unparse(srid, fmt): 

492 return _WRITE_FORMATS[fmt].format(srid) 

493 

494 

495## 

496 

497 

498_obj_cache: dict = { 

499 WGS84.srid: WGS84, 

500 WEBMERCATOR.url: WEBMERCATOR, 

501} 

502 

503_pyproj_cache: dict = {} 

504 

505_transformer_cache: dict = {}