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.28'

37_bounds__ = ' bounds '

38_forward_ = 'forward'

39_outside__ = 'outside '

40_region4hBG = Bounds4Tuple(48.5, 1.0, 52.5, 7.0, name='hBG 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"

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

52 _raiser = False

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

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

57 @kwarg datum: Conic's datum to use (C{pygeodesy.Datums}, 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}='Be*LBG'} (C{str}).

61 '''

62 if name: # typename(self)

63 self.name = name

64 if datum:

65 c = self._conic.toDatum(datum)

66 if self._conic is not c:

67 self._conic = c

68# E = self.datum.ellipsoid

69# if not E.isOblate:

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

71 if raiser:

72 self.raiser = True

73# T = self.datum.transform

74# if not T.isunity:

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

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: Method C{isinside} of L{bounds4<_BeLBGbase.bounds4>},

268 L{region4<_BeLBGbase.region4>} and L{bounds4(asLb)

269 <Lb4Tuple.isinside>}, L{region4(asLB)<Lb4Tuple.isinside>}

270 '''

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

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

273 _region4hBG)

274

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

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

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

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

279 _NAN = True

280 elif t < self._bounds2:

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

282 else:

283 _NAN = False

284 return lat, lon, _NAN, raiser, name

285

286 def _outsidError(self, *lat_lon):

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

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

289

290 @property

291 def raiser(self):

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

293 '''

294 return self._raiser

295

296 @raiser.setter # PYCHOK setter!

297 def raiser(self, raiser):

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

299 '''

300 self._raiser = bool(raiser)

301

302 def region4(self, asLb=False):

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

304

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

306 easting and northing, otherwise in geodetic lat- and

307 longitudes (C{bool}).

308

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

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

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

312

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

314 southern half of the Netherlands.

315 '''

316 return self._region4Lb if asLb else _region4hBG

317

318 @property_RO

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

320 n = _region4hBG.name.split()

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

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

323 self.__class__._region4Lb = r

324 return r

325

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

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

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

329

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

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

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

333 to ignore C{hBGh} interpolation.

334

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

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

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

338

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

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

341

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

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

344

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

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

347 '''

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

349 if _NAN:

350 lat = lon = height = NAN

351 else:

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

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

354

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

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

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

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

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

360 return lat, lon, height

361

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

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

364

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

366

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

368 '''

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

370

371 @property_RO

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

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

374 '''

375 lat, lon, H = self._Uccle3

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

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

378 self.__class__.Uccle = u

379 return u

380

381 @property_ROver

382 def _Uccle3(self):

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

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

385

386

387class Be08LBG(_BeLBGbase):

388 '''Belgian Lambert 2008 C{pygeodesy.Conics.Be08Lb} converter.

389 '''

390 _conic = 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} converter.

399 '''

400 _conic = Conics.Be72Lb

401

402 @property_ROver

403 def _bounds2(self):

404 return _bounds2(49.30, 2.31, self)

405

406

407class Be72NLBG(_BeLBGbase):

408 '''Belgian Lambert 1972N C{pygeodesy.Conics.Be72NLb} converter.

409 '''

410 _conic = Conics.Be72NLb

411

412 @property_ROver

413 def _bounds2(self):

414 return _bounds2(49.31, 2.15, self)

415

416

417class Be72RLBG(_BeLBGbase):

418 '''Belgian Lambert 1972R C{pygeodesy.Conics.Be72RLb} converter.

419 '''

420 _conic = Conics.Be72RLb

421

422 @property_ROver

423 def _bounds2(self):

424 return _bounds2(49.21, 2.14, self)

425

426

427class Be50LBG(_BeLBGbase):

428 '''Belgian Lambert 1950 C{pygeodesy.Conics.Be50Lb} converter.

429 '''

430 _conic = Conics.Be50Lb

431

432 @property_ROver

433 def _bounds2(self):

434 return _bounds2(49.31, 5.26, self)

435

436

437class _Bounds2Tuple(_NamedTuple):

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

439 '''

440 _Names_ = Bounds4Tuple._Names_[:2]

441 _Units_ = Bounds4Tuple._Units_[:2]

442

443

444class _Lb2Tuple(_NamedTuple):

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

446 '''

447 _Names_ = Lb4Tuple._Names_[:2]

448 _Units_ = Lb4Tuple._Units_[:2]

449

450

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

452 c_lonI, f_lonI, lonN_f):

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

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

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

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

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

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

459

460

461def _bounds2(latS, lonW, beLBG):

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

463 n = typename(beLBG) + _bounds__

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

465

466

467def _c_f_N_f3(*deg_SWD):

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

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

470 N = _degN(*deg_SWD)

471 # assert N >= 0, N

472 f = floor(N)

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

474

475

476def _degN(deg, degSW, degD):

477 # return C{deg} Normalized

478 return (deg - degSW) / degD

479

480

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

482

483if _FOR_DOCS:

484 # for epydoc to include __doc__ for all classes

485 for B in _Be5LBGs:

486 B.bounds4 = _BeLBGbase.bounds4

487 B.conic = _BeLBGbase.conic

488 B.datum = _BeLBGbase.datum

489 B.forward = _BeLBGbase.forward

490 B.hBGh = _BeLBGbase.hBGh

491 B.hBGh3 = _BeLBGbase.hBGh3

492 B.isinside = _BeLBGbase.isinside

493 B.region4 = _BeLBGbase.region4

494 B.reverse = _BeLBGbase.reverse

495

496__all__ += _ALL_DOCS(_BeLBGbase)

497__all__ += _all_OTHER(Conics, *_Be5LBGs)

498del _ALL_DOCS, _all_OTHER, _Be5LBGs

499

500# **) MIT License

501#

503#

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

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

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

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

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

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

510#

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

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

513#

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

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

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

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

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

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

520# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pybelbg / belbgs.py: 96%

214 statements