Coverage for pyrdnap/rdnap2018.py: 92%

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, LatLonN3Tuple, 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, # basics, "consterns"

26 typename, Bounds4Tuple, LatLonDatum3Tuple, 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.07.31'

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{LatLonN3Tuple}C{(lat, lon, N)} with the 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 r = self._reverse(RDx, RDy, 0, raiser=False)

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

261

262 def _rdNAPh(self, lat, lon):

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

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

265 if _isinside(lat, lon):

266 c_f_N_f6 = _RD._c_f_N_f6(lat, lon)

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

268 return Meter(NAPh=N)

269 return NAN # c0 2.5.1e+

270

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

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

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

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

275 self._rdNAPh(lat1, lon1)

276

277 @deprecated_property_RO

278 def region(self): # PYCHOK no cover

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

280 return self._region4()

281

282 def region4(self, asRD=False):

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

284

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

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

287

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

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

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

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

292 '''

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

294

295 @property_ROver

296 def _region4RD(self):

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

298 n = _region4.name

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

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

301 return d

302

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

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

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

306 '''

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

308 if _NAN:

309 h = lat = lon = NAN

310 else:

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

312 return RDNAP7Tuple(RDx, RDy, H,

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

314

315 def _reverse2(self, lat, lon):

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

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

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

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

320

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

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

323 phiClamC = _oblique2spherical(RDx, RDy)

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

325

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

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

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

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

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

331 return lat, lon, h

332

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

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

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

336

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

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

339 '''

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

341 if _NAN:

342 lat = lon = NAN

343 else:

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

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

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

347

348 @property_RO

349 def reverseDatum(self):

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

351 '''

352 return self._datum # sae as .forwardDatum

353

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

355 return self._notOverloaded(inverse=inverse)

356

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

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

359

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

361

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

363 '''

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

365

366 @property_RO

367 def variant(self): # PYCHOK no cover

368 return self._notOverloaded()

369

370

371class RDNAP2018v1(_RDNAPbase):

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

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

374 '''

375 if _FOR_DOCS:

376 __init__ = _RDNAPbase.__init__

377

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

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

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

381 height C{H}.

382

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

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

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

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

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

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

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

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

391

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

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

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

395 C{RD} region.

396

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

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

399

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

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

402 '''

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

404

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

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

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

408 lat_lon = self._forward2(*lat_lon)

409 return self._inside2(raiser, *lat_lon)

410

411 if _FOR_DOCS:

412 forward3 = _RDNAPbase.forward3

413 isinside = _RDNAPbase.isinside

414 isinsideRD = _RDNAPbase.isinsideRD

415

416 @property_ROver

417 def _rdgrid(self):

418 try:

419 from pyrdnap import v1grid

420 except Exception as x:

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

422 return v1grid

423

424 if _FOR_DOCS:

425 rdNAPh = _RDNAPbase.rdNAPh

426 region4 = _RDNAPbase.region4

427

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

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

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

431 height.

432

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

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

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

436 to ignore C{NAPh} interpolation.

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

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

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

440 (C{str}).

441

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

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

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

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

446

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

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

449

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

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

452 '''

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

454

455 if _FOR_DOCS:

456 reverse3 = _RDNAPbase.reverse3

457

458 def similarity(self, inverse=False):

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

460

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

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

463 '''

464 return _RD._xRD2ETRS if inverse else _RD._xETRS2RD

465

466 @property_ROver

467 def variant(self):

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

469 '''

470 return 1

471

472

473class RDNAP2018v2(_RDNAPbase):

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

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

476 '''

477 if _FOR_DOCS:

478 __init__ = _RDNAPbase.__init__

479

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

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

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

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

484

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

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

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

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

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

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

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

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

493

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

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

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

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

498

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

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

501

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

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

504 '''

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

506 try: # force outside exception

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

508 except RDNAPError as x:

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

510 raise # reraise

511 d = self.forwardDatum

512 n = name.get(_name_, _forward_)

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

514 return r

515

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

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

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

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

520 return self._inside2(*raiser_lat_lon)

521

522 if _FOR_DOCS:

523 forward3 = _RDNAPbase.forward3

524 isinside = _RDNAPbase.isinside

525 isinsideRD = _RDNAPbase.isinsideRD

526

527 @property_ROver

528 def _rdgrid(self):

529 try:

530 from pyrdnap import v2grid

531 except Exception as x:

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

533 return v2grid

534

535 if _FOR_DOCS:

536 rdNAPh = _RDNAPbase.rdNAPh

537 region4 = _RDNAPbase.region4

538

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

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

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

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

543

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

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

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

547 to ignore C{NAPh} interpolation.

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

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

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

551 (C{str}).

552

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

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

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

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

557

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

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

560

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

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

563 '''

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

565 try: # force outside exception

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

567 except RDNAPError as x:

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

569 raise # reraise

570 d = self.reverseDatum

571 n = name.get(_name_, _reverse_)

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

573 return r

574

575 if _FOR_DOCS:

576 reverse3 = _RDNAPbase.reverse3

577

578 def similarity(self, inverse=False):

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

580 '''

581 return None if inverse else None

582

583 @property_ROver

584 def variant(self):

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

586 '''

587 return 2

588

589

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

591 # equiv to math.atan2 iff x0 is y

592 if x > 0:

593 r = atan(y / x)

594 elif x < 0:

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

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

597# r = NAN

598 else:

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

600 return r

601

602

603def _atan_exp(w): # 2.4.1c

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

605

606

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

608 c_lonI, f_lonI, lonN_f):

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

610 # assert isinstance(v_grid, _V_grid), v_grid

611 ne = v_grid(c_latI, c_lonI)

612 nw = v_grid(c_latI, f_lonI)

613 se = v_grid(f_latI, c_lonI)

614 sw = v_grid(f_latI, f_lonI)

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

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

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

618

619

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

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

622 r = hypot(x, y)

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

624# return NAN, NAN

625 if r > _TOL_M:

626 a = E.a * E.e2

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

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

629 s = sin(phi_)

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

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

632 if fabs(phi - phi_) < _TOL_R:

633 break

634 phi_ = phi

635 else:

636 phi = copysign(PI_2, z)

637 lam = _atan3(y, x, y)

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

639

640

641def _ellipsoidal2spherical(lat, lon): # 2.4.1

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

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

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

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

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

647 phiC = _atan_exp(w)

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

649 return phiC, lamC # -Capital 𝛷, 𝛬

650

651

652def _eq0(r, r0=_0_0):

653 return fabs(r - r0) < _TOL_R

654

655

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

657# return fabs(d - d0) < _TOL_D

658

659

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

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

662 y, x = sincos2d(lon)

663 z, c = sincos2d(lat)

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

665 H = _isNAN0(h)

666 c *= n + H

667 x *= c

668 y *= c

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

670 return x, y, z

671

672

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

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

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

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

677 # C{True} otherwise.

678 S, W, N, E = region4

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

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

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

682

683

684def _LatLon3(lat, lon):

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

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

687

688

689def _ne0(r, r0=_0_0):

690 return fabs(r - r0) > _TOL_R

691

692

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

694# return fabs(d - d0) > _TOL_D

695

696

697def _oblique2spherical(x, y): # 3.1.1

698 # inverse oblique stereographic conformal projection from

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

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

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

702 x -= A0.X0

703 y -= A0.Y0

704 r = hypot(x, y)

705 if r > _TOL_M: # x and y

706 s0, c0 = A0.sincos2PHI0C

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

708 ca = sp * y / r

709 xN = cp * c0 - ca * s0

710 yN = sp * x / r

711 zN = cp * s0 + ca * c0

712 phiC = asin(zN)

713# elif _isNAN(r):

714# return NAN, NAN

715 else:

716 _, xN = A0.sincos2PHI0C

717 yN = _0_0

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

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

720 return phiC, lamC # -Capital 𝛷, 𝛬

721

722

723def _RDxRDy3(RDx, RDy):

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

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

726

727

728def _spherical2ellipsoidal(phiC, lamC): # 3.1.2

729 # inverse Gauss conformal projection from

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

731 phi = phiC

732 if PI_2 > phi > -PI_2:

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

734# w = A0.log_tan(phi)

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

736 phi_ = phi

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

738 if fabs(phi - phi_) < _TOL_R:

739 break

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

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

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

743

744

745def _spherical2oblique(phiC, lamC): # 2.4.2

746 # oblique stereographic conformal projection

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

748 x = A0.X0 # 2.4.2g

749 y = A0.Y0 # 2.4.2h

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

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

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

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

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

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

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

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

758 if EPS0 < sp_22 < EPS1:

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

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

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

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

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

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

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

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

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

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

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

770 pass

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

772 x = y = NAN

773# else:

774# raise RDNAPError((phiC, lamC))

775 return x, y

776

777

778__all__ += _ALL_DOCS(_RDNAPbase)

779__all__ += _all_OTHER(RDNAP2018v1, RDNAP2018v2, RD4Tuple)

780del _ALL_DOCS, _all_OTHER

781

782# **) MIT License

783#

785#

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

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

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

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

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

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

792#

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

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

795#

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

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

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

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

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

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

802# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pyrdnap / rdnap2018.py: 92%

301 statements