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 INF, NAN, PI_2, PI_4, R_M, _0_0, _0_5, _1_0, _2_0, \ 

17 _90_0, _isfinite, isnear0, _1_over, _over, \ 

18 _EPS2e4 as _TOL 

19from pygeodesy.datums import _ellipsoidal_datum, _WGS84 

20from pygeodesy.errors import _ValueError, _xkwds 

21from pygeodesy.fmath import fdot, fdot_ 

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

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

24 _x_, _y_, _0_, _1_, _2_ 

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

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

27from pygeodesy.namedTuples import LatLon2Tuple, Vector2Tuple 

28from pygeodesy.props import Property_RO, property_doc_, property_RO 

29from pygeodesy.streprs import Fmt, pairs, unstr 

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

31 Lon, Phid, Radius, Scalar 

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

33 

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

35 sin, sinh, sqrt, tan, tanh 

36 

37__all__ = _ALL_LAZY.hom 

38__version__ = '26.10.21' 

39 

40 

41class _HotineObliqueBase(_NamedBase): 

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

43 ''' 

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

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

46 _args = () # __init__ args 

47 _g0 = _zC = NAN # rectified-skew, azimuth 

48 _k0 = _1_0 # central scale factor 

49 _lat0 = _0_0 # origin latitude (degrees) NAN 

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

51 

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

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

54 if name: 

55 self.name = name 

56 if k0 != 1: 

57 self._k0 = Scalar(k0=k0) 

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

59 if lat0 is not None: 

60 self._lat0 = Lat(lat0=lat0) 

61 self._a0 = Phid(self._lat0) 

62 sa0, ca0 = self._sca0 

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

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

65 if x0: 

66 self._x0 = Meter(x0=x0) 

67 if y0: 

68 self._y0 = Meter(y0=y0) 

69 

70 @Property_RO 

71 def _args_kwds(self): 

72 kwds = self._kwds() 

73 if self.name: 

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

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

76 

77 @Property_RO 

78 def azimuthC(self): 

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

80 ''' 

81 z = self._zC 

82 if _isfinite(z): 

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

84 return z 

85 

86 def _Error(self, **kwds): 

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

88 

89 @Property_RO 

90 def gamma0(self): 

91 '''Get the rectified-skew angle (C{degrees}) or C{NAN}. 

92 ''' 

93 g = self._g0 

94 if _isfinite(g): 

95 g = Degrees(gamma0=degrees(g)) 

96 return g 

97 

98 @property_RO 

99 def k0(self): 

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

101 ''' 

102 return self._k0 

103 

104 def _kwds(self, **kwds): 

105 if self.k0 != 1: 

106 kwds.update(k0=self.k0) 

107 if self._lat0: 

108 kwds.update(lat0=self._lat0) 

109 if self._x0: 

110 kwds.update(x0=self._x0) 

111 if self._y0: 

112 kwds.update(y0=self._y0) 

113 return kwds 

114 

115 @property_RO 

116 def latlon0(self): 

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

118 ''' 

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

120 

121 @Property_RO 

122 def _lon0(self): 

123 return degrees(self._b0) 

124 

125 def _philam2(self, lat_latlon, lon): 

126 lat, lon = \ 

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

128 ( Lat(lat_latlon), Lon(lon)) 

129 return Phid(lat), Lamd(lon) 

130 

131 def _philam4(self, latlon1, latlon2): 

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

133 t = self._args 

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

135 self._args = t, self._args 

136 return a1, b1, a2, b2 

137 

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

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

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

141 self._args += z, 

142 self._zC = z = radians(z) 

143 self._b0 = b # 𝜆C 

144 self._a0 = a # 𝜑C 

145 return a, b, z 

146 

147 @Property_RO 

148 def _sca0(self): # 𝜑0 

149 return sincos2(self._a0) 

150 

151 @Property_RO 

152 def _scb0(self): # 𝜆0 

153 return sincos2(self._b0) 

154 

155 @Property_RO 

156 def _scap(self): # 𝜑p 

157 return sincos2(self._ap) 

158 

159 @Property_RO 

160 def _scbp(self): # 𝜆p 

161 return sincos2(self._bp) 

162 

163 @Property_RO 

164 def _scg0(self): # 𝛾0 

165 return sincos2(self._g0) 

166 

167 @Property_RO 

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

169 return sincos2(self._zC) 

170 

171 def toStr(self, **kwds): 

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

173 ''' 

174 t = self._args_kwds 

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

176 

177 @property_RO 

178 def xy0(self): 

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

180 ''' 

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

182 

183 

184class _HotineOblique123(_HotineObliqueBase): 

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

186 ''' 

187 _Rk0 = \ 

188 _radius = R_M 

189 

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

191 self._ap = ap 

192 self._bp = bp = _PI(bp) 

193 self._b0 = _PI(bp + PI_2) 

194 

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

196 

197 if radius != _HotineOblique123._radius or \ 

198 self.k0 != _HotineObliqueBase._k0: 

199 self._radius = R = Radius(radius) 

200 self._Rk0 = R * self.k0 

201 

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

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

204 

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

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

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

208 is scalar, ignored otherwise. 

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

210 

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

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

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

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

215 ''' 

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

217 

218 d = Lamd(lon) - self._b0 

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

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

221 

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

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

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

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

226 

227 k = _1_0 - A**2 

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

229 x = u + self._x0 

230 y = v + self._y0 

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

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

233 

234 def _kwds(self, **unused): 

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

236 

237 @Property_RO 

238 def latlonP(self): 

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

240 ''' 

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

242 

243 @property_RO 

244 def radius(self): 

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

246 ''' 

247 return self._radius 

248 

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

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

251 

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

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

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

255 

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

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

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

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

260 ''' 

261 t = _1_over(self._Rk0) 

262 u = x - self._x0 

263 v = y - self._y0 

264 e = u * t 

265 n = v * t 

266 

267 sx, cx = sincos2(e) 

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

269 sp, cp = self._scap 

270 

271 a = degrees(asin1(fdot_(sp, ty, cp, sx / cy))) if cy else INF 

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

273 return Hotine8Tuple(x, y, u, v, 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 ca1_sa2 = ca1 * sa2 

299 sa1_ca2 = sa1 * ca2 

300 bp = atan2(fdot_(ca1_sa2, cb1, -sa1_ca2, cb2), 

301 fdot_(sa1_ca2, sb2, -ca1_sa2, 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 _xyuvIOGP = False 

376 

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

378 

379 if earth is not _HotineObliqueAB._earth: 

380 self._earth = _ellipsoidal_datum(earth) 

381 

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

383 

384 @Property_RO 

385 def _ABD21EFt0(self): # 7-tuple 

386 sa0, ca0 = self._sca0 

387 a0 = self._a0 

388 E = self.ellipsoid 

389 e2 = E.e2not 

390 r = _1_0 - e2 

391 B = ca0**4 * e2 / r 

392 B = sqrt(_1_0 + B) 

393 s = _1_0 - sa0**2 * e2 

394 r = sqrt(r) * B 

395 A = r * E.a * self.k0 / s 

396 D = F = r / (sqrt(s) * ca0) 

397 D2 = D**2 

398 if D2 > 1: 

399 D21 = sqrt(D2 - _1_0) 

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

401 else: # PYCHOK no cover 

402 D21 = _0_0 

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

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

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

406 

407 @property_RO 

408 def ABmap(self): 

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

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

411 ''' 

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

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

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

415 self._a0, self._b0, self._g0, t0, 

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

417 

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

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

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

421# d.update(locals__) 

422# d.update(kwds) 

423# n = self.methodname(where) 

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

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

426 

427 @property_RO 

428 def earth(self): 

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

430 ''' 

431 return self._earth 

432 

433 @Property_RO 

434 def ellipsoid(self): 

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

436 ''' 

437 return self.earth.ellipsoid 

438 

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

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

441 

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

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

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

445 is scalar, ignored otherwise. 

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

447 

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

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

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

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

452 ''' 

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

454 

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

456 

457 a = Phid(lat) 

458 b = Lamd(lon) 

459 Bd = _PI(b - self._b0) * B 

460 E = self.ellipsoid 

461 sg0, cg0 = self._scg0 

462 

463 sa, ca, V, cBd = sincos2_(a, Bd) # phi 

464 if ca: # PYLINT W0125 

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

466 S, T = _ST_Q(Q) 

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

468 

469 v = _d_s(_1_0, U) 

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

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

472 elif ca: # PYLINT W0125 

473 pass # 0.0 or INF 

474 else: # PYCHOK no cover lat +/-90 

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

476 v = log(tan(PI_4 + v)) * A_B 

477 

478 if cBd: # PYLINT W0125 

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

480 u = atan2(u, cBd) * A_B # Snyder 9-31 XXX atan1 

481 else: # PYCHOK no cover 

482 u = A * Bd # Note after 9-31 

483 

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

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

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

487 else: # PYCHOK no cover 

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

489 v = _0_0 

490 

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

492 if x > 0: 

493 k = E.a * ca * cBd # Snyder 9-32 

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

495 else: # PYCHOK no cover 

496 k = NAN if x < 0 else _0_0 

497 

498 s, c = self._scg0 if self.xyuvIOGP else self._sczC # PYLINT W0125 

499 x = fdot_(v, c, u, s, start=self._x0) 

500 y = fdot_(u, c, -v, s, start=self._y0) 

501 

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

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

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

505 

506 def _kwds(self, **unused): 

507 return _HotineObliqueBase._kwds(self, earth=self.ellipsoid) 

508 

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

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

511 

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

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

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

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

516 

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

518 

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

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

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

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

523 ''' 

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

525 

526 e, n = (x - self._x0), (y - self._y0) 

527 s, c = self._scg0 if self.xyuvIOGP else self._sczC # PYLINT W0125 

528 v, u = (fdot_(e, c, -n, s), # 9-40 

529 fdot_(n, c, e, s)) # 9-41 

530 u += uC # IOGP uC 

531 V, cBu_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 = _over(E, sqrt(t)) 

540 t = pow(t, _1_over(B)) 

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

542 n = Int(triPhi=triPhi) 

543 e = self.ellipsoid.e 

544 if e and n > 0: 

545 e_2 = e * _0_5 

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

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

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

549 d = fabs(a - X) 

550 if d < _TOL: 

551 i += 1 

552 break 

553 X = a 

554 else: # PYCHOK no cover 

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

556 e = Fmt.no_convergence(d, _TOL) 

557 raise HotineError(e, txt=t) 

558 

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

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

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

562 

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

564 b = atan2(b, cBu_A) 

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

566 

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

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

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

570 

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

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

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

574 

575 @Property_RO 

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

577 e2 = self.ellipsoid.e2not 

578 e4 = e2 * e2 

579 e6 = e4 * e2 

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

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

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

583 e6 * e2 * (4279 / 161280)) 

584 

585 @property_doc_(' the IOGP I{xy-uv} C{gamma0} transform instead of C{azimuthC} (C{bool})') 

586 def xyuvIOGP(self): 

587 return self._xyuvIOGP 

588 

589 @xyuvIOGP.setter # PYCHOK setter! 

590 def xyuvIOGP(self, IOGP): 

591 self._xyuvIOGP = bool(IOGP) 

592 

593 

594class HotineObliqueA(_HotineObliqueAB): 

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

596 ''' 

597 

598 def __init__(self, latlon1, latlon2, earth=_WGS84, **k0_lat0_x0_y0_name): 

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

600 

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

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

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

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

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

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, **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, **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 k0_lat0_x0_y0_name: Central scale factor C{B{k0}=1.0} (C{scalar}), origin 

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

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

670 ''' 

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

672 _HotineObliqueAB.__init__(self, earth, **k0_lat0_x0_y0_name) 

673 szC, czC = self._sczC 

674 

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

676 A_B = _over(A, B) 

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

678 

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

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

681 

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

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

684 (atan2(D21, czC) * A_B) 

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

686 

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

688 

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

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

691 

692 @property_RO 

693 def latlonC(self): 

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

695 ''' 

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

697 

698 

699class HotineError(_ValueError): 

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

701 ''' 

702 pass # PYLINT W0107 

703 

704 

705class Hotine8Tuple(_NamedTuple): 

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

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

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

709 

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

711 ''' 

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

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

714 

715 

716class HotineAB14Tuple(_NamedTuple): 

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

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

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

720 ''' 

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

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

723 

724 

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

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

727 

728 

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

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

731 

732 

733def _ST_Q(Q): 

734 _Q = _1_over(Q) 

735 S = (Q - _Q) * _0_5 

736 T = (Q + _Q) * _0_5 

737 return S, T 

738 

739 

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

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

742 t = _d_s(_1_0, e * s)