Coverage for pygeodesy/hom.py: 99%

1 

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

3 

4u'''A pure Python implementation of 2 ellipsoidal and 3 spherical variants of I{Hotine}'s 

5U{Oblique Mercator<https://WikiPedia.org/wiki/Oblique_Mercator_projection>} projection from 

6U{Snyder Chapter 9<https://pubs.USGS.gov/pp/1395/report.pdf>} page 66ff, on-line U{Neacsu 

7<https://Neacsu.net/geodesy/snyder/3-cylindrical/sect_9>} and U{IOGP section 3.2.4 

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

9 

10@see: U{Oblique Mercator Usage<https://Neacsu.net/geodesy/snyder/3-cylindrical/sect_9/#usage>}. 

11''' 

12# make sure int/int division yields float quotient 

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

14 

15from pygeodesy.basics import _copysign, _ALL_LAZY, typename 

16from pygeodesy.constants import NAN, PI_2, R_M, _0_0, _0_5, _1_0, _2_0, _90_0, \ 

17 _EPS2e4 as _TOL, _isfinite, isnear0, _1_over, _over 

18from pygeodesy.datums import _ellipsoidal_datum, _WGS84 

19from pygeodesy.errors import _ValueError, _xkwds 

20from pygeodesy.fmath import fdot, fdot_ 

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

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

23 _x_, _y_, _0_, _1_, _2_ 

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

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

26from pygeodesy.namedTuples import LatLon2Tuple, Vector2Tuple 

27from pygeodesy.props import Property_RO, property_doc_, property_RO 

28from pygeodesy.streprs import Fmt, pairs, unstr 

29from pygeodesy.units import Azimuth, Degrees, Int, Meter, Lamd, Lat, _LatLon2, Lon, \ 

30 Phid, Radius, Scalar 

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

32 

33from math import cos, cosh, degrees, e as _base_e, fabs, log, radians, \ 

34 sin, sinh, sqrt, tan, tanh 

35 

36__all__ = _ALL_LAZY.hom 

37__version__ = '26.10.24' 

38 

39 

40class _HotineObliqueBase(_NamedBase): 

41 '''(INTERNAL) Base class for all C{HotineOblique*} variants. 

42 ''' 

43 _a0 = _b0 = _0_0 # origin (phi0, lam0) 

44 _ap = _bp = _0_0 # pole (phip, lamp) 

45 _args = () # __init__ args 

46 _g0 = _zC = NAN # RSO angle, azimuth 

47 _k0 = _1_0 # central scale factor 

48 _lat0 = _0_0 # origin latitude (degrees) NAN 

49 _x0 = _y0 = _0_0 # false east-, northing 

50 

51 def __init__(self, unused, k0=_1_0, lat0=None, 

52 x0=_0_0, y0=_0_0, **name): 

53 if name: 

54 self.name = name 

55 if k0 != 1: 

56 self._k0 = Scalar(k0=k0) 

57# self._Rk0 = self.R * self.k0 

58 if lat0 is not None: 

59 self._lat0 = Lat(lat0=lat0) 

60 self._a0 = Phid(self._lat0) 

61 sa0, ca0 = self._sca0 

62 if isnear0(sa0) or isnear0(ca0): 

63 raise self._Error(lat0=self._lat0) 

64 if x0: 

65 self._x0 = Meter(x0=x0) 

66 if y0: 

67 self._y0 = Meter(y0=y0) 

68 

69 @Property_RO 

70 def _args_kwds(self): 

71 kwds = self._kwds() 

72 if self.name: 

73 kwds = _xkwds(kwds, name=self.name) 

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

75 

76 @Property_RO 

77 def azimuthC(self): 

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

79 ''' 

80 z = self._zC 

81 if _isfinite(z): 

82 z = Azimuth(azimuth=degrees(z), _180=True) # Snyder page 275, 𝛼c 

83 return z 

84 

85 def _Error(self, **kwds): 

86 return HotineError(txt=str(self), **kwds) 

87 

88 @Property_RO 

89 def gamma0(self): 

90 '''Get the Rectified Skew Orthomorphic (RSO) angle (C{degrees}) or C{NAN}. 

91 ''' 

92 g = self._g0 

93 if _isfinite(g): 

94 g = Degrees(gamma0=degrees(g)) 

95 return g 

96 

97 @property_RO 

98 def k0(self): 

99 '''Get the scale factor at origin (C{scalar}). 

100 ''' 

101 return self._k0 

102 

103 def _kwds(self, **kwds): 

104 if self.k0 != 1: 

105 kwds.update(k0=self.k0) 

106 if self._lat0: 

107 kwds.update(lat0=self._lat0) 

108 if self._x0: 

109 kwds.update(x0=self._x0) 

110 if self._y0: 

111 kwds.update(y0=self._y0) 

112 return kwds 

113 

114 @property_RO 

115 def latlon0(self): 

116 '''Get the origin, other pole (L{LatLon2Tuple}). 

117 ''' 

118 return LatLon2Tuple(self._lat0, self._lon0, name=_0_) 

119 

120 @Property_RO 

121 def _lon0(self): 

122 return degrees(self._b0) 

123 

124 def _philam2(self, lat_latlon, lon): 

125 lat, lon = \ 

126 self._args = _LatLon2(lat_latlon, Error=HotineError) if lon is None else ( 

127 Lat(lat_latlon), Lon(lon)) 

128 return Phid(lat), Lamd(lon) 

129 

130 def _philam4(self, latlon1, latlon2): 

131 a1, b1 = self._philam2(latlon1, None) 

132 t = self._args 

133 a2, b2 = self._philam2(latlon2, None) 

134 self._args = t, self._args 

135 return a1, b1, a2, b2 

136 

137 def _philamalp3(self, latlon, lon, azimuth): 

138 a, b = self._philam2(latlon, lon) 

139 z = Azimuth(azimuth) # == % 360 

140 self._args += z, 

141 self._zC = z = radians(z) 

142 self._b0 = b # 𝜆C 

143 self._a0 = a # 𝜑C 

144 return a, b, z 

145 

146 @Property_RO 

147 def _sca0(self): # 𝜑0 

148 return sincos2(self._a0) 

149 

150# @Property_RO 

151# def _scb0(self): # 𝜆0 

152# return sincos2(self._b0) 

153 

154 @Property_RO 

155 def _scap(self): # 𝜑p 

156 return sincos2(self._ap) 

157 

158# @Property_RO 

159# def _scbp(self): # 𝜆p 

160# return sincos2(self._bp) 

161 

162 @Property_RO 

163 def _scg0(self): # 𝛾0 

164 return sincos2(self._g0) 

165 

166 @Property_RO 

167 def _sczC(self): # 𝛼C or 𝛽C 

168 return sincos2(self._zC) 

169 

170 def toStr(self, **kwds): 

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

172 ''' 

173 t = self._args_kwds 

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

175 

176 @property_RO 

177 def xy0(self): 

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

179 ''' 

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

181 

182 

183class _HotineOblique123(_HotineObliqueBase): 

184 '''(INTERNAL) Spherical C{HotineOblique1}, C{-2} and C{-3} base class. 

185 ''' 

186 _Rk0 = \ 

187 _radius = R_M 

188 

189 def __init__(self, ap, bp, radius, **k0_lat0_x0_y0_name): 

190 self._ap = ap 

191 self._bp = bp = _PI(bp) 

192 self._b0 = _PI(bp + PI_2) 

193 

194 _HotineObliqueBase.__init__(self, None, **k0_lat0_x0_y0_name) 

195 

196 if radius != _HotineOblique123._radius or \ 

197 self.k0 != _HotineObliqueBase._k0: 

198 self._radius = R = Radius(radius) 

199 self._Rk0 = R * self.k0 

200 

201 def forward(self, lat_latlon, lon=None, **name): 

202 '''Project a geodetic lat- and longitude to Hotine east- and northing. 

203 

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

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

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

207 is scalar, ignored otherwise. 

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

209 

210 @return: L{Hotine8Tuple}C{(x, y, u, v, lat, lon, scale, hom)} with falsed 

211 east- and northing C{x} and C{y} and "unrectified" C{u} and C{v} 

212 in C{meter}, geodetic C{lat}- and C{lon}gitude in degrees, C{scale} 

213 factor (aka C{k}) and this C{HotineOblique*} Mercator instance C{hom}. 

214 ''' 

215 lat, lon = _LatLon2(lat_latlon, lon, HotineError) 

216 

217 d = Lamd(lon) - self._b0 

218 sa, ca, sd, cd = sincos2_(Phid(lat), d) 

219 sp, cp = self._scap # phi-p 

220 

221 A = fdot_(sp, sa, -cp, ca * sd) 

222 u = fdot_(cp, sa, sp, ca * sd) # / ca 

223 u = atan2( u, ca * cd) * self._Rk0 

224 v = log(_d_s(_1_0, -A)) * self._Rk0 * _0_5 

225 

226 k = _1_0 - A**2 

227 k = (self.k0 / sqrt(k)) if k > 0 else NAN 

228 x = u + self._x0 

229 y = v + self._y0 

230 return Hotine8Tuple(x, y, u, v, lat, lon, k, self, 

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

232 

233 def _kwds(self, **unused): 

234 return _HotineObliqueBase._kwds(self, radius=self.radius) 

235 

236 @Property_RO 

237 def latlonP(self): 

238 '''Get the main pole (L{LatLon2Tuple}). 

239 ''' 

240 return LatLon2Tuple(degrees(self._ap), degrees(self._bp), name='P') 

241 

242 @property_RO 

243 def radius(self): 

244 '''Get the the earth radius (C{meter}, conventionally). 

245 ''' 

246 return self._radius 

247 

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

249 '''Project an Hotine east- and northing to geodetic lat- and longitude. 

250 

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

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

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

254 

255 @return: L{Hotine8Tuple}C{(x, y, u, v, lat, lon, scale, hom)} with geodetic 

256 C{lat}- and C{lon}gitude in degrees, C{scale} factor (aka C{k}), 

257 falsed east- and northing C{x} and C{y} and "unrectified" C(u) and 

258 C{v} in C{meter} and this C{HotineOblique*} Mercator instance C{hom}. 

259 ''' 

260 t = _1_over(self._Rk0) 

261 u = x - self._x0 

262 v = y - self._y0 

263 e = u * t 

264 n = v * t 

265 

266 sx, cx = sincos2(e) 

267 sy, cy, ty = sinh(n), cosh(n), tanh(n) 

268 sp, cp = self._scap 

269 

270 t = _over(sx, cy) 

271 a = asin1( fdot_(sp, ty, cp, t)) 

272 b = _PI(atan2(fdot_(sp, sx, -cp, sy), cx) + self._b0) 

273 return Hotine8Tuple(x, y, u, v, degrees(a), degrees(b), self.k0, self, 

274 name=_namethod(self.reverse, **name)) 

275 

276 

277class HotineOblique1(_HotineOblique123): 

278 '''Spherical variant 1: Hotine Oblique Mercator projection given 

279 C{two points}. 

280 ''' 

281 

282 def __init__(self, latlon1, latlon2, radius=R_M, **k0_lat0_x0_y0_name): 

283 '''New I{spherical} L{HotineOblique1} Mercator, variant 1. 

284 

285 @arg latlon1: Center line point (C{LatLon}, C{2-tuple (lat, lon)} or 

286 a comma-separated I{"latDMS, lonDMS"} C{str}). 

287 @arg latlon2: Other center line point (C{LatLon}, C{2-tuple (lat, lon)} 

288 or a comma separated I{"latDMS, lonDMS"} C{str}). 

289 @kwarg radius: Earth radius (C{meter}, conventionally). 

290 @kwarg k0_lat0_x0_y0_name: Central scale factor C{B{k0}=1.0} (C{scalar}), origin 

291 latitude C{B{lat0}=0.0} (C{degrees}), false easting C{B{x0}=0.0} and 

292 false northing C{B{y0}=0.0} (C{meter}) and C{B{name}=NN} (C{str}). 

293 ''' 

294 a1, b1, a2, b2 = self._philam4(latlon1, latlon2) 

295 

296 sa1, ca1, sb1, cb1 = sincos2_(a1, b1) 

297 sa2, ca2, sb2, cb2 = sincos2_(a2, b2) 

298 ca1sa2 = ca1 * sa2 

299 ca2sa1 = ca2 * sa1 

300 bp = atan2(fdot_(ca1sa2, cb1, -ca2sa1, cb2), 

301 fdot_(ca2sa1, sb2, -ca1sa2, sb1)) 

302 ap = atan1(-cos(bp - b1), tan(a1)) 

303 _HotineOblique123.__init__(self, ap, bp, radius, **k0_lat0_x0_y0_name) 

304 

305 @property_RO 

306 def latlon1(self): 

307 '''Get the 1st center line point (L{LatLon2Tuple}). 

308 ''' 

309 return LatLon2Tuple(*self._args[0], name=_1_) 

310 

311 @property_RO 

312 def latlon2(self): 

313 '''Get the 2nd center line point (L{LatLon2Tuple}). 

314 ''' 

315 return LatLon2Tuple(*self._args[1], name=_2_) 

316 

317 

318class HotineOblique2(_HotineOblique123): 

319 '''Spherical variant 2: Hotine Oblique Mercator projection given the C{central point} 

320 and C{azimuth}. 

321 ''' 

322 

323 def __init__(self, lat_latlon, lon=None, azimuth=_0_0, radius=R_M, **k0_lat0_x0_y0_name): 

324 '''New spherical L{HotineOblique2} Mercator, variant 2. 

325 

326 @arg lat_latlon: Center latitude (scalar C{degrees90} or any C{LatLon}, 2-tuple 

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

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

329 =ignored otherwise. 

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

331 @kwarg radius: Earth radius (C{meter}, conventionally). 

332 @kwarg k0_lat0_x0_y0_name: Central scale factor C{B{k0}=1.0} (C{scalar}), origin 

333 latitude C{B{lat0}=0.0} (C{degrees}), false easting C{B{x0}=0.0} and 

334 false northing C{B{y0}=0.0} (C{meter}) and C{B{name}=NN} (C{str}). 

335 ''' 

336 aC, bp, zC = self._philamalp3(lat_latlon, lon, azimuth) 

337 sa, ca, sz, cz = sincos2_(aC, zC) # == ._sca0, ._sczC 

338 

339 ap = asin1( sz * ca) 

340 bp += atan2(-cz, -sa * sz) # PYLINT E1130 

341 _HotineOblique123.__init__(self, ap, bp, radius, **k0_lat0_x0_y0_name) 

342 

343 @property_RO 

344 def latlonC(self): 

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

346 ''' 

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

348 

349 

350class HotineOblique3(_HotineOblique123): 

351 '''Spherical variant 3: Hotine Oblique Mercator projection given the C{main pole}. 

352 ''' 

353 

354 def __init__(self, lat_latlon, lon=None, radius=R_M, **k0_lat0_x0_y0_name): 

355 '''New spherical L{HotineOblique3} Mercator, variant 3. 

356 

357 @arg lat_latlon: Main pole latitude (scalar C{degrees90} or any C{LatLon}, 

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

359 @arg lon: Main pole longitude (C{degrees180}), required if B{C{lat_latlon}} is 

360 scalar, ignored otherwise. 

361 @kwarg radius: Earth radius (C{meter}, conventionally). 

362 @kwarg k0_lat0_x0_y0_name: Central scale factor C{B{k0}=1.0} (C{scalar}), origin 

363 latitude C{B{lat0}=0.0} (C{degrees}), false easting C{B{x0}=0.0} and 

364 false northing C{B{y0}=0.0} (C{meter}) and C{B{name}=NN} (C{str}). 

365 ''' 

366 ap, bp = self._philam2(lat_latlon, lon) 

367 _HotineOblique123.__init__(self, ap, bp, radius, **k0_lat0_x0_y0_name) 

368 

369 

370class _HotineObliqueAB(_HotineObliqueBase): 

371 '''(INTERNAL) Ellipsoidal C{HotineObliqueA} and C{-B} base class. 

372 ''' 

373 _ABEG_BuCz = () # 7-tuple 

374 _earth = _WGS84 

375 _useRSO = False 

376 

377 def __init__(self, earth, useRSO, **k0_lat0_x0_y0_name): 

378 

379 if earth is not _HotineObliqueAB._earth: 

380 self._earth = _ellipsoidal_datum(earth) 

381 if useRSO: 

382 self.useRSO = useRSO 

383 

384 _HotineObliqueBase.__init__(self, None, **k0_lat0_x0_y0_name) 

385 

386 @Property_RO 

387 def _ABD21EFt0(self): # 8-tuple 

388 sa0, ca0 = self._sca0 

389 E = self.ellipsoid 

390 e2 = E.e2not 

391 r = _1_0 - e2 

392 B = _over(ca0**4 * e2, r) 

393 B = sqrt(_1_0 + B) 

394 s = _1_0 - sa0**2 * e2 

395 r = sqrt(r) * B 

396 A = _over(r * E.a * self.k0, s) 

397 a0 = self._a0 

398 D = F = _over(r, sqrt(s) * ca0) 

399 D2 = D**2 

400 if D2 > 1: 

401 D21 = sqrt(D2 - _1_0) 

402 F += _copysign(D21, a0) # Snyder F, NOT! 

403 else: # PYCHOK no cover 

404 D21 = _0_0 

405 t0 = _tesB(E.e, sa0, 1, a0) # 9-13, not 9-13a! 

406 E = pow(t0, B) * F # Snyder E, 9-15 & 9-36 IOGP H 

407 return A, B, D, D2, D21, E, F, t0 

408 

409 @property_RO 

410 def ABmap(self): 

411 '''Get 14 map parameters C{(A, B, D, D2, E, G, uC, zC, 

412 a0, b0, g0, t0, x0, y0)} (L{HotineAB14Tuple}). 

413 ''' 

414 A, B, D, D2, _, E, _, t0 = self._ABD21EFt0 

415 _, _, _, G, _, uC, _ = self._ABEG_BuCz # PYLINT W0632 

416 return HotineAB14Tuple(A, B, D, D2, E, G, uC, self._zC, 

417 self._a0, self._b0, self._g0, t0, 

418 self._x0, self._y0, name=typename(self)) 

419 

420# def _debug(self, where, locals__, **kwds): 

421# d = dict(a0=self._a0, b0=self._b0, latlon0=self.latlon0, 

422# g0=self._g0, k0=self.k0, zC=self._zC) 

423# d.update(locals__) 

424# d.update(kwds) 

425# n = self.methodname(where) 

426# s = _COMMA_ + _NL_ + _SPACE_ * (len(n) + 1) 

427# print(n, s.join(pairs(d))) 

428 

429 @property_RO 

430 def earth(self): 

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

432 ''' 

433 return self._earth 

434 

435 @Property_RO 

436 def ellipsoid(self): 

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

438 ''' 

439 return self.earth.ellipsoid 

440 

441 def forward(self, lat_latlon, lon=None, **name): 

442 '''Project a geodetic lat- and longitude to Hotine east- and northing. 

443 

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

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

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

447 is scalar, ignored otherwise. 

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

449 

450 @return: L{Hotine8Tuple}C{(x, y, u, v, lat, lon, scale, hom)} with falsed east- 

451 and northing C{x} and C{y} and "unrectified" C{u} and C{v} in C{meter}, 

452 geodetic C{lat}- and C{lon}gitude in degrees, C{scale} factor (aka 

453 C{k}) and this C{HotineOblique*} Mercator instance C{hom}. 

454 ''' 

455 lat, lon = _LatLon2(lat_latlon, lon, HotineError) 

456 

457 A, B, Q, _, A_B, uC, zCE = self._ABEG_BuCz # PYLINT W0632 

458 

459 a = Phid(lat) 

460 b = Lamd(lon) 

461 dB = _PI(b - self._b0) * B 

462 E = self.ellipsoid 

463 sg0, cg0 = self._scg0 

464 

465 sa, ca, V, cdB = sincos2_(a, dB) # phi 

466 if ca: # PYLINT W0125 

467 Q *= _tesB(E.e, sa, -B, a) # PYLINT E1120 Snyder 9-25 ... 

468 S, T = _ST_Q(Q) 

469 U = _over(fdot_(S, sg0, -V, cg0), T) # PYLINT E1130 

470 

471 v = _d_s(_1_0, U) 

472 if _isfinite(v) and v > 0: # U +/-1 

473 v = log(v) * A_B * _0_5 # ... Snyder 9-30 

474# elif ca: # PYLINT W0125 

475# pass # 0.0 or INF 

476# else: # lat +/-90 

477# v = _copysign(self._g0 * _0_5, -sa) # PYLINT E1130 

478# v = log(tan(PI_4 + v)) * A_B 

479 

480 if cdB: # PYLINT W0125 

481 u = fdot_(S, cg0, V, sg0) 

482 u = atan2(u, cdB) * A_B # Snyder 9-31 XXX atan1 

483 else: # PYCHOK no cover 

484 u = A * dB # Note after 9-31 

485 

486 if uC or zCE: # variant B, IOGP page 62 

487 u -= (uC if b < self._b0 else # ABS(uC) * SIGN(a0) * SIGN(b0 - b) 

488 (-uC if b > self._b0 else u)) if zCE else uC 

489 else: # PYCHOK no cover 

490 u = a * A_B # lat +/-90 idem 

491 v = _0_0 

492 

493 x = _1_0 - sa**2 * E.e2not 

494 if x > 0: 

495 k = E.a * ca * cdB # Snyder 9-32 

496 k = Scalar(k=cos(u / A_B) * sqrt(x) * A / k) if k and A_B else NAN 

497 else: # PYCHOK no cover 

498 k = NAN if x < 0 else _0_0 

499 x, y = self._rectify2(u, v, self._x0, self._y0) 

500 

501# self._debug(self.forward, locals(), t=_tesB(E.e, sa, 1, a)) 

502 return Hotine8Tuple(x, y, u, v, lat, lon, k, self, 

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

504 

505 def _kwds(self, **unused): 

506 return _HotineObliqueBase._kwds(self, earth=self.ellipsoid, useRSO=self.useRSO) 

507 

508 def _rectify2(self, e, n, uC, v0): 

509 s, c = self._scg0 if self.useRSO else self._sczC # PYLINT W0125 

510 return (fdot_(n, c, e, s, start=uC), # 9-40 

511 fdot_(e, c, -n, s, start=v0)) # 9-41 

512 

513 def reverse(self, x, y, triPhi=0, **name): 

514 '''Project an Hotine east- and northing to geodetic lat- and longitude. 

515 

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

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

518 @kwarg triPhi: Number of C{phi} iterations (C{int}), set C{B{triPhi}=0} 

519 to use the Horner-like C{e2-sin(X2)} series. 

520 

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

522 

523 @return: L{Hotine8Tuple}C{(x, y, u, v, lat, lon, scale, hom)} with geodetic 

524 C{lat}- and C{lon}gitude in degrees, C{scale} factor (aka C{k}), 

525 falsed east- and northing C{x} and C{y} and "unrectified" C{u} and 

526 C{v} in C{meter} and this C{HotineOblique*} Mercator instance C{hom}. 

527 ''' 

528 _, B, E, _, A_B, uC, _ = self._ABEG_BuCz # PYLINT W0632 

529 

530 u, v = self._rectify2(x - self._x0, y - self._y0, uC, _0_0) 

531 V, cuB_A = sincos2(_over(u, A_B)) 

532 sg0, cg0 = self._scg0 

533 

534 Q = pow(_base_e, _over(-v, A_B)) 

535 S, T = _ST_Q(Q) 

536 U = _over(fdot_(V, cg0, S, sg0), T) 

537 t, i = _d_s(_1_0, -U), None 

538 if _isfinite(t) and t > 0: 

539 t = pow(E / sqrt(t), _1_over(B)) 

540 X = _phit(t) # Snyder chi, 7-13 

541 n = Int(triPhi=triPhi) 

542 e = self.ellipsoid.e 

543 if e and n > 0: 

544 e_2 = e * _0_5 

545 for i in range(n): # 4..5 

546 a = _d_s(_1_0, sin(X) * e) 

547 a = _phit(pow(a, e_2) * t) # 7-9 

548 d = fabs(a - X) 

549 if d < _TOL: 

550 i += 1 

551 break 

552 X = a 

553 else: # PYCHOK no cover 

554 t = unstr(self.reverse, x, y, triPhi=n) 

555 raise HotineError(Fmt.no_convergence(d, _TOL), txt=t) 

556 

557 else: # Horner-like e2-sin(X2) series, 3-5 

558 sX4 = [sin(X * _) for _ in range(2, 9, 2)] 

559 a = fdot(self._Xe5, X, *sX4) 

560 

561 b = fdot_(S, cg0, -V, sg0) 

562 b = atan2(b, cuB_A) 

563 b = _PI(self._b0 - _over(b, B)) 

564 

565 lat, lon = degrees(a), degrees(b) 

566 else: # PYCHOK no cover a = +/-PI_2, b = self._b0 

567 lat, lon = _copysign(_90_0, U), self._lon0 

568 

569# self._debug(self.reverse, locals()) 

570 return Hotine8Tuple(x, y, u, v, lat, lon, self.k0, self, iteration=i, 

571 name=_namethod(self.reverse, **name)) 

572 

573 @property_doc_(' apply Rectified Skew Orthomorphic (RSO) angle C{gamma0}, otherwise Oblique Mercator C{azimuthC} (C{bool}).') 

574 def useRSO(self): 

575 return self._useRSO 

576 

577 @useRSO.setter # PYCHOK setter! 

578 def useRSO(self, RSO): 

579 self._useRSO = bool(RSO) 

580 

581 @Property_RO 

582 def _Xe5(self): # Snyder page 75, 3-5 

583 e2 = self.ellipsoid.e2not 

584 e4 = e2 * e2 

585 e6 = e4 * e2 

586 return (_1_0, e2 * (_0_5 + e2 * (5 / 24 + e2 * ( 1 / 12 + e2 * ( 13 / 360)))), 

587 e4 * ( 7 / 48 + e2 * (29 / 240 + e2 * ( 811 / 11520))), 

588 e6 * ( 7 / 120 + e2 * ( 81 / 1120)), 

589 e6 * e2 * (4279 / 161280)) 

590 

591 

592class HotineObliqueA(_HotineObliqueAB): 

593 '''Ellipsoidal variant A: Hotine Oblique Mercator projection given C{two points}. 

594 ''' 

595 

596 def __init__(self, latlon1, latlon2, earth=_WGS84, useRSO=False, **k0_lat0_x0_y0_name): 

597 '''New ellipsoidal L{HotineObliqueA} Mercator, variant A. 

598 

599 @arg latlon1: Center line point (C{LatLon}, 2-tuple C{(lat, lon)} or 

600 comma-separated I{"latDMS, lonDMS"} C{str}). 

601 @arg latlon2: Other center line point (C{LatLon}, 2-tuple C{(lat, lon)} or 

602 comma-separated I{"latDMS, lonDMS"} C{str}). 

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

604 @kwarg useRSO: If C{True} apply Rectified Skew Orthomorphic (RSO) angle C{gamma0}, 

605 otherwise Oblique Mercator C{azimuthC} (C{bool}). 

606 @kwarg k0_lat0_x0_y0_name: Central scale factor C{B{k0}=1.0} (C{scalar}), origin 

607 latitude C{B{lat0}=0.0} (C{degrees}), false easting C{B{x0}=0.0} and 

608 false northing C{B{y0}=0.0} (C{meter}) and C{B{name}=NN} (C{str}). 

609 ''' 

610 a1, b1, a2, b2 = self._philam4(latlon1, latlon2) 

611 _HotineObliqueAB.__init__(self, earth, useRSO, **k0_lat0_x0_y0_name) 

612 

613 (lat1, _), (lat2, _) = self._args 

614 if isnear0(a1 - a2): # phi1 == phi2 

615 raise self._Error(lat1=lat1, lat2=lat2) 

616 sa1, ca1, sa2, ca2 = sincos2_(a1, a2) 

617 if isnear0(sa1) or isnear0(ca1): # phi1 +/-PI_2 or 0 

618 raise self._Error(lat1=lat1) 

619 elif isnear0(ca2 + PI_2): # phi2 == -PI_2 

620 raise self._Error(lat2=lat2) 

621 

622 A, B, D, _, _, E, _, _ = self._ABD21EFt0 

623 e = self.ellipsoid.e 

624 H = _tesB(e, sa1, B, a1) # Snyder H 

625 L = _tesB(e, sa2, B, a2) 

626 G = (_over(E, H) - _over(H, E)) * _0_5 

627 J = _d_s(E**2, L * H) 

628 P = _d_s(L, H) 

629 

630 b0 = atan1(tan(B * _PI(b1 - b2) * _0_5) * J, P) # atan2? 

631 b0 = (b1 + b2) * _0_5 - _over(b0, B) 

632 g0 = atan1(sin(B * _PI(b1 - b0)), G) # atan2? 

633 zC = asin1(sin(g0) * D) 

634 uC = _0_0 

635 

636 self._b0, self._g0, self._zC = b0, g0, zC 

637 self._ABEG_BuCz = A, B, E, G, _over(A, B), uC, False 

638# self._debug(HotineObliqueA.__init__, locals(), gamma0=self.gamma0, zCd=self.azimuth) 

639 

640 @property_RO 

641 def latlon1(self): 

642 '''Get the 1st center line point (L{LatLon2Tuple}). 

643 ''' 

644 return LatLon2Tuple(*self._args[0], name=_1_) 

645 

646 @property_RO 

647 def latlon2(self): 

648 '''Get the 2nd center line point (L{LatLon2Tuple}). 

649 ''' 

650 return LatLon2Tuple(*self._args[1], name=_2_) 

651 

652 

653class HotineObliqueB(_HotineObliqueAB): 

654 '''Ellipsoidal variant B: Hotine Oblique Mercator projection given the C{central point} 

655 and C{azimuth}. 

656 ''' 

657 

658 def __init__(self, lat_latlon, lon=None, azimuth=_0_0, earth=_WGS84, useRSO=False, **k0_lat0_x0_y0_name): 

659 '''New ellipsoidal L{HotineObliqueB} Mercator, variant B. 

660 

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

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

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

664 ignored otherwise. 

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

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

667 @kwarg useRSO: If C{True} apply Rectified Skew Orthomorphic (RSO) angle C{gamma0}, 

668 otherwise Oblique Mercator C{azimuthC} (C{bool}). 

669 @kwarg k0_lat0_x0_y0_name: Central scale factor C{B{k0}=1.0} (C{scalar}), origin 

670 latitude C{B{lat0}=0.0} (C{degrees}), false easting C{B{x0}=0.0} and 

671 false northing C{B{y0}=0.0} (C{meter}) and C{B{name}=NN} (C{str}). 

672 ''' 

673 aC, bC, zC = self._philamalp3(lat_latlon, lon, azimuth) 

674 _HotineObliqueAB.__init__(self, earth, useRSO, **k0_lat0_x0_y0_name) 

675 szC, czC = self._sczC 

676 

677 A, B, D, _, D21, E, F, _ = self._ABD21EFt0 

678 A_B = _over(A, B) 

679 G = (F - _1_over(F)) * _0_5 

680 

681 g0 = asin1(_over(szC, D)) # == IOGP 

682 b0 = bC - _over(asin1(tan(g0) * G), B) # == IOGP 

683 

684 zCE = isnear0(zC - PI_2) # IOGP page 61, zC == 90? 

685 uC = (A * (bC - b0)) if zCE else ( 

686 atan2(D21, czC) * A_B) 

687 uC = _copysign(uC, aC) # == ABS(uC) * SIGN(𝜑C) 

688 

689 self._b0, self._g0 = _PI(b0), g0 

690 

691 self._ABEG_BuCz = A, B, E, G, A_B, uC, zCE 

692# self._debug(HotineObliqueB.__init__, locals(), gamma0=self.gamma0) 

693 

694 @property_RO 

695 def latlonC(self): 

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

697 ''' 

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

699 

700 

701class HotineError(_ValueError): 

702 '''Error raised for C{HotineOblique*} issues. 

703 ''' 

704 pass # PYLINT W0107 

705 

706 

707class Hotine8Tuple(_NamedTuple): 

708 '''8-Tuple C{(x, y, u, v, lat, lon, scale, hom)} with (falsed) east- and northing C{x} 

709 and C{y} and "unrectified" C{u} and C{v} in C{meter}, geodetic C{lat}- and C{lon}gitude 

710 in degrees, C{scale} factor (aka C{k}) and C{HotineOblique*} Mercator instance C{hom}. 

711 

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

713 ''' 

714 _Names_ = (_x_, _y_, 'u', 'v', _lat_, _lon_, _scale_, 'hom') 

715 _Units_ = ( Meter, Meter, Meter, Meter, Lat, Lon, _Pass, _Pass) 

716 

717 

718class HotineAB14Tuple(_NamedTuple): 

719 '''14-Tuple C{(A, B, D, D2, E, G, uC, zC, a0, b0, g0, t0, x0, y0)} of ellipsoidal map 

720 parameters I{only} for variants L{HotineObliqueA} and L{HotineObliqueB} with 

721 angles C{zC}, C{a0}, C{b0} and C{g0} in C{radians}. 

722 ''' 

723 _Names_ = tuple('A,B,D,D2,E,G,uC,zC,a0,b0,g0,t0,x0,y0'.split(_COMMA_)) 

724 _Units_ = (Scalar,) * 7 + (_Pass,) * 4 + (Scalar,) + (Meter,) * 2 

725 

726 

727def _d_s(a, b): # diff / sum 

728 return _over(a - b, a + b) 

729 

730 

731def _phit(t): # Snyder page 75, X 7-13 

732 return PI_2 - atan1(t) * _2_0 # 7-13 

733 

734 

735def _ST_Q(Q): 

736 _Q = _1_over(Q) 

737 S = (Q - _Q) * _0_5 

738 T = (Q + _Q) * _0_5 

739 return S, T 

740 

741 

742def _tesB(e, s, B, *a): # t(e, s) power B! 

743 if a and B: # ϕ, Snyder page 71, 9-13 

744 t = _d_s(_1_0, e * s) 

745 t = pow( t, -e * _0_5) 

746 t *= tan(_0_5 * (PI_2 - a[0]))