Coverage for pyrdnap/rdnap2018.py: 91%

231 statements  

« prev     ^ index     » next       coverage.py v7.10.7, created at 2026-05-08 18:11 -0400

1 

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

3 

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

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

6Each provide a C{forward} method to convert geodetic lat-/longitudes and heights to C{RD} 

7coodinates and C{NAP} heights and a C{reverse} method for converting the other way. 

8 

9The C{forward} and C{reverse} results of L{RDNAP2018v1} meet the C{RDNAPTRANS(tm)2018_v220627} 

10self-validation requirements of C{0.000000010 degrees} and C{0.0010 meter} for tests inside 

11the C{RD} region, see B{C{Note below}}. Class L{RDNAP2018v2} does not and is not required to. 

12 

13The original C{RDNAPTRANS(tm)2018_v220627} grid files for both variants are I{not included} 

14in C{PyPRDNAP} and C{pyrdnap} due to the size of those files. Instead, the grid files for each 

15variant I{include only} the C{lat_corr_}, C{lon_corr_} and C{_NAP_quasi_geoid_height_...} columns, 

16extracted from the original grid files with leading and trailing zeros removed and formatted as 

17row-order matrices. 

18 

19@note: L{RDNAP2018v1}, C{PyRDNAP} and C{pyrdnap} have B{not been formally validated} and are 

20 B{not certified} to carry the trademark name C{RDNAPTRANS(tm)}. 

21''' 

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

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

24 

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

26from pyrdnap.v_grids import _V_grid, _v_gridz_import 

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

28 _isNAN, _earth_datum, 

29 _ALL_DOCS, _ALL_OTHER, _FOR_DOCS, 

30 _NamedBase, notOverloaded) 

31from pygeodesy import (EPS0, EPS1, NAN, PI_2, PI, PI2, # "consterns" 

32 Datum, Datums, Ellipsoid, # datums, ellipsoids 

33 property_RO, property_ROnce, # props 

34 Lamd, Phid, # units 

35 sincos2, sincos2d) # utily 

36 

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

38 fabs, floor, hypot, radians, sin, sqrt 

39 

40__all__ = () 

41__version__ = '26.05.08' 

42 

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

44_TOL_M = 1e-6 # meter 

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

46_TRIPS = 16 # 5..6 sufficient 

47 

48 

49class _RDNAPbase(_NamedBase): 

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

51 ''' 

52 _datum = None # forward, v1 reverse Datum, lazily 

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

54 _raiser = False 

55 

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

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

58 

59 @kwarg a_ellipsoid: An ellipsoid (L{Ellipsoid}), a datum (L{Datum}) or the 

60 ellipsoid's equatorial radius (C{scalar}, conventionally 

61 in C{meter}), see B{C{f}}. Default C{Datums.GRS80}. 

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

63 is specified as C{scalar}, ignored otherwise. 

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

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

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

67 

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

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

70 ''' 

71 if a_ellipsoid is f is None: 

72 self._datum = Datums.GRS80 

73 else: 

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

75 E = self._datum.ellipsoid 

76 if not E.isOblate: 

77 raise RDNAPError('not oblate %r' % (E,)) 

78 self._EETRS = E 

79 if raiser: 

80 T = self._datum.transform 

81 if not T.isunity: 

82 raise RDNAPError('not unity %r' % (T,)) 

83 self._raiser = True 

84 if name: 

85 self.name = name 

86 

87 def forward(self, lat, lon, height=0, raiser=None, name='forward'): 

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

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

90 

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

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

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

94 to ignore C{NAPh} interpolation. 

95 @kwarg raiser: If C{True} raise an L{RDNAPError} if B{C{lat}} or 

96 B{C{lon}} is outside the C{RD} region (C{bool}), 

97 if C{False} don't, overriding property C{raiser}. 

98 @kwarg name: Optional C{B{name}='forward'} (C{str}). 

99 

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

101 with local C{RDx}, C{RDy} coordinates and C{NAPh} height. 

102 ''' 

103 lat0, lon0 = \ 

104 lat_, lon_ = self._forwardXform(lat, lon, raiser) 

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

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

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

108 fabs(lonc - lon_) < _TOL_D: 

109 break 

110 lat_, lon_ = latc, lonc 

111 

112 philamC = _ellipsoidal2spherical(latc, lonc) 

113 RDx, RDy = _spherical2oblique(*philamC) 

114 NAPh = NAN if _isNAN(height) else (height - self.rdNAPh(lat, lon)) # 2.5.2 

115 return RDNAP7Tuple(RDx, RDy, NAPh, lat, lon, height, 

116 self.forwardDatum, name=name) 

117 

118 @property_RO 

119 def forwardDatum(self): 

120 '''Get the geodetic C{forward} datum (L{Datum}). 

121 ''' 

122 return self._datum 

123 

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

125 # default and variant 2: no datum Xform 

126 if (raiser or (raiser is None and self._raiser)) and \ 

127 not _RD.isinside(lat, lon): 

128 raise self._outsidError(lat, lon) 

129 return lat, lon 

130 

131 _forwardXform = _inside2 # no datum Xform 

132 

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

134 '''Is C{(B{lat}, B{lon})} inside the C{RD} region (C{bool})? 

135 

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

137 ''' 

138 return _RD.isinside(lat, lon, eps) 

139 

140 def _outsidError(self, *lat_lon): 

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

142 return RDNAPError('%r outside %s' % (lat_lon, self.region)) 

143 

144 @property_RO 

145 def _rdgrid(self): 

146 raise notOverloaded(self) 

147 

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

149 # return the RD-corrected C{(lat, lon)} 

150 if _RD.isinside(lat, lon): 

151 c_f_N_f6 = _RD.c_f_N_f6(lat, lon) 

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

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

154 

155 if lat0 is lon0 is None: # reverse 

156 lat += lat_corr 

157 lon += lon_corr 

158 else: # forward 

159 lat = lat0 - lat_corr 

160 lon = lon0 - lon_corr 

161 return lat, lon # NAN, NAN? 

162 

163 def rdNAPh(self, lat, lon, raiser=False): # 2.5.1 and 3.5 

164 '''Interpolate the C{NAP_h} quasi-geoid-height 

165 I{within} the C{RD} region. 

166 

167 @arg lat: Latitude (C{degrees} GRS80 (ETRS89), geodetic). 

168 @arg lon: Longitude (C{degrees} GRS80 (ETRS89), geodetic). 

169 @kwarg raiser: If C{True} raise an L{RDNAPError} if B{C{lat}} or 

170 B{C{lon}} is outside the C{RD} region (C{bool}), 

171 otherwise don't abd return C{NAN}. 

172 

173 @return: C{NAPh} (C{meter}) or C{NAN} if B{C{lat}} 

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

175 ''' 

176 if _RD.isinside(lat, lon): 

177 c_f_N_f6 = _RD.c_f_N_f6(lat, lon) 

178 NAP_h = _bilinear(self._rdgrid._NAP_h, *c_f_N_f6) 

179 return NAP_h 

180 elif raiser: 

181 raise self._outsidError(lat, lon) 

182 return NAN # c0 2.5.1e+ 

183 

184 @property_RO 

185 def region(self): 

186 '''Get the C{RD} region as L{RDregion4Tuple}C{(S, W, N, E)}, all C{GRS80 (ETRS89) degrees}. 

187 ''' 

188 return _RD.region 

189 

190 def _reverse(self, RDx, RDy, NAPh, raiser, name='reverse'): 

191 '''(INTERNAL) Convert local C{(B{RDx}, B{RDy})} coordinates and 

192 B{C{NAPh}} quasi-geoid-height to geodetic GRS80 (ETRS89) or 

193 Bessel1841 (RD-Bessel) C{lat}, C{lon} and C{height}. 

194 ''' 

195 philamC = _oblique2spherical(RDx, RDy) 

196 lat, lon = _spherical2ellipsoidal(*philamC) 

197 

198 lat, lon = self._rdlatlon2(lat, lon) 

199 lat, lon = self._reverseXform(lat, lon, raiser) 

200 h = NAN if _isNAN(NAPh) else (NAPh + self.rdNAPh(lat, lon)) 

201 return RDNAP7Tuple(RDx, RDy, NAPh, lat, lon, h, 

202 self.reverseDatum, name=name) 

203 

204 @property_RO 

205 def reverseDatum(self): 

206 '''Get the geodetic C{reverse} datum (L{Datum}), GRS80 or Bessel1841. 

207 ''' 

208 return {1: self._datum, 

209 2: A0.D0}.get(self.variant) 

210 

211 _reverseXform = _inside2 # no datum Xform 

212 

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

214 '''Return this C{RDNAP2018#v} instance as a string. 

215 

216 @kwarg prec: Precision, number of decimal digits (0..9). 

217 

218 @return: This C{RDNAP2018#v} (C{str}). 

219 ''' 

220 return self.attrs('name', 'variant', 'forwardDatum', prec=prec) # _ellipsoid_, _name__ 

221 

222 @property_RO 

223 def variant(self): 

224 raise None 

225 

226 

227class RDNAP2018v1(_RDNAPbase): 

228 '''Transformer implementing RD NAP 2018 C{variant 1}. 

229 ''' 

230 if _FOR_DOCS: 

231 __init__ = _RDNAPbase.__init__ 

232 forward = _RDNAPbase.forward 

233 

234 def _forwardXform(self, lat, lon, raiser): 

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

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

237 x, y, z = _RD._xETRS2RD.transform(x, y, z) 

238 lat, lon = _cartesian2geodetic(x, y, z, A0.E0) 

239 return self._inside2(lat, lon, raiser) 

240 

241 @property_ROnce 

242 def _rdgrid(self): 

243 try: 

244 from pyrdnap import v1grid 

245 except (AttributeError, ImportError): 

246 v1grid = _v_gridz_import(self.variant) 

247 return v1grid 

248 

249 def reverse(self, RDx, RDy, NAPh=0, raiser=None, **name): # RDNAP to GRS80 (ETRS89) 

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

251 to geodetic B{GRS80 (ETRS89)} C{(lat, lon, height)}. 

252 

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

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

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

256 or C{NAN} to ignore C{NAPh} interpolation. 

257 @kwarg raiser: If C{True} raise an L{RDNAPError} if C{lat} or 

258 C{lon} is outside the C{RD} region (C{bool}), if 

259 C{False} don't, overriding property C{raiser}. 

260 @kwarg name: Optional C{B{name}='reverse'} (C{str}). 

261 

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

263 datum)} with geodetic C{lat} and C{lon}, C{height} 

264 and C{datum} B{GRS80 (ETRS89)}. 

265 ''' 

266 return self._reverse(RDx, RDy, NAPh, raiser, **name) 

267 

268 def _reverseXform(self, lat, lon, raiser): 

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

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

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

272 lat, lon = _cartesian2geodetic(x, y, z, self._EETRS) 

273 return self._inside2(lat, lon, raiser) 

274 

275 def similarity(self, inverse=False): 

276 '''Get the forward or reverse datum transform (C{Similarity}). 

277 ''' 

278 return _RD._xRD2ETRS if inverse else _RD._xETRS2RD 

279 

280 @property_ROnce 

281 def variant(self): 

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

283 ''' 

284 return 1 

285 

286 

287class RDNAP2018v2(_RDNAPbase): 

288 '''Transformer implementing RD NAP 2018 C{variant 2}. 

289 

290 @note: Method L{RDNAP2018v2.reverse} returns geodetic B{Bessel1841 

291 (RD-Bessel)} and B{not GRS80 (ETRS89)} lat-/longitudes. 

292 ''' 

293 if _FOR_DOCS: 

294 __init__ = _RDNAPbase.__init__ 

295 forward = _RDNAPbase.forward 

296 

297 @property_ROnce 

298 def _rdgrid(self): 

299 try: 

300 from pyrdnap import v2grid 

301 except (AttributeError, ImportError): 

302 v2grid = _v_gridz_import(self.variant) 

303 return v2grid 

304 

305 def reverse(self, RDx, RDy, NAPh=0, raiser=None, **name): # RDNAP to RD-Bessel 

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

307 to geodetic B{Bessel1841 (RD-Bessel)} C{(lat, lon, height)}. 

308 

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

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

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

312 or C{NAN} to ignore C{NAPh} interpolation. 

313 @kwarg raiser: If C{True} raise an L{RDNAPError} if C{lat} or 

314 C{lon} is outside the C{RD} region (C{bool}), if 

315 C{False} don't, overriding property C{raiser}. 

316 @kwarg name: Optional C{B{name}='reverse'} (C{str}). 

317 

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

319 datum)} with geodetic C{lat} and C{lon}, C{height} 

320 and C{datum} B{Bessel1841 (RD-Bessel)}. 

321 ''' 

322 return self._reverse(RDx, RDy, NAPh, raiser, **name) 

323 

324 def similarity(self, **unused): # PYCHOK signature 

325 '''Get the forward or reverse datum transform, always C{None}. 

326 ''' 

327 return None 

328 

329 @property_ROnce 

330 def variant(self): 

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

332 ''' 

333 return 2 

334 

335 

336class RDNAPError(ValueError): 

337 '''Error raised for C{RD} and C{NAP} issues. 

338 ''' 

339 pass 

340 

341 

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

343 # equiv to math.atan2 iff x0 is y 

344 if x > 0: 

345 r = atan(y / x) 

346 elif x < 0: 

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

348 else: 

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

350 return r 

351 

352 

353def _atan_exp(w): # 2.4.1c 

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

355 

356 

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

358 c_lonI, f_lonI, lonN_f): 

359 # interpolate a lat_corr_, lon_corr_ or NAP_h_... 

360 assert isinstance(v_grid, _V_grid) 

361 nw = v_grid(c_latI, f_lonI) 

362 ne = v_grid(c_latI, c_lonI) 

363 sw = v_grid(f_latI, f_lonI) 

364 se = v_grid(f_latI, c_lonI) 

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

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

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

368 

369 

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

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

372 r = hypot(x, y) 

373 if r > _TOL_M: 

374 a = E.a * E.e2 

375 phi_ = atan(z / r) 

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

377 s = sin(phi_) 

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

379 phi = atan((z + s) / r) 

380 if fabs(phi - phi_) < _TOL_R: 

381 break 

382 phi_ = phi 

383 else: 

384 phi = copysign(PI_2, z) 

385 lam = _atan3(y, x, y) 

386 return degrees(phi), degrees(lam) 

387 

388 

389def _ellipsoidal2spherical(lat, lon): # 2.4.1 

390 # convert geodetic C{(lat, lon)} to spherical C{(𝛷, 𝛬)} 

391 phiC = phi = Phid(lat) 

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

393 q = A0.ln_tan(phi) - A0.ln_e_2(phi) 

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

395 phiC = _atan_exp(w) 

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

397 return phiC, lamC # -Capital 

398 

399 

400def _eq0(r, r0=_0_0): 

401 return fabs(r - r0) < _TOL_R 

402 

403 

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

405# return fabs(d - d0) < _TOL_D 

406 

407 

408def _geodetic2cartesian(lat, lon, E, h0=0): # 2.2.1 

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

410 y, x = sincos2d(lon) 

411 z, c = sincos2d(lat) 

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

413 c *= n + h0 

414 x *= c 

415 y *= c 

416 z *= n * (_1_0 - E.e2) + h0 

417 return x, y, z 

418 

419 

420def _ne0(r, r0=_0_0): 

421 return fabs(r - r0) > _TOL_R 

422 

423 

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

425# return fabs(d - d0) > _TOL_D 

426 

427 

428def _oblique2spherical(x, y): # 3.1.1 

429 # inverse oblique stereographic conformal projection 

430 # from 2-D C{(x, y)} to spherical C{(𝛷, 𝛬)} 

431 x -= A0.X0 

432 y -= A0.Y0 

433 r = hypot(x, y) 

434 if r > _TOL_M: # x and y 

435 s0, c0 = A0.sincos2PHI0C 

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

437 ca = sp * y / r 

438 xN = cp * c0 - ca * s0 

439 yN = sp * x / r 

440 zN = ca * c0 + cp * s0 

441 phiC = asin(zN) 

442 else: 

443 _, xN = A0.sincos2PHI0C 

444 yN = _0_0 

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

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

447 return phiC, lamC # -Capital 

448 

449 

450def _spherical2ellipsoidal(phiC, lamC): # 3.1.2 

451 # inverse Gauss conformal projection from 

452 # spherical C{(𝛷, 𝛬)} to geodetic C{(phi, lam)} 

453 phi = phiC 

454 if PI_2 > phi > -PI_2: 

455 q = (A0.ln_tan(phi) - A0.M0) / A0.N0 

456# w = A0.ln_tan(phi) 

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

458 phi_ = phi 

459 phi = _atan_exp(A0.ln_e_2(phi) + q) 

460 if fabs(phi - phi_) < _TOL_R: 

461 break 

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

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

464 return degrees(phi), degrees(lam) 

465 

466 

467def _spherical2oblique(phiC, lamC): # 2.4.2 

468 # oblique stereographic conformal projection 

469 # from spherical C{(𝛷, 𝛬)} to 2-D C{(x, y)} 

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

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

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

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

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

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

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

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

478 if EPS0 < sp_22 < EPS1: 

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

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

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

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

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

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

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

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

487 x = q * (c * sin(b)) 

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

489 else: 

490 x = y = _0_0 # NAN? 

491 x += A0.X0 

492 y += A0.Y0 

493 elif _eq0(phiC, A0.PHI0C) and _eq0(lamC, A0.LAM0C): 

494 x = A0.X0 # x0 2.4.2g 

495 y = A0.Y0 # y0 2.4.2h 

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

497 x = y = NAN 

498# else: 

499# raise RDNAPError(str((phiC, lamC))) 

500 return x, y 

501 

502 

503__all__ += _ALL_DOCS(_RDNAPbase) 

504__all__ += _ALL_OTHER(RDNAP2018v1, RDNAP2018v2, RDNAPError, 

505 Datum, Ellipsoid) # passed along from PyGeodesy 

506 

507 

508# **) MIT License 

509# 

510# Copyright (C) 2026-2026 -- mrJean1 at Gmail -- All Rights Reserved. 

511# 

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

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

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

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

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

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

518# 

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

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

521# 

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

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

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

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

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

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

528# OTHER DEALINGS IN THE SOFTWARE.