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. Each provide a C{forward} method to convert geodetic lat-/longitudes and height

7to local C{RD} coodinates and C{NAP} heights and a C{reverse} method for converting vice-versa.

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

10validated to meet the C{RDNAPTRANS(tm)2018_v220627} requirements.

12Likewise for the L{RDNAP2018v2.forward} and L{.reverse<RDNAP2018v2.reverse>} results.

13'''

14# make sure int/int division yields float quotient, see .basics

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

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

18from pyrdnap.v_grids import _v_grid # _V_grid

19from pyrdnap.__pygeodesy import (_0_0, _0_5, _1_0, _2_0,

20 _isNAN, _isNAN0, _earth_datum, _xkwds_pop2,

21 _name_, _ALL_DOCS, _all_OTHER, _FOR_DOCS,

22 _NamedBase, RDNAPError)

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

24 typename, LatLonDatum3Tuple, RD4Tuple, # namedTuples

25 deprecated_property_RO, property_RO, property_ROnce, # props

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

27 sincos2, sincos2d) # utily

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

30 fabs, floor, hypot, radians, sin, sqrt

32__all__ = ()

33__version__ = '26.07.09'

35_forward_ = 'forward'

36_outside__ = 'outside '

37_region4 = _RD._region4

38_reverse_ = 'reverse'

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

40_TOL_M = 1e-6 # meter

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

42_TRIPS = 16 # 5..6 sufficient

45class _RDNAPbase(_NamedBase):

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

47 '''

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

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

50 _raiser = False

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

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

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

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

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

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

59 specified as C{scalar}, ignored otherwise.

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

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

62 @kwarg name: Optional name (C{str}).

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

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

66 '''

67 if a_ellipsoid is f is None:

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

69 else:

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

71 self._EETRS = E = self._datum.ellipsoid

72 if not E.isOblate:

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

74 if raiser: # PYCHOK no cover

75 T = self._datum.transform

76 if not T.isunity:

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

78 self.raiser = True

79 if name:

80 self.name = name

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

83 '''(INTERNAL) Convert geodetic C{(B{lat}, B{lon})} and B{C{height}}

84 to local C{(RDx, RDy)} coordinates and C{NAPh} quasi-geoid-height.

85 '''

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

87 if _NAN:

88 RDx = RDy = NAPh = NAN

89 else:

90 RDx, RDy, NAPh = self._forward3(raiser, lat, lon, height)

91 return RDNAP7Tuple(RDx, RDy, NAPh,

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

94 def _forward2(self, lat, lon):

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

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

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

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

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

101 return self._notOverloaded(*args)

102

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

104 # C{_forward} core, returning C{(RDx, RDy, NAPh)}

105 lat0, lon0 = \

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

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

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

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

110 fabs(lonc - lon_) < _TOL_D:

111 break

112 lat_, lon_ = latc, lonc

113

114 phiClamC = _ellipsoidal2spherical(latc, lonc)

115 RDx, RDy = _spherical2oblique(*phiClamC)

116 NAPh = NAN if height is None or _isNAN(height) else (

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

118 return RDx, RDy, NAPh

119

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

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

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

123

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

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

126 '''

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

128 if _NAN:

129 lat = lon = NAN

130 else:

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

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

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

134

135 @property_RO

136 def forwardDatum(self):

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

138 '''

139 return self._datum

140

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

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

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

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

145 not _isinside(lat, lon): # _region4

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

147 return lat, lon

148

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

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

151

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

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

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

155

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

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

158 '''

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

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

161

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

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

164

165 @arg RDx: X coordinate (C{meter}, local).

166 @arg RDy: Y coordinate (C{meter}, local).

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

168

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

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

171 '''

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

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

174

175 def _outsidError(self, *lat_lon):

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

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

178 return E

179

180 @property

181 def raiser(self):

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

183 '''

184 return self._raiser

185

186 @raiser.setter # PYCHOK setter!

187 def raiser(self, raiser):

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

189 '''

190 self._raiser = bool(raiser)

191

192 @property_RO

193 def _rdgrid(self): # PYCHOK no cover

194 return self._notOverloaded()

195

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

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

198 if _isinside(lat, lon):

199 c_f_N_f6 = _RD._c_f_N_f6(lat, lon)

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

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

202

203 if lat0 is lon0 is None: # reverse

204 lat += lat_corr

205 lon += lon_corr

206 else: # forward

207 lat = lat0 - lat_corr

208 lon = lon0 - lon_corr

209 return lat, lon # NAN, NAN?

210

211 def rdNAPh(self, lat, lon, height=0): # 2.5.1 and 3.5

212 '''Interpolate the C{NAPh} quasi-geoid-height for a point

213 C{(lat, lon)} I{within} the C{RD} region.

214

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

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

217 @kwarg height: Optional geoid height (C{meter}, conventionally).

218

219 @return: C{NAPh} quasi-geoid-height (C{meter}) or C{NAN} if

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

221 '''

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

223 if _NAN:

224 h = NAN

225 else:

226 h = self._rdNAPh(lat, lon)

227 if not _isNAN(h):

228 h = Meter(height=height) - h

229 return Meter(NAPh=h)

230

231 def _rdNAPh(self, lat, lon):

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

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

234 if _isinside(lat, lon):

235 c_f_N_f6 = _RD._c_f_N_f6(lat, lon)

236 return _bilinear(self._rdgrid._NAP_h, *c_f_N_f6)

237 return NAN # c0 2.5.1e+

238

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

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

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

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

243 self._rdNAPh(lat1, lon1)

244

245 @deprecated_property_RO

246 def region(self): # PYCHOK no cover

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

248 return self._region4()

249

250 def region4(self, asRD=False):

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

252

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

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

255

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

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

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

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

260 '''

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

262

263 @property_ROnce

264 def _region4RD(self):

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

266 n = _region4.name

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

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

269 return d

270

271 def _reverse(self, RDx, RDy, NAPh, raiser=None, name=_reverse_):

272 '''(INTERNAL) Convert local C{(B{RDx}, B{RDy})} and B{C{NAPh}}

273 quasi-geoid-height to geodetic C{lat}, C{lon} and C{height}.

274 '''

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

276 if _NAN:

277 h = lat = lon = NAN

278 else:

279 lat, lon, h = self._reverse3(raiser, RDx, RDy, NAPh)

280 return RDNAP7Tuple(RDx, RDy, NAPh,

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

282

283 def _reverse2(self, lat, lon):

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

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

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

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

288

289 def _reverse3(self, raiser, RDx, RDy, NAPh): # in .__main__

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

291 phiClamC = _oblique2spherical(RDx, RDy)

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

293

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

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

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

297 h = NAN if NAPh is None or _isNAN(NAPh) else (

298 NAPh + self._rdNAPh_v(lat, lon, *latclonc))

299 return lat, lon, h

300

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

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

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

304

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

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

307 '''

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

309 if _NAN:

310 lat = lon = NAN

311 else:

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

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

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

315

316 @property_RO

317 def reverseDatum(self):

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

319 '''

320 return self._datum # sae as .forwardDatum

321

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

323 return self._notOverloaded(inverse=inverse)

324

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

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

327

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

329

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

331 '''

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

333

334 @property_RO

335 def variant(self): # PYCHOK no cover

336 return self._notOverloaded()

337

338

339class RDNAP2018v1(_RDNAPbase):

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

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

342 '''

343 if _FOR_DOCS:

344 __init__ = _RDNAPbase.__init__

345

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

347 '''Convert GRS80 (ETRS98) geodetic C{(B{lat}, B{lon})} and B{C{height}}

348 to local C{RDx}, C{RDy} coordinates and C{NAPh} quasi-geoid-height.

349

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

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

352 @kwarg height: Height, optional (C{meter} above geoid) or C{NAN}

353 to ignore C{NAPh} interpolation.

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

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

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

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

358

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

360 with local C{RDx}, C{RDy} coordinates and C{NAPh} height, all

361 in C{meter} or with C{height} is C{NAN} if C{lat} or C{lon} is

362 outside the C{RD} region.

363

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

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

366 '''

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

368

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

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

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

372 lat_lon = self._forward2(*lat_lon)

373 return self._inside2(raiser, *lat_lon)

374

375 if _FOR_DOCS:

376 forward3 = _RDNAPbase.forward3

377 isinside = _RDNAPbase.isinside

378 isinsideRD = _RDNAPbase.isinsideRD

379

380 @property_ROnce

381 def _rdgrid(self):

382 try:

383 from pyrdnap import v1grid

384 except Exception as x:

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

386 return v1grid

387

388 if _FOR_DOCS:

389 rdNAPh = _RDNAPbase.rdNAPh

390 region4 = _RDNAPbase.region4

391

392 def reverse(self, RDx, RDy, NAPh=0, **raiser_name):

393 '''Convert a local C{(B{RDx}, B{RDy})} point and B{C{NAPh}} height to

394 GRS80 (ETRS89) geodetic lat-, longitude and height, B{by default}.

395

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

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

398 @kwarg NAPh: C{NAP} quasi-geoid-height (C{meter}, conventionally) or

399 C{NAN} to ignore C{NAPh} interpolation.

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

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

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

403 (C{str}).

404

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

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

407 C{height} in C{meter} or C{NAN} if C{lat} or C{lon} is outside

408 the C{RD} region.

409

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

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

412 '''

413 return self._reverse(RDx, RDy, NAPh, **raiser_name)

414

415 if _FOR_DOCS:

416 reverse3 = _RDNAPbase.reverse3

417

418 def similarity(self, inverse=False):

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

420

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

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

423 '''

424 return _RD._xRD2ETRS if inverse else _RD._xETRS2RD

425

426 @property_ROnce

427 def variant(self):

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

429 '''

430 return 1

431

432

433class RDNAP2018v2(_RDNAPbase):

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

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

436 '''

437 if _FOR_DOCS:

438 __init__ = _RDNAPbase.__init__

439

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

441 '''Convert GRS80 (ETRS98) geodetic C{(B{lat}, B{lon})} and B{C{height}}

442 to local C{RDx, RDy} coordinates and C{NAPh} quasi-geoid-height,

443 provided the point is not outside the C{RD} region.

444

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

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

447 @kwarg height: Height, optional (C{meter} above geoid) or C{NAN}

448 to ignore C{NAPh} interpolation.

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

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

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

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

453

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

455 with local C{RDx}, C{RDy} coordinates and C{NAPh} height, all

456 in C{meter} or all C{NAN} if C{lat} or C{lon} is outside the

457 C{RD} region.

458

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

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

461 '''

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

463 try: # force outside exception

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

465 except RDNAPError as x:

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

467 raise # reraise

468 d = self.forwardDatum

469 n = name.get(_name_, _forward_)

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

471 return r

472

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

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

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

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

477 return self._inside2(*raiser_lat_lon)

478

479 if _FOR_DOCS:

480 forward3 = _RDNAPbase.forward3

481 isinside = _RDNAPbase.isinside

482 isinsideRD = _RDNAPbase.isinsideRD

483

484 @property_ROnce

485 def _rdgrid(self):

486 try:

487 from pyrdnap import v2grid

488 except Exception as x:

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

490 return v2grid

491

492 if _FOR_DOCS:

493 rdNAPh = _RDNAPbase.rdNAPh

494 region4 = _RDNAPbase.region4

495

496 def reverse(self, RDx, RDy, NAPh=0, **raiser_name):

497 '''Convert a local C{(B{RDx}, B{RDy})} point and B{C{NAPh}} height to

498 GRS80 (ETRS89) geodetic lat-, longitude and height, provided the

499 point is not outside the C{RD} region.

500

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

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

503 @kwarg NAPh: C{NAP} quasi-geoid-height (C{meter}, conventionally) or

504 C{NAN} to ignore C{NAPh} interpolation.

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

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

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

508 (C{str}).

509

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

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

512 C{height} in C{meter} or with C{lat}, C{lon} and C{height}

513 all C{NAN} if outside the C{RD} region.

514

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

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

517 '''

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

519 try: # force outside exception

520 r = self._reverse(RDx, RDy, NAPh, raiser=True, **name)

521 except RDNAPError as x:

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

523 raise # reraise

524 d = self.reverseDatum

525 n = name.get(_name_, _reverse_)

526 r = RDNAP7Tuple(RDx, RDy, NAPh, NAN, NAN, NAN, d, name=n)

527 return r

528

529 if _FOR_DOCS:

530 reverse3 = _RDNAPbase.reverse3

531

532 def similarity(self, inverse=False):

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

534 '''

535 return None if inverse else None

536

537 @property_ROnce

538 def variant(self):

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

540 '''

541 return 2

542

543

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

545 # equiv to math.atan2 iff x0 is y

546 if x > 0:

547 r = atan(y / x)

548 elif x < 0:

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

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

551# r = NAN

552 else:

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

554 return r

555

556

557def _atan_exp(w): # 2.4.1c

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

559

560

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

562 c_lonI, f_lonI, lonN_f):

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

564 # assert isinstance(v_grid, _V_grid), v_grid

565 ne = v_grid(c_latI, c_lonI)

566 nw = v_grid(c_latI, f_lonI)

567 se = v_grid(f_latI, c_lonI)

568 sw = v_grid(f_latI, f_lonI)

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

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

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

572

573

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

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

576 r = hypot(x, y)

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

578# return NAN, NAN

579 if r > _TOL_M:

580 a = E.a * E.e2

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

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

583 s = sin(phi_)

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

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

586 if fabs(phi - phi_) < _TOL_R:

587 break

588 phi_ = phi

589 else:

590 phi = copysign(PI_2, z)

591 lam = _atan3(y, x, y)

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

593

594

595def _ellipsoidal2spherical(lat, lon): # 2.4.1

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

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

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

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

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

601 phiC = _atan_exp(w)

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

603 return phiC, lamC # -Capital 𝛷, 𝛬

604

605

606def _eq0(r, r0=_0_0):

607 return fabs(r - r0) < _TOL_R

608

609

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

611# return fabs(d - d0) < _TOL_D

612

613

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

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

616 y, x = sincos2d(lon)

617 z, c = sincos2d(lat)

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

619 H = _isNAN0(h)

620 c *= n + H

621 x *= c

622 y *= c

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

624 return x, y, z

625

626

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

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

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

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

631 # C{True} otherwise.

632 S, W, N, E = region4

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

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

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

636

637

638def _LatLon3(lat, lon):

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

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

641

642

643def _ne0(r, r0=_0_0):

644 return fabs(r - r0) > _TOL_R

645

646

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

648# return fabs(d - d0) > _TOL_D

649

650

651def _oblique2spherical(x, y): # 3.1.1

652 # inverse oblique stereographic conformal projection from

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

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

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

656 x -= A0.X0

657 y -= A0.Y0

658 r = hypot(x, y)

659 if r > _TOL_M: # x and y

660 s0, c0 = A0.sincos2PHI0C

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

662 ca = sp * y / r

663 xN = cp * c0 - ca * s0

664 yN = sp * x / r

665 zN = cp * s0 + ca * c0

666 phiC = asin(zN)

667# elif _isNAN(r):

668# return NAN, NAN

669 else:

670 _, xN = A0.sincos2PHI0C

671 yN = _0_0

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

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

674 return phiC, lamC # -Capital 𝛷, 𝛬

675

676

677def _RDxRDy3(RDx, RDy):

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

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

680

681

682def _spherical2ellipsoidal(phiC, lamC): # 3.1.2

683 # inverse Gauss conformal projection from

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

685 phi = phiC

686 if PI_2 > phi > -PI_2:

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

688# w = A0.log_tan(phi)

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

690 phi_ = phi

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

692 if fabs(phi - phi_) < _TOL_R:

693 break

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

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

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

697

698

699def _spherical2oblique(phiC, lamC): # 2.4.2

700 # oblique stereographic conformal projection

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

702 x = A0.X0 # 2.4.2g

703 y = A0.Y0 # 2.4.2h

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

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

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

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

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

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

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

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

712 if EPS0 < sp_22 < EPS1:

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

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

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

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

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

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

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

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

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

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

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

724 pass

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

726 x = y = NAN

727# else:

728# raise RDNAPError((phiC, lamC))

729 return x, y

730

731

732__all__ += _ALL_DOCS(_RDNAPbase)

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

734del _ALL_DOCS, _all_OTHER

735

736# **) MIT License

737#

739#

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

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

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

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

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

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

746#

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

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

749#

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

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

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

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

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

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

756# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pyrdnap / rdnap2018.py: 94%

293 statements