Coverage for pyrdnap/rd0.py: 84%

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

4u'''(INTERNAL) C{RD} and C{RD0} constants and C{RDNAP7Tuple} and C{RDregion} classes.

5'''

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

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

9from pyrdnap.v_grids import _v_assert

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

11 _isNAN, _xinstanceof, _xsubclassof,

12 _LLEB, _xkwds,

13 _COMMASPACE_, _datum_, _E_, _lat_, _lon_,

14 _height_, _N_, _S_, _UNDER_, _W_,

15 _ALL_OTHER, _Pass, _NamedTuple)

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

17 Datum, Datums, Similarity, # datums

18 LatLon2Tuple, PhiLam2Tuple, Vector2Tuple, Vector3Tuple, # namedTuples

19 Property_RO, property_ROnce, # props

20 pairs, # streprs

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

22 sincos2) # utily

24from math import atan, ceil, fabs, floor, log, sin, sqrt, tan

26__all__ = ()

27__version__ = '26.05.14'

30class _RDbase(object):

31 '''(INTERNAL) Base.

32 '''

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

34 # return updated dict C{d}

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

36 if _pred(n):

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

38 return d

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

41 # return this C{_RDx} as string

42 d = self._toDict() # PYCHOK OK

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

44 return _COMMASPACE_(*t)

47class _RD(_RDbase):

48 '''(INTERNAL) Limits, constants for RDNAP2018 (ASCII.txt).

49 '''

50 LAT_INC = Lat(LAT_INC= 0.0125) # degrees, all

51 LAT_MAX = Lat(LAT_MAX=56.0)

52 LAT_MIN = Lat(LAT_MIN=50.0)

53 LON_INC = Lon(LON_INC= 0.02)

54 LON_MAX = Lon(LON_MAX= 8.0)

55 LON_MIN = Lon(LON_MIN=_2_0)

56 N_LAT_LON = ((LAT_MAX - LAT_MIN) / LAT_INC + _1_0, # 2.3.2g n-phi

57 (LON_MAX - LON_MIN) / LON_INC + _1_0) # 2.3.2g n-lambda

59 def __init__(self):

60 _v_assert(map2(int, self.N_LAT_LON))

62 def c_f_N_f6(self, lat, lon):

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

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

65 return _c_f_N_f3(lat, self.LAT_MIN, self.LAT_INC) + \

66 _c_f_N_f3(lon, self.LON_MIN, self.LON_INC)

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

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

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

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

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

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

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

76 @property_ROnce

77 def region(self):

78 t = RDregion4Tuple(self.LAT_MIN, self.LON_MIN,

79 self.LAT_MAX, self.LON_MAX, name='RD region ')

80 assert t.S < t.N and t.W < t.E, t.name

81 return t

83 def _toDict(self):

84 def _p(n): # lambda

85 return n.replace(_UNDER_, NN).isupper()

87 return self._preDict(_p, _xETRS2RD=self._xETRS2RD,

88 _xRD2ETRS=self._xRD2ETRS)

90 @property_ROnce

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

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

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

95 @property_ROnce

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

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

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

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

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

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

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

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

105 # LAT_INC=0.0125, LAT_MAX=56.0, LAT_MIN=50.0, LON_INC=0.02, LON_MAX=8.0, LON_MIN=2.0,

106 # N_LAT_LON=(481.0, 301.0)

107

108_RD = _RD() # PYCHOK singleton

109

110

111class _RD0(_RDbase):

112 '''(INTERNAL) C{RD} Amersfoort, NL references.

113

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

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

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

117 '''

118 H0 = Meter(H0=_0_0) # E0 height

119 H0_ETRS = Meter(H0_ETRS=43.0)

120 K0 = 0.9999079 # scale factor

121 LAT0 = Lat(LAT0='52 9 22.178N') # 52.15616055+°

122 LON0 = Lon(LON0=' 5 23 15.5E') # 5.387638888+°

123 LAM0C = \

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

125 PHI0 = Phid(PHI0=LAT0) # 𝜑0 0.910296727, 𝛷0 below

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

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

128

129# @property_ROnce

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

131# s, _ = self.sincos2PHI0

132# w = self._w1(s)

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

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

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

136

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

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

139# s, _ = sincos2d(lat)

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

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

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

143# return asin(s), r

144

145 @property_ROnce

146 def D0(self): # lazily

147 return Datums.Bessel1841

148

149 @property_ROnce

150 def E0(self): # lazily

151 return self.D0.ellipsoid

152

153 def ln_e_2(self, phi):

154 e = self.E0.e

155 p = e * sin(phi)

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

157

158 def ln_tan(self, phi):

159 return log(tan((phi + PI_2) * _0_5))

160

161 @property_ROnce

162 def M0(self): # 2.4.1+ m

163 q0 = self.ln_tan(self.PHI0) - self.ln_e_2(self.PHI0)

164 w0 = self.ln_tan(self.PHI0C)

165 return w0 - self.N0 * q0

166

167 @property_ROnce

168 def N0(self): # 2.4.1+ n, sphere

169 E = self.E0

170 _, c = self.sincos2PHI0

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

172

173 @property_ROnce

174 def PHI0C(self): # 2.4.1+

175 # get 𝛷0, Amersfoort latitude on sphere

176 rM, rN = self._rMN2

177 return Phi(PHI0C=atan((sqrt(rM) / sqrt(rN)) * tan(self.PHI0)))

178

179 @property_ROnce

180 def R(self): # radius conformal sphere

181 rM, rN = self._rMN2

182 return sqrt(rM * rN)

183

184 @property_ROnce

185 def RK2(self): # 2.4.1+

186 return self.R * self.K0 * _2_0

187

188 @property_ROnce

189 def _rMN2(self): # 2.4.1+

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

191 E = self.E0

192 s, _ = self.sincos2PHI0

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

194 rN = E.a / sqrt(s)

195 rM = E.e21 * rN / s

196 return rM, rN

197

198 @property_ROnce

199 def sincos2PHI0(self):

200 return sincos2(self.PHI0)

201

202 @property_ROnce

203 def sincos2PHI0C(self):

204 return sincos2(self.PHI0C)

205

206 def _toDict(self):

207 def _p(n): # lambda

208 return n.endswith('0') or n.endswith('0C') # _0_

209

210 return self._preDict(_p, H0_ETRS=self.H0_ETRS, R=self.R,

211 RK2=self.RK2, _rMN2=self._rMN2)

212

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

214# w1 = NAN

215# if _1_0 > sphi > _N_1_0:

216# e = self.E0.e

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

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

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

220# return w1

221

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

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

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

225 # H0=0.0, H0_ETRS=43.0, K0=0.9999079, LAM0=0.094032038, LAM0C=0.094032038, LAT0=52.156160556, LON0=5.387638889,

226 # M0=0.003773954, N0=1.000475857, PHI0=0.910296727, PHI0C=0.909684757, R=6382644.571035412, RK2=12764113.458940839,

227 # sincos2PHI0=(0.7896858198001045, 0.6135114554811807), sincos2PHI0C=(0.7893102212553742, 0.6139946047171687),

228 # X0=155000.0, Y0=463000.0, _rMN2=(6374588.709792872, 6390710.612840701)

229

230_RD0 = _RD0() # PYCHOK singleton, in .test

231

232

233class RDNAP7Tuple(_NamedTuple): # in .v_self

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

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

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

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

238

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

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

241 L{RDNAP2018v2.reverse}.

242 '''

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

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

245

246 @property_ROnce

247 def _datum_index(self):

248 return self._Names_.index(_datum_)

249

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

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

252

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

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

255

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

257 C{datum=B{datum}}.

258 '''

259 def _diff(a, b):

260 return fabs(a - b)

261

262 _xinstanceof(RDNAP7Tuple, other=other)

263 d = self._datum_index

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

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

266

267 @Property_RO

268 def lam(self):

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

270 '''

271 return Lamd(self.lon) # PYCHOK lon

272

273 @Property_RO

274 def latlon(self):

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

276 '''

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

278

279 @Property_RO

280 def latlonheight(self):

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

282 '''

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

284

285 @Property_RO

286 def latlonheightdatum(self):

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

288 '''

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

290

291 @Property_RO

292 def phi(self):

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

294 '''

295 return Phid(self.lat) # PYCHOK lat

296

297 @Property_RO

298 def philam(self):

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

300 '''

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

302

303 @Property_RO

304 def philamheight(self):

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

306 '''

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

308

309 @Property_RO

310 def philamheightdatum(self):

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

312 '''

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

314

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

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

317

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

319

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

321 '''

322 _xinstanceof(Datum, datum2=datum2)

323 h = self.height # PYCHOK preserve height NAN, because Ecef._forward ...

324 g = self.toLatLon(_LLEB).toDatum(datum2) # ... treats height NAN as 0

325 return self.dup(lat=g.lat, lon=g.lon, datum=g.datum,

326 height=h if _isNAN(h) else g.height, **name)

327

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

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

330

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

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

333

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

335

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

337 '''

338 _xsubclassof(_LLEB, LatLon=LatLon)

339 h = self.height # PYCHOK treat height NAN as 0, like Ecef._forward

340 kwds = _xkwds(LatLon_kwds, name=self.name, height=_0_0 if _isNAN(h) else h)

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

342

343 @Property_RO

344 def xy(self):

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

346 '''

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

348

349 @Property_RO

350 def xyz(self):

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

352 '''

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

354

355

356class RDregion4Tuple(_NamedTuple):

357 '''4-Tuple C{(S, W, N, E)} with C{RD} region in C{GRS80 (ETRS89) degrees}.

358 '''

359 _Names_ = (_S_, _W_, _N_, _E_)

360 _Units_ = ( Lat, Lon, Lat, Lon)

361

362 @Property_RO

363 def SW(self):

364 '''Get the C{SW} corner as (L{LatLon2Tuple}C{(lat, lon)}).

365 '''

366 return LatLon2Tuple(self.S, self.W, name=self.name)

367

368 @Property_RO

369 def NE(self):

370 '''Get the C{NE} corner as (L{LatLon2Tuple}C{(lat, lon)}).

371 '''

372 return LatLon2Tuple(self.N, self.E, name=self.name)

373

374

375def _c_f_N_f3(deg, deg_MIN, deg_INC):

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

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

378 N = (deg - deg_MIN) / deg_INC

379 f = floor(N)

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

381

382

383__all__ += _ALL_OTHER(RDNAP7Tuple, RDregion4Tuple,

384 # passed along from PyGeodesy

385 LatLon2Tuple, PhiLam2Tuple, Similarity, Vector2Tuple, Vector3Tuple)

386

387# **) MIT License

388#

390#

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

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

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

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

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

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

397#

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

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

400#

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

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

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

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

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

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

407# OTHER DEALINGS IN THE SOFTWARE.