Coverage for pygeodesy/lom.py: 99%

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

4u'''A pure Python implementation of I{Laborde}'s U{Oblique Mercator

5<https://WikiPedia.org/wiki/Oblique_Mercator_projection>} projection U{IOGP section

63.2.4.2<https://www.IOGP.org/wp-content/uploads/2019/09/373-07-02.pdf>} page 64ff.

7'''

8# make sure int/int division yields float quotient

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

11from pygeodesy.basics import _copysign, _ALL_LAZY

12from pygeodesy.constants import PI_2, _0_0, _0_5, _1_0, _2_0, _3_0, _over, \

13 _EPS2e4 as _TOL

14from pygeodesy.datums import _ellipsoidal_datum, _WGS84

15from pygeodesy.errors import _ValueError

16from pygeodesy.fmath import fdot_, hypot

17# from pygeodesy.internals import typename # from .basics

18from pygeodesy.interns import _C_, _COMMASPACE_, _lat_, _lon_, _scale_, \

19 _V_, _x_, _y_, _0_

20# from pygeodesy.lazily import _ALL_LAZY # from .basics

21from pygeodesy.named import _NamedBase, _NamedTuple, _namethod, _Pass

22from pygeodesy.namedTuples import LatLon2Tuple, PhiLam2Tuple, Vector2Tuple

23from pygeodesy.props import Property_RO, property_RO, property_ROver

24from pygeodesy.streprs import Fmt, pairs, unstr

25from pygeodesy.units import Azimuth, Meter, Lamd, Lat, _LatLon2, Lon, \

26 Phid, Scalar

27from pygeodesy.utily import asin1, atan1, sincos2, sincos2_, wrapPI as _PI

29from math import atan, degrees, e as _base_e, fabs, log, radians, \

30 sin, sqrt, tan

32__all__ = _ALL_LAZY.lom

33__version__ = '26.10.20'

35_TRIPS = 17 # per iteration loop

38class LabordeOblique(_NamedBase):

39 '''Ellipsoidal Laborde Oblique Mercator projection given the C{central point}

40 and C{azimuth}.

41 '''

42 _aC = _0_0 # phiC, radians

43 _args = () # 3-tuple

44 _aS = _0_0 # phiS, radians

45 _bC = _0_0 # lamC, radians

46 _BCGR = () # 4-tuple

47 _d = _0_0 # iteration delta

48 _earth = _WGS84

49 _kC = _1_0 # central scale

50 _x0 = _0_0 # false easting

51 _y0 = _0_0 # false northing

52 _zC = _0_0 # central azimuth, radians

54 def __init__(self, lat_latlon, lon=None, azimuth=_0_0, earth=_WGS84, kC=_1_0, x0=_0_0, y0=_0_0, **name):

55 '''New ellipsoidal L{LabordeOblique} Mercator.

57 @arg lat_latlon: Central latitude (scalar C{degrees90} or C{LatLon}, 2-tuple

58 C{(lat, lon)} or comma-separated I{"latDMS, lonDMS"} C{str}).

59 @arg lon: Central longitude (C{degrees180}), required if B{C{lat_latlon}} is scalar,

60 ignored otherwise.

61 @kwarg azimuth: Azimuth of the center line (C{degrees360} from North).

62 @kwarg earth: Earth model (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).

63 @kwarg kC: Central scale factor (C{scalar}).

64 @kwarg x0: False easting (C{meter}).

65 @kwarg y0: False northing (C{meter}).

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

68 @see: U{Roundtrip error<ttps://docs.Carmenta.com/pages/labordeobliquemercatorprojection.html>}

69 '''

70 self._args = lat, lon = _LatLon2(lat_latlon, lon, LabordeError)

71 self._aC, self._bC = Phid(lat), Lamd(lon)

73 zC = Azimuth(azimuth)

74 self._args += zC,

75 if zC:

76 self._zC = zC = radians(zC)

77 szC, czC, s, c = sincos2_(zC, zC * _2_0)

78 G = complex(_1_0 - c, s) / 12.0

79 del s, c

80 else: # PYCHOK no cover

81 czC = _1_0

82 szC = G = zC = _0_0

83 if earth is not LabordeOblique._earth:

84 self._earth = _ellipsoidal_datum(earth)

86 if name:

87 self.name = name

88 if kC != 1:

89 self._kC = Scalar(kC=kC)

90 if x0:

91 self._x0 = Meter(x0=x0)

92 if y0:

93 self._y0 = Meter(y0=y0)

95 E = self.ellipsoid

96 r = _1_0 - E.e2not

98 B = _over(czC**4 * E.e2not, r)

99 B = sqrt(B + _1_0)

100

101 aS = asin1(_over(szC, B))

102 C = log(_tan_2(aS)) - _log_2(zC, E.e) * B

103

104 R = _1_0 - szC**2 * E.e2not

105 R = _over(sqrt(r) * E.a * kC, R)

106

107 self._aS = aS

108 self._BCGR = B, C, G, R

109

110 @Property_RO

111 def _args_kwds(self):

112 kwds = dict(earth=self.ellipsoid)

113 if self.kC != 1:

114 kwds.update(kC=self.kC)

115 if self._x0:

116 kwds.update(x0=self._x0)

117 if self._y0:

118 kwds.update(y0=self._y0)

119 if self.name:

120 kwds.update(name=self.name)

121 return unstr(None, *self._args, **kwds)

122

123 @Property_RO

124 def azimuthC(self):

125 '''Get the center line azimuth (C{degrees180}) or C{NAN}.

126 '''

127 return Azimuth(azimuth=degrees(self._zC), _180=True)

128

129 @property_RO

130 def earth(self):

131 '''Get the earth model (L{Datum}).

132 '''

133 return self._earth

134

135 @Property_RO

136 def ellipsoid(self):

137 '''Get the earth' ellipsoid (L{Ellipsoid}).

138 '''

139 return self.earth.ellipsoid

140

141 def forward(self, lat_latlon, lon, **name):

142 '''Project a geodetic lat- and longitude in C{degrees} to Laborde east- and northing.

143

144 @arg lat_latlon: Latitude of the point (scalar C{degrees90} or C{LatLon},

145 2-tuple C{(lat, lon)} or comma-separated I{"latDMS, lonDMS"} C{str}).

146 @arg lon: Longitude of the point (C{degrees180}), required if B{C{lat_latlon}}

147 is scalar, ignored otherwise.

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

149

150 @return: L{Laborde8Tuple}C{(x, y, U, V, lat, lon, scale, lom)} with falsed

151 east- and northing C{x} and C{y} in C{meter}, "unrectified" C{U} and

152 C{V}, geodetic C{lat}- and C{lon}gitude in degrees, C{scale} factor

153 (aka C{k}) and this C{LabordeOblique*} Mercator instance C{lom}.

154 '''

155 lat, lon = _LatLon2(lat_latlon, lon, LabordeError)

156 a, b = Phid(lat), Lamd(lon)

157

158 B, C, G, R = self._BCGR

159 e = self.ellipsoid.e

160

161 L = (b - self._bC) * B

162 q = _log_2(a, e) * B + C

163 _, P = _phiq2(q)

164 L, P, U, V = _LPUV4(L, P, True, *self._scaS)

165

166 H = complex(-L, P)

167 H += pow(H, 3) * G

168 x = H.imag * R + self._x0

169 y = H.real * R + self._y0

170 return Laborde8Tuple(x, y, U, V, lat, lon, self.kC, self,

171 name=_namethod(self.forward, **name))

172

173 @property_RO

174 def kC(self):

175 '''Get the scale factor at center (C{scalar}).

176 '''

177 return self._kC

178

179 @property_RO

180 def latlonC(self):

181 '''Get the center point (L{LatLon2Tuple}).

182 '''

183 return LatLon2Tuple(*self._args[:2], name=_C_)

184

185 @property_ROver

186 def latlonParis(self):

187 '''Get Paris Observatory lat- and Paris Meridian longitude (L{LatLon2Tuple}).

188 '''

189 # <https://Wikipedia.org/wiki/Paris_Observatory> <https://WikiPedia.org/wiki/Paris_meridian>

190 return LatLon2Tuple(Lat('48 48 18.36'), Lon('2 20 14.025'), name='Paris')

191

192 @property_RO

193 def philamC(self):

194 '''Get the center point (L{PhiLam2Tuple}).

195 '''

196 return _philam2(self.latlonC)

197

198 @property_ROver

199 def philamParis(self):

200 '''Get Paris Observatory lat- and Paris Meridian longitude (L{PhiLam2Tuple}).

201 '''

202 return _philam2(self.latlonParis)

203

204 def reverse(self, x, y, **name):

205 '''Project an Laborde east- and northing to geodetic lat- and longitude.

206

207 @arg x: Easting of the point (C{meter}).

208 @arg y: Northing of the point (C{meter}).

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

210

211 @return: L{Laborde8Tuple}C{(x, y, U, V, lat, lon, scale, lom)} with geodetic

212 C{lat}- and C{lon}gitude in degrees, "unrectified" C{U} and C{V},

213 C{scale} factor (aka C{k}), falsed east- and northing C{x} and C{y}

214 and this C{LabordeOblique*} Mercator instance C{lom}.

215 '''

216 B, C, G, R = self._BCGR

217

218 H0 = complex(y - self._y0,

219 x - self._x0) / R

220 P = pow( H0, 3) * G

221 H_ = _over(H0, P + H0)

222 for i in self._TRIPS(0, x, y):

223 L = pow(H_, 2) * G

224 P = pow(H_, 3) * G

225 H = _over(P * _2_0 + H0,

226 L * _3_0 + _1_0)

227 P += H

228 P -= H0

229 if self._TOL(P.real):

230 break

231 H_ = H

232 P = pow(_base_e, H.imag)

233 P = atan(P) * _2_0 - PI_2

234 L, P, U, V = _LPUV4(-H.real, P, False, *self._scaS)

235

236 b = _PI(_over(L, B) + self._bC)

237 e = self.ellipsoid.e

238 q, a_ = _phiq2(_over(P - C, B))

239 for i in self._TRIPS(i, x, y): # PYLINT W0631

240 a = _1_se(-sin(a_), e) * q

241 a = _phiq(a)

242 if self._TOL(a - a_):

243 break

244 a_ = a

245

246 return Laborde8Tuple(x, y, U, V, degrees(a), degrees(b), self.kC, self,

247 iteration=i, name=_namethod(self.reverse, **name))

248

249 @Property_RO

250 def _scaS(self):

251 return sincos2(self._aS)

252

253# @Property_RO

254# def _sczC(self):

255# return sincos2(self._zC)

256

257 def toStr(self, **kwds):

258 '''Get this C{HotineOblique*} instance as C{str}.

259 '''

260 t = self._args_kwds

261 return _COMMASPACE_(t, *pairs(kwds))

262

263 def _TOL(self, d):

264 self._d = d = fabs(d)

265 return d < _TOL

266

267 def _TRIPS(self, n, *xy):

268 for i in range(n + 1, n + _TRIPS):

269 yield i

270 else: # PYCHOK no cover

271 t = unstr(self.reverse, *xy)

272 e = Fmt.no_convergence(self._d, _TOL)

273 raise LabordeError(e, txt=t)

274

275 @property_RO

276 def xy0(self):

277 '''Get the false east- and northing (L{Vector2Tuple}, C{meter}).

278 '''

279 return Vector2Tuple(self._x0, self._y0, name=_0_)

280

281

282class LabordeError(_ValueError):

283 '''Error raised for C{LabordeOblique*} issues.

284 '''

285 pass # PYLINT W0107

286

287

288class Laborde8Tuple(_NamedTuple):

289 '''8-Tuple C{(x, y, U, V, lat, lon, scale, lom)} with (falsed) east- and northing

290 C{x} and C{y} in C{meter}, "unrectified" C{U} and C{V}, geodetic C{lat}- and

291 C{lon}gitude in degrees, C{scale} factor (aka C{k}) and C{LabordeOblique*}

292 Mercator instance C{lom}.

293

294 @note: Except C{hom}, any attribute may be C{NAN} or non-finite for singulaties.

295 '''

296 _Names_ = (_x_, _y_, 'U', _V_, _lat_, _lon_, _scale_, 'lom')

297 _Units_ = ( Meter, Meter, Scalar, Scalar, Lat, Lon, Scalar, _Pass)

298

299

300def _log_2(a, e):

301 r = _1_se(sin(a), e)

302 r *= _tan_2(a)

303 return log(r)

304

305

306def _LPUV4(L, P, forward, saS, caS):

307 sL, cL, sP, cP = sincos2_(L, P)

308 V = cP * cL

309 W = sP if forward else (cP * sL)

310 U = fdot_(V, caS, W, saS)

311 V = fdot_(V, saS, -W, caS)

312 if forward:

313 W = cP * sL

314 else:

315 W = V

316 V = sP

317 d = hypot(U, V)

318 if d:

319 L = atan1(V, U + d) * _2_0

320 P = atan1(W, d)

321 else: # PYCHOK no cover

322 L = _0_0

323 P = _copysign(PI_2, W)

324 return L, log(_tan_2(P)), U, V

325

326

327def _philam2(t):

328 return PhiLam2Tuple(radians(t.lat), radians(t.lon), name=t.name)

329

330

331def _phiq(q):

332 return atan(q) * _2_0 - PI_2

333

334

335def _phiq2(q):

336 q = pow(_base_e, q)

337 return q, _phiq(q)

338

339

340def _1_se(s, e): # diff / sum

341 s *= e

342 r = _over(_1_0 - s, _1_0 + s)

343 return pow(r, e * _0_5)

344

345

346def _tan_2(a):

347 # tan(x/2) == sin(x) / (cos(x) + 1)

348 # == (1 - cos(x)) / sin(x)

349 return tan((PI_2 + a) * _0_5)

350

351# **) MIT License

352#

354#

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

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

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

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

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

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

361#

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

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

364#

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

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

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

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

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

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

371# OTHER DEALINGS IN THE SOFTWARE.