Coverage for pyrdnap/rd0.py: 87%

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, # "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, floor, log, sin, sqrt, tan

26__all__ = ()

27__version__ = '26.05.08'

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(tuple(map(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 geodetic B{GRS80 (ETRS89)} 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_RO

247 def lam(self):

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

249 '''

250 return Lamd(self.lon) # PYCHOK lon

251

252 @Property_RO

253 def latlon(self):

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

255 '''

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

257

258 @Property_RO

259 def latlonheight(self):

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

261 '''

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

263

264 @Property_RO

265 def latlonheightdatum(self):

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

267 '''

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

269

270 @Property_RO

271 def phi(self):

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

273 '''

274 return Phid(self.lat) # PYCHOK lat

275

276 @Property_RO

277 def philam(self):

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

279 '''

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

281

282 @Property_RO

283 def philamheight(self):

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

285 '''

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

287

288 @Property_RO

289 def philamheightdatum(self):

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

291 '''

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

293

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

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

296

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

298

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

300 '''

301 _xinstanceof(Datum, datum2=datum2)

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

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

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

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

306

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

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

309

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

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

312

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

314

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

316 '''

317 _xsubclassof(_LLEB, LatLon=LatLon)

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

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

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

321

322 @Property_RO

323 def xy(self):

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

325 '''

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

327

328 @Property_RO

329 def xyz(self):

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

331 '''

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

333

334

335class RDregion4Tuple(_NamedTuple):

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

337 '''

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

339 _Units_ = ( Lat, Lon, Lat, Lon)

340

341 @Property_RO

342 def SW(self):

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

344 '''

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

346

347 @Property_RO

348 def NE(self):

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

350 '''

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

352

353

354def _c_f_N_f3(deg, deg_MIN, deg_INC):

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

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

357 N = (deg - deg_MIN) / deg_INC

358 f = floor(N)

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

360

361

362__all__ += _ALL_OTHER(RDNAP7Tuple, RDregion4Tuple,

363 # passed along from PyGeodesy

364 LatLon2Tuple, PhiLam2Tuple, Similarity, Vector2Tuple, Vector3Tuple)

365

366# **) MIT License

367#

369#

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

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

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

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

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

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

376#

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

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

379#

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

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

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

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

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

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

386# OTHER DEALINGS IN THE SOFTWARE.