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, transforming from and to

11ETRS89 (GRS80) geodetic points.

13The L{RDNAP2018v2.forward} and L{.reverse<RDNAP2018v2.reverse>} results have been formally

14validated to meet the C{RDNAPTRANS(tm)2018_v220627} requirements, transforming from ETRS89

15(GRS80) and to RD-Bessel (Bessel1841) geodetic points.

16'''

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

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

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

21from pyrdnap.v_grids import _v_grid # _V_grid

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

23 _isNAN, _isNAN0, _earth_datum,

24 _ALL_DOCS, _all_OTHER, _FOR_DOCS,

25 _NamedBase, RDNAPError)

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

27 LatLonDatum3Tuple, # namedTuples

28 deprecated_property_RO, property_RO, property_ROnce, # props

29 Lamd, Lat, Lon, Phid, # units

30 sincos2, sincos2d) # utily

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

33 fabs, floor, hypot, radians, sin, sqrt

35__all__ = ()

36__version__ = '26.06.18'

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

39_TOL_M = 1e-6 # meter

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

41_TRIPS = 16 # 5..6 sufficient

44class _RDNAPbase(_NamedBase): # in .rd0._RD.regionB

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

46 '''

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

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

49 _raiser = False

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

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

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

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

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

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

58 specified as C{scalar}, ignored otherwise.

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

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

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

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

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

65 '''

66 if a_ellipsoid is f is None:

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

68 else:

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

70 self._EETRS = E = self._datum.ellipsoid

71 if not E.isOblate:

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

73 if raiser: # PYCHOK no cover

74 T = self._datum.transform

75 if not T.isunity:

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

77 self._raiser = True

78 if name:

79 self.name = name

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

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

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

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

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

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

88 to ignore C{NAPh} interpolation.

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

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

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

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

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

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

96 '''

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

98 lat0, lon0 = \

99 lat_, lon_ = self._forwardX2(raiser, lat, lon)

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

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

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

103 fabs(lonc - lon_) < _TOL_D:

104 break

105 lat_, lon_ = latc, lonc

106

107 phiClamC = _ellipsoidal2spherical(latc, lonc)

108 RDx, RDy = _spherical2oblique(*phiClamC)

109 NAPh = NAN if _isNAN(height) else (height - # NOT lat0, lon0

110 self._rdNAPh_v(lat, lon, latc, lonc)) # 2.5.2

111 return RDNAP7Tuple(RDx, RDy, NAPh,

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

113

114 def forward3(self, lat, lon, name='forward3'):

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

116 (RD-Bessel) as specified by C{RDNAPTRANS(tm)2018_v220627}.

117

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

119 C{lat} and C{lon} all Bessel1841 (RD-Bessel).

120 '''

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

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

123 lat, lon = _cartesian2geodetic(x, y, z, A0.E0) # pseudo

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

125

126 @property_RO

127 def forwardDatum(self):

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

129 '''

130 return self._datum

131

132 def _forwardX2(self, *args): # PYCHOK no cover

133 return self._notOverloaded(*args)

134

135 def _inside2(self, raiser, lat, lon, **asRD):

136 # if RD-Bessel C{(lat, lon)} is not inside C{RD} region

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

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

139 not _RD.isinside(lat, lon, **asRD):

140 raise self._outsidError(lat, lon, **asRD)

141 return lat, lon

142

143 def isinside(self, lat, lon, asRD=True, eps=0):

144 '''Is geodetic C{(B{lat}, B{lon})} inside the C{RD} or C{ETRS} region?

145

146 @kwarg asRD: Use C{B{asRD}=False} for the C{ETRS} region and in case

147 C{(B{lat}, B{lon})} are ETRS89 (GRS80), not Bessel1841

148 (RD_Bessel) (C{bool}).

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

150

151 @return: C{True} if inside or on, C{False} otherwise.

152 '''

153 return _RD.isinside(Lat(lat), Lon(lon), asRD, eps)

154

155 def _outsidError(self, *lat_lon, **asRD):

156 # format an RDNAPError for C{lat_lon} outside RD or ETRS region

157 r = self.region4(**asRD).toRepr()

158 return RDNAPError(lat_lon, txt='outside ' + r)

159

160 @property_RO

161 def _rdgrid(self): # PYCHOK no cover

162 return self._notOverloaded()

163

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

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

166 if _RD.isinside(lat, lon):

167 c_f_N_f6 = _RD._c_f_N_f6(lat, lon)

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

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

170

171 if lat0 is lon0 is None: # reverse

172 lat += lat_corr

173 lon += lon_corr

174 else: # forward

175 lat = lat0 - lat_corr

176 lon = lon0 - lon_corr

177 return lat, lon # NAN, NAN?

178

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

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

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

182

183 @arg lat: Latitude (C{degrees}, RD-Bessel geodetic).

184 @arg lon: Longitude (C{degrees}, RD-Bessel geodetic).

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

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

187 otherwise don't and return C{NAN}.

188

189 @return: C{NAPh} (C{meter}) or C{NAN} if C{B{raiser} is False}

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

191 '''

192 return self._rdNAPh(Lat(lat), Lon(lon), raiser)

193

194 def _rdNAPh(self, lat, lon, raiser):

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

196 if _RD.isinside(lat, lon): # asRD=True, eps=0

197 c_f_N_f6 = _RD._c_f_N_f6(lat, lon)

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

199 elif raiser:

200 raise self._outsidError(lat, lon) # asRD=True

201 return NAN # c0 2.5.1e+

202

203 def _rdNAPh_v(self, lat1, lon1, lat2, lon2): # asRD=True

204 '''(INTERNAL) Interpolate C{NAPh} at ETRS C{lat1, lon1} for variant 1 or

205 at RD-corrected or inverse-projected C{lat2, lon2} for variant 2.

206 '''

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

208 self._rdNAPh(lat1, lon1, False)

209

210 @deprecated_property_RO

211 def region(self): # PYCHOK no cover

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

213 return self._region4RD

214

215 def region4(self, asRD=True): # in .rd0._RD

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

217

218 @kwarg asRD: Use C{B{asRD}=False} to get the C{ETRS} (ETRS89) instead of the

219 C{RD} (RD-Bessel) region (C{bool}).

220

221 @return: A L{Bounds4Tuple}C{(latS, lonW, latN, lonE)} with C{RD-Bessel}

222 (Bessel1841) or C{ETRS} (ETRS89) geodetic lat- and longtudes.

223 '''

224 return self._region4RD if asRD else self._region4ETRS

225

226 @property_ROnce

227 def _region4ETRS(self): # as ETRS (ETRS89) L{Bounds4Tuple}

228 S, W, N, E = r = self._region4RD

229 s, w, _ = self.reverse3(S, W)

230 n, e, _ = self.reverse3(N, E)

231 _ETRS = r.name.replace('RD', 'ETRS')

232 return r.classof(s, w, n, e, name=_ETRS) # r.dup(latS=s, ...)

233

234 @property_ROnce

235 def _region4RD(self): # as RD-Bessel L{Bounds4Tuple}

236 return _RD._region4RD

237

238 def _reverse(self, RDx, RDy, NAPh, asRD=False, raiser=None, name='reverse', asETRS=None):

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

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

241 to RD-Bessel C{B{asRD}=True} or ETRS C{B{asRB}=False} (or use

242 C{B{asETRS}=True} respectively C{False} overriding C{B{asRD}}).

243 '''

244 phiClamC = _oblique2spherical(RDx, RDy)

245 latlon = _spherical2ellipsoidal(*phiClamC)

246

247 latc, lonc = self._rdlatlon2(*latlon)

248 lat, lon, d = self._reverseX3(raiser, latc, lonc)

249 h = NAN if _isNAN(NAPh) else (NAPh +

250 self._rdNAPh_v(lat, lon, *latlon))

251

252 if (asRD if asETRS is None else (not asETRS)):

253 lat, lon, d = latc, lonc, A0.D0

254 return RDNAP7Tuple(RDx, RDy, NAPh,

255 lat, lon, h, d, name=name)

256

257 def reverse3(self, lat, lon, name='reverse3'):

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

259 GRS80 (ETRS98) as specified by C{RDNAPTRANS(tm)2018_v220627}.

260

261 @return: A L{LatLon3Tuple}C{(lat, lon, datum)} with C{datum} and

262 C{lat} and C{lon} all GRS80 (ETRS89).

263 '''

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

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

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

267 return LatLonDatum3Tuple(lat, lon, self.forwardDatum, name=name)

268

269 @property_RO

270 def reverseDatum(self):

271 '''Get the I{default} C{reverse} datum (L{Datum}), GRS80 or Bessel1841.

272 '''

273 return {1: self._datum, # self.forwardDatum

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

275

276 def _reverseX3(self, *raiser_lat_lon):

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

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

279 lat, lon = self._inside2(*raiser_lat_lon)

280 return self.reverse3(lat, lon)

281

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

283 return self._notOverloaded(inverse=inverse)

284

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

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

287

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

289

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

291 '''

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

293

294 @property_RO

295 def variant(self): # PYCHOK no cover

296 return self._notOverloaded()

297

298

299class RDNAP2018v1(_RDNAPbase):

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

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

302

303 @note: Method L{RDNAP2018v2.reverse} returns B{by default validated GRS80

304 (ETRS89)} geodetic lat-, longitude and datum.

305 '''

306 if _FOR_DOCS:

307 __init__ = _RDNAPbase.__init__

308 forward = _RDNAPbase.forward

309 forward3 = _RDNAPbase.forward3

310

311 def _forwardX2(self, raiser, *lat_lon): # PYCHOK signature

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

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

314 lat, lon, _ = self.forward3(*lat_lon)

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

316

317 if _FOR_DOCS:

318 isinside = _RDNAPbase.isinside

319 rdNAPh = _RDNAPbase.rdNAPh

320 region4 = _RDNAPbase.region4

321

322 @property_ROnce

323 def _rdgrid(self):

324 try:

325 from pyrdnap import v1grid

326 except Exception as x:

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

328 return v1grid

329

330 def reverse(self, RDx, RDy, NAPh=0, asRD=False, **raiser_name):

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

332 B{GRS80 (ETRS89)} geodetic C{(lat, lon, height)} B{by default}.

333

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

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

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

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

338 @kwarg asRD: Use C{B{asRD}=True} to return (non-validated) Bessel1841

339 (RD-Bessel) instead of (validated) GRS80 (ETRS89) geodetic

340 lat- and longitudes (C{bool}).

341 @kwarg raiser_name: Like the C{forward} method, C{B{raiser}=None}

342 (C{bool}) and optional C{B{name}='reverse'} (C{str}).

343

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

345 with geodetic C{lat} and C{lon}, C{height} and C{datum}

346 B{GRS80 (ETRS89)} or C{Bessel1841 (RD-Bessel)}.

347

348 @note: L{RDNAP2018v1.reverse} has been validated only for default

349 C{B{asRD}=False} per C{RDNAPTRANS(tm)2018_v220627}.

350 '''

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

352

353 if _FOR_DOCS:

354 reverse3 = _RDNAPbase.reverse3

355

356 def similarity(self, inverse=False):

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

358

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

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

361 '''

362 return _RD._xRD2ETRS if inverse else _RD._xETRS2RD

363

364 @property_ROnce

365 def variant(self):

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

367 '''

368 return 1

369

370

371class RDNAP2018v2(_RDNAPbase):

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

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

374

375 @note: Method L{RDNAP2018v2.reverse} returns B{by default validated

376 Bessel1841 (RD-Bessel)} geodetic lat-, longitude and datum.

377 '''

378 if _FOR_DOCS:

379 __init__ = _RDNAPbase.__init__

380 forward = _RDNAPbase.forward

381 forward3 = _RDNAPbase.forward3

382

383 def _forwardX2(self, *raiser_lat_lon):

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

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

386 return self._inside2(*raiser_lat_lon) # asRD=False?

387

388 @property_ROnce

389 def _rdgrid(self):

390 try:

391 from pyrdnap import v2grid

392 except Exception as x:

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

394 return v2grid

395

396 if _FOR_DOCS:

397 isinside = _RDNAPbase.isinside

398 rdNAPh = _RDNAPbase.rdNAPh

399 region4 = _RDNAPbase.region4

400

401 def reverse(self, RDx, RDy, NAPh=0, asRD=True, **raiser_name):

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

403 B{Bessel1841 (RD-Bessel)} geodetic C{(lat, lon, height)} B{by default}.

404

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

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

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

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

409 @kwarg asRD: Use C{B{asRD}=False} to return (non-validated) GRS80

410 (ETRS89) instead of (validated) Bessel1841 (RD-Bessel)

411 geodetic lat- and longitudes (C{bool}).

412 @kwarg raiser_name: Like the C{forward} method, C{B{raiser}=None}

413 (C{bool}) and optional C{B{name}='reverse'} (C{str}).

414

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

416 with geodetic C{lat} and C{lon}, C{height} and C{datum}

417 B{Bessel1841 (RD-Bessel)} or C{GRS80 (ETRS89)}.

418

419 @note: L{RDNAP2018v2.reverse} has been validated only for default

420 C{B{asRD}=True} per C{RDNAPTRANS(tm)2018_v220627}.

421 '''

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

423

424 if _FOR_DOCS:

425 reverse3 = _RDNAPbase.reverse3

426

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

428 '''Get the similarity transform, always C{None}.

429 '''

430 return None

431

432 @property_ROnce

433 def variant(self):

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

435 '''

436 return 2

437

438

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

440 # equiv to math.atan2 iff x0 is y

441 if x > 0:

442 r = atan(y / x)

443 elif x < 0:

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

445 else:

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

447 return r

448

449

450def _atan_exp(w): # 2.4.1c

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

452

453

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

455 c_lonI, f_lonI, lonN_f):

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

457 # assert isinstance(v_grid, _V_grid), v_grid

458 ne = v_grid(c_latI, c_lonI)

459 nw = v_grid(c_latI, f_lonI)

460 se = v_grid(f_latI, c_lonI)

461 sw = v_grid(f_latI, f_lonI)

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

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

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

465

466

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

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

469 r = hypot(x, y)

470 if r > _TOL_M:

471 a = E.a * E.e2

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

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

474 s = sin(phi_)

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

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

477 if fabs(phi - phi_) < _TOL_R:

478 break

479 phi_ = phi

480 else:

481 phi = copysign(PI_2, z)

482 lam = _atan3(y, x, y)

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

484

485

486def _ellipsoidal2spherical(lat, lon): # 2.4.1

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

488 phiC = phi = Phid(lat)

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

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

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

492 phiC = _atan_exp(w)

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

494 return phiC, lamC # -Capital 𝛷, 𝛬

495

496

497def _eq0(r, r0=_0_0):

498 return fabs(r - r0) < _TOL_R

499

500

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

502# return fabs(d - d0) < _TOL_D

503

504

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

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

507 y, x = sincos2d(lon)

508 z, c = sincos2d(lat)

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

510 H = _isNAN0(h)

511 c *= n + H

512 x *= c

513 y *= c

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

515 return x, y, z

516

517

518def _ne0(r, r0=_0_0):

519 return fabs(r - r0) > _TOL_R

520

521

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

523# return fabs(d - d0) > _TOL_D

524

525

526def _oblique2spherical(x, y): # 3.1.1

527 # inverse oblique stereographic conformal projection from

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

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

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

531 x -= A0.X0

532 y -= A0.Y0

533 r = hypot(x, y)

534 if r > _TOL_M: # x and y

535 s0, c0 = A0.sincos2PHI0C

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

537 ca = sp * y / r

538 xN = cp * c0 - ca * s0

539 yN = sp * x / r

540 zN = cp * s0 + ca * c0

541 phiC = asin(zN)

542 else:

543 _, xN = A0.sincos2PHI0C

544 yN = _0_0

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

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

547 return phiC, lamC # -Capital 𝛷, 𝛬

548

549

550def _spherical2ellipsoidal(phiC, lamC): # 3.1.2

551 # inverse Gauss conformal projection from

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

553 phi = phiC

554 if PI_2 > phi > -PI_2:

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

556# w = A0.log_tan(phi)

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

558 phi_ = phi

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

560 if fabs(phi - phi_) < _TOL_R:

561 break

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

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

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

565

566

567def _spherical2oblique(phiC, lamC): # 2.4.2

568 # oblique stereographic conformal projection

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

570 x = A0.X0 # 2.4.2g

571 y = A0.Y0 # 2.4.2h

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

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

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

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

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

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

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

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

580 if EPS0 < sp_22 < EPS1:

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

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

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

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

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

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

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

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

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

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

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

592 pass

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

594 x = y = NAN

595# else:

596# raise RDNAPError((phiC, lamC))

597 return x, y

598

599

600__all__ += _ALL_DOCS(_RDNAPbase)

601__all__ += _all_OTHER(RDNAP2018v1, RDNAP2018v2)

602del _ALL_DOCS, _all_OTHER

603

604# **) MIT License

605#

607#

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

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

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

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

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

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

614#

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

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

617#

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

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

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

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

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

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

624# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pyrdnap / rdnap2018.py: 94%

241 statements