Coverage for pybelbg/belbgs.py: 96%

2# -*- coding: utf-8 -*-

4u'''Main classes L{Be08LBG}, L{Be72LBG}, L{Be72NLBG}, L{Be72RLBG} and L{Be50LBG} implementing

5the C{Belgian Lambert 2018}, C{-1972}, C{-1972N}, C{-1972R} respectively C{-1950} conic projection

6and quasi-geoid heights with bilinear interpolation from C{Belgian hybrid quasi-geoid} U{hBG18

7<https://DOI.org/10.5880/isg.2018.003>}.

9Each of the 5 classes provides a C{forward} method to transform a geodetic lat-, longitude and

10(ellipsoidal) height to local easting, northing and (orthometric) height and a C{reverse} method

11for converting local to geodetic coordinates and (orthometric to ellipsoidal) height.

13All classes use the same hybrid quasi-geoid C{hBG18} C{region4} but slighly different C{bounds4}

14for valid lat- and longitudes. Heights for points outside C{region4} are not interpolated and

15are C{NAN} or throw a L{BeLBGError} exception. Likewise, lat-, longitudes I{below} C{bounds4}

16and easting, northing I{below} C{bounds4(asLb)} are considered invalid and replaced with C{NAN}

17or raise a L{BeLBGError}.

18'''

19# make sure int/int division yields float quotient in Py2-

20from __future__ import division as _; del _ # noqa: E702 ;

22from pybelbg.__pygeodesy import (BeLBGError, BeLBG7Tuple, LatLonN3Tuple, Lb4Tuple,

23 _1_0, _3600_0, _isNAN, _name_,

24 _ALL_DOCS, _all_OTHER, _FOR_DOCS,

25 _isinside, _NamedBase, _NamedTuple)

26from pybelbg.__pygeodesy import _COMMA_, _SPACE_ # PYCHOK used!

27from pygeodesy import (typename, NAN, NN, # "consterns"

28 Conics, Bounds4Tuple, # lcc, namedTuples

29 property_RO, property_ROver, # props

30 Degrees, Easting, Height, Lat, Lon, Northing) # units

32from math import ceil, floor

34__all__ = ()

35__version__ = '26.07.31'

37_bounds__ = ' bounds '

38_forward_ = 'forward'

39_outside__ = 'outside '

40_region4hBG = Bounds4Tuple(48.5, 1.0, 52.5, 7.0, name='hBG18 region ') # <https://EPSG.io/4937>

41_reverse_ = 'reverse'

42_Uccle_ = 'Uccle' # PYCHOK == _BelBGbase.Uccle.name

45class _BeLBGbase(_NamedBase):

46 '''(INTERNAL) C{Be*LBG} base class.

47 '''

48 _bounds2 = None # overloaded

49 _conic = None # overloaded

50 _latD = Lat(latD=36 / _3600_0) # 36" in degrees

51 _lonD = Lon(lonD=54 / _3600_0) # 54" in degrees

52 _raiser = False

54 def __init__(self, datum=None, raiser=False, **name):

55 '''New C{Be*LBG} transformer instance, optionally with a different C{conic}'s datum.

57 @kwarg datum: Conic's datum to use (C{pygeodesy.Datums}, I{ellipsoidal} only).

58 @kwarg raiser: If C{True} raise a L{BeLBGError} for lat- or longitudes

59 outside L{region4} or below L{bounds4} (C{bool}).

60 @kwarg name: Optional name C{B{name}=NN} (C{str}).

61 '''

62 if datum:

63 c = self._conic.toDatum(datum)

64 if self._conic is not c:

65 self._conic = c

66# E = self.datum.ellipsoid

67# if not E.isOblate:

68# raise BeLBGError(repr(E), txt='not oblate')

69 if raiser:

70 self.raiser = True

71# T = self.datum.transform

72# if not T.isunity:

73# raise BeLBGError(repr(T), txt='not unity')

74 if name:

75 self.name = name # or typename(self)

77 def _as4Lb(self, t4, name=NN):

78 # return C{t4} as L{Lb4Tuple}

79 S, W, N, E = t4

80 s, w, _ = self._forward3(False, S, W, None)

81 n, e, _ = self._forward3(False, N, E, None)

82 return Lb4Tuple(s, w, n, e, name=name or t4.name)

84 def _belowError(self, raiser, coords, bounds2):

85 # throw a below bounds2 exception if requested

86 if raiser or (raiser is None and self.raiser):

87 raise BeLBGError(coords, txt=_outside__ + bounds2.toRepr())

88 return True

90 @property_RO

91 def _bounds2Lb(self): # overwrite __class__._bounds2Lb

92 # get lower-left of C{bounds4} as C{_Lb2Tuple}

93 b = self.bounds4(True)

94 b = _Lb2Tuple(b.minE, b.minN, name=b.name)

95 self.__class__._bounds2Lb = b

96 return b

98 def bounds4(self, asLb=False):

99 '''Get the South, West, North and East bounds of this C{Lambert} conic projection.

100

101 @kwarg asLb: Use C{B{asLb}=True} for the bounds in local C{Lambert} easting and

102 northing, otherwise in geodetic lat- and longitudes (C{bool}).

103

104 @return: A L{Bounds4Tuple}C{(latS, lonW, latN, lonE)} with (WGS84) geodetic lat-

105 and longitudes in C{degrees} or an L{Lb4Tuple}C{(minE, minN, maxE, maxN)}

106 in C{meter}.

107

108 @note: The C{bounds4} cover Belgium, Belgium's U{EEZ<http://MarineRegions.org/mrgid/3293>},

109 the Netherlands, the Netherlands' U{EEZ<http://MarineRegions.org/mrgid/5668>}

110 and Luxemburg.

111 '''

112 return self._bounds4Lb if asLb else self._bounds4

113

114 @property_ROver

115 def _bounds4BeNeLux(self):

116 # BeNeLux covering BE, BE EEZ, NL, NL EEZ and LU

117 def _b4(swne): # <http://MarineRegions.org/mrgid/XXXX>

118 return Bounds4Tuple(*swne.split(_COMMA_)).toUnits()

119

120 b = _b4('49 26 50.6N, 5 44 5.6E, 50 10 59.9N, 6 31 49.2E') # Lux <2233>

121 b = _b4('50 45 5.8N, 3 21 31.3E, 53 33 38.5N, 7 13 37.9E').union(b) # Ne <15>

122 b = _b4('51 19 48.6N, 2 32 21.6E, 55 45 54.0N, 7 12 37.0E').union(b) # Ne EEZ <5668>

123 b = _b4('49 29 50.3N, 2 32 47.8E, 51 30 15.1N, 6 24 27.0E').union(b) # Be <14>

124 b = _b4('51 2 23.1N, 2 14 18.0E, 51 52 34.0N, 4 24 7.3E').union(b) # Be EEZ <3293>

125 # BeNeLux bounds (latS=49.447389, lonW=2.238333, latN=55.765, lonE=7.227194)

126 return b.toUnits(name='BeNeLux' + _bounds__)

127

128 @property_RO

129 def _bounds4(self): # overwrite class._bounds4

130 # get C{bounds4} with lower-left adjusted upward

131 s, w = self._bounds2

132 S, W, N, E = self._bounds4BeNeLux

133 n = typename(self) + _bounds__

134 b = Bounds4Tuple(max(s, S), max(w, W), N, E, name=n)

135 self.__class__._bounds4 = b

136 return b

137

138 @property_RO

139 def _bounds4Lb(self): # overwrite class._bounds4Lb

140 # get C{bounds4} as C{Lb4Tuple}

141 b = self._as4Lb(self._bounds4)

142 self.__class__._bounds4Lb = b

143 return b

144

145 def _c_f_N_f6(self, lat, lon):

146 # return (int(ceil), int(floor), Normalized less floor) of C{lat}) + \

147 # (int(ceil), int(floor), Normalized less floor) of C{lon})

148 S, W, _, _ = _region4hBG

149 return _c_f_N_f3(lat, S, self._latD) + \

150 _c_f_N_f3(lon, W, self._lonD)

151

152 @property_RO

153 def conic(self):

154 '''Get the C{Lambert} conic (C{pygeodesy.Conic}).

155 '''

156 return self._conic

157

158 @property_RO

159 def datum(self):

160 '''Get the C{Lambert} conic's datum (C{pygeodesy.Datum}).

161 '''

162 return self.conic.datum

163

164 def _EasNor5(self, e, n, raiser=None, name=_reverse_):

165 # return e, n, ... if non-NAN and not below bounds2Lb

166 e, n = t = Easting(e), Northing(n)

167 if _isNAN(e) or _isNAN(n):

168 _NAN = True

169 elif t < self._bounds2Lb:

170 _NAN = self._belowError(raiser, t, self._bounds2Lb)

171 else:

172 _NAN = False

173 return e, n, _NAN, raiser, name

174

175 def forward(self, lat, lon, height=0, **raiser_name):

176 '''Convert geodetic C{(B{lat}, B{lon})} and (ellipsoidal) B{C{height}}

177 to C{easting}, C{northing} and (orthometric) height C{H}.

178

179 @arg lat: Latitude (C{degrees}, geodetic).

180 @arg lon: Longitude (C{degrees}, geodetic).

181 @kwarg height: The (ellipsoidal) height (C{meter}, conventionally) or

182 C{None} to ignore C{hBGh} interpolation.

183

184 @return: A L{BeLBG7Tuple}C{(easting, northing, H, lat, lon, height, beLBG)}

185 with C{easting}, C{northing} and (orthometric) height C{H} in

186 C{meter} or C{NAN} and C{beLBG} is this C{Be*LBG} instance.

187

188 @raise BeLBGError: If the point is outside the C{BG} region and property

189 C{raiser is True} or keyword argument C{B{raiser}=True}.

190

191 @note: C{H}, C{easting} and C{northing} will all be C{NAN} if B{C{lat}} or

192 B{C{lon}} is below this converter's L{bounds4}.

193

194 @note: Orthometric height C{(H = h - N)} euals ellipsoidal height C{h}

195 less (hybrid quasi-) geoid height C{N}.

196 '''

197 lat, lon, _NAN, raiser, name = self._LatLon5(lat, lon, **raiser_name)

198 if _NAN:

199 e = n = H = NAN

200 else:

201 e, n, H = self._forward3(raiser, lat, lon, height)

202 return BeLBG7Tuple(e, n, H, lat, lon, height, self, name=name)

203

204 def _forward3(self, raiser, lat, lon, height): # in .__main__

205 # C{forward} core, returning C{(easting, northing, H)}

206 H = NAN if height is None or _isNAN(height) else (

207 Height(height) - self._hBGh(lat, lon, raiser))

208 e, n, _ = self._conic.forward3(lat, lon)

209 return e, n, H

210

211 def hBGh(self, lat, lon):

212 '''Interpolate the hybrid quasi-geoid C{hBG} height for a geodetic point.

213

214 @arg lat: Latitude (C{degrees}, geodetic).

215 @arg lon: Longitude (C{degrees}, geodetic).

216

217 @return: Hybrid quasi-geoid C{hBG} height C{N} (C{meter}) or C{NAN}

218 if B{C{lat}} or B{C{lon}} is outside L{region4}.

219 '''

220 lat, lon, _NAN, _, _ = self._LatLon5(lat, lon, False)

221 return NAN if _NAN else self._hBGh(lat, lon)

222

223 def _hBGh(self, lat, lon, raiser=False):

224 # interpolate C{N} at C{(lat, lon)} or C{NAN} if

225 # outside or ... if _isNAN(lat) or _isNAN(lon)

226 if _isinside(lat, lon, 0, _region4hBG):

227 c_f_N_f6 = self._c_f_N_f6(lat, lon)

228 N = _bilinear(self._hBG18, *c_f_N_f6)

229 N = Height(N=N)

230 elif raiser or (raiser is None and self._raiser):

231 raise self._outsidError(lat, lon)

232 else:

233 N = NAN

234 return N

235

236 def hBGh3(self, easting, northing):

237 '''Interpolate the hybrid quasi-geoid C{hBG} height for a local point.

238

239 @arg easting: Easting (C{meter}, local).

240 @arg northing: Northing (C{meter}, local).

241

242 @return: L{LatLonN3Tuple}C{(lat, lon, N)} with the (hybrid

243 quasi-) geoid C{hBG} height C{N} in C{meter} or

244 C{NAN} if C{lat} or C{lon} is outside C{region4}.

245 '''

246 r = self.reverse(easting, northing, H=0, raiser=False)

247 return LatLonN3Tuple(r.lat, r.lon, r.height, name=self.name)

248

249 @property_ROver

250 def _hBG18(self): # load the hBG18 geoid, I{once}

251 from pybelbg.hBG18 import _hBG18 as hBG

252 S, W, N, E = _region4hBG

253 assert int(_degN(N, S, self._latD) + _1_0) == len(hBG)

254 assert int(_degN(E, W, self._lonD) + _1_0) == len(hBG)

255 return hBG

256

257 def isinside(self, lat, lon, eps=0):

258 '''Is geodetic C{(B{lat}, B{lon})} inside the C{hBG} region?

259

260 @arg lat: Latitude (C{degrees}, geodetic).

261 @arg lon: Longitude (C{degrees}, geodetic).

262 @kwarg eps: Over-/undersize the C{hBG} region (C{degrees}).

263

264 @return: C{None} if B{C{lat}} or B{C{lon}} is NAN, C{False}

265 if outside the C{hBG} region, C{True} otherwise.

266

267 @see: Methods C{Bounds4Tuple.isinside} and L{Lb4Tuple.isinside}.

268 '''

269 lat, lon, _NAN, _, _ = self._LatLon5(lat, lon, False)

270 return None if _NAN else _isinside(lat, lon, Degrees(eps=eps),

271 _region4hBG)

272

273 def _LatLon5(self, lat, lon, raiser=None, name=_forward_):

274 # return lat, lon, ... if non-NAN and not below bounds2

275 lat, lon = t = Lat(lat), Lon(lon)

276 if _isNAN(lat) or _isNAN(lon):

277 _NAN = True

278 elif t < self._bounds2:

279 _NAN = self._belowError(raiser, t, self._bounds2)

280 else:

281 _NAN = False

282 return lat, lon, _NAN, raiser, name

283

284 def _outsidError(self, *lat_lon):

285 # format an BeLBGError for C{lat_lon} outside the C{hBG} region

286 return BeLBGError(lat_lon, txt=_outside__ + _region4hBG.toRepr())

287

288 @property

289 def raiser(self):

290 '''Is an C{BeLBGError} thrown for points outside the C{hBG} region?

291 '''

292 return self._raiser

293

294 @raiser.setter # PYCHOK setter!

295 def raiser(self, raiser):

296 '''Use C{True} to throw an C{BeLBGError} for points outside the C{hBG} region.

297 '''

298 self._raiser = bool(raiser)

299

300 def region4(self, asLb=False):

301 '''Get the South, West, North and East bounds of the C{hBG} region.

302

303 @kwarg asLb: Use C{B{asLb}=True} for the bounds in local C{Lambert}

304 easting and northing, otherwise in geodetic lat- and

305 longitudes (C{bool}).

306

307 @return: A L{Bounds4Tuple}C{(latS, lonW, latN, lonE)} with (WGS84)

308 geodetic lat- and longitudes in C{degrees} or an

309 L{Lb4Tuple}C{(minE, minN, maxE, maxN)} in C{meter}.

310

311 @note: The C{hBG} region covers all of Belgium, Luxemburg and the

312 southern half of the Netherlands.

313 '''

314 return self._region4Lb if asLb else _region4hBG

315

316 @property_RO

317 def _region4Lb(self): # overwrite class._region4Lb

318 n = _region4hBG.name.split()

319 n = _SPACE_(typename(self), *n[1:])

320 r = self._as4Lb(_region4hBG, name=n)

321 self.__class__._region4Lb = r

322 return r

323

324 def reverse(self, easting, northing, H=0, **raiser_name):

325 '''Convert local B{C{easting}}, B{C{northing}} and (orthometric) height

326 B{C{H}} to geodetic C{lat-}, C{longitude} and (ellipsoidal) C{height}.

327

328 @arg easting: Easting (C{meter}, local).

329 @arg northing: Northing (C{meter}, local).

330 @kwarg H: The (orthometric) height (C{meter}, conventionally) or C{None}

331 to ignore C{hBGh} interpolation.

332

333 @return: A L{BeLBG7Tuple}C{(easting, northing, H, lat, lon, height, beLBG)}

334 with geodetic C{lat} and C{lon} and (ellipsoidal) C{height} in

335 C{meter} or C{NAN} and C{beLBG} is this C{Be*LBG} instance.

336

337 @raise BeLBGError: If the point is outside the C{hBG} region and property

338 C{raiser is True} or keyword argument C{B{raiser}=True}.

339

340 @note: All C{lon}, C{lat} and C{height} will be C{NAN} if B{C{easting}} or

341 B{C{northing}} is below this converter's C{bounds4(asLb)}.

342

343 @note: Ellipsoidal height C{(h = H + N)} equals orthometric height C{H}

344 plus (hybrid quasi-) geoid height C{N}.

345 '''

346 e, n, _NAN, raiser, name = self._EasNor5(easting, northing, **raiser_name)

347 if _NAN:

348 lat = lon = height = NAN

349 else:

350 lat, lon, height = self._reverse3(raiser, e, n, H)

351 return BeLBG7Tuple(e, n, H, lat, lon, height, self, name=name)

352

353 def _reverse3(self, raiser, e, n, H): # in .__main__

354 # C{reverse} core, returning C{(lat, lon, height)}

355 lat, lon, _, _= self._conic.reverse4(e, n)

356 height = NAN if H is None or _isNAN(H) else (

357 Height(H=H) + self._hBGh(lat, lon, raiser))

358 return lat, lon, height

359

360 def toStr(self, prec=9, **unused): # PYCHOK signature

361 '''Return this C{Be*LBG} instance as a string.

362

363 @kwarg prec: Precision, number of decimal digits (C{int}, 0..9).

364

365 @return: This C{Be*LBG} (C{str}).

366 '''

367 return self.attrs(_name_, 'conic', 'raiser', prec=prec) # _datum_, _Uccle_

368

369 @property_RO

370 def Uccle(self): # overwrite class.Uccle

371 '''Get C{Uccle<https://ROBinfo.OMA.BE/en/astro-info/geographical-coordinates-of-our-sites>} (aka Ukkel) as L{BeLBG7Tuple}.

372 '''

373 lat, lon, H = self._Uccle3

374 h = self.hBGh(lat, lon) + H # height=147.815887 Be08LBG

375 u = self.forward(lat, lon, height=h, name=_Uccle_) # Be08LBG easting=649250.118675, northing=665257.86004

376 self.__class__.Uccle = u

377 return u

378

379 @property_ROver

380 def _Uccle3(self):

381 # # lat=50.797778, lon=4.358111, H=104.9

382 return Lat('50 47 52N'), Lon('4 21 29.2E'), Height(H=104.9)

383

384

385class Be08LBG(_BeLBGbase):

386 '''Belgian Lambert 2008 C{pygeodesy.Conics.Be08Lb} transformer.

387 '''

388 @property_ROver

389 def _conic(self):

390 return Conics.Be08Lb

391

392 @property_ROver

393 def _bounds2(self):

394 return _bounds2(44.77, -3.82, self)

395

396

397class Be72LBG(_BeLBGbase):

398 '''Belgian Lambert 1972 C{pygeodesy.Conics.Be72Lb} transformer.

399 '''

400 @property_ROver

401 def _conic(self):

402 return Conics.Be72Lb

403

404 @property_ROver

405 def _bounds2(self):

406 return _bounds2(49.30, 2.31, self)

407

408

409class Be72NLBG(_BeLBGbase):

410 '''Belgian Lambert 1972N C{pygeodesy.Conics.Be72NLb} transformer.

411 '''

412 @property_ROver

413 def _conic(self):

414 return Conics.Be72NLb

415

416 @property_ROver

417 def _bounds2(self):

418 return _bounds2(49.31, 2.15, self)

419

420

421class Be72RLBG(_BeLBGbase):

422 '''Belgian Lambert 1972R C{pygeodesy.Conics.Be72RLb} transformer.

423 '''

424 @property_ROver

425 def _conic(self):

426 return Conics.Be72RLb

427

428 @property_ROver

429 def _bounds2(self):

430 return _bounds2(49.21, 2.14, self)

431

432

433class Be50LBG(_BeLBGbase):

434 '''Belgian Lambert 1950 C{pygeodesy.Conics.Be50Lb} transformer.

435 '''

436 @property_ROver

437 def _conic(self):

438 return Conics.Be50Lb

439

440 @property_ROver

441 def _bounds2(self):

442 return _bounds2(49.31, 5.26, self)

443

444

445class _Bounds2Tuple(_NamedTuple):

446 '''2-Tuple C{(latS, lonW)} lower-left corner.

447 '''

448 _Names_ = Bounds4Tuple._Names_[:2]

449 _Units_ = Bounds4Tuple._Units_[:2]

450

451

452class _Lb2Tuple(_NamedTuple):

453 '''2-Tuple C{(minE, minN)} lower-left corner.

454 '''

455 _Names_ = Lb4Tuple._Names_[:2]

456 _Units_ = Lb4Tuple._Units_[:2]

457

458

459def _bilinear(hBG, c_latI, f_latI, latN_f,

460 c_lonI, f_lonI, lonN_f):

461 # interpolate hybrid quasi-geoid C{hBG} height

462 ne, nw = hBG(c_latI, c_lonI, f_lonI)

463 se, sw = hBG(f_latI, c_lonI, f_lonI)

464 lonN_f1 = _1_0 - lonN_f # == 1 - (lonN - f_lonN)

465 return (ne * lonN_f + nw * lonN_f1) * latN_f + \

466 (se * lonN_f + sw * lonN_f1) * (_1_0 - latN_f)

467

468

469def _bounds2(latS, lonW, beLBG):

470 # return the C{beLBG} lower-left bounds

471 n = typename(beLBG) + _bounds__

472 return _Bounds2Tuple(latS, lonW, name=n)

473

474

475def _c_f_N_f3(*deg_SWD):

476 # return int(ceil) and int(floor) of Normalized

477 # and (Normalized less floor) of C{deg} degrees

478 N = _degN(*deg_SWD)

479 # assert N >= 0, N

480 f = floor(N)

481 return int(ceil(N)), int(f), (N - f)

482

483

484def _degN(deg, degSW, degD):

485 # return C{deg} Normalized

486 return (deg - degSW) / degD

487

488

489_Be5LBGs = Be08LBG, Be72LBG, Be72NLBG, Be72RLBG, Be50LBG

490

491if _FOR_DOCS: # force epydoc to document all ...

492 for B in _Be5LBGs: # ... public methods

493 B.__init__ = _BeLBGbase.__init__

494 B.bounds4 = _BeLBGbase.bounds4

495 B.conic = _BeLBGbase.conic

496 B.datum = _BeLBGbase.datum

497 B.forward = _BeLBGbase.forward

498 B.hBGh = _BeLBGbase.hBGh

499 B.hBGh3 = _BeLBGbase.hBGh3

500 B.isinside = _BeLBGbase.isinside

501 B.region4 = _BeLBGbase.region4

502 B.reverse = _BeLBGbase.reverse

503

504__all__ += _ALL_DOCS(_BeLBGbase)

505__all__ += _all_OTHER(Conics, *_Be5LBGs)

506del _ALL_DOCS, _all_OTHER, _Be5LBGs

507

508# **) MIT License

509#

511#

512# Permission is hereby granted, free of charge, to any person obtaining a

513# copy of this software and associated documentation files (the "Software"),

514# to deal in the Software without restriction, including without limitation

515# the rights to use, copy, modify, merge, publish, distribute, sublicense,

516# and/or sell copies of the Software, and to permit persons to whom the

517# Software is furnished to do so, subject to the following conditions:

518#

519# The above copyright notice and this permission notice shall be included

520# in all copies or substantial portions of the Software.

521#

522# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

523# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

524# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

525# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

526# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

527# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

528# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pybelbg / belbgs.py: 96%

224 statements