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 C{.reverse} results have been formally validated to meet the

10C{RDNAPTRANS(tm)2018_v220627} requirements, transforming from and to ETRS89 (GRS80) points.

12The L{RDNAP2018v2.forward} and C{.reverse} results have been formally validated to meet the

13C{RDNAPTRANS(tm)2018_v220627} requirements, transforming from ETRS89 (GRS80) and to RD-Bessel

14(Bessel1841) points.

15'''

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

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

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

20from pyrdnap.v_grids import RDNAPError, _V_grid, _v_gridz_import

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

22 _isNAN, _isNAN0, _earth_datum,

23 _ALL_DOCS, _ALL_OTHER, _FOR_DOCS,

24 _NamedBase)

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

26 Datum, Ellipsoid, LatLonDatum3Tuple, # datums, ellipsoids

27 deprecated_property_RO, property_RO, property_ROnce, # props

28 Lamd, Lat, Lon, 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.06.14'

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

38_TOL_M = 1e-6 # meter

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

40_TRIPS = 16 # 5..6 sufficient

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

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

45 '''

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

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

48 _raiser = False

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

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

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

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

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

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

57 specified as C{scalar}, ignored otherwise.

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

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

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

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

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

64 '''

65 if a_ellipsoid is f is None:

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

67 else:

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

69 self._EETRS = E = self._datum.ellipsoid

70 if not E.isOblate:

71 raise RDNAPError('not oblate: %r' % (E,))

72 if raiser: # PYCHOK no cover

73 T = self._datum.transform

74 if not T.isunity:

75 raise RDNAPError('not unity: %r' % (T,))

76 self._raiser = True

77 if name:

78 self.name = name

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

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

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

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

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

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

87 to ignore C{NAPh} interpolation.

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

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

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

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

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

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

95 '''

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

97 lat0, lon0 = \

98 lat_, lon_ = self._forwardXform2(raiser, lat, lon)

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

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

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

102 fabs(lonc - lon_) < _TOL_D:

103 break

104 lat_, lon_ = latc, lonc

105

106 phiClamC = _ellipsoidal2spherical(latc, lonc)

107 RDx, RDy = _spherical2oblique(*phiClamC)

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

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

110 return RDNAP7Tuple(RDx, RDy, NAPh,

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

112

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

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

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

116

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

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

119 '''

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

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

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

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

124

125 @property_RO

126 def forwardDatum(self):

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

128 '''

129 return self._datum

130

131 def _forwardXform2(self, *args): # PYCHOK no cover

132 return self._notOverloaded(*args)

133

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

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

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

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

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

139 raise self._outsidError(lat, lon)

140 return lat, lon

141

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

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

144

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

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

147 (RD_Bessel) (C{bool}).

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

149 '''

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

151

152 def _outsidError(self, *lat_lon):

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

154 return RDNAPError('%r outside %r' % (lat_lon, self.region4()))

155

156 @property_RO

157 def _rdgrid(self): # PYCHOK no cover

158 return self._notOverloaded()

159

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

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

162 if _RD.isinside(lat, lon):

163 c_f_N_f6 = _RD._c_f_N_f6(lat, lon)

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

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

166

167 if lat0 is lon0 is None: # reverse

168 lat += lat_corr

169 lon += lon_corr

170 else: # forward

171 lat = lat0 - lat_corr

172 lon = lon0 - lon_corr

173 return lat, lon # NAN, NAN?

174

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

176 '''Interpolate the C{NAPh} quasi-geoid-height I{within} the C{RD} region.

177

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

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

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

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

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

183

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

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

186 '''

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

188

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

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

191 if _RD.isinside(lat, lon): # eps=0

192 c_f_N_f6 = _RD._c_f_N_f6(lat, lon)

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

194 elif raiser:

195 raise self._outsidError(lat, lon)

196 return NAN # c0 2.5.1e+

197

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

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

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

201 '''

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

203 self._rdNAPh(lat1, lon1, False)

204

205 @deprecated_property_RO

206 def region(self): # PYCHOK no cover

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

208 return self._region4RD

209

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

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

212

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

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

215

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

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

218 '''

219 return self._region4RD if asRD else self._region4ETRS

220

221 @property_ROnce

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

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

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

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

226 _ETRS_ = r.name.replace('RD', 'ETRS')

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

228

229 @property_ROnce

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

231 return _RD._region4RD

232

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

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

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

236 as RD-Bessel C{B{asRD}=True} or ETRS C{B{asRB}=False} or use

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

238 '''

239 phiClamC = _oblique2spherical(RDx, RDy)

240 latlon = _spherical2ellipsoidal(*phiClamC)

241

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

243 lat, lon, d = self._reverseXform3(raiser, latc, lonc)

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

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

246

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

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

249 return RDNAP7Tuple(RDx, RDy, NAPh,

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

251

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

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

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

255

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

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

258 '''

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

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

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

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

263

264 @property_RO

265 def reverseDatum(self):

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

267 '''

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

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

270

271 def _reverseXform3(self, *raiser_lat_lon):

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

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

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

275 return self.reverse3(lat, lon)

276

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

278 return self._notOverloaded(inverse=inverse)

279

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

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

282

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

284

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

286 '''

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

288

289 @property_RO

290 def variant(self): # PYCHOK no cover

291 return self._notOverloaded()

292

293

294class RDNAP2018v1(_RDNAPbase):

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

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

297

298 @note: Method L{RDNAP2018v2.reverse} returns B{by default GRS80 (ETRS89)}

299 validated, geodetic lat- and longitudes and datum.

300 '''

301 if _FOR_DOCS:

302 __init__ = _RDNAPbase.__init__

303 forward = _RDNAPbase.forward

304 forward3 = _RDNAPbase.forward3

305

306 def _forwardXform2(self, raiser, *lat_lon): # PYCHOK signature

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

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

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

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

311

312 if _FOR_DOCS:

313 isinside = _RDNAPbase.isinside

314 rdNAPh = _RDNAPbase.rdNAPh

315 region4 = _RDNAPbase.region4

316

317 @property_ROnce

318 def _rdgrid(self):

319 try:

320 from pyrdnap import v1grid

321 except (AttributeError, ImportError, RDNAPError):

322 v1grid = _v_gridz_import(self.variant)

323 return v1grid

324

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

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

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

328

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

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

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

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

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

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

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

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

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

338

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

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

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

342

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

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

345 '''

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

347

348 if _FOR_DOCS:

349 reverse3 = _RDNAPbase.reverse3

350

351 def similarity(self, inverse=False):

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

353

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

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

356 '''

357 return _RD._xRD2ETRS if inverse else _RD._xETRS2RD

358

359 @property_ROnce

360 def variant(self):

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

362 '''

363 return 1

364

365

366class RDNAP2018v2(_RDNAPbase):

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

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

369

370 @note: Method L{RDNAP2018v2.reverse} returns B{by default Bessel1841

371 (RD-Bessel)} validated, geodetic lat- and longitudes and datum.

372 '''

373 if _FOR_DOCS:

374 __init__ = _RDNAPbase.__init__

375 forward = _RDNAPbase.forward

376 forward3 = _RDNAPbase.forward3

377

378 def _forwardXform2(self, *raiser_lat_lon):

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

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

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

382

383 @property_ROnce

384 def _rdgrid(self):

385 try:

386 from pyrdnap import v2grid

387 except (AttributeError, ImportError, RDNAPError):

388 v2grid = _v_gridz_import(self.variant)

389 return v2grid

390

391 if _FOR_DOCS:

392 isinside = _RDNAPbase.isinside

393 rdNAPh = _RDNAPbase.rdNAPh

394 region4 = _RDNAPbase.region4

395

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

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

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

399

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

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

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

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

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

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

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

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

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

409

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

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

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

413

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

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

416 '''

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

418

419 if _FOR_DOCS:

420 reverse3 = _RDNAPbase.reverse3

421

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

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

424 '''

425 return None

426

427 @property_ROnce

428 def variant(self):

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

430 '''

431 return 2

432

433

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

435 # equiv to math.atan2 iff x0 is y

436 if x > 0:

437 r = atan(y / x)

438 elif x < 0:

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

440 else:

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

442 return r

443

444

445def _atan_exp(w): # 2.4.1c

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

447

448

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

450 c_lonI, f_lonI, lonN_f):

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

452 assert isinstance(v_grid, _V_grid), v_grid

453 ne = v_grid(c_latI, c_lonI)

454 nw = v_grid(c_latI, f_lonI)

455 se = v_grid(f_latI, c_lonI)

456 sw = v_grid(f_latI, f_lonI)

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

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

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

460

461

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

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

464 r = hypot(x, y)

465 if r > _TOL_M:

466 a = E.a * E.e2

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

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

469 s = sin(phi_)

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

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

472 if fabs(phi - phi_) < _TOL_R:

473 break

474 phi_ = phi

475 else:

476 phi = copysign(PI_2, z)

477 lam = _atan3(y, x, y)

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

479

480

481def _ellipsoidal2spherical(lat, lon): # 2.4.1

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

483 phiC = phi = Phid(lat)

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

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

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

487 phiC = _atan_exp(w)

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

489 return phiC, lamC # -Capital 𝛷, 𝛬

490

491

492def _eq0(r, r0=_0_0):

493 return fabs(r - r0) < _TOL_R

494

495

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

497# return fabs(d - d0) < _TOL_D

498

499

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

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

502 y, x = sincos2d(lon)

503 z, c = sincos2d(lat)

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

505 H = _isNAN0(h)

506 c *= n + H

507 x *= c

508 y *= c

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

510 return x, y, z

511

512

513def _ne0(r, r0=_0_0):

514 return fabs(r - r0) > _TOL_R

515

516

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

518# return fabs(d - d0) > _TOL_D

519

520

521def _oblique2spherical(x, y): # 3.1.1

522 # inverse oblique stereographic conformal projection from

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

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

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

526 x -= A0.X0

527 y -= A0.Y0

528 r = hypot(x, y)

529 if r > _TOL_M: # x and y

530 s0, c0 = A0.sincos2PHI0C

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

532 ca = sp * y / r

533 xN = cp * c0 - ca * s0

534 yN = sp * x / r

535 zN = cp * s0 + ca * c0

536 phiC = asin(zN)

537 else:

538 _, xN = A0.sincos2PHI0C

539 yN = _0_0

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

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

542 return phiC, lamC # -Capital 𝛷, 𝛬

543

544

545def _spherical2ellipsoidal(phiC, lamC): # 3.1.2

546 # inverse Gauss conformal projection from

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

548 phi = phiC

549 if PI_2 > phi > -PI_2:

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

551# w = A0.log_tan(phi)

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

553 phi_ = phi

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

555 if fabs(phi - phi_) < _TOL_R:

556 break

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

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

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

560

561

562def _spherical2oblique(phiC, lamC): # 2.4.2

563 # oblique stereographic conformal projection

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

565 x = A0.X0 # 2.4.2g

566 y = A0.Y0 # 2.4.2h

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

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

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

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

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

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

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

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

575 if EPS0 < sp_22 < EPS1:

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

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

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

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

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

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

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

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

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

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

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

587 pass

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

589 x = y = NAN

590# else:

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

592 return x, y

593

594

595__all__ += _ALL_DOCS(_RDNAPbase)

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

597 Datum, Ellipsoid, LatLonDatum3Tuple) # passed along from PyGeodesy

598

599# **) MIT License

600#

602#

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

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

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

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

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

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

609#

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

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

612#

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

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

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

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

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

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

619# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pyrdnap / rdnap2018.py: 94%

240 statements