Coverage for pyrdnap/rd0.py: 99%

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

4u'''(INTERNAL) RijksDriehoeksmeting C{_RD} and reference C{_RD0}

5constants and classes C{RDNAP7Tuple} and C{LqRD}.

6'''

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

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

10from pyrdnap.v_grids import _v_assert

11from pyrdnap.__pygeodesy import (_0_0, _0_5, _1_0, _2_0, # PYCHOK used!

12 _isNAN, _isNAN0, _xinstanceof, _xsubclassof,

13 _LLEB, _xkwds,

14 _COMMASPACE_, _datum_, _lat_, _lon_, _height_,

15 _ALL_OTHER, _FOR_DOCS, _Pass, _NamedTuple)

16from pygeodesy import (NAN, NN, map1, map2, # basics, "consterns"

17 Datum, Datums, Similarity, # datums

18 Bounds4Tuple, LatLon2Tuple, PhiLam2Tuple, # namedTuples

19 Vector2Tuple, Vector3Tuple, LqRD as _LqRD, # ltp

20 Property_RO, property_ROnce, # props

21 pairs, # streprs

22 Height, Lamd, Lat, Lon, Meter, Phi, Phid, # units

23 sincos2, tanPI_2_2) # utily

25from math import atan2, ceil, fabs, floor, log, sin, sqrt

27__all__ = ()

28__version__ = '26.05.24'

30_LQRD = _LqRD() # get Amersfoort, bounds, etc. (deleted below)

33def _c_f_N_f3(*deg_SW_D):

34 # return int(ceil) and int(floor) of Normalized

35 # and (Normalized less floor) of C{deg} degrees

36 N = _degN(*deg_SW_D)

37 # assert N >= 0, N

38 f = floor(N)

39 return int(ceil(N)), int(f), (N - f)

42def _degN(deg, degSW, degD):

43 # return C{deg} Normalized

44 return (deg - degSW) / degD

47class _RDbase(object):

48 '''(INTERNAL) Base.

49 '''

50 def _preDict(self, _pred, **d):

51 # return updated dict C{d}

52 for n in self.__class__.__dict__.keys():

53 if _pred(n):

54 d[n] = getattr(self, n)

55 return d

57 def toStr(self, prec=9, **fmt_ints):

58 # return this C{_RDx} as string

59 d = self._toDict() # PYCHOK OK

60 t = pairs(d, prec=prec, **fmt_ints)

61 return _COMMASPACE_(*t)

64class _RD(_RDbase):

65 '''(INTERNAL) Bounds, constants for RDNAP2018 (ASCII.txt).

66 '''

67 latD = Lat(latD=0.0125) # degrees, all

68 lonD = Lon(lonD=0.02)

70 def __init__(self):

71 S, W, N, E = self.region

72 nlat = _degN(N, S, self.latD) + _1_0 # 2.3.2g n-phi

73 nlon = _degN(E, W, self.lonD) + _1_0 # 2.3.2g n-lambda

74 _v_assert(map1(int, nlat, nlon))

76 def _c_f_N_f6(self, lat, lon):

77 # return (int(ceil), int(floor), Normalized less floor) of C{lat}) + \

78 # (int(ceil), int(floor), Normalized less floor) of C{lon})

79 return _c_f_N_f3(lat, self.region.latS, self.latD) + \

80 _c_f_N_f3(lon, self.region.lonW, self.lonD)

82 def isinside(self, lat, lon, eps=0): # eps=_TOL_D, 0 or -_TOLD_D

83 # is C{(lat, lon)} inside the this C{RD} region, optionally

84 # over-/undersized by positive respectively negative C{eps}?

85 S, W, N, E = self.region

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

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

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

90 region = _LQRD.region # as L{Bounds4Tuple}

92 def _toDict(self):

93 def _p(n): # lambda

94 return n.endswith('D') or n.endswith('S')

96 return self._preDict(_p, region=self.region)

98 @property_ROnce

99 def _xETRS2RD(self): # transform ETRS (GRS80) to RD-Bessel

100 return Similarity(tx=-565.7346, ty=-50.4058, tz=-465.2895, s=-4.07242,

101 rx=-1.91513, ry=1.60365, rz=-9.09546, name='_xETRS2RD')

102

103 @property_ROnce

104 def _xRD2ETRS(self): # transform RD-Bessel to ETRS (GRS80)

105 return Similarity(tx=565.7381, ty=50.4018, tz=465.2904, s=4.07244,

106 rx=1.91514, ry=-1.60363, rz=9.09546, name='_xRD2ETRS')

107

108 # % python -c "import pyrdnap; print(pyrdnap.rd0._RD.toStr())"

109 # _xETRS2RD=Similarity(name='_xETRS2RD', tx=-565.73, ty=-50.406, tz=-465.29, s=-4.0724,

110 # rx=-1.9151, ry=1.6037, rz=-9.0955),

111 # _xRD2ETRS=Similarity(name='_xRD2ETRS', tx=565.74, ty=50.402, tz=465.29, s=4.0724,

112 # rx=1.9151, ry=-1.6036, rz=9.0955),

113 # latD=0.0125, lonD=0.02, region=RD region (latS=50.0, lonW=2.0, latN=56.0, lonE=8.0)

114

115_RD = _RD() # PYCHOK singleton, in .test/testRndTrips

116

117

118class _RD0(_RDbase):

119 '''(INTERNAL) C{RD} Amersfoort, NL constancts for RDNAP2018 (ASCII.txt).

120

121 @see: U{EPSG:9809<https://EPSG.io/9809-method>}, U{"Oblique Stereographic"

122 <https://PROJ.org/en/stable/operations/projections/sterea.html>} and

123 <http://geotiff.maptools.org/proj_list/oblique_stereographic.html>

124 '''

125 H0 = Meter(H0 =_LQRD.height0) # Amersfoort.height0 0.0 m

126 H0_ETRS = Meter(H0_ETRS=_LQRD.height0_ETRS) # 43.0 m

127 K0 = 0.9999079 # 2.4.1 scale factor

128 LAT0 = Lat(LAT0=_LQRD.Amersfoort.lat) # '52 9 22.178N' == 52.15616055+°

129 LON0 = Lon(LON0=_LQRD.Amersfoort.lon) # ' 5 23 15.5E' == 5.387638888+°

130 LAM0C = \

131 LAM0 = Lamd(LAM0=LON0) # 𝜆0, 𝛬0 = 𝜆0 on sphere 0.094032038

132 PHI0 = Phid(PHI0=LAT0) # 𝜑0 0.910296727, PHI0C 𝛷0 set below

133 X0 = Meter(X0=155000.0) # false Easting 155029.784?

134 Y0 = Meter(Y0=463000.0) # false Norting 463109.889?

135

136# @property_ROnce

137# def C0(self): # c, sphere

138# s, _ = self.sincos2PHI0

139# w = self._w1(s)

140# c = (w - _1_0) / (w + _1_0)

141# return (((self.N0 + s) * (_1_0 - c)) /

142# ((self.N0 - s) * (_1_0 + c)))

143

144# def chilam(self, lat, lon): # EPSG:9809

145# # return 2-tuple (chi, lam), conformal in radians

146# s, _ = sincos2d(lat)

147# w2 = self._w1(s) * self.C0

148# s = (w2 - _1_0) / (w2 + _1_0)

149# r = radians(lon - self.LON0) * self.N0

150# return asin(s), r

151

152 @property_ROnce

153 def D0(self): # lazily

154 return Datums.Bessel1841

155

156 @property_ROnce

157 def E0(self): # lazily

158 return self.D0.ellipsoid

159

160 def log_e_2(self, phi):

161 e = self.E0.e

162 p = e * sin(phi)

163 return log((_1_0 + p) / (_1_0 - p)) * (e * _0_5)

164

165 def log_tan(self, phi):

166 return log(tanPI_2_2(phi)) # tan((phi + PI/2) / 2)

167

168 @property_ROnce

169 def M0(self): # 2.4.1 p 15 m

170 return self.W0 - self.N0 * self.Q0

171

172 @property_ROnce

173 def N0(self): # 2.4.1 p 15 n, sphere

174 E = self.E0

175 _, c = self.sincos2PHI0

176 return sqrt(c**4 * E.e2 / E.e21 + _1_0)

177

178 @property_ROnce

179 def PHI0C(self): # 2.4.1 p 15 𝛷0

180 # get 𝛷0, Amersfoort latitude on sphere

181 m, n = self.Rmn2

182 s, c = self.sincos2PHI0

183 return Phi(PHI0C=atan2(m * s, n * c)) # atan((m / n) * tan(PHI0))

184

185 @property_ROnce

186 def Q0(self): # 2.4.1 p 15 q0

187 return self.log_tan(self.PHI0) - self.log_e_2(self.PHI0)

188

189 @property_ROnce

190 def R(self): # 2.4.1 p 15 R, radius conformal sphere

191 m, n = self.Rmn2

192 return m * n

193

194 @property_ROnce

195 def RK2(self): # 2.4.2

196 return self.R * self.K0 * _2_0

197

198 @property_ROnce

199 def Rmn2(self): # 2.4.1 p 15 sqrt(RsubM), sqrt(RsubN)

200 # get 2-tuple (RHO0, NU0) EPSG:9809

201 E = self.E0

202 s, _ = self.sincos2PHI0

203 s = _1_0 - s**2 * E.e2

204 # assert s > 0

205 N = E.a / sqrt(s)

206 # assert N > 0

207 M = E.e21 * N / s

208 # assert M > 0

209 return map1(sqrt, M, N) # sqrt!

210

211 @property_ROnce

212 def sincos2PHI0(self): # 𝜑0

213 return sincos2(self.PHI0)

214

215 @property_ROnce

216 def sincos2PHI0C(self): # 𝛷0

217 return sincos2(self.PHI0C)

218

219 def _toDict(self):

220 def _p(n): # lambda

221 return n.endswith('0') or n.startswith('R') or \

222 n.endswith('0C') # _0_

223

224 return self._preDict(_p, H0_ETRS=self.H0_ETRS)

225

226 @property_ROnce

227 def W0(self): # 2.4.1 p 15 w0

228 return self.log_tan(self.PHI0C) # 𝛷0

229

230# def _w1(self, sphi): # EPSG:9809

231# w1 = NAN

232# if _1_0 > sphi > _N_1_0:

233# e = self.E0.e

234# S = (_1_0 + sphi) / (_1_0 - sphi)

235# T = (_1_0 - sphi * e) / (_1_0 + sphi * e)

236# w1 = pow(pow(T, e) * S, self.N0)

237# return w1

238

239 # % python -c "import pyrdnap; print(pyrdnap.rd0._RD0.toStr())"

240 # D0=Datum(name='Bessel1841', ellipsoid=Ellipsoids.Bessel1841, transform=Transforms.Bessel1841),

241 # E0=Ellipsoid(name='Bessel1841', a=6377397.155, f=0.00334277, f_=299.1528128, b=6356078.962818),

242 # H0=0.0, H0_ETRS=43.0, K0=0.9999079, LAM0=0.094032038, LAM0C=0.094032038,

243 # LAT0=52.156160556, LON0=5.387638889, M0=0.003773954, N0=1.000475857,

244 # PHI0=0.910296727, PHI0C=0.909684757, Q0=1.06531844,

245 # R=6382644.571035411, RK2=12764113.458940838, Rmn2=(2524.794785679199, 2527.9854850929623),

246 # sincos2PHI0=(0.7896858198001045, 0.6135114554811807), sincos2PHI0C=(0.7893102212553742, 0.6139946047171686),

247 # W0=1.069599332, X0=155000.0, Y0=463000.0

248

249_RD0 = _RD0() # PYCHOK singleton, in .test/testRndTrips

250

251

252class RDNAP7Tuple(_NamedTuple): # in .v_self

253 '''7-Tuple C{(RDx, RDy, NAPh, lat, lon, height, datum)} with I{local} C{RDx},

254 C{RDy} and C{NAPh} quasi-geoid_height, geodetic C{lat}, C{lon}, C{height}

255 and C{datum} with C{lat} and C{lon} in C{degrees} and C{RDx}, C{RDy}, C{NAPh}

256 and C{height} in C{meter}, conventionally.

257

258 @note: The C{lat} and {lon} are B{GRS80 (ETRS89)} geodetic coordinates from

259 L{RDNAP2018v1.reverse} but B{Bessel1841 (RD-Bessel)} when returned from

260 L{RDNAP2018v2.reverse}.

261 '''

262 _Names_ = ('RDx', 'RDy', 'NAPh', _lat_, _lon_, _height_, _datum_)

263 _Units_ = ( Meter, Meter, Meter, Lat, Lon, Height, _Pass)

264

265 @property_ROnce

266 def _datum_index(self):

267 return self._Names_.index(_datum_)

268

269 def diff(self, other, datum=None, **name):

270 '''Return the difference between this and an C{other} C{RDNAP7Tuple}.

271

272 @kwarg datum: Optional difference C{B{datum}=None} (C{Latum}).

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

274

275 @return: An L{RDNAP7Tuple} with the C{_diff} for each item, but

276 C{datum=B{datum}}.

277 '''

278 def _diff(a, b):

279 return fabs(a - b)

280

281 _xinstanceof(RDNAP7Tuple, other=other)

282 d = self._datum_index

283 t = map2(_diff, self[:d], other[:d])

284 return RDNAP7Tuple(t + (datum,), **name)

285

286 @Property_RO

287 def lam(self):

288 '''Get the longitude (B{C{radians}}).

289 '''

290 return Lamd(self.lon) # PYCHOK lon

291

292 @Property_RO

293 def latlon(self):

294 '''Get the lat-, longitude in C{degrees} (L{LatLon2Tuple}C{(lat, lon)}).

295 '''

296 return LatLon2Tuple(self.lat, self.lon, name=self.name)

297

298 @Property_RO

299 def latlonheight(self):

300 '''Get the lat-, longitude in C{degrees} and height (L{LatLon3Tuple}C{(lat, lon, height)}).

301 '''

302 return self.latlon.to3Tuple(self.height)

303

304 @Property_RO

305 def latlonheightdatum(self):

306 '''Get the lat-, longitude in C{degrees} with height and datum (L{LatLon4Tuple}C{(lat, lon, height, datum)}).

307 '''

308 return self.latlonheight.to4Tuple(self.datum)

309

310 @Property_RO

311 def phi(self):

312 '''Get the latitude (B{C{radians}}).

313 '''

314 return Phid(self.lat) # PYCHOK lat

315

316 @Property_RO

317 def philam(self):

318 '''Get the lat- and longitude in C{radians} (L{PhiLam2Tuple}C{(phi, lam)}).

319 '''

320 return PhiLam2Tuple(self.phi, self.lam, name=self.name) # PYCHOK lam, phi

321

322 @Property_RO

323 def philamheight(self):

324 '''Get the lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}).

325 '''

326 return self.philam.to3Tuple(self.height) # PYCHOK height

327

328 @Property_RO

329 def philamheightdatum(self):

330 '''Get the lat-, longitude in C{radians} with height and datum (L{PhiLamn4Tuple}C{(phi, lam, height, datum)}).

331 '''

332 return self.philamheight.to4Tuple(self.datum)

333

334 def toDatum(self, datum2, **name):

335 '''Convert this C{lat}, C{lon} and C{height} to B{C{datum2}}.

336

337 @arg datum2: Datum to convert I{to} (L{Datum}).

338

339 @return: An L{RDNAP7Tuple} with converted C{lat}, C{lon} and C{height}.

340 '''

341 _xinstanceof(Datum, datum2=datum2)

342 g = self.toLatLon(_LLEB).toDatum(datum2)

343 h = NAN if _isNAN(self.height) else g.height # PYCHOK preserve height NAN

344 return self.dup(lat=g.lat, lon=g.lon, datum=g.datum, height=h,

345 **_xkwds(name, name=self.name))

346

347 def toLatLon(self, LatLon, **LatLon_kwds):

348 '''Return this C{lat}, C{lon}, C{datum} and C{height} as B{C{LatLon}}.

349

350 @arg LatLon: An ellipsodial C{LatLon} class (C{pygeodesy.ellipsoidal*}).

351 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword arguments.

352

353 @return: An B{C{LatLon}} instance.

354

355 @raise TypeError: B{C{LatLon}} not ellipsoidal or an other issue.

356 '''

357 _xsubclassof(_LLEB, LatLon=LatLon)

358 h = _isNAN0(self.height) # PYCHOK height

359 kwds = _xkwds(LatLon_kwds, name=self.name, height=h)

360 return LatLon(self.lat, self.lon, datum=self.datum, **kwds) # PYCHOK datum

361

362 @Property_RO

363 def xy(self):

364 '''Get the I{local} C{(RDx, RDy)} coordinates (L{Vector2Tuple}C{(x, y)}).

365 '''

366 return Vector2Tuple(self.RDx, self.RDy, name=self.name)

367

368 @Property_RO

369 def xyz(self):

370 '''Get the I{local} C{(RDx, RDy, NAPh)} coordinates and height (L{Vector3Tuple}C{(x, y, z)}).

371 '''

372 return Vector3Tuple(self.RDx, self.RDy, self.NAPh, name=self.name)

373

374

375class LqRD(_LqRD):

376 '''Like U{pygeodesy.LqRD<https://mrJean1.GitHub.io/PyGeodesy/docs/pygeodesy.ltp.LqRD-class.html>}

377 but with methods C{forward} and C{reverse} returning an L{RDNAP7Tuple} with C{NAPh} replaced

378 by I{local} C{z}, the perpendicular distance to the local tangent plane (LTP).

379

380 This C{quasi-RD} transformer B{does not} implement any U{RD NAP<https://www.NSGI.NL/

381 coordinatenstelsels-en-transformaties/coordinatentransformaties/rdnap-etrs89-rdnaptrans>}

382 specification and B{does not} provide I{Netherlands}' C{B{N}ormaal B{A}msterdams B{P}eil

383 (NAP)} quasi-geodetic-height.

384 '''

385 if _FOR_DOCS:

386 __init__ = _LqRD.__init__

387

388 def forward(self, lat_latlonh, lon=None, height=0, **name): # PYCHOK signature

389 '''Convert I{geodetic} C{(lat, lon, height)} to I{local} C{quasi-RD (x, y, z)}.

390

391 @arg lat_latlonh: C{Scalar} (geodetic) latitude (C{degrees}) or a I{local}

392 C{quasi-RD} L{RDNAP7Tuple}.

393 @kwarg lon: C{Scalar} (geodetic) longitude (C{degrees}) iff B{C{lat_latlonh}}

394 is C{scalar}, ignored otherwise.

395 @kwarg height: Optional height (C{meter}, conventionally) perpendicular to and

396 above (or below) the ellipsoid's surface, iff B{C{lat_latlonh}}

397 is C{scalar}, ignored otherwise.

398 @kwarg name: Optional C{B{name}=NN} (C{str}).

399

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

401 C{NAPh} set to I{local} C{z}.

402

403 @see: B{pygeodesy.LqRD.forward} for more information.

404 '''

405 t = _LqRD.forward(self, lat_latlonh, lon=lon, height=height)

406 return LqRD._l9t2r7t(t, **name)

407

408 def reverse(self, x_xyz, y=None, z=None, **name): # PYCHOK signature

409 '''Convert I{local} C{quasi-RD (x, y, z)} to I{geodetic} C{(lat, lon, height)}.

410

411 @arg x_xyz: Local C{quasi-RD x} coordinate (C{scalar}) or a I{local}

412 C{quasi-RD} L{RDNAP7Tuple}.

413 @kwarg y: Local C{quasi-RD y} coordinate (C{meter}) iff B{C{x_xyz}} is

414 C{scalar}, ignored otherwise.

415 @kwarg z: Local C{z} coordinate (C{meter}) iff B{C{x_xyz}} is C{scalar},

416 ignored otherwise.

417 @kwarg name: Optional C{B{name}=NN} (C{str}).

418

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

420 with C{NAPh} set to I{local} B{C{z}}.

421

422 @see: B{pygeodesy.LqRD.reverse} for more information.

423 '''

424 t = _LqRD.reverse(self, x_xyz, y=y, z=z)

425 return LqRD._l9t2r7t(t, **name)

426

427 @staticmethod

428 def _l9t2r7t(t, name=NN, **unused): # M=False

429 return RDNAP7Tuple(t.x, t.y, t.z, # NAPh = t.z

430 t.lat, t.lon, t.height, t.ecef.datum, name=name or t.name)

431

432

433__all__ += _ALL_OTHER(LqRD, RDNAP7Tuple, Bounds4Tuple, # passed along from PyGeodesy

434 LatLon2Tuple, PhiLam2Tuple, Similarity, Vector2Tuple, Vector3Tuple)

435del _LQRD

436

437# **) MIT License

438#

440#

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

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

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

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

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

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

447#

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

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

450#

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

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

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

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

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

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

457# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pyrdnap / rd0.py: 99%

181 statements