Coverage for pyrdnap/rdnap2018.py: 94%

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

4u'''Main classes L{RDNAP2018v1} and L{RDNAP2018v2} follow C{variant 1} respectively C{variant

52} of the U{RDNAPTRANS(tm)2018_v220627<https://formulieren.kadaster.nl/aanvragen_rdnaptrans>}

6specification.

8Each class provides a C{forward} method to transform geodetic lat-, longitudes and (ellipsoidal)

9height to local C{RD} coodinates and (orthometric) height and a C{reverse} method for converting

10local to geodetic coordinates and (orthometric to ellipsoidal) height.

12The L{RDNAP2018v1.forward} and L{.reverse<RDNAP2018v1.reverse>} results have been formally

13validated to meet the C{RDNAPTRANS(tm)2018_v220627} requirements. Likewise for the results

14of L{RDNAP2018v2.forward} and L{.reverse<RDNAP2018v2.reverse>}.

15'''

16# make sure int/int division yields float quotient in Py2

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

19from pyrdnap.rd0 import _RD, _RD0 as A0, RDNAP7Tuple

20from pyrdnap.v_grids import _v_grid # _V_grid

21from pyrdnap.__pygeodesy import (_0_0, _0_5, _1_0, _2_0, _xinstanceof,

22 _isNAN, _isNAN0, _earth_datum, _xkwds_pop2,

23 _name_, _ALL_DOCS, _ALL_OTHER, _FOR_DOCS,

24 _NamedBase, RDNAPError)

25from pygeodesy import (map1, EPS0, EPS1, NAN, PI_2, PI, PI2, typename, # basics, "consterns"

26 Bounds4Tuple, LatLonDatum3Tuple, LatLonNgeoid3Tuple, RD4Tuple, # namedTuples

27 deprecated_property_RO, property_RO, property_ROver, # props

28 Degrees, Lamd, Lat, Lon, Meter, Phid, # units

29 sincos2, sincos2d) # utily

31from math import asin, atan, copysign, degrees, exp, \

32 fabs, floor, hypot, radians, sin, sqrt

34__all__ = ()

35__version__ = '26.08.18'

37_forward_ = 'forward'

38_outside__ = 'outside '

39_region4 = _RD._region4

40_reverse_ = 'reverse'

41_TOL_D = 1e-9 # degrees 2.3.3f+

42_TOL_M = 1e-6 # meter

43_TOL_R = radians(_TOL_D) # 2e-11

44_TRIPS = 16 # 5..6 sufficient

47class _RDNAPbase(_NamedBase):

48 '''(INTERNAL) L{RDNAP2018v1}C{/-v2} base class.

49 '''

50 _datum = None # forward, v1 reverse Datum, lazily (GRS80)

51 _EETRS = None # forward, v1 reverse Ellipsoid, lazily

52 _raiser = False

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

55 '''New C{RDNAP2018v1} or C{-v2} instance.

57 @kwarg a_ellipsoid: An ellipsoid (L{Ellipsoid}) or the ellipsoid's equatorial

58 radius (C{scalar}, conventionally in C{meter}), see B{C{f}}

59 or a datum (L{Datum}). Default C{Datums.GRS80} for ETRS89.

60 @kwarg f: The flattening of the ellipsoid (C{scalar}) if B{C{a_ellipsoid}} is

61 specified as C{scalar}, ignored otherwise.

62 @kwarg raiser: If C{True} raise an L{RDNAPError} for lat-/longitudes outside

63 the C{RD} region (C{bool}).

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

66 @raise RDNAPError: Ellipsoid (or datum) is not oblate (i.e. is spherical or

67 prolate) or the datum's C{transform} is not C{unity}.

68 '''

69 if a_ellipsoid is f is None:

70 self._datum = A0.D80 # GRS80 (ETRS89)

71 else:

72 _earth_datum(self, a_ellipsoid, f, **name) # sets self._datum

73 self._EETRS = E = self._datum.ellipsoid

74 if not E.isOblate:

75 raise RDNAPError(repr(E), txt='not oblate')

76 if raiser: # PYCHOK no cover

77 T = self._datum.transform

78 if not T.isunity:

79 raise RDNAPError(repr(T), txt='not unity')

80 self.raiser = True

81 if name:

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

84 def _asRD(self, b4): # bounds of C{r} as C{RD4Tuple}

85 _xinstanceof(Bounds4Tuple, b4=b4)

86 S, W, N, E = b4

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

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

89 # assert b.x < t.x and b.y < t.y

90 return RD4Tuple(s, w, n, e, name=b4.name)

92 def bounds4(self, asRD=False):

93 '''Get the South, West, North and East bounds of the Netherlands' U{EEZ

94 <http://MarineRegions.org/mrgid/5668>} and U{EPSG:28992<https://EPSG.io/28992>}.

96 @kwarg asRd: Use C{B{asRD}=True} for the bounds in C{meter}, otherwise in

97 C{degrees} (C{bool}).

99 @return: A L{Bounds4Tuple}C{(latS, lonW, latN, lonE)} with lat- and longitudes

100 in C{degrees} or an L{RD4Tuple}C{(minRDx, minRDy, maxRDx, maxRDy)} with

101 the C{quasi-RD} bounds in C{meter}.

102

103 @see: U{EEZ<https://NL.WikiPedia.org/wiki/Nederlandse_Exclusieve_Economische_Zone>}

104 '''

105 b = _RD._bounds4

106 return self._asRD(b) if asRD else b

107

108 def _forward(self, lat, lon, height=0, raiser=None, name=_forward_):

109 '''(INTERNAL) Convert geodetic C{(lat, lon)} and ellipsoidal C{height}

110 to local C{(RDx, RDy)} coordinates and orthometric C{H}.

111 '''

112 lat, lon, _NAN = _LatLon3(lat, lon)

113 if _NAN:

114 RDx = RDy = H = NAN

115 else:

116 RDx, RDy, H = self._forward3(raiser, lat, lon, height)

117 return RDNAP7Tuple(RDx, RDy, H,

118 lat, lon, height, self.forwardDatum, name=name)

119

120 def _forward2(self, lat, lon):

121 # datum-transform C{(lat, lon)} from ETRS to RD-Bessel

122 x, y, z = _geodetic2cartesian(lat, lon, A0.H0_ETRS, self._EETRS)

123 x, y, z = _RD._xETRS2RD.transform(x, y, z) # pseudo

124 return _cartesian2geodetic(x, y, z, A0.E0) # pseudo

125

126 def _forward2x(self, *args): # PYCHOK no cover

127 return self._notOverloaded(*args)

128

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

130 # C{_forward} core, returning C{(RDx, RDy, H)}

131 lat0, lon0 = \

132 lat_, lon_ = self._forward2x(raiser, lat, lon)

133 for _ in range(_TRIPS): # 2.3.3a-f, 1..2

134 latc, lonc = self._rdlatlon2(lat_, lon_, lat0, lon0)

135 if fabs(latc - lat_) < _TOL_D and \

136 fabs(lonc - lon_) < _TOL_D:

137 break

138 lat_, lon_ = latc, lonc

139

140 phiClamC = _ellipsoidal2spherical(latc, lonc)

141 RDx, RDy = _spherical2oblique(*phiClamC)

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

143 height - self._rdNAPh_v(lat, lon, latc, lonc))

144 return RDx, RDy, H

145

146 def forward3(self, lat, lon, **name):

147 '''Datum-transform C{(B{lat}, B{lon})} from GRS80 (ETRS98) to Bessel1841

148 (RD-Bessel) using only C{RDNAPTRANS(tm)2018_v220627}'s similarity.

149

150 @return: A L{LatLonDatum3Tuple}C{(lat, lon, datum)} with C{lat},

151 C{lon} and C{datum} all Bessel1841 (RD-Bessel).

152 '''

153 lat, lon, _NAN = _LatLon3(lat, lon)

154 if _NAN:

155 lat = lon = NAN

156 else:

157 lat, lon = self._forward2(lat, lon)

158 n = name.get(_name_, typename(_RDNAPbase.forward3))

159 return LatLonDatum3Tuple(lat, lon, A0.D0, name=n)

160

161 @property_RO

162 def forwardDatum(self):

163 '''Get the C{forward} datum (L{Datum}, default GRS80).

164 '''

165 return self._datum

166

167 def _inside2(self, raiser, lat, lon):

168 # if RD-Bessel or ETRS C{(lat, lon)} is outside the C{RD}

169 # region raise an error if C{raiser} or self.raiser is True

170 if (raiser or (raiser is None and self.raiser)) and \

171 not _isinside(lat, lon): # _region4

172 raise self._outsidError(lat, lon) # _region4

173 return lat, lon

174

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

176 '''Is geodetic C{(B{lat}, B{lon})} inside the C{RD B{region4}}?

177

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

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

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

181

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

183 if outside the C{RD} region, C{True} otherwise.

184 '''

185 lat, lon, _NAN = _LatLon3(lat, lon)

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

187

188 def isinsideRD(self, RDx, RDy, eps=0):

189 '''Is local C{(B{RDx}, B{RDy})} inside the C{RD B{region4}}?

190

191 @arg RDx: Local C{RD} X (C{meter}, conventionally).

192 @arg RDy: Local C{RD} Y (C{meter}, conventionally).

193 @kwarg eps: Over-/undersize the C{RD} region (C{meter}).

194

195 @return: C{None} if B{C{RDx}} or B{C{RDy}} is NAN, C{False}

196 if outside the C{RD} region, C{True} otherwise.

197 '''

198 x, y, _NAN = _RDxRDy3(RDx, RDy)

199 return None if _NAN else _isinside(x, y, Meter(eps=eps), self.region4(True))

200

201 def _outsidError(self, *lat_lon):

202 # format an RDNAPError for C{lat_lon} outside the RD region

203 E = RDNAPError(lat_lon, txt=_outside__ + _region4.toRepr())

204 return E

205

206 @property

207 def raiser(self):

208 '''Do points outside the C{RD} region cause an C{RDNAPError}?

209 '''

210 return self._raiser

211

212 @raiser.setter # PYCHOK setter!

213 def raiser(self, raiser):

214 '''Use C{True} to throw an C{RDNAPError} for points outside the C{RD} region.

215 '''

216 self._raiser = bool(raiser)

217

218 @property_RO

219 def _rdgrid(self): # PYCHOK no cover

220 return self._notOverloaded()

221

222 def _rdlatlon2(self, lat, lon, lat0=None, lon0=None): # 2.3.2

223 # return the RD-corrected C{(lat, lon)} if inside

224 if _isinside(lat, lon):

225 c_f_N_f6 = _RD._c_f_N_f6(lat, lon)

226 lat_corr = _bilinear(self._rdgrid._lat_corr, *c_f_N_f6)

227 lon_corr = _bilinear(self._rdgrid._lon_corr, *c_f_N_f6)

228

229 if lat0 is lon0 is None: # reverse

230 lat += lat_corr

231 lon += lon_corr

232 else: # forward

233 lat = lat0 - lat_corr

234 lon = lon0 - lon_corr

235 return lat, lon # NAN, NAN?

236

237 def rdNAPh(self, lat, lon): # 2.5.1 and 3.5

238 '''Interpolate the quasi-geoid C{NAPh} height for a geodetic point.

239

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

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

242

243 @return: Quasi-geoid C{NAPh} height C{N} (C{meter}) or C{NAN} if

244 B{C{lat}} or B{C{lon}} is outside the C{RD} region.

245 '''

246 lat, lon, _NAN = _LatLon3(lat, lon)

247 return NAN if _NAN else self._rdNAPh(lat, lon)

248

249 def rdNAPh3(self, RDx, RDy):

250 '''Interpolate the quasi-geoid C{NAPh} height for a local point.

251

252 @arg RDx: Local C{RD} X (C{meter}, conventionally).

253 @arg RDy: Local C{RD} Y (C{meter}, conventionally).

254

255 @return: L{LatLonNgeoid3Tuple}C{(lat, lon, N)} with quasi-geoid

256 C{NAPh} height C{N} in C{meter} or C{NAN} if C{lat} or

257 C{lon} is outside C{RD} region.

258 '''

259 return self._reverse(RDx, RDy, 0, raiser=False).latlonNgeoid

260

261 def _rdNAPh(self, lat, lon):

262 # return C{NAPh} at C{(lat, lon)} or C{NAN} if

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

264 if _isinside(lat, lon):

265 c_f_N_f6 = _RD._c_f_N_f6(lat, lon)

266 N = _bilinear(self._rdgrid._NAP_h, *c_f_N_f6)

267 return Meter(NAPh=N)

268 return NAN # c0 2.5.1e+

269

270 def _rdNAPh_v(self, lat1, lon1, lat2, lon2):

271 # interpolate C{NAPh} at ETRS C{lat1, lon1} for variant 1 or at

272 # RD-corrected or inverse-projected C{lat2, lon2} for variant 2

273 return self._rdNAPh(lat2, lon2) if self.variant == 2 else \

274 self._rdNAPh(lat1, lon1)

275

276 @deprecated_property_RO

277 def region(self): # PYCHOK no cover

278 '''DEPRECATED on 2026.06.12, use method L{region4()<_RDNAPbase.region4>}.'''

279 return self._region4()

280

281 def region4(self, asRD=False):

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

283

284 @kwarg asRd: Use C{B{asRD}=True} for the bounds in C{RD meter},

285 otherwise C{degrees} (C{bool}).

286

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

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

289 L{RD4Tuple}C{(minRDx, minRDy, maxRDx, maxRDy)} with

290 the bounds in C{meter}, truncated to C{millimeter}.

291 '''

292 return self._region4RD[self.variant] if asRD else _region4

293

294 @property_ROver

295 def _region4RD(self):

296 # C{RD} regions in C{meter}, see .__main__._RD4Tuple

297 n = _region4.name

298 d = {1: RD4Tuple(-87853.981, 228817.837, 318159.693, 894090.744, name=n),

299 2: RD4Tuple(-87776.807, 228895.002, 317993.007, 893924.047, name=n)}

300 return d

301

302 def _reverse(self, RDx, RDy, H, raiser=None, name=_reverse_):

303 '''(INTERNAL) Convert local C{(RDx, RDy)} and orthometric height

304 C{H} to geodetic C{lat}, C{lon} and ellipsoidal C{height}.

305 '''

306 RDx, RDy, _NAN = _RDxRDy3(RDx, RDy)

307 if _NAN:

308 h = lat = lon = NAN

309 else:

310 lat, lon, h = self._reverse3(raiser, RDx, RDy, H)

311 return RDNAP7Tuple(RDx, RDy, H,

312 lat, lon, h, self.reverseDatum, name=name)

313

314 def _reverse2(self, lat, lon):

315 # datum-transform C{(lat, lon)} from RD-Bessel to ETRS

316 x, y, z = _geodetic2cartesian(lat, lon, A0.H0, A0.E0)

317 x, y, z = _RD._xRD2ETRS.transform(x, y, z)

318 return _cartesian2geodetic(x, y, z, self._EETRS)

319

320 def _reverse3(self, raiser, RDx, RDy, H): # in .__main__

321 # C{_reverse} core, returning C{(lat, lon, height)}

322 phiClamC = _oblique2spherical(RDx, RDy)

323 latlon = _spherical2ellipsoidal(*phiClamC) # RD-Bessel

324

325 latlon = self._inside2(raiser, *latlon)

326 latclonc = self._rdlatlon2(*latlon) # RD-corrected

327 lat, lon = self._reverse2(*latclonc)

328 h = NAN if H is None or _isNAN(H) else (

329 H + self._rdNAPh_v(lat, lon, *latclonc))

330 return lat, lon, h

331

332 def reverse3(self, lat, lon, **name):

333 '''Datum-transform C{(B{lat}, B{lon})} from Bessel1841 (RD-Bessel) to

334 GRS80 (ETRS98) using only C{RDNAPTRANS(tm)2018_v220627}'s similarity.

335

336 @return: A L{LatLonDatum3Tuple}C{(lat, lon, datum)} with C{lat},

337 C{lon} and C{datum} all GRS80 (ETRS89).

338 '''

339 lat, lon, _NAN = _LatLon3(lat, lon)

340 if _NAN:

341 lat = lon = NAN

342 else:

343 lat, lon = self._reverse2(lat, lon)

344 n = name.get(_name_, typename(_RDNAPbase.reverse3))

345 return LatLonDatum3Tuple(lat, lon, self.reverseDatum, name=n)

346

347 @property_RO

348 def reverseDatum(self):

349 '''Get the C{reverse} datum (L{Datum}, default GRS80).

350 '''

351 return self._datum # sae as .forwardDatum

352

353 def similarity(self, inverse=None): # PYCHOK no cover

354 return self._notOverloaded(inverse=inverse)

355

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

357 '''Return this C{RDNAP20181v1} or C{-v2} instance as a string.

358

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

360

361 @return: This C{RDNAP2018v1} or C{-v2} (C{str}).

362 '''

363 return self.attrs(_name_, 'variant', 'forwardDatum', prec=prec) # _ellipsoid_

364

365 @property_RO

366 def variant(self): # PYCHOK no cover

367 return self._notOverloaded()

368

369

370class RDNAP2018v1(_RDNAPbase):

371 '''Transformer implementing C{variant 1} of the U{RDNAPTRANS(tm)2018_v220627

372 <https://formulieren.kadaster.nl/aanvragen_rdnaptrans>} specification.

373 '''

374 if _FOR_DOCS:

375 __init__ = _RDNAPbase.__init__

376

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

378 '''Convert GRS80 (ETRS98) geodetic C{(B{lat}, B{lon})} and (ellipsoidal)

379 B{C{height}} to local C{RDx}, C{RDy} coordinates and (orthometric)

380 height C{H}.

381

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

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

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

385 C{None} to ignore C{NAPh} interpolation.

386 @kwarg raiser_name: Use C{B{raiser}=True} to raise an L{RDNAPError}

387 if B{C{lat}} or B{C{lon}} is outside the C{RD} region,

388 overriding property C{raiser} (C{bool}) and optional

389 C{B{name}='forward'} (C{str}).

390

391 @return: An L{RDNAP7Tuple}C{(RDx, RDy, H, lat, lon, height, datum)}

392 with local C{RDx}, C{RDy} coordinates and (orthometric) height

393 C{H} in C{meter} or C{NAN} if C{lat} or C{lon} is outside the

394 C{RD} region.

395

396 @raise RDNAPError: If the point is outside the C{RD} region and property

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

398

399 @note: Orthometric height C{(H = h - NAPh)} equals ellipsoidal height C{h}

400 less the quasi-geoid height C{NAPh}.

401 '''

402 return self._forward(lat, lon, height, **raiser_name)

403

404 def _forward2x(self, raiser, *lat_lon): # PYCHOK signature

405 # datum-transform C{(lat, lon)} from ETRS89 to RD-Bessel

406 # and raise an C{RDNAPError} if outside the C{RD} region

407 lat_lon = self._forward2(*lat_lon)

408 return self._inside2(raiser, *lat_lon)

409

410 if _FOR_DOCS:

411 forward3 = _RDNAPbase.forward3

412 isinside = _RDNAPbase.isinside

413 isinsideRD = _RDNAPbase.isinsideRD

414

415 @property_ROver

416 def _rdgrid(self):

417 try:

418 from pyrdnap import v1grid

419 except Exception as x:

420 raise RDNAPError(_v_grid(1), cause=x)

421 return v1grid

422

423 if _FOR_DOCS:

424 rdNAPh = _RDNAPbase.rdNAPh

425 region4 = _RDNAPbase.region4

426

427 def reverse(self, RDx, RDy, H=0, **raiser_name):

428 '''Convert a local C{(B{RDx}, B{RDy})} point and (orthometric) height

429 B{C{H}} to GRS80 (ETRS89) geodetic lat-, longitude and (ellipsoidal)

430 height.

431

432 @arg RDx: Local C{RD} X (C{meter}, conventionally).

433 @arg RDy: Local C{RD} Y (C{meter}, conventionally).

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

435 to ignore C{NAPh} interpolation.

436 @kwarg raiser_name: Use C{B{raiser}=True} to raise an L{RDNAPError}

437 for points outside the C{RD} region, overriding property

438 C{raiser} (C{bool}) and an optional C{B{name}='reverse'}

439 (C{str}).

440

441 @return: An L{RDNAP7Tuple}C{(RDx, RDy, H, lat, lon, height, datum)}

442 with geodetic C{lat}, C{lon} and C{datum} GRS80 (ETRS89) and

443 (ellipsoidal) C{height} in C{meter} or C{NAN} if C{lat} or

444 C{lon} is outside the C{RD} region.

445

446 @raise RDNAPError: If the point is outside the C{RD} region and property

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

448

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

450 plus the quasi-geoid height C{NAPh}.

451 '''

452 return self._reverse(RDx, RDy, H, **raiser_name)

453

454 if _FOR_DOCS:

455 reverse3 = _RDNAPbase.reverse3

456

457 def similarity(self, inverse=False):

458 '''Get the similarity transform (C{Similarity}).

459

460 @kwarg inverse: Use C{True} for the C{reverse} or C{False}

461 for the C{forward} transform (C{bool}).

462 '''

463 return _RD._xRD2ETRS if inverse else _RD._xETRS2RD

464

465 @property_ROver

466 def variant(self):

467 '''Get this C{RDNAP2018}'s variant (C{int}).

468 '''

469 return 1

470

471

472class RDNAP2018v2(_RDNAPbase):

473 '''Transformer implementing C{variant 2} of the U{RDNAPTRANS(tm)2018_v220627

474 <https://formulieren.kadaster.nl/aanvragen_rdnaptrans>} specification.

475 '''

476 if _FOR_DOCS:

477 __init__ = _RDNAPbase.__init__

478

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

480 '''Convert GRS80 (ETRS98) geodetic C{(B{lat}, B{lon})} and (ellipsoidal)

481 B{C{height}} to local C{RDx, RDy} coordinates and (orthometric) height

482 C{H}, provided the point is not outside the C{RD} region.

483

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

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

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

487 C{None} to ignore C{NAPh} interpolation.

488 @kwarg raiser_name: Use C{B{raiser}=True} to raise an L{RDNAPError}

489 if B{C{lat}} or B{C{lon}} is outside the C{RD} region,

490 overriding property C{raiser} (C{bool}) and an optional

491 C{B{name}='forward'} (C{str}).

492

493 @return: An L{RDNAP7Tuple}C{(RDx, RDy, H, lat, lon, height, datum)}

494 with local C{RDx}, C{RDy} coordinates and (orthometric) height

495 C{H} in C{meter}. C{RDx}, C{RDy} and C{H} are all C{NAN} if

496 C{lat} or C{lon} is outside the C{RD} region.

497

498 @raise RDNAPError: If the point is outside the C{RD} region and property

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

500

501 @note: Orthometric height C{(H = h - NAPh)} equals ellipsoidal height C{h}

502 less the quasi-geoid height C{NAPh}.

503 '''

504 raiser, name = _xkwds_pop2(raiser_name, raiser=self.raiser)

505 try: # force outside exception

506 r = self._forward(lat, lon, height, raiser=True, **name)

507 except RDNAPError as x:

508 if raiser or _outside__ not in str(x):

509 raise # reraise

510 d = self.forwardDatum

511 n = name.get(_name_, _forward_)

512 r = RDNAP7Tuple(NAN, NAN, NAN, lat, lon, height, d, name=n)

513 return r

514

515 def _forward2x(self, *raiser_lat_lon): # 2.3.4

516 # NO datum-transform C{(lat, lon)} to RD-Bessel, but

517 # raise an C{RDNAPError} if outside the C{RD} region

518 # (using the ETRS as RD-Bessel lat- and longitudes)

519 return self._inside2(*raiser_lat_lon)

520

521 if _FOR_DOCS:

522 forward3 = _RDNAPbase.forward3

523 isinside = _RDNAPbase.isinside

524 isinsideRD = _RDNAPbase.isinsideRD

525

526 @property_ROver

527 def _rdgrid(self):

528 try:

529 from pyrdnap import v2grid

530 except Exception as x:

531 raise RDNAPError(_v_grid(2), cause=x)

532 return v2grid

533

534 if _FOR_DOCS:

535 rdNAPh = _RDNAPbase.rdNAPh

536 region4 = _RDNAPbase.region4

537

538 def reverse(self, RDx, RDy, H=0, **raiser_name):

539 '''Convert a local C{(B{RDx}, B{RDy})} point and (orthometric) height

540 B{C{H}} to GRS80 (ETRS89) geodetic lat-, longitude and (ellispoidal)

541 height, provided the point is not outside the C{RD} region.

542

543 @arg RDx: Local C{RD} X (C{meter}, conventionally).

544 @arg RDy: Local C{RD} Y (C{meter}, conventionally).

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

546 to ignore C{NAPh} interpolation.

547 @kwarg raiser_name: Use C{B{raiser}=True} to raise an L{RDNAPError}

548 for points outside the C{RD} region, overriding property

549 C{raiser} (C{bool}) and an optional C{B{name}='reverse'}

550 (C{str}).

551

552 @return: An L{RDNAP7Tuple}C{(RDx, RDy, H, lat, lon, height, datum)}

553 with geodetic C{lat}, C{lon} and C{datum} GRS80 (ETRS89) and

554 (ellipsoidal) C{height} in C{meter}. Outside the C{RD} region

555 all C{lat}, C{lon} and C{height} are C{NAN}.

556

557 @raise RDNAPError: If the point is outside the C{RD} region and property

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

559

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

561 plus the quasi-geoid height C{NAPh}.

562 '''

563 raiser, name = _xkwds_pop2(raiser_name, raiser=self.raiser)

564 try: # force outside exception

565 r = self._reverse(RDx, RDy, H, raiser=True, **name)

566 except RDNAPError as x:

567 if raiser or _outside__ not in str(x):

568 raise # reraise

569 d = self.reverseDatum

570 n = name.get(_name_, _reverse_)

571 r = RDNAP7Tuple(RDx, RDy, H, NAN, NAN, NAN, d, name=n)

572 return r

573

574 if _FOR_DOCS:

575 reverse3 = _RDNAPbase.reverse3

576

577 def similarity(self, inverse=False):

578 '''Get the similarity transform (C{None}, always).

579 '''

580 return None if inverse else None

581

582 @property_ROver

583 def variant(self):

584 '''Get this C{RDNAP2018}'s variant (C{int}).

585 '''

586 return 2

587

588

589def _atan3(y, x, x0): # 2.2.3e and 3.1.1i

590 # equiv to math.atan2 iff x0 is y

591 if x > 0:

592 r = atan(y / x)

593 elif x < 0:

594 r = atan(y / x) + copysign(PI, x0)

595# elif _isNAN(x) or _isNAN(y) or _isNAN(x0):

596# r = NAN

597 else:

598 r = copysign(PI_2, x0) if x0 else _0_0

599 return r

600

601

602def _atan_exp(w): # 2.4.1c

603 return atan(exp(w)) * _2_0 - PI_2

604

605

606def _bilinear(v_grid, c_latI, f_latI, latN_f, # 2.3.1f and g

607 c_lonI, f_lonI, lonN_f):

608 # interpolate a lat_corr_, lon_corr_ or NAP_h...

609 # assert isinstance(v_grid, _V_grid), v_grid

610 ne = v_grid(c_latI, c_lonI)

611# if c_latI == f_latI and c_lonI == f_lonI:

612# return ne

613 nw = v_grid(c_latI, f_lonI)

614 se = v_grid(f_latI, c_lonI)

615 sw = v_grid(f_latI, f_lonI)

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

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

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

619

620

621def _cartesian2geodetic(x, y, z, E): # 2.2.3 == EcefUPC.reverse?

622 # convert cartesian C{(x, y, z)} to C{E}-geodetic C{(lat, lon)}

623 r = hypot(x, y)

624# if _isNAN(r) or _isNAN(z):

625# return NAN, NAN

626 if r > _TOL_M:

627 a = E.a * E.e2

628 phi_ = atan(z / r) # atan2(z, r)

629 for _ in range(_TRIPS): # 4..6

630 s = sin(phi_)

631 s *= a / sqrt(_1_0 - s**2 * E.e2)

632 phi = atan((z + s) / r) # atan2(z + s, r)

633 if fabs(phi - phi_) < _TOL_R:

634 break

635 phi_ = phi

636 else:

637 phi = copysign(PI_2, z)

638 lam = _atan3(y, x, y)

639 return map1(degrees, phi, lam) # lat, lon

640

641

642def _ellipsoidal2spherical(lat, lon): # 2.4.1

643 # convert RD-Bessel C{(lat, lon)} to spherical C{(𝛷, 𝛬)}

644 phiC = phi = Phid(lat) # clip=90

645 if PI_2 > phi > -PI_2: # 2.4.1c

646 q = A0.log_tan(phi) - A0.log_e_2(phi)

647 w = A0.N0 * q + A0.M0 # 2.4.1b

648 phiC = _atan_exp(w)

649 lamC = (Lamd(lon) - A0.LAM0) * A0.N0 + A0.LAM0C # 2.4.1d

650 return phiC, lamC # -Capital 𝛷, 𝛬

651

652

653def _eq0(r, r0=_0_0):

654 return fabs(r - r0) < _TOL_R

655

656

657# def _eq0d(d, d0=_0_0):

658# return fabs(d - d0) < _TOL_D

659

660

661def _geodetic2cartesian(lat, lon, h, E): # 2.2.1

662 # convert C{E}-geodetic C{(lat, lon)} to cartesian C{(x, y, z)}

663 y, x = sincos2d(lon)

664 z, c = sincos2d(lat)

665 n = E.a / sqrt(_1_0 - z**2 * E.e2)

666 H = _isNAN0(h)

667 c *= n + H

668 x *= c

669 y *= c

670 z *= n * (_1_0 - E.e2) + H

671 return x, y, z

672

673

674def _isinside(lat, lon, eps=0, region4=_region4):

675 # is C{(lat, lon)} inside C{region4}, optionally over- or

676 # undersized by positive respectively negative C{eps}?

677 # returns: C{False} if C{lat} or C{lon} outside or NAN,

678 # C{True} otherwise.

679 S, W, N, E = region4

680 return ((S - lat) <= eps and (lat - N) <= eps and

681 (W - lon) <= eps and (lon - E) <= eps) if eps else \

682 (S <= lat <= N and W <= lon <= E)

683

684

685def _LatLon3(lat, lon):

686 lat, lon = Lat(lat), Lon(lon)

687 return lat, lon, (_isNAN(lon) or _isNAN(lat))

688

689

690def _ne0(r, r0=_0_0):

691 return fabs(r - r0) > _TOL_R

692

693

694# def _ne0d(d, d0=_0_0):

695# return fabs(d - d0) > _TOL_D

696

697

698def _oblique2spherical(x, y): # 3.1.1

699 # inverse oblique stereographic conformal projection from

700 # C{RD (x, y)} to spherical C{(𝛷, 𝛬)}, see C++ function

701 # sterea_e_inverse in U{Proj/src/projections/sterea.cpp

702 # <https://Proj.org/en/stable/operations/projections/sterea.html>}

703 x -= A0.X0

704 y -= A0.Y0

705 r = hypot(x, y)

706 if r > _TOL_M: # x and y

707 s0, c0 = A0.sincos2PHI0C

708 sp, cp = sincos2(atan(r / A0.RK2) * _2_0) # psi atan2(r, A0.RK2)

709 ca = sp * y / r

710 xN = cp * c0 - ca * s0

711 yN = sp * x / r

712 zN = cp * s0 + ca * c0

713 phiC = asin(zN)

714# elif _isNAN(r):

715# return NAN, NAN

716 else:

717 _, xN = A0.sincos2PHI0C

718 yN = _0_0

719 phiC = A0.PHI0C # asin(sin(PHI0C))

720 lamC = _atan3(yN, xN, x) + A0.LAM0C

721 return phiC, lamC # -Capital 𝛷, 𝛬

722

723

724def _RDxRDy3(RDx, RDy):

725 x, y = map1(Meter, RDx, RDy)

726 return x, y, (_isNAN(x) or _isNAN(y))

727

728

729def _spherical2ellipsoidal(phiC, lamC): # 3.1.2

730 # inverse Gauss conformal projection from

731 # spherical C{(𝛷, 𝛬)} to RD-Bessel C{(lat, lon)}

732 phi = phiC

733 if PI_2 > phi > -PI_2:

734 q = (A0.log_tan(phi) - A0.M0) / A0.N0

735# w = A0.log_tan(phi)

736 for _ in range(_TRIPS): # 3..6

737 phi_ = phi

738 phi = _atan_exp(A0.log_e_2(phi) + q)

739 if fabs(phi - phi_) < _TOL_R:

740 break

741 lam = (lamC - A0.LAM0C) / A0.N0 + A0.LAM0

742 lam += floor((PI - lam) / PI2) * PI2

743 return map1(degrees, phi, lam) # lat, lon

744

745

746def _spherical2oblique(phiC, lamC): # 2.4.2

747 # oblique stereographic conformal projection

748 # from spherical C{(𝛷, 𝛬)} to C{RD (x, y)}

749 x = A0.X0 # 2.4.2g

750 y = A0.Y0 # 2.4.2h

751 a = phiC - A0.PHI0C # 𝛷 - 𝛷0

752 b = lamC - A0.LAM0C # 𝛬 - 𝛬0

753 if (_ne0(a) or _ne0(b)) and (_ne0(phiC, -A0.PHI0C) or

754 _ne0(lamC, -A0.LAM0C + PI)):

755 s0, c0 = A0.sincos2PHI0C # sin(𝛷0), cos(𝛷0)

756 s, c = sincos2(phiC) # sin(𝛷), cos(𝛷)

757 sp_22 = sin(a * _0_5)**2 + \

758 sin(b * _0_5)**2 * c * c0 # sin(𝜓/2)**2

759 if EPS0 < sp_22 < EPS1:

760 # r = 2kR * tan(𝜓/2)

761 # q = r / (sin(𝜓/2) * cos(𝜓/2) * 2)

762 # = 2kR * sin(𝜓/2) / (sin(𝜓/2) * cos(𝜓/2)**2 * 2)

763 # = 2kR / (cos(𝜓/2)**2 * 2)

764 # = 2kR / ((1 - sin(𝜓/2)**2) * 2)

765 # = 2kR / (2 - sin(𝜓/2)**2 * 2)

766 t = sp_22 * _2_0 # 0 < t < 2

767 q = A0.RK2 / (_2_0 - t)

768 x += q * (c * sin(b))

769 y += q * (s - s0 + s0 * t) / c0

770 elif _eq0(a) and _eq0(b):

771 pass

772 else: # if _eq0(phiC, -A0.PHI0C) and _eq0(lamC, A0.LAM0C - PI):

773 x = y = NAN

774# else:

775# raise RDNAPError((phiC, lamC))

776 return x, y

777

778

779__all__ += _ALL_DOCS(_RDNAPbase)

780__all__ += _ALL_OTHER(RDNAP2018v1, RDNAP2018v2, # passed along from PyGeodesy

781 Bounds4Tuple, LatLonNgeoid3Tuple, RD4Tuple)

782del _ALL_DOCS, _ALL_OTHER

783

784# **) MIT License

785#

787#

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

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

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

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

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

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

794#

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

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

797#

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

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

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

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

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

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

804# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pyrdnap / rdnap2018.py: 94%

300 statements