Coverage for pychlv/chlvs.py: 97%

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

4u'''Classes L{ChLV}, L{ChLVa}, L{ChLVe} and L{ChLVr} for conversion between

5(WGS84 or Bessel1841) geodetic lat-, longitude and height and U{SwissTopo

6<https://www.SwissTopo.admin.CH/en/the-swiss-coordinates-system>} I{LandesVermessung

71995 (LV95)} and I{-1903 (LV03)} local projections.

9Classes L{ChLVa} and L{ChLVe} are C{SwissTopo} I{approximations}, class L{ChLVr}

10implements the C{SwissTopo} I{rigorous formulas} and class L{ChLV} uses U{EcefKarney

11<https://MrJean1.GitHub.io/PyGeodesy/docs/pygeodesy.ecef-module.html>}'s Earth-Centered,

12Earth-Fixed (ECEF) methods.

13'''

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

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

17from pychlv.__pygeodesy import (ChLVError, _args_kwds_names, _isin, _over,

18 _0_01, _0_5, _1_0, _2_0, _60_0, _100_0, _3600_0,

19 _WGS84, _llhn4, _xyzn4,

20 _NotImplementedError, _xkwds, _floor,

21 _DOT_, _not_falsed_, _all_OTHER)

22from pychlv.chlvTuples import ChLVEN2Tuple, ChLV9Tuple, ChLVYX2Tuple, ChLVyx2Tuple

23from pygeodesy import (EcefKarney, EPS, Fdot_, fdot_, Fhorner, fsumf_, Lamd,

24 LatLon4Tuple, LocalCartesian, Ltp, map1, Phid, PI_2,

25 property_RO, property_ROver, sincos2_, tanPI_2_2,

26 typename, unstr) # Ellipsoids, sincos2, toDMS

28from math import asin, atan, atan2, degrees, exp, fabs, log, sin, sqrt

30__all__ = ()

31__version__ = '26.07.17'

33_height0_ = 'height0'

34_TRIPS = 16 # 6..7 sufficient

37class _ChLV(object):

38 '''(INTERNAL) Base class for C{ChLV*} classes.

39 '''

40 _03_falsing = ChLVyx2Tuple(0.6e6, 0.2e6)

41# _92_falsing = ChLVYX2Tuple(2.0e6, 1.0e6) # _95_ - _03_

42 _95_falsing = ChLVEN2Tuple(2.6e6, 1.2e6)

43 _9t0 = None # from LocalCartesian

45 def _ChLV9Tuple(self, fw, M, name, *Y_X_h_lat_lon_h):

46 '''(INTERNAL) Helper for C{ChLV*.forward} and C{.reverse}.

47 '''

48 if bool(M): # PYCHOK no cover

49 m = self.forward if fw else self.reverse # PYCHOK attr

50 n = _DOT_(*map1(typename, type(self), m))

51 raise _NotImplementedError(unstr(n, M=M), txt=None)

52 t = Y_X_h_lat_lon_h + (self, self._9t0, None)

53 return ChLV9Tuple(t, name=name)

55 @property_RO

56 def ellipsoid(self):

57 '''Get the ellipsoid (C{WGS84}).

58 '''

59 return _WGS84.ellipsoid

61 def forward(self, latlonh, lon=None, height=0, M=None, **name): # PYCHOK no cover

62 '''Convert (WGS84) geodetic to I{Swiss} projection coordinates.

64 @arg latlonh: Either a C{LatLon}, L{Ltp} or C{scalar} (geodetic) latitude (C{degrees}).

65 @kwarg lon: Optional, C{scalar} (geodetic) longitude (C{degrees}) iff B{C{latlonh}} is

66 C{scalar}, ignored otherwise.

67 @kwarg height: Optional, height, vertically above (or below) the surface of the ellipsoid

68 (C{meter}) iff B{C{latlonh}} and B{C{lon}} are C{scalar}, ignored otherwise.

69 @kwarg M: If C{True}, return the I{concatenated} rotation U{EcefMatrix<https://MrJean1.GitHub.io/

70 PyGeodesy/docs/pygeodesy.ecef-module.html>} iff available and for C{ChLV} only,

71 C{None} otherwise (C{bool}).

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

74 @return: A L{ChLV9Tuple}C{(Y, X, h_, lat, lon, height, ltp, ecef, M)} with the unfalsed

75 I{Swiss Y, X} coordinates, I{Swiss h_} height, the given I{geodetic} C{lat},

76 C{lon} and C{height}, C{ecef} (U{Ecef9Tuple<https://MrJean1.GitHub.io/PyGeodesy/

77 docs/pygeodesy.ecef-module.html>}) at I{Bern, Ch}, rotation matrix C{M} and

78 C{ltp} this C{ChLV}, C{ChLVa} or C{ChLVe} instance.

80 @raise ChLVError: Invalid or non-C{scalar} B{C{latlonh}}, B{C{lon}} or B{C{height}}.

81 '''

82 LocalCartesian._notOverloaded(self, latlonh, lon=lon, height=height, M=M, **name)

84 def reverse(self, enh_, n=None, h_=0, M=None, **name): # PYCHOK no cover

85 '''Convert I{Swiss} projection to (WGS84) geodetic coordinates.

87 @arg enh_: A Swiss projection (L{ChLV9Tuple}) or the C{scalar}, falsed I{Swiss E_LV95}

88 or I{y_LV03} easting (C{meter}).

89 @kwarg n: Falsed I{Swiss N_LV95} or I{x_LV03} northing (C{meter}) iff B{C{enh_}} is

90 C{scalar}, ignored otherwise.

91 @kwarg h_: I{Swiss h'} height (C{meter}) iff B{C{enh_}} and B{C{n}} are C{scalar},

92 ignored otherwise.

93 @kwarg M: If C{True}, return the I{concatenated} rotation U{EcefMatrix<https://MrJean1.GitHub.io/

94 PyGeodesy/docs/pygeodesy.ecef-module.html>} iff available and for C{ChLV} only,

95 C{None} otherwise (C{bool}).

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

98 @return: A L{ChLV9Tuple}C{(Y, X, h_, lat, lon, height, ltp, ecef, M)} with the unfalsed

99 I{Swiss Y, X} coordinates, I{Swiss h_} height, the given I{geodetic} C{lat},

100 C{lon} and C{height}, C{ecef} (U{Ecef9Tuple<https://MrJean1.GitHub.io/PyGeodesy/

101 docs/pygeodesy.ecef-module.html>}) at I{Bern, Ch}, rotation matrix C{M} and

102 C{ltp} this C{ChLV}, C{ChLVa} or C{ChLVe} instance.

103

104 @raise ChLVError: Invalid or non-C{scalar} B{C{enh_}}, B{C{n}} or B{C{h_}}.

105 '''

106 LocalCartesian._notOverloaded(self, enh_, n=n, h_=h_, M=M, **name)

107

108 @staticmethod

109 def _deg2ab(deg, sLL):

110 # convert degrees to arc-seconds

111 def _dms(ds, p, q, swap):

112 d = _floor(ds)

113 t = (ds - d) * p

114 m = _floor(t)

115 s = (t - m) * p

116 if swap:

117 d, s = s, d

118 return d + (m + s * q) * q

119

120 s = _dms(deg, _60_0, _0_01, False) # deg2sexag

121 s = _dms( s, _100_0, _60_0, True) # sexag2asec

122 return (s - sLL) / ChLV._s_ab

123

124 @property_ROver

125 def _enh_n_h(self): # PYCHOK no cover

126 # Get C{ChLV*.reverse} args[1:4] names, I{once}

127 t = _args_kwds_names(_ChLV.reverse)[1:4]

128 for C in (ChLV, ChLVa, ChLVe, ChLVr):

129 assert _args_kwds_names(C.reverse)[1:4] == t

130 return t # overwrite property_ROver

131

132 @staticmethod

133 def _falsing2(LV95):

134 '''(INTERNAL) Get the C{LV95} or C{LV03} falsing.

135 '''

136 return _ChLV._95_falsing if _isin(LV95, True, 95) else (

137 _ChLV._03_falsing if _isin(LV95, False, 3) else ChLVYX2Tuple(0, 0))

138

139 @staticmethod

140 def _llh2abh_3(lat, lon, h):

141 '''(INTERNAL) Helper for C{ChLV*.forward}.

142 '''

143 a = ChLV._deg2ab(lat, ChLV._sLat) # phi', lat_aux

144 b = ChLV._deg2ab(lon, ChLV._sLon) # lam', lng_aux

145 h_ = fsumf_(h, -ChLV.Bern.height, 2.73 * b, 6.94 * a)

146 return a, b, h_

147

148 @staticmethod

149 def _YXh_2abh3(Y, X, h_):

150 '''(INTERNAL) Helper for C{ChLV*.reverse}.

151 '''

152 def _YX2ab(YX):

153 return YX * ChLV._ab_m

154

155 a, b = map1(_YX2ab, Y, X)

156 h = fsumf_(h_, ChLV.Bern.height, -12.6 * a, -22.64 * b)

157 return a, b, h

158

159 @staticmethod

160 def _YXh_n4(enh_, n, h_, **name):

161 '''(INTERNAL) Helper for C{ChLV*.reverse}.

162 '''

163 Y, X, h_, name = _xyzn4(enh_, n, h_, (ChLV9Tuple,), Error=ChLVError,

164 _xyz_y_z_names=_ChLV._enh_n_h, **name)

165 if isinstance(enh_, ChLV9Tuple):

166 Y, X = enh_.Y, enh_.X

167 else: # isscalar(enh_)

168 Y, X = ChLV.unfalse2(Y, X) # PYCHOK ChLVYX2Tuple

169 return Y, X, h_, name

170

171

172class ChLV(_ChLV, Ltp):

173 '''Conversion between (WGS84) geodetic and I{Swiss} projection coordinates using U{PyGeodesy

174 EcefKarney<https://MrJean1.GitHub.io/PyGeodesy/docs/pygeodesy.ecef-module.html>} ECEF methods.

175

176 @see: U{Swiss projection formulas<https://www.SwissTopo.admin.CH/en/maps-data-online/

177 calculation-services.html>}, page 7ff, U{NAVREF<https://www.SwissTopo.admin.CH/en/

178 maps-data-online/calculation-services/navref.html>}, U{REFRAME<https://www.SwissTopo.admin.CH/

179 en/maps-data-online/calculation-services/reframe.html>} and U{SwissTopo Scripts GPS WGS84

180 <-> LV03<https://GitHub.com/ValentinMinder/Swisstopo-WGS84-LV03>}.

181 '''

182 _9Tuple = ChLV9Tuple

183

184 _ab_d = 0.36 # a, b units per degree, ...

185 _ab_m = 1.0e-6 # ... per meter and ...

186 _ab_M = _1_0 # ... per 1,000 Km or 1 Mm

187 _s_d = _3600_0 # arc-seconds per degree ...

188 _s_ab = _s_d / _ab_d # ... and per a, b unit

189 _sLat = 169028.66 # Bern, Ch in ...

190 _sLon = 26782.5 # ... arc-seconds ...

191 # lat, lon, height == 46°57'08.66", 7°26'22.50", 49.55m ("new" 46°57'07.89", 7°26'22.335")

192 Bern = LatLon4Tuple(_sLat / _s_d, _sLon / _s_d, 49.55, _WGS84, name='Bern')

193

194 def __init__(self, latlonh0=Bern, **other_Ltp_kwds):

195 '''New L{ChLV} I{ECEF-based} converter, centered at I{Bern, Ch}.

196

197 @kwarg latlonh0: The I{geodetic} origin and height, overriding C{Bern, Ch}.

198 @kwarg other_Ltp_kwds: Optional, other L{Ltp.__init__} keyword arguments.

199

200 @see: L{Ltp.__init__} for more information.

201 '''

202 Ltp.__init__(self, latlonh0, **_xkwds(other_Ltp_kwds, ecef=None, name=ChLV.Bern.name))

203

204 def forward(self, latlonh, lon=None, height=0, M=None, **name): # PYCHOK unused M

205 # overloaded for the _ChLV.forward.__doc__

206 return Ltp.forward(self, latlonh, lon=lon, height=height, M=M, **name)

207

208 def reverse(self, enh_, n=None, h_=0, M=None, **name): # PYCHOK signature

209 # overloaded for the _ChLV.reverse.__doc__

210 Y, X, h_, n = _ChLV._YXh_n4(enh_, n, h_, **name)

211 return Ltp.reverse(self, Y, X, h_, M=M, name=n)

212

213 @staticmethod

214 def false2(Y, X, LV95=True, **name):

215 '''Add the I{Swiss LV95} or I{LV03} falsing.

216

217 @arg Y: Unfalsed I{Swiss Y} easting (C{meter}).

218 @arg X: Unfalsed I{Swiss X} northing (C{meter}).

219 @kwarg LV95: If C{True}, add C{LV95} falsing, if C{False} add

220 C{LV03} falsing, otherwise leave unfalsed.

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

222

223 @return: A L{ChLVEN2Tuple}C{(E_LV95, N_LV95)} or a

224 L{ChLVyx2Tuple}C{(y_LV03, x_LV03)} with falsed B{C{Y}}

225 and B{C{X}}, otherwise a L{ChLVYX2Tuple}C{(Y, X)}

226 with B{C{Y}} and B{C{X}} as-is.

227 '''

228 e, n = t = _ChLV._falsing2(LV95)

229 return t.classof(e + Y, n + X, **name)

230

231 @staticmethod

232 def isLV03(e, n):

233 '''Is C{(B{e}, B{n})} a valid I{Swiss LV03} projection?

234

235 @arg e: Falsed (or unfalsed) I{Swiss} easting (C{meter}).

236 @arg n: Falsed (or unfalsed) I{Swiss} northing (C{meter}).

237

238 @return: C{True} if C{(B{e}, B{n})} is a valid, falsed I{Swiss

239 LV03}, projection C{False} otherwise.

240 '''

241 # @see: U{Map<https://www.SwissTopo.admin.CH/en/knowledge-facts/

242 # surveying-geodesy/reference-frames/local/lv95.html>}

243 return 400.0e3 < e < 900.0e3 and 40.0e3 < n < 400.0e3

244

245 @staticmethod

246 def isLV95(e, n, raiser=True):

247 '''Is C{(B{e}, B{n})} a valid I{Swiss LV95} or I{LV03} projection?

248

249 @arg e: Falsed (or unfalsed) I{Swiss} easting (C{meter}).

250 @arg n: Falsed (or unfalsed) I{Swiss} northing (C{meter}).

251 @kwarg raiser: If C{True}, throw a L{ChLVError} if B{C{e}} and

252 B{C{n}} are invalid I{Swiss LV95} nor I{LV03}.

253

254 @return: C{True} or C{False} if C{(B{e}, B{n})} is a valid I{Swiss

255 LV95} respectively I{LV03} projection, C{None} otherwise.

256 '''

257 if ChLV.isLV03(e, n):

258 return False

259 elif ChLV.isLV03(e - 2.0e6, n - 1.0e6): # _92_falsing = _95_ - _03_

260 return True

261 elif raiser: # PYCHOK no cover

262 raise ChLVError(unstr(ChLV.isLV95, e=e, n=n), txt=_not_falsed_)

263 return None

264

265 @staticmethod

266 def unfalse2(e, n, LV95=None, **name):

267 '''Remove the I{Swiss LV95} or I{LV03} falsing.

268

269 @arg e: Falsed I{Swiss E_LV95} or I{y_LV03} easting (C{meter}).

270 @arg n: Falsed I{Swiss N_LV95} or I{x_LV03} northing (C{meter}).

271 @kwarg LV95: If C{True}, remove I{LV95} falsing, if C{False} remove

272 I{LV03} falsing, otherwise use method C{isLV95(B{e}, B{n})}.

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

274

275 @return: A L{ChLVYX2Tuple}C{(Y, X)} with the unfalsed B{C{e}}

276 respectively B{C{n}}.

277 '''

278 Y, X = _ChLV._falsing2(ChLV.isLV95(e, n) if LV95 is None else LV95)

279 return ChLVYX2Tuple(e - Y, n - X, **name)

280

281

282class ChLVa(_ChLV, LocalCartesian):

283 '''Conversion between (WGS84) geodetic and I{Swiss} projection coordinates using the U{SwissTopo

284 Approximate<https://www.SwissTopo.admin.ch/dam/en/sd-web/KLRCX9XIdXDu/ch1903wgs84-EN.pdf>}

285 formulas.

286

287 @see: U{Services<https://www.SwissTopo.admin.CH/en/maps-data-online/calculation-services.html>}

288 page 13 and U{older references<https://GitHub.com/alphasldiallo/Swisstopo-WGS84-LV03>}.

289 '''

290 def __init__(self, name=ChLV.Bern.name):

291 '''New L{ChLVa} I{Approximate SwissTopo} converter, centered at I{Bern, Ch}.

292

293 @kwarg name: Optional C{B{name}=Bern.name} (C{str}).

294 '''

295 LocalCartesian.__init__(self, latlonh0=ChLV.Bern, name=name)

296

297 def forward(self, latlonh, lon=None, height=0, M=None, **name):

298 # overloaded for the _ChLV.forward.__doc__

299 lat, lon, h, n = _llhn4(latlonh, lon, height, Error=ChLVError, **name)

300 a, b, h_ = _ChLV._llh2abh_3(lat, lon, h)

301 a2, b2 = a**2, b**2

302

303 Y = fdot_(211455.93, b,

304 -10938.51, b * a,

305 -0.36, b * a2,

306 -44.54, b * b2, start=72.37) # + 600_000

307 X = fdot_(308807.95, a,

308 3745.25, b2,

309 76.63, a2,

310 -194.56, b2 * a,

311 119.79, a2 * a, start=147.07) # + 200_000

312 return self._ChLV9Tuple(True, M, n, Y, X, h_, lat, lon, h)

313

314 def reverse(self, enh_, n=None, h_=0, M=None, **name): # PYCHOK signature

315 # overloaded for the _ChLV.reverse.__doc__

316 Y, X, h_, n = _ChLV._YXh_n4(enh_, n, h_, **name)

317 a, b, h = _ChLV._YXh_2abh3(Y, X, h_)

318 ab_d, a2, b2 = ChLV._ab_d, a**2, b**2

319

320 lat = Fdot_(3.238272, b,

321 -0.270978, a2,

322 -0.002528, b2,

323 -0.0447, a2 * b,

324 -0.014, b2 * b, start=16.9023892).fover(ab_d)

325 lon = Fdot_(4.728982, a,

326 0.791484, a * b,

327 0.1306, a * b2,

328 -0.0436, a * a2, start=2.6779094).fover(ab_d)

329 return self._ChLV9Tuple(False, M, n, Y, X, h_, lat, lon, h)

330

331

332class ChLVe(_ChLV, LocalCartesian):

333 '''Conversion between (WGS84) geodetic and I{Swiss} projection coordinates using the U{SwissTopo

334 Ellipsoidal approximate<https://www.swisstopo.admin.ch/dam/en/sd-web/c46Fz-MHIc3u/refsys-EN.pdf>}

335 formulas and U{Bolliger, J.<https://eMuseum.GGGS.CH/literatur-lv/liste-Dateien/1967_Bolliger_a.pdf>}

336 pp 148-151 (also U{GGGS<https://eMuseum.GGGS.CH/literatur-lv/liste.htm>}).

337

338 @note: Methods L{ChLVe.forward} and L{ChLVe.reverse} have an additional keyword argument

339 C{B{gamma}=False} to approximate the I{meridian convergence}. If C{B{gamma}=True}

340 a 2-tuple C{(t, gamma)} is returned with C{t} the usual result (C{ChLV9Tuple}) and

341 C{gamma}, the I{meridian convergence} (decimal C{degrees}). To convert C{gamma}

342 to C{grades} or C{gons}, use function L{pygeodesy.degrees2grades}.

343

344 @see: U{Services<https://www.SwissTopo.admin.CH/en/maps-data-online/calculation-services.html>}

345 page 10-11 and U{older references<https://GitHub.com/alphasldiallo/Swisstopo-WGS84-LV03>}.

346 '''

347 def __init__(self, name=ChLV.Bern.name):

348 '''New L{ChLVe} I{Ellipsoidal SwissTopo} converter, centered at I{Bern, Ch}.

349

350 @kwarg name: Optional C{B{name}=Bern.name} (C{str}).

351 '''

352 LocalCartesian.__init__(self, latlonh0=ChLV.Bern, name=name)

353

354 def forward(self, latlonh, lon=None, height=0, M=None, gamma=False, **name): # PYCHOK gamma

355 # overloaded for the _ChLV.forward.__doc__

356 lat, lon, h, n = _llhn4(latlonh, lon, height, Error=ChLVError, **name)

357 a, b, h_ = _ChLV._llh2abh_3(lat, lon, h)

358 ab_M, z, _H = ChLV._ab_M, 0, Fhorner

359

360 B1 = _H(a, 211428.533991, -10939.608605, -2.658213, -8.539078, -0.00345, -0.007992)

361 B3 = _H(a, -44.232717, 4.291740, -0.309883, 0.013924)

362 B5 = _H(a, 0.019784, -0.004277)

363 Y = _H(b, z, B1, z, B3, z, B5).fover(ab_M) # 1,000 Km!

364

365 B0 = _H(a, z, 308770.746371, 75.028131, 120.435227, 0.009488, 0.070332, -0.00001)

366 B2 = _H(a, 3745.408911, -193.792705, 4.340858, -0.376174, 0.004053)

367 B4 = _H(a, -0.734684, 0.144466, -0.011842)

368 X = _H(b, B0, z, B2, z, B4, z, 0.000488).fover(ab_M) # 1,000 Km!

369

370 t = self._ChLV9Tuple(True, M, n, Y, X, h_, lat, lon, h)

371 if gamma:

372 U1 = _H(a, 2255515.207166, 2642.456961, 1.284180, 2.577486, 0.001165)

373 U3 = _H(a, -412.991934, 64.106344, -2.679566, 0.123833)

374 U5 = _H(a, 0.204129, -0.037725)

375 g = _H(b, z, U1, z, U3, z, U5).fover(ChLV._ab_m) # * ChLV._ab_d degrees?

376 t = t, g

377 return t

378

379 def reverse(self, enh_, n=None, h_=0, M=None, gamma=False, **name): # PYCHOK gamma

380 # overloaded for the _ChLV.reverse.__doc__

381 Y, X, h_, n = _ChLV._YXh_n4(enh_, n, h_, **name)

382 a, b, h = _ChLV._YXh_2abh3(Y, X, h_)

383 s_d, _H, z = ChLV._s_d, Fhorner, 0

384

385 A0 = _H(b, ChLV._sLat, 32386.4877666, -25.486822, -132.457771, 0.48747, 0.81305, -0.0069)

386 A2 = _H(b, -2713.537919, -450.442705, -75.53194, -14.63049, -2.7604)

387 A4 = _H(b, 24.42786, 13.20703, 4.7476)

388 lat = _H(a, A0, z, A2, z, A4, z, -0.4249).fover(s_d)

389

390 A1 = _H(b, 47297.3056722, 7925.714783, 1328.129667, 255.02202, 48.17474, 9.0243)

391 A3 = _H(b, -442.709889, -255.02202, -96.34947, -30.0808)

392 A5 = _H(b, 9.63495, 9.0243)

393 lon = _H(a, ChLV._sLon, A1, z, A3, z, A5).fover(s_d)

394 # == (ChLV._sLon + a * (A1 + a**2 * (A3 + a**2 * A5))) / s_d

395

396 t = self._ChLV9Tuple(False, M, n, Y, X, h_, lat, lon, h)

397 if gamma:

398 U1 = _H(b, 106679.792202, 17876.57022, 4306.5241, 794.87772, 148.1545, 27.8725)

399 U3 = _H(b, -1435.508, -794.8777, -296.309, -92.908)

400 U5 = _H(b, 29.631, 27.873)

401 g = _H(a, z, U1, z, U3, z, U5).fover(ChLV._s_ab) # degrees

402 t = t, g

403 return t

404

405

406class ChLVr(ChLV, LocalCartesian):

407 '''Conversion between (Bessel1841) geodetic and I{Swiss} projection coordinates using the

408 U{SwissTopo Rigorous<https://www.SwissTopo.admin.CH/dam/en/sd-web/c46Fz-MHIc3u/refsys-EN.pdf>}

409 formulas 3.2 and 3.3, December 2018, pp 20-21.

410 '''

411 _a0 = _b0 = _K0 = \

412 _lam0 = _phi0 = _R0 = None

413

414 def __init__(self, name=ChLV.Bern.name):

415 '''New L{ChLVr} I{Rigorous SwissTopo} converter, centered at I{Bern, Ch}.

416

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

418 '''

419 LocalCartesian.__init__(self, latlonh0=ChLV.Bern, name=name)

420 assert self._E0

421

422 def _auxiliaries(self): # in .test/testChLVs.py

423 from pygeodesy import toDMS

424 d = type(self).__dict__

425 d = dict((n, v) for n, v in d.items() if n.endswith('0'))

426 D = degrees(self._b0)

427 d.update(_b0DMS=toDMS(D, prec=8))

428 return d

429

430 @property_ROver

431 def _E0(self): # set the C{ChLVr} auxiliaries, I{once}

432 from pygeodesy import Ellipsoids, sincos2

433

434 C = ChLVr

435 E = Ellipsoids.Bessel1841 # e=0.0816968312228869 e2=0.0066743722318609 b=6356078.962818

436

437 C._phi0 = phi0 = Phid(self.lat0) # 0.8194740686761218 == 46.9524055555+° == 46°57'08.66"

438 C._lam0 = Lamd(self.lon0) # 0.1298452241431616 == 7.4395833333+° == 7°26'22.50"

439 s0, c0 = sincos2(phi0)

440

441 C._a0 = sqrt(c0**4 * E.e2 / E.e21 + _1_0) # 1.000729138430515

442 C._b0 = asin(s0 / C._a0) # 0.8186943585816748 == 46.9077314579° == 46°54'27.83324843"

443 C._K0 = (_log_E(phi0, E) - _log_tan(phi0)) * C._a0 + _log_tan(C._b0) # 0.003066732377846

444 C._R0 = _over(E.b, _1_0 - s0**2 * E.e2) # 6378815.9036478335, E.a * sqrt(1 - E.e2) == E.b

445 return E # overwrite property_ROver

446

447 @property_RO

448 def ellipsoid(self):

449 '''Get the ellipsoid (C{Bessel1841}).

450 '''

451 return self._E0

452

453 def forward(self, latlonh, lon=None, height=0, M=None, **name):

454 '''Convert (Bessel1841) geodetic to I{Swiss} (Mercator) projection coordinates.

455

456 @arg latlonh: Either a C{LatLon}, L{Ltp} or C{scalar} (geodetic) latitude (C{degrees}).

457 @kwarg lon: Optional, C{scalar} (geodetic) longitude (C{degrees}) iff B{C{latlonh}} is

458 C{scalar}, ignored otherwise.

459 @kwarg height: Optional, height, vertically above (or below) the surface of the ellipsoid

460 (C{meter}) iff B{C{latlonh}} and B{C{lon}} are C{scalar}, ignored otherwise.

461 @kwarg M: Ignored, passed as-is (C{None}).

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

463

464 @return: A L{ChLV9Tuple}C{(Y, X, h_, lat, lon, height, ltp, ecef, M)} with the unfalsed

465 I{Swiss Y, X} coordinates, I{Swiss h_} height, the given I{geodetic} C{lat},

466 C{lon} and C{height}, C{ecef} (U{Ecef9Tuple<https://MrJean1.GitHub.io/PyGeodesy/

467 docs/pygeodesy.ecef-module.html>}) at I{Bern, Ch} and C{ltp} this C{ChLVr} instance.

468

469 @raise ChLVError: Invalid or non-C{scalar} B{C{latlonh}}, B{C{lon}} or B{C{height}}.

470 '''

471 lat, lon, h, n = _llhn4(latlonh, lon, height, Error=ChLVError, **name)

472 _, _, h_ = _ChLV._llh2abh_3(lat, lon, h)

473 phi, lam = Phid(lat), Lamd(lon)

474

475 S = (_log_tan(phi) - _log_E(phi, self._E0)) * self._a0 + self._K0

476 b = _atan_exp(S) # spherical (b, l)

477 B = (lam - self._lam0) * self._a0

478

479 sb, cb, sB, cB, sb0, cb0 = sincos2_(b, B, self._b0)

480

481 Y = _atan4(sB, cb0 * cB, sb0 * sb, cb) * self._R0

482 X = _log_R(sb * cb0 - cB * sb0 * cb, self._R0)

483

484 return self._ChLV9Tuple(True, M, n, Y, X, h_, lat, lon, h)

485

486 def reverse(self, enh_, n=None, h_=0, M=None, **name):

487 '''Convert I{Swiss} (Mercator) projection to (Bessel1841) geodetic coordinates.

488

489 @arg enh_: A Swiss projection (L{ChLV9Tuple}) or the C{scalar}, falsed I{Swiss E_LV95}

490 or I{y_LV03} easting (C{meter}).

491 @kwarg n: Falsed I{Swiss N_LV95} or I{x_LV03} northing (C{meter}) iff B{C{enh_}} is

492 C{scalar}, ignored otherwise.

493 @kwarg h_: I{Swiss h'} height (C{meter}) iff B{C{enh_}} and B{C{n}} are C{scalar},

494 ignored otherwise.

495 @kwarg M: Ignored, passed as-is (C{None}).

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

497

498 @return: A L{ChLV9Tuple}C{(Y, X, h_, lat, lon, height, ltp, ecef, M)} with the unfalsed

499 I{Swiss Y, X} coordinates, I{Swiss h_} height, the given I{geodetic} C{lat},

500 C{lon} and C{height}, C{ecef} (U{Ecef9Tuple<https://MrJean1.GitHub.io/PyGeodesy/

501 docs/pygeodesy.ecef-module.html>}) at I{Bern, Ch} and C{ltp} this C{ChLVr} instance.

502

503 @raise ChLVError: Invalid or non-C{scalar} B{C{enh_}}, B{C{n}} or B{C{h_}}.

504 '''

505 Y, X, h_, n = _ChLV._YXh_n4(enh_, n, h_, **name)

506 _, _, h = _ChLV._YXh_2abh3(Y, X, h_)

507

508 B_ = Y / self._R0 # spherical (b_, l_)

509 b_ = _atan_exp(X / self._R0)

510 sb_, cb_, sB_, cB_, sb0, cb0 = sincos2_(b_, B_, self._b0)

511

512 b = asin(sb_ * cb0 + cB_ * sb0 * cb_) # b, l

513 B = _atan4(sB_, cb0 * cB_, -sb0 * sb_, cb_)

514

515 e = self._E0.e

516 s = (_log_tan(b) - self._K0) / self._a0

517 for _ in range(_TRIPS): # 6..7 trips

518 r = asin(sin(b) * e)

519 S = _log_tan(r) * e + s

520 phi = _atan_exp(S)

521 if fabs(phi - b) < EPS:

522 break

523 b = phi

524

525 lat, lon = map1(degrees, phi, B / self._a0 + self._lam0)

526 return self._ChLV9Tuple(False, M, n, Y, X, h_, lat, lon, h)

527

528

529def _atan4(a, b, c, d): # avoid ZeroDivisionError

530 return atan2(a, b + c / d) if d else atan2(a * d, b * d + c)

531

532

533def _atan_exp(x):

534 return atan(exp(x)) * _2_0 - PI_2

535

536

537def _log_E(phi, E):

538 return _log_R(sin(phi) * E.e, E.e)

539

540

541def _log_R(sb_, R):

542 return log(_over(_1_0 + sb_, _1_0 - sb_)) * (R * _0_5)

543

544

545def _log_tan(phi):

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

547

548

549del _ChLV()._enh_n_h # PYCHOK set/check the _ChLV._enh_n_h names

550

551__all__ += _all_OTHER(ChLV, ChLVa, ChLVe, ChLVr, # passed along from PyGeodesy

552 EcefKarney, LatLon4Tuple, LocalCartesian, Ltp)

553

554# **) MIT License

555#

557#

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

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

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

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

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

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

564#

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

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

567#

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

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

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

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

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

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

574# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pychlv / chlvs.py: 97%

220 statements