1

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

3

4 u'''Formulary of basic geodesy functions and approximations.

5'''

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

7 from __future__ import division as _; del _ # PYCHOK semicolon

8

9# from pygeodesy.cartesianBase import CartesianBase # _MODS

10 from pygeodesy.constants import EPS, EPS0, EPS1, PI, PI2, PI3, PI_2, R_M, \

11 _0_0s, float0_, isnon0, remainder, _umod_PI2, \

12 _0_0, _0_125, _0_25, _0_5, _1_0, _2_0, _4_0, \

13 _32_0, _90_0, _180_0, _360_0

14 from pygeodesy.datums import Datum, Ellipsoid, _ellipsoidal_datum, \

15 _mean_radius, _spherical_datum, _WGS84, _EWGS84

16# from pygeodesy.ellipsoids import Ellipsoid, _EWGS84 # from .datums

17 from pygeodesy.errors import IntersectionError, LimitError, limiterrors, \

18 _TypeError, _ValueError, _xattr, _xError, \

19 _xcallable,_xkwds, _xkwds_pop2

20 from pygeodesy.fmath import euclid, hypot, hypot2, sqrt0

21 from pygeodesy.fsums import fsumf_, Fmt, unstr

22# from pygeodesy.internals import _dunder_nameof # from .named

23 from pygeodesy.interns import _delta_, _distant_, _inside_, _SPACE_, _too_

24 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

25 from pygeodesy.named import _name__, _name2__, _NamedTuple, _xnamed, \

26 _dunder_nameof

27 from pygeodesy.namedTuples import Bearing2Tuple, Distance4Tuple, LatLon2Tuple, \

28 Intersection3Tuple, PhiLam2Tuple, Vector3Tuple

29# from pygeodesy.streprs import Fmt, unstr # from .fsums

30# from pygeodesy.triaxials import _hartzell3 # _MODS

31 from pygeodesy.units import _isHeight, _isRadius, Bearing, Degrees_, Distance, \

32 Distance_, Height, Lamd, Lat, Lon, Meter_, Phid, \

33 Radians, Radians_, Radius, Radius_, Scalar, _100km

34 from pygeodesy.utily import acos1, atan2b, atan2d, degrees2m, _loneg, m2degrees, \

35 tan_2, sincos2, sincos2_, sincos2d_, _Wrap

36# from pygeodesy.vector3d import _otherV3d # _MODS

37# from pygeodesy.vector3dBase import _xyz_y_z3 # _MODS

38# from pygeodesy import ellipsoidalExact, ellipsoidalKarney, vector3d, \

39# sphericalNvector, sphericalTrigonometry # _MODS

40

41 from contextlib import contextmanager

42 from math import asin, atan, atan2, cos, degrees, fabs, radians, sin, sqrt # pow

43

44__all__ = _ALL_LAZY.formy

45__version__ = '24.07.29'

46

47_RADIANS2 = (PI / _180_0)**2 # degrees- to radians-squared

48_ratio_ = 'ratio'

49_xline_ = 'xline'

50

51

52 def _anti2(a, b, n_2, n, n2):

53 '''(INTERNAL) Helper for C{antipode} and C{antipode_}.

54 '''

55 r = remainder(a, n) if fabs(a) > n_2 else a

56 if r == a:

57 r = -r

58 b += n

59 if fabs(b) > n:

60 b = remainder(b, n2)

61 return float0_(r, b)

62

63

64 def antipode(lat, lon, **name):

65 '''Return the antipode, the point diametrically opposite to a given

66 point in C{degrees}.

67

68 @arg lat: Latitude (C{degrees}).

69 @arg lon: Longitude (C{degrees}).

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

71

72 @return: A L{LatLon2Tuple}C{(lat, lon)}.

73

74 @see: Functions L{antipode_} and L{normal} and U{Geosphere

75 <https://CRAN.R-Project.org/web/packages/geosphere/geosphere.pdf>}.

76 '''

77 return LatLon2Tuple(*_anti2(lat, lon, _90_0, _180_0, _360_0), **name)

78

79

80 def antipode_(phi, lam, **name):

81 '''Return the antipode, the point diametrically opposite to a given

82 point in C{radians}.

83

84 @arg phi: Latitude (C{radians}).

85 @arg lam: Longitude (C{radians}).

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

87

88 @return: A L{PhiLam2Tuple}C{(phi, lam)}.

89

90 @see: Functions L{antipode} and L{normal_} and U{Geosphere

91 <https://CRAN.R-Project.org/web/packages/geosphere/geosphere.pdf>}.

92 '''

93 return PhiLam2Tuple(*_anti2(phi, lam, PI_2, PI, PI2), **name)

94

95

96 def bearing(lat1, lon1, lat2, lon2, **final_wrap):

97 '''Compute the initial or final bearing (forward or reverse azimuth)

98 between two (spherical) points.

99

100 @arg lat1: Start latitude (C{degrees}).

101 @arg lon1: Start longitude (C{degrees}).

102 @arg lat2: End latitude (C{degrees}).

103 @arg lon2: End longitude (C{degrees}).

104 @kwarg final_wrap: Optional keyword arguments for function

105 L{pygeodesy.bearing_}.

106

107 @return: Initial or final bearing (compass C{degrees360}) or zero if

108 both points coincide.

109 '''

110 r = bearing_(Phid(lat1=lat1), Lamd(lon1=lon1),

111 Phid(lat2=lat2), Lamd(lon2=lon2), **final_wrap)

112 return degrees(r)

113

114

115 def bearing_(phi1, lam1, phi2, lam2, final=False, wrap=False):

116 '''Compute the initial or final bearing (forward or reverse azimuth) between

117 two (spherical) points.

118

119 @arg phi1: Start latitude (C{radians}).

120 @arg lam1: Start longitude (C{radians}).

121 @arg phi2: End latitude (C{radians}).

122 @arg lam2: End longitude (C{radians}).

123 @kwarg final: If C{True}, return the final, otherwise the initial bearing

124 (C{bool}).

125 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{phi2}} and

126 B{C{lam2}} (C{bool}).

127

128 @return: Initial or final bearing (compass C{radiansPI2}) or zero if both

129 are coincident.

130

131 @see: U{Bearing<https://www.Movable-Type.co.UK/scripts/latlong.html>}, U{Course

132 between two points<https://www.EdWilliams.org/avform147.htm#Crs>} and

133 U{Bearing Between Two Points<https://web.Archive.org/web/20020630205931/

134 https://MathForum.org/library/drmath/view/55417.html>}.

135 '''

136 db, phi2, lam2 = _Wrap.philam3(lam1, phi2, lam2, wrap)

137 if final: # swap plus PI

138 phi1, lam1, phi2, lam2, db = phi2, lam2, phi1, lam1, -db

139 r = PI3

140 else:

141 r = PI2

142 sa1, ca1, sa2, ca2, sdb, cdb = sincos2_(phi1, phi2, db)

143

144 x = ca1 * sa2 - sa1 * ca2 * cdb

145 y = sdb * ca2

146 return _umod_PI2(atan2(y, x) + r) # .utily.wrapPI2

147

148

149 def _bearingTo2(p1, p2, wrap=False): # for points.ispolar, sphericalTrigonometry.areaOf

150 '''(INTERNAL) Compute initial and final bearing.

151 '''

152 try: # for LatLon_ and ellipsoidal LatLon

153 return p1.bearingTo2(p2, wrap=wrap)

154 except AttributeError:

155 pass

156 # XXX spherical version, OK for ellipsoidal ispolar?

157 t = p1.philam + p2.philam

158 return Bearing2Tuple(degrees(bearing_(*t, final=False, wrap=wrap)),

159 degrees(bearing_(*t, final=True, wrap=wrap)),

160 name__=_bearingTo2)

161

162

163 def compassAngle(lat1, lon1, lat2, lon2, adjust=True, wrap=False):

164 '''Return the angle from North for the direction vector M{(lon2 - lon1,

165 lat2 - lat1)} between two points.

166

167 Suitable only for short, not near-polar vectors up to a few hundred

168 Km or Miles. Use function L{pygeodesy.bearing} for longer vectors.

169

170 @arg lat1: From latitude (C{degrees}).

171 @arg lon1: From longitude (C{degrees}).

172 @arg lat2: To latitude (C{degrees}).

173 @arg lon2: To longitude (C{degrees}).

174 @kwarg adjust: Adjust the longitudinal delta by the cosine of the

175 mean latitude (C{bool}).

176 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}}

177 and B{C{lon2}} (C{bool}).

178

179 @return: Compass angle from North (C{degrees360}).

180

181 @note: Courtesy of Martin Schultz.

182

183 @see: U{Local, flat earth approximation

184 <https://www.EdWilliams.org/avform.htm#flat>}.

185 '''

186 d_lon, lat2, lon2 = _Wrap.latlon3(lon1, lat2, lon2, wrap)

187 if adjust: # scale delta lon

188 d_lon *= _scale_deg(lat1, lat2)

189 return atan2b(d_lon, lat2 - lat1)

190

191

192 def cosineAndoyerLambert(lat1, lon1, lat2, lon2, datum=_WGS84, wrap=False):

193 '''Compute the distance between two (ellipsoidal) points using the U{Andoyer-Lambert

194 <https://books.google.com/books?id=x2UiAQAAIAAJ>} correction of the U{Law of

195 Cosines<https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>} formula.

196

197 @arg lat1: Start latitude (C{degrees}).

198 @arg lon1: Start longitude (C{degrees}).

199 @arg lat2: End latitude (C{degrees}).

200 @arg lon2: End longitude (C{degrees}).

201 @kwarg datum: Datum (L{Datum}) or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or

202 L{a_f2Tuple}) to use.

203 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}} and

204 B{C{lon2}} (C{bool}).

205

206 @return: Distance (C{meter}, same units as the B{C{datum}}'s ellipsoid axes or

207 C{radians} if C{B{datum} is None}).

208

209 @raise TypeError: Invalid B{C{datum}}.

210

211 @see: Functions L{cosineAndoyerLambert_}, L{cosineForsytheAndoyerLambert},

212 L{cosineLaw}, L{equirectangular}, L{euclidean}, L{flatLocal}/L{hubeny},

213 L{flatPolar}, L{haversine}, L{thomas} and L{vincentys} and method

214 L{Ellipsoid.distance2}.

215 '''

216 return _dE(cosineAndoyerLambert_, datum, wrap, lat1, lon1, lat2, lon2)

217

218

219 def cosineAndoyerLambert_(phi2, phi1, lam21, datum=_WGS84):

220 '''Compute the I{angular} distance between two (ellipsoidal) points using the U{Andoyer-Lambert

221 <https://books.google.com/books?id=x2UiAQAAIAAJ>} correction of the U{Law of

222 Cosines<https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>} formula.

223

224 @arg phi2: End latitude (C{radians}).

225 @arg phi1: Start latitude (C{radians}).

226 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).

227 @kwarg datum: Datum (L{Datum}) or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or

228 L{a_f2Tuple}) to use.

229

230 @return: Angular distance (C{radians}).

231

232 @raise TypeError: Invalid B{C{datum}}.

233

234 @see: Functions L{cosineAndoyerLambert}, L{cosineForsytheAndoyerLambert_},

235 L{cosineLaw_}, L{euclidean_}, L{flatLocal_}/L{hubeny_}, L{flatPolar_},

236 L{haversine_}, L{thomas_} and L{vincentys_} and U{Geodesy-PHP

237 <https://GitHub.com/jtejido/geodesy-php/blob/master/src/Geodesy/Distance/

238 AndoyerLambert.php>}.

239 '''

240 s2, c2, s1, c1, r, c21 = _sincosa6(phi2, phi1, lam21)

241 if isnon0(c1) and isnon0(c2):

242 E = _ellipsoidal(datum, cosineAndoyerLambert_)

243 if E.f: # ellipsoidal

244 r2 = atan2(E.b_a * s2, c2)

245 r1 = atan2(E.b_a * s1, c1)

246 s2, c2, s1, c1 = sincos2_(r2, r1)

247 r = acos1(s1 * s2 + c1 * c2 * c21)

248 if r:

249 sr, _, sr_2, cr_2 = sincos2_(r, r * _0_5)

250 if isnon0(sr_2) and isnon0(cr_2):

251 s = (sr + r) * ((s1 - s2) / sr_2)**2

252 c = (sr - r) * ((s1 + s2) / cr_2)**2

253 r += (c - s) * E.f * _0_125

254 return r

255

256

257 def cosineForsytheAndoyerLambert(lat1, lon1, lat2, lon2, datum=_WGS84, wrap=False):

258 '''Compute the distance between two (ellipsoidal) points using the U{Forsythe-Andoyer-Lambert

259 <https://www2.UNB.Ca/gge/Pubs/TR77.pdf>} correction of the U{Law of Cosines

260 <https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>} formula.

261

262 @arg lat1: Start latitude (C{degrees}).

263 @arg lon1: Start longitude (C{degrees}).

264 @arg lat2: End latitude (C{degrees}).

265 @arg lon2: End longitude (C{degrees}).

266 @kwarg datum: Datum (L{Datum}) or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or

267 L{a_f2Tuple}) to use.

268 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}} and

269 B{C{lon2}} (C{bool}).

270

271 @return: Distance (C{meter}, same units as the B{C{datum}}'s ellipsoid axes or

272 C{radians} if C{B{datum} is None}).

273

274 @raise TypeError: Invalid B{C{datum}}.

275

276 @see: Functions L{cosineForsytheAndoyerLambert_}, L{cosineAndoyerLambert},

277 L{cosineLaw}, L{equirectangular}, L{euclidean}, L{flatLocal}/L{hubeny},

278 L{flatPolar}, L{haversine}, L{thomas} and L{vincentys} and method

279 L{Ellipsoid.distance2}.

280 '''

281 return _dE(cosineForsytheAndoyerLambert_, datum, wrap, lat1, lon1, lat2, lon2)

282

283

284 def cosineForsytheAndoyerLambert_(phi2, phi1, lam21, datum=_WGS84):

285 '''Compute the I{angular} distance between two (ellipsoidal) points using the

286 U{Forsythe-Andoyer-Lambert<https://www2.UNB.Ca/gge/Pubs/TR77.pdf>} correction of

287 the U{Law of Cosines<https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>}

288 formula.

289

290 @arg phi2: End latitude (C{radians}).

291 @arg phi1: Start latitude (C{radians}).

292 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).

293 @kwarg datum: Datum (L{Datum}) or ellipsoid to use (L{Ellipsoid},

294 L{Ellipsoid2} or L{a_f2Tuple}).

295

296 @return: Angular distance (C{radians}).

297

298 @raise TypeError: Invalid B{C{datum}}.

299

300 @see: Functions L{cosineForsytheAndoyerLambert}, L{cosineAndoyerLambert_},

301 L{cosineLaw_}, L{euclidean_}, L{flatLocal_}/L{hubeny_}, L{flatPolar_},

302 L{haversine_}, L{thomas_} and L{vincentys_} and U{Geodesy-PHP

303 <https://GitHub.com/jtejido/geodesy-php/blob/master/src/Geodesy/

304 Distance/ForsytheCorrection.php>}.

305 '''

306 s2, c2, s1, c1, r, _ = _sincosa6(phi2, phi1, lam21)

307 if r and isnon0(c1) and isnon0(c2):

308 E = _ellipsoidal(datum, cosineForsytheAndoyerLambert_)

309 if E.f: # ellipsoidal

310 sr, cr, s2r, _ = sincos2_(r, r * 2)

311 if isnon0(sr) and fabs(cr) < EPS1:

312 s = (s1 + s2)**2 / (1 + cr)

313 t = (s1 - s2)**2 / (1 - cr)

314 x = s + t

315 y = s - t

316

317 s = 8 * r**2 / sr

318 a = 64 * r + s * cr * 2 # 16 * r**2 / tan(r)

319 d = 48 * sr + s # 8 * r**2 / tan(r)

320 b = -2 * d

321 e = 30 * s2r

322 c = fsumf_(30 * r, e * _0_5, s * cr) # 8 * r**2 / tan(r)

323 t = fsumf_( a * x, e * y**2, b * y, -c * x**2, d * x * y)

324

325 r += fsumf_(-r * x, 3 * y * sr, t * E.f / _32_0) * E.f * _0_25

326 return r

327

328

329 def cosineLaw(lat1, lon1, lat2, lon2, radius=R_M, wrap=False):

330 '''Compute the distance between two points using the U{spherical Law of Cosines

331 <https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>} formula.

332

333 @arg lat1: Start latitude (C{degrees}).

334 @arg lon1: Start longitude (C{degrees}).

335 @arg lat2: End latitude (C{degrees}).

336 @arg lon2: End longitude (C{degrees}).

337 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})

338 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or

339 L{a_f2Tuple}) to use.

340 @kwarg wrap: If C{True}, wrap or I{normalize} and B{C{lat2}}

341 and B{C{lon2}} (C{bool}).

342

343 @return: Distance (C{meter}, same units as B{C{radius}} or the

344 ellipsoid or datum axes).

345

346 @raise TypeError: Invalid B{C{radius}}.

347

348 @see: Functions L{cosineLaw_}, L{cosineAndoyerLambert},

349 L{cosineForsytheAndoyerLambert}, L{equirectangular}, L{euclidean},

350 L{flatLocal}/L{hubeny}, L{flatPolar}, L{haversine}, L{thomas} and

351 L{vincentys} and method L{Ellipsoid.distance2}.

352

353 @note: See note at function L{vincentys_}.

354 '''

355 return _dS(cosineLaw_, radius, wrap, lat1, lon1, lat2, lon2)

356

357

358 def cosineLaw_(phi2, phi1, lam21):

359 '''Compute the I{angular} distance between two points using the U{spherical Law of

360 Cosines<https://www.Movable-Type.co.UK/scripts/latlong.html#cosine-law>} formula.

361

362 @arg phi2: End latitude (C{radians}).

363 @arg phi1: Start latitude (C{radians}).

364 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).

365

366 @return: Angular distance (C{radians}).

367

368 @see: Functions L{cosineLaw}, L{cosineAndoyerLambert_},

369 L{cosineForsytheAndoyerLambert_}, L{euclidean_},

370 L{flatLocal_}/L{hubeny_}, L{flatPolar_}, L{haversine_},

371 L{thomas_} and L{vincentys_}.

372

373 @note: See note at function L{vincentys_}.

374 '''

375 return _sincosa6(phi2, phi1, lam21)[4]

376

377

378 def _d3(wrap, lat1, lon1, lat2, lon2):

379 '''(INTERNAL) Helper for _dE, _dS and _eA.

380 '''

381 if wrap:

382 d_lon, lat2, _ = _Wrap.latlon3(lon1, lat2, lon2, wrap)

383 return radians(lat2), Phid(lat1=lat1), radians(d_lon)

384 else: # for backward compaibility

385 return Phid(lat2=lat2), Phid(lat1=lat1), Phid(d_lon=lon2 - lon1)

386

387

388 def _dE(func_, earth, *wrap_lls):

389 '''(INTERNAL) Helper for ellipsoidal distances.

390 '''

391 E = _ellipsoidal(earth, func_)

392 r = func_(*_d3(*wrap_lls), datum=E)

393 return r * E.a

394

395

396 def _dS(func_, radius, *wrap_lls, **adjust):

397 '''(INTERNAL) Helper for spherical distances.

398 '''

399 r = func_(*_d3(*wrap_lls), **adjust)

400 if radius is not R_M:

401 _, lat1, _, lat2, _ = wrap_lls

402 radius = _mean_radius(radius, lat1, lat2)

403 return r * radius

404

405

406 def _eA(excess_, radius, *wrap_lls):

407 '''(INTERNAL) Helper for spherical excess or area.

408 '''

409 r = excess_(*_d3(*wrap_lls))

410 if radius:

411 _, lat1, _, lat2, _ = wrap_lls

412 r *= _mean_radius(radius, lat1, lat2)**2

413 return r

414

415

416 def _ellipsoidal(earth, where):

417 '''(INTERNAL) Helper for distances.

418 '''

419 return _EWGS84 if earth in (_WGS84, _EWGS84) else (

420 earth if isinstance(earth, Ellipsoid) else

421 (earth if isinstance(earth, Datum) else # PYCHOK indent

422 _ellipsoidal_datum(earth, name__=where)).ellipsoid)

423

424

425 def equirectangular(lat1, lon1, lat2, lon2, radius=R_M, **adjust_limit_wrap):

426 '''Compute the distance between two points using the U{Equirectangular Approximation

427 / Projection<https://www.Movable-Type.co.UK/scripts/latlong.html#equirectangular>}.

428

429 @arg lat1: Start latitude (C{degrees}).

430 @arg lon1: Start longitude (C{degrees}).

431 @arg lat2: End latitude (C{degrees}).

432 @arg lon2: End longitude (C{degrees}).

433 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum}) or ellipsoid

434 (L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).

435 @kwarg adjust_limit_wrap: Optional keyword arguments for function L{equirectangular4}.

436

437 @return: Distance (C{meter}, same units as B{C{radius}} or the ellipsoid or datum axes).

438

439 @raise TypeError: Invalid B{C{radius}}.

440

441 @see: Function L{equirectangular4} for more details, the available B{C{options}},

442 errors, restrictions and other, approximate or accurate distance functions.

443 '''

444 d = sqrt(equirectangular4(Lat(lat1=lat1), Lon(lon1=lon1),

445 Lat(lat2=lat2), Lon(lon2=lon2),

446 **adjust_limit_wrap).distance2) # PYCHOK 4 vs 2-3

447 return degrees2m(d, radius=_mean_radius(radius, lat1, lat2))

448

449

450 def _equirectangular(lat1, lon1, lat2, lon2, **adjust_limit_wrap):

451 '''(INTERNAL) Helper for the L{frechet._FrechetMeterRadians}

452 and L{hausdorff._HausdorffMeterRedians} classes.

453 '''

454 return equirectangular4(lat1, lon1, lat2, lon2, **adjust_limit_wrap).distance2 * _RADIANS2

455

456

457 def equirectangular4(lat1, lon1, lat2, lon2, adjust=True, limit=45, wrap=False):

458 '''Compute the distance between two points using the U{Equirectangular Approximation

459 / Projection<https://www.Movable-Type.co.UK/scripts/latlong.html#equirectangular>}.

460

461 This approximation is valid for short distance of several hundred Km or Miles, see

462 the B{C{limit}} keyword argument and L{LimitError}.

463

464 @arg lat1: Start latitude (C{degrees}).

465 @arg lon1: Start longitude (C{degrees}).

466 @arg lat2: End latitude (C{degrees}).

467 @arg lon2: End longitude (C{degrees}).

468 @kwarg adjust: Adjust the wrapped, unrolled longitudinal delta by the cosine of the mean

469 latitude (C{bool}).

470 @kwarg limit: Optional limit for lat- and longitudinal deltas (C{degrees}) or C{None}

471 or C{0} for unlimited.

472 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}} and B{C{lon2}}

473 (C{bool}).

474

475 @return: A L{Distance4Tuple}C{(distance2, delta_lat, delta_lon, unroll_lon2)}

476 in C{degrees squared}.

477

478 @raise LimitError: If the lat- and/or longitudinal delta exceeds the B{C{-limit..limit}}

479 range and L{limiterrors<pygeodesy.limiterrors>} is C{True}.

480

481 @see: U{Local, flat earth approximation<https://www.EdWilliams.org/avform.htm#flat>},

482 functions L{equirectangular}, L{cosineAndoyerLambert}, L{cosineForsytheAndoyerLambert},

483 L{cosineLaw}, L{euclidean}, L{flatLocal}/L{hubeny}, L{flatPolar}, L{haversine},

484 L{thomas} and L{vincentys} and methods L{Ellipsoid.distance2}, C{LatLon.distanceTo*}

485 and C{LatLon.equirectangularTo}.

486 '''

487 d_lon, lat2, ulon2 = _Wrap.latlon3(lon1, lat2, lon2, wrap)

488 d_lat = lat2 - lat1

489

490 if limit and limit > 0 and limiterrors():

491 d = max(fabs(d_lat), fabs(d_lon))

492 if d > limit:

493 t = _SPACE_(_delta_, Fmt.PAREN_g(d), Fmt.exceeds_limit(limit))

494 s = unstr(equirectangular4, lat1, lon1, lat2, lon2,

495 limit=limit, wrap=wrap)

496 raise LimitError(s, txt=t)

497

498 if adjust: # scale delta lon

499 d_lon *= _scale_deg(lat1, lat2)

500

501 d2 = hypot2(d_lat, d_lon) # degrees squared!

502 return Distance4Tuple(d2, d_lat, d_lon, ulon2 - lon2)

503

504

505 def euclidean(lat1, lon1, lat2, lon2, radius=R_M, adjust=True, wrap=False):

506 '''Approximate the C{Euclidean} distance between two (spherical) points.

507

508 @arg lat1: Start latitude (C{degrees}).

509 @arg lon1: Start longitude (C{degrees}).

510 @arg lat2: End latitude (C{degrees}).

511 @arg lon2: End longitude (C{degrees}).

512 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})

513 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or

514 L{a_f2Tuple}) to use.

515 @kwarg adjust: Adjust the longitudinal delta by the cosine of

516 the mean latitude (C{bool}).

517 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}}

518 and B{C{lon2}} (C{bool}).

519

520 @return: Distance (C{meter}, same units as B{C{radius}} or the

521 ellipsoid or datum axes).

522

523 @raise TypeError: Invalid B{C{radius}}.

524

525 @see: U{Distance between two (spherical) points

526 <https://www.EdWilliams.org/avform.htm#Dist>}, functions L{euclid},

527 L{euclidean_}, L{cosineAndoyerLambert}, L{cosineForsytheAndoyerLambert},

528 L{cosineLaw}, L{equirectangular}, L{flatLocal}/L{hubeny}, L{flatPolar},

529 L{haversine}, L{thomas} and L{vincentys} and methods L{Ellipsoid.distance2},

530 C{LatLon.distanceTo*} and C{LatLon.equirectangularTo}.

531 '''

532 return _dS(euclidean_, radius, wrap, lat1, lon1, lat2, lon2, adjust=adjust)

533

534

535 def euclidean_(phi2, phi1, lam21, adjust=True):

536 '''Approximate the I{angular} C{Euclidean} distance between two (spherical) points.

537

538 @arg phi2: End latitude (C{radians}).

539 @arg phi1: Start latitude (C{radians}).

540 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).

541 @kwarg adjust: Adjust the longitudinal delta by the cosine

542 of the mean latitude (C{bool}).

543

544 @return: Angular distance (C{radians}).

545

546 @see: Functions L{euclid}, L{euclidean}, L{cosineAndoyerLambert_},

547 L{cosineForsytheAndoyerLambert_}, L{cosineLaw_},

548 L{flatLocal_}/L{hubeny_}, L{flatPolar_}, L{haversine_},

549 L{thomas_} and L{vincentys_}.

550 '''

551 if adjust:

552 lam21 *= _scale_rad(phi2, phi1)

553 return euclid(phi2 - phi1, lam21)

554

555

556 def excessAbc_(A, b, c):

557 '''Compute the I{spherical excess} C{E} of a (spherical) triangle from two sides

558 and the included (small) angle.

559

560 @arg A: An interior triangle angle (C{radians}).

561 @arg b: Frist adjacent triangle side (C{radians}).

562 @arg c: Second adjacent triangle side (C{radians}).

563

564 @return: Spherical excess (C{radians}).

565

566 @raise UnitError: Invalid B{C{A}}, B{C{b}} or B{C{c}}.

567

568 @see: Functions L{excessGirard_}, L{excessLHuilier_} and U{Spherical

569 trigonometry<https://WikiPedia.org/wiki/Spherical_trigonometry>}.

570 '''

571 A = Radians_(A=A)

572 b = Radians_(b=b) * _0_5

573 c = Radians_(c=c) * _0_5

574

575 sA, cA, sb, cb, sc, cc = sincos2_(A, b, c)

576 return atan2(sA * sb * sc, cb * cc + cA * sb * sc) * _2_0

577

578

579 def excessCagnoli_(a, b, c):

580 '''Compute the I{spherical excess} C{E} of a (spherical) triangle using U{Cagnoli's

581 <https://Zenodo.org/record/35392>} (D.34) formula.

582

583 @arg a: First triangle side (C{radians}).

584 @arg b: Second triangle side (C{radians}).

585 @arg c: Third triangle side (C{radians}).

586

587 @return: Spherical excess (C{radians}).

588

589 @raise UnitError: Invalid B{C{a}}, B{C{b}} or B{C{c}}.

590

591 @see: Function L{excessLHuilier_} and U{Spherical trigonometry

592 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.

593 '''

594 a = Radians_(a=a)

595 b = Radians_(b=b)

596 c = Radians_(c=c)

597

598 s = fsumf_(a, b, c) * _0_5

599 _s = sin

600 r = _s(s) * _s(s - a) * _s(s - b) * _s(s - c)

601 c = cos(a * _0_5) * cos(b * _0_5) * cos(c * _0_5)

602 r = asin(sqrt(r) * _0_5 / c) if c and r > 0 else _0_0

603 return Radians(Cagnoli=r * _2_0)

604

605

606 def excessGirard_(A, B, C):

607 '''Compute the I{spherical excess} C{E} of a (spherical) triangle using U{Girard's

608 <https://MathWorld.Wolfram.com/GirardsSphericalExcessFormula.html>} formula.

609

610 @arg A: First interior triangle angle (C{radians}).

611 @arg B: Second interior triangle angle (C{radians}).

612 @arg C: Third interior triangle angle (C{radians}).

613

614 @return: Spherical excess (C{radians}).

615

616 @raise UnitError: Invalid B{C{A}}, B{C{B}} or B{C{C}}.

617

618 @see: Function L{excessLHuilier_} and U{Spherical trigonometry

619 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.

620 '''

621 return Radians(Girard=fsumf_(Radians_(A=A),

622 Radians_(B=B),

623 Radians_(C=C), -PI))

624

625

626 def excessLHuilier_(a, b, c):

627 '''Compute the I{spherical excess} C{E} of a (spherical) triangle using U{L'Huilier's

628 <https://MathWorld.Wolfram.com/LHuiliersTheorem.html>}'s Theorem.

629

630 @arg a: First triangle side (C{radians}).

631 @arg b: Second triangle side (C{radians}).

632 @arg c: Third triangle side (C{radians}).

633

634 @return: Spherical excess (C{radians}).

635

636 @raise UnitError: Invalid B{C{a}}, B{C{b}} or B{C{c}}.

637

638 @see: Function L{excessCagnoli_}, L{excessGirard_} and U{Spherical

639 trigonometry<https://WikiPedia.org/wiki/Spherical_trigonometry>}.

640 '''

641 a = Radians_(a=a)

642 b = Radians_(b=b)

643 c = Radians_(c=c)

644

645 s = fsumf_(a, b, c) * _0_5

646 _t = tan_2

647 r = _t(s) * _t(s - a) * _t(s - b) * _t(s - c)

648 r = atan(sqrt(r)) if r > 0 else _0_0

649 return Radians(LHuilier=r * _4_0)

650

651

652 def excessKarney(lat1, lon1, lat2, lon2, radius=R_M, wrap=False):

653 '''Compute the surface area of a (spherical) quadrilateral bounded by a

654 segment of a great circle, two meridians and the equator using U{Karney's

655 <https://MathOverflow.net/questions/97711/the-area-of-spherical-polygons>}

656 method.

657

658 @arg lat1: Start latitude (C{degrees}).

659 @arg lon1: Start longitude (C{degrees}).

660 @arg lat2: End latitude (C{degrees}).

661 @arg lon2: End longitude (C{degrees}).

662 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})

663 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or

664 L{a_f2Tuple}) or C{None}.

665 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll

666 B{C{lat2}} and B{C{lon2}} (C{bool}).

667

668 @return: Surface area, I{signed} (I{square} C{meter} or the same units as

669 B{C{radius}} I{squared}) or the I{spherical excess} (C{radians})

670 if C{B{radius}=0} or C{None}.

671

672 @raise TypeError: Invalid B{C{radius}}.

673

674 @raise UnitError: Invalid B{C{lat2}} or B{C{lat1}}.

675

676 @raise ValueError: Semi-circular longitudinal delta.

677

678 @see: Functions L{excessKarney_} and L{excessQuad}.

679 '''

680 return _eA(excessKarney_, radius, wrap, lat1, lon1, lat2, lon2)

681

682

683 def excessKarney_(phi2, phi1, lam21):

684 '''Compute the I{spherical excess} C{E} of a (spherical) quadrilateral bounded by

685 a segment of a great circle, two meridians and the equator using U{Karney's

686 <https://MathOverflow.net/questions/97711/the-area-of-spherical-polygons>}

687 method.

688

689 @arg phi2: End latitude (C{radians}).

690 @arg phi1: Start latitude (C{radians}).

691 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).

692

693 @return: Spherical excess, I{signed} (C{radians}).

694

695 @raise ValueError: Semi-circular longitudinal delta B{C{lam21}}.

696

697 @see: Function L{excessKarney} and U{Area of a spherical polygon

698 <https://MathOverflow.net/questions/97711/the-area-of-spherical-polygons>}.

699 '''

700 # from: Veness <https://www.Movable-Type.co.UK/scripts/latlong.html> Area

701 # method due to Karney: for each edge of the polygon,

702 #

703 # tan(Δλ / 2) · (tan(φ1 / 2) + tan(φ2 / 2))

704 # tan(E / 2) = -----------------------------------------

705 # 1 + tan(φ1 / 2) · tan(φ2 / 2)

706 #

707 # where E is the spherical excess of the trapezium obtained by extending

708 # the edge to the equator-circle vector for each edge (see also ***).

709 _t = tan_2

710 t2 = _t(phi2)

711 t1 = _t(phi1)

712 t = _t(lam21, lam21=None)

713 return Radians(Karney=atan2(t * (t1 + t2),

714 _1_0 + (t1 * t2)) * _2_0)

715

716

717# ***) Original post no longer available, following is a copy of the main part

718# <http://OSGeo-org.1560.x6.Nabble.com/Area-of-a-spherical-polygon-td3841625.html>

719#

720# The area of a polygon on a (unit) sphere is given by the spherical excess

721#

722# A = 2 * pi - sum(exterior angles)

723#

724# However this is badly conditioned if the polygon is small. In this case, use

725#

726# A = sum(S12{i, i+1}) over the edges of the polygon

727#

728# where S12 is the area of the quadrilateral bounded by an edge of the polygon,

729# two meridians and the equator, i.e. with vertices (phi1, lambda1), (phi2,

730# lambda2), (0, lambda1) and (0, lambda2). S12 is given by

731#

732# tan(S12 / 2) = tan(lambda21 / 2) * (tan(phi1 / 2) + tan(phi2 / 2)) /

733# (tan(phi1 / 2) * tan(phi2 / 2) + 1)

734#

735# = tan(lambda21 / 2) * tanh((Lamb(phi1) + Lamb(phi2)) / 2)

736#

737# where lambda21 = lambda2 - lambda1 and Lamb(x) is the Lambertian (or the

738# inverse Gudermannian) function

739#

740# Lambertian(x) = asinh(tan(x)) = atanh(sin(x)) = 2 * atanh(tan(x / 2))

741#

742# Notes: The formula for S12 is exact, except that...

743# - it is indeterminate if an edge is a semi-circle

744# - the formula for A applies only if the polygon does not include a pole

745# (if it does, then add +/- 2 * pi to the result)

746# - in the limit of small phi and lambda, S12 reduces to the trapezoidal

747# formula, S12 = (lambda2 - lambda1) * (phi1 + phi2) / 2

748# - I derived this result from the equation for the area of a spherical

749# triangle in terms of two edges and the included angle given by, e.g.

750# U{Todhunter, I. - Spherical Trigonometry (1871), Sec. 103, Eq. (2)

751# <http://Books.Google.com/books?id=3uBHAAAAIAAJ&pg=PA71>}

752# - I would be interested to know if this formula for S12 is already known

753# - Charles Karney

754

755

756 def excessQuad(lat1, lon1, lat2, lon2, radius=R_M, wrap=False):

757 '''Compute the surface area of a (spherical) quadrilateral bounded by a segment

758 of a great circle, two meridians and the equator.

759

760 @arg lat1: Start latitude (C{degrees}).

761 @arg lon1: Start longitude (C{degrees}).

762 @arg lat2: End latitude (C{degrees}).

763 @arg lon2: End longitude (C{degrees}).

764 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})

765 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or

766 L{a_f2Tuple}) or C{None}.

767 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll

768 B{C{lat2}} and B{C{lon2}} (C{bool}).

769

770 @return: Surface area, I{signed} (I{square} C{meter} or the same units as

771 B{C{radius}} I{squared}) or the I{spherical excess} (C{radians})

772 if C{B{radius}=0} or C{None}.

773

774 @raise TypeError: Invalid B{C{radius}}.

775

776 @raise UnitError: Invalid B{C{lat2}} or B{C{lat1}}.

777

778 @see: Function L{excessQuad_} and L{excessKarney}.

779 '''

780 return _eA(excessQuad_, radius, wrap, lat1, lon1, lat2, lon2)

781

782

783 def excessQuad_(phi2, phi1, lam21):

784 '''Compute the I{spherical excess} C{E} of a (spherical) quadrilateral bounded

785 by a segment of a great circle, two meridians and the equator.

786

787 @arg phi2: End latitude (C{radians}).

788 @arg phi1: Start latitude (C{radians}).

789 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).

790

791 @return: Spherical excess, I{signed} (C{radians}).

792

793 @see: Function L{excessQuad} and U{Spherical trigonometry

794 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.

795 '''

796 s = sin((phi2 + phi1) * _0_5)

797 c = cos((phi2 - phi1) * _0_5)

798 return Radians(Quad=atan2(tan_2(lam21) * s, c) * _2_0)

799

800

801 def flatLocal(lat1, lon1, lat2, lon2, datum=_WGS84, scaled=True, wrap=False):

802 '''Compute the distance between two (ellipsoidal) points using

803 the U{ellipsoidal Earth to plane projection<https://WikiPedia.org/

804 wiki/Geographical_distance#Ellipsoidal_Earth_projected_to_a_plane>}

805 aka U{Hubeny<https://www.OVG.AT/de/vgi/files/pdf/3781/>} formula.

806

807 @arg lat1: Start latitude (C{degrees}).

808 @arg lon1: Start longitude (C{degrees}).

809 @arg lat2: End latitude (C{degrees}).

810 @arg lon2: End longitude (C{degrees}).

811 @kwarg datum: Datum (L{Datum}) or ellipsoid (L{Ellipsoid},

812 L{Ellipsoid2} or L{a_f2Tuple}) to use.

813 @kwarg scaled: Scale prime_vertical by C{cos(B{phi})} (C{bool}),

814 see method L{pygeodesy.Ellipsoid.roc2_}.

815 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll

816 B{C{lat2}} and B{C{lon2}} (C{bool}).

817

818 @return: Distance (C{meter}, same units as the B{C{datum}}'s

819 ellipsoid axes).

820

821 @raise TypeError: Invalid B{C{datum}}.

822

823 @note: The meridional and prime_vertical radii of curvature

824 are taken and scaled at the mean of both latitude.

825

826 @see: Functions L{flatLocal_} or L{hubeny_}, L{cosineLaw}, L{flatPolar},

827 L{cosineAndoyerLambert}, L{cosineForsytheAndoyerLambert},

828 L{equirectangular}, L{euclidean}, L{haversine}, L{thomas},

829 L{vincentys}, method L{Ellipsoid.distance2} and U{local, flat

830 earth approximation<https://www.EdWilliams.org/avform.htm#flat>}.

831 '''

832 E = _ellipsoidal(datum, flatLocal)

833 return E._hubeny_2(*_d3(wrap, lat1, lon1, lat2, lon2),

834 scaled=scaled, squared=False) * E.a

835

836 hubeny = flatLocal # PYCHOK for Karl Hubeny

837

838

839 def flatLocal_(phi2, phi1, lam21, datum=_WGS84, scaled=True):

840 '''Compute the I{angular} distance between two (ellipsoidal) points using

841 the U{ellipsoidal Earth to plane projection<https://WikiPedia.org/

842 wiki/Geographical_distance#Ellipsoidal_Earth_projected_to_a_plane>}

843 aka U{Hubeny<https://www.OVG.AT/de/vgi/files/pdf/3781/>} formula.

844

845 @arg phi2: End latitude (C{radians}).

846 @arg phi1: Start latitude (C{radians}).

847 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).

848 @kwarg datum: Datum (L{Datum}) or ellipsoid (L{Ellipsoid},

849 L{Ellipsoid2} or L{a_f2Tuple}) to use.

850 @kwarg scaled: Scale prime_vertical by C{cos(B{phi})} (C{bool}),

851 see method L{pygeodesy.Ellipsoid.roc2_}.

852

853 @return: Angular distance (C{radians}).

854

855 @raise TypeError: Invalid B{C{datum}}.

856

857 @note: The meridional and prime_vertical radii of curvature

858 are taken and scaled I{at the mean of both latitude}.

859

860 @see: Functions L{flatLocal} or L{hubeny}, L{cosineAndoyerLambert_},

861 L{cosineForsytheAndoyerLambert_}, L{cosineLaw_}, L{flatPolar_},

862 L{euclidean_}, L{haversine_}, L{thomas_} and L{vincentys_} and

863 U{local, flat earth approximation

864 <https://www.EdWilliams.org/avform.htm#flat>}.

865 '''

866 E = _ellipsoidal(datum, flatLocal_)

867 return E._hubeny_2(phi2, phi1, lam21, scaled=scaled, squared=False)

868

869 hubeny_ = flatLocal_ # PYCHOK for Karl Hubeny

870

871

872 def flatPolar(lat1, lon1, lat2, lon2, radius=R_M, wrap=False):

873 '''Compute the distance between two (spherical) points using

874 the U{polar coordinate flat-Earth <https://WikiPedia.org/wiki/

875 Geographical_distance#Polar_coordinate_flat-Earth_formula>}

876 formula.

877

878 @arg lat1: Start latitude (C{degrees}).

879 @arg lon1: Start longitude (C{degrees}).

880 @arg lat2: End latitude (C{degrees}).

881 @arg lon2: End longitude (C{degrees}).

882 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})

883 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or

884 L{a_f2Tuple}) to use.

885 @kwarg wrap: If C{True}, wrap or I{normalize} and B{C{lat2}}

886 and B{C{lon2}} (C{bool}).

887

888 @return: Distance (C{meter}, same units as B{C{radius}} or the

889 ellipsoid or datum axes).

890

891 @raise TypeError: Invalid B{C{radius}}.

892

893 @see: Functions L{flatPolar_}, L{cosineAndoyerLambert},

894 L{cosineForsytheAndoyerLambert},L{cosineLaw},

895 L{flatLocal}/L{hubeny}, L{equirectangular},

896 L{euclidean}, L{haversine}, L{thomas} and

897 L{vincentys}.

898 '''

899 return _dS(flatPolar_, radius, wrap, lat1, lon1, lat2, lon2)

900

901

902 def flatPolar_(phi2, phi1, lam21):

903 '''Compute the I{angular} distance between two (spherical) points

904 using the U{polar coordinate flat-Earth<https://WikiPedia.org/wiki/

905 Geographical_distance#Polar_coordinate_flat-Earth_formula>}

906 formula.

907

908 @arg phi2: End latitude (C{radians}).

909 @arg phi1: Start latitude (C{radians}).

910 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).

911

912 @return: Angular distance (C{radians}).

913

914 @see: Functions L{flatPolar}, L{cosineAndoyerLambert_},

915 L{cosineForsytheAndoyerLambert_}, L{cosineLaw_},

916 L{euclidean_}, L{flatLocal_}/L{hubeny_}, L{haversine_},

917 L{thomas_} and L{vincentys_}.

918 '''

919 a = fabs(PI_2 - phi1) # co-latitude

920 b = fabs(PI_2 - phi2) # co-latitude

921 if a < b:

922 a, b = b, a

923 if a < EPS0:

924 a = _0_0

925 elif b > 0:

926 b = b / a # /= chokes PyChecker

927 c = b * cos(lam21) * _2_0

928 c = fsumf_(_1_0, b**2, -fabs(c))

929 a *= sqrt0(c)

930 return a

931

932

933 def _hartzell(pov, los, earth, **kwds):

934 '''(INTERNAL) Helper for C{CartesianBase.hartzell} and C{LatLonBase.hartzell}.

935 '''

936 if earth is None:

937 earth = pov.datum

938 else:

939 earth = _spherical_datum(earth, name__=hartzell)

940 pov = pov.toDatum(earth)

941 h = pov.height

942 if h < 0: # EPS0

943 t = _SPACE_(Fmt.PARENSPACED(height=h), _inside_)

944 raise IntersectionError(pov=pov, earth=earth, txt=t)

945 return hartzell(pov, los=los, earth=earth, **kwds) if h > 0 else pov # EPS0

946

947

948 def hartzell(pov, los=False, earth=_WGS84, **name_LatLon_and_kwds):

949 '''Compute the intersection of the earth's surface and a Line-Of-Sight from

950 a Point-Of-View in space.

951

952 @arg pov: Point-Of-View outside the earth (C{LatLon}, C{Cartesian},

953 L{Ecef9Tuple} or L{Vector3d}).

954 @kwarg los: Line-Of-Sight, I{direction} to earth (L{Los}, L{Vector3d}),

955 C{True} for the I{normal, plumb} onto the surface or C{False}

956 or C{None} to point to the center of the earth.

957 @kwarg earth: The earth model (L{Datum}, L{Ellipsoid}, L{Ellipsoid2},

958 L{a_f2Tuple} or a C{scalar} earth radius in C{meter}).

959 @kwarg name_LatLon_and_kwds: Optional, overriding C{B{name}="hartzell"}

960 (C{str}), class C{B{LatLon}=None} to return the intersection

961 plus additional C{LatLon} keyword arguments, include the

962 B{C{datum}} if different and to convert from B{C{earth}}.

963

964 @return: The intersection (L{Vector3d}, B{C{pov}}'s C{cartesian type} or the

965 given B{C{LatLon}} instance) with attribute C{height} set to the

966 distance to the B{C{pov}}.

967

968 @raise IntersectionError: Invalid B{C{pov}} or B{C{pov}} inside the earth or

969 invalid B{C{los}} or B{C{los}} points outside or

970 away from the earth.

971

972 @raise TypeError: Invalid B{C{earth}}, C{ellipsoid} or C{datum}.

973

974 @see: Class L{Los}, functions L{tyr3d} and L{hartzell4} and methods

975 L{Ellipsoid.hartzell4} and any C{Cartesian.hartzell} and C{LatLon.hartzell}.

976 '''

977 n, LatLon_and_kwds = _name2__(name_LatLon_and_kwds, name__=hartzell)

978 try:

979 D = _spherical_datum(earth, name__=hartzell)

980 r, h, i = _MODS.triaxials._hartzell3(pov, los, D.ellipsoid._triaxial)

981

982 C = _MODS.cartesianBase.CartesianBase

983 if LatLon_and_kwds:

984 c = C(r, datum=D)

985 r = c.toLatLon(**_xkwds(LatLon_and_kwds, height=h))

986 elif isinstance(r, C):

987 r.height = h

988 if i:

989 r._iteration = i

990 except Exception as x:

991 raise IntersectionError(pov=pov, los=los, earth=earth, cause=x,

992 **LatLon_and_kwds)

993 return _xnamed(r, n) if n else r

994

995

996 def haversine(lat1, lon1, lat2, lon2, radius=R_M, wrap=False):

997 '''Compute the distance between two (spherical) points using the

998 U{Haversine<https://www.Movable-Type.co.UK/scripts/latlong.html>}

999 formula.

1000

1001 @arg lat1: Start latitude (C{degrees}).

1002 @arg lon1: Start longitude (C{degrees}).

1003 @arg lat2: End latitude (C{degrees}).

1004 @arg lon2: End longitude (C{degrees}).

1005 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})

1006 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or

1007 L{a_f2Tuple}) to use.

1008 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll

1009 B{C{lat2}} and B{C{lon2}} (C{bool}).

1010

1011 @return: Distance (C{meter}, same units as B{C{radius}}).

1012

1013 @raise TypeError: Invalid B{C{radius}}.

1014

1015 @see: U{Distance between two (spherical) points

1016 <https://www.EdWilliams.org/avform.htm#Dist>}, functions

1017 L{cosineLaw}, L{cosineAndoyerLambert}, L{cosineForsytheAndoyerLambert},

1018 L{equirectangular}, L{euclidean}, L{flatLocal}/L{hubeny}, L{flatPolar},

1019 L{thomas} and L{vincentys} and methods L{Ellipsoid.distance2},

1020 C{LatLon.distanceTo*} and C{LatLon.equirectangularTo}.

1021

1022 @note: See note at function L{vincentys_}.

1023 '''

1024 return _dS(haversine_, radius, wrap, lat1, lon1, lat2, lon2)

1025

1026

1027 def haversine_(phi2, phi1, lam21):

1028 '''Compute the I{angular} distance between two (spherical) points

1029 using the U{Haversine<https://www.Movable-Type.co.UK/scripts/latlong.html>}

1030 formula.

1031

1032 @arg phi2: End latitude (C{radians}).

1033 @arg phi1: Start latitude (C{radians}).

1034 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).

1035

1036 @return: Angular distance (C{radians}).

1037

1038 @see: Functions L{haversine}, L{cosineAndoyerLambert_},

1039 L{cosineForsytheAndoyerLambert_}, L{cosineLaw_},

1040 L{euclidean_}, L{flatLocal_}/L{hubeny_}, L{flatPolar_},

1041 L{thomas_} and L{vincentys_}.

1042

1043 @note: See note at function L{vincentys_}.

1044 '''

1045 def _hsin(rad):

1046 return sin(rad * _0_5)**2

1047

1048 h = _hsin(phi2 - phi1) + cos(phi1) * cos(phi2) * _hsin(lam21) # haversine

1049 return atan2(sqrt0(h), sqrt0(_1_0 - h)) * _2_0 # == asin(sqrt(h)) * 2

1050

1051

1052 def heightOf(angle, distance, radius=R_M):

1053 '''Determine the height above the (spherical) earth' surface after

1054 traveling along a straight line at a given tilt.

1055

1056 @arg angle: Tilt angle above horizontal (C{degrees}).

1057 @arg distance: Distance along the line (C{meter} or same units as

1058 B{C{radius}}).

1059 @kwarg radius: Optional mean earth radius (C{meter}).

1060

1061 @return: Height (C{meter}, same units as B{C{distance}} and B{C{radius}}).

1062

1063 @raise ValueError: Invalid B{C{angle}}, B{C{distance}} or B{C{radius}}.

1064

1065 @see: U{MultiDop geog_lib.GeogBeamHt<https://GitHub.com/NASA/MultiDop>}

1066 (U{Shapiro et al. 2009, JTECH

1067 <https://Journals.AMetSoc.org/doi/abs/10.1175/2009JTECHA1256.1>}

1068 and U{Potvin et al. 2012, JTECH

1069 <https://Journals.AMetSoc.org/doi/abs/10.1175/JTECH-D-11-00019.1>}).

1070 '''

1071 r = h = Radius(radius)

1072 d = fabs(Distance(distance))

1073 if d > h:

1074 d, h = h, d

1075

1076 if d > EPS0: # and h > EPS0

1077 d = d / h # /= h chokes PyChecker

1078 s = sin(Phid(angle=angle, clip=_180_0))

1079 s = fsumf_(_1_0, s * d * _2_0, d**2)

1080 if s > 0:

1081 return h * sqrt(s) - r

1082

1083 raise _ValueError(angle=angle, distance=distance, radius=radius)

1084

1085

1086 def heightOrthometric(h_ll, N):

1087 '''Get the I{orthometric} height B{H}, the height above the geoid, earth surface.

1088

1089 @arg h_ll: The height above the ellipsoid (C{meter}) or an I{ellipsoidal}

1090 location (C{LatLon} with a C{height} or C{h} attribute).

1091 @arg N: The I{geoid} height (C{meter}), the height of the geoid above the

1092 ellipsoid at the same B{C{h_ll}} location.

1093

1094 @return: I{Orthometric} height C{B{H} = B{h} - B{N}} (C{meter}, same units

1095 as B{C{h}} and B{C{N}}).

1096

1097 @see: U{Ellipsoid, Geoid, and Othometric Heights<https://www.NGS.NOAA.gov/

1098 GEOID/PRESENTATIONS/2007_02_24_CCPS/Roman_A_PLSC2007notes.pdf>}, page

1099 6 and module L{pygeodesy.geoids}.

1100 '''

1101 h = h_ll if _isHeight(h_ll) else _xattr(h_ll, height=_xattr(h_ll, h=0))

1102 return Height(H=Height(h=h) - Height(N=N))

1103

1104

1105 def horizon(height, radius=R_M, refraction=False):

1106 '''Determine the distance to the horizon from a given altitude above the

1107 (spherical) earth.

1108

1109 @arg height: Altitude (C{meter} or same units as B{C{radius}}).

1110 @kwarg radius: Optional mean earth radius (C{meter}).

1111 @kwarg refraction: Consider atmospheric refraction (C{bool}).

1112

1113 @return: Distance (C{meter}, same units as B{C{height}} and B{C{radius}}).

1114

1115 @raise ValueError: Invalid B{C{height}} or B{C{radius}}.

1116

1117 @see: U{Distance to horizon<https://www.EdWilliams.org/avform.htm#Horizon>}.

1118 '''

1119 h, r = Height(height), Radius(radius)

1120 if min(h, r) < 0:

1121 raise _ValueError(height=height, radius=radius)

1122

1123 d2 = ((r * 2.415750694528) if refraction else # 2.0 / 0.8279

1124 fsumf_(r, r, h)) * h

1125 return sqrt0(d2)

1126

1127

1128 class _idllmn6(object): # see also .geodesicw._wargs, .latlonBase._toCartesian3, .vector2d._numpy

1129 '''(INTERNAL) Helper for C{intersection2} and C{intersections2}.

1130 '''

1131 @contextmanager # <https://www.Python.org/dev/peps/pep-0343/> Examples

1132 def __call__(self, datum, lat1, lon1, lat2, lon2, small, wrap, s, **kwds):

1133 try:

1134 if wrap:

1135 _, lat2, lon2 = _Wrap.latlon3(lon1, lat2, lon2, wrap)

1136 kwds = _xkwds(kwds, wrap=wrap) # for _xError

1137 m = small if small is _100km else Meter_(small=small)

1138 n = _dunder_nameof(intersections2 if s else intersection2)

1139 if datum is None or euclidean(lat1, lon1, lat2, lon2) < m:

1140 d, m = None, _MODS.vector3d

1141 _i = m._intersects2 if s else m._intersect3d3

1142 elif _isRadius(datum) and datum < 0 and not s:

1143 d = _spherical_datum(-datum, name=n)

1144 m = _MODS.sphericalNvector

1145 _i = m.intersection

1146 else:

1147 d = _spherical_datum(datum, name=n)

1148 if d.isSpherical:

1149 m = _MODS.sphericalTrigonometry

1150 _i = m._intersects2 if s else m._intersect

1151 elif d.isEllipsoidal:

1152 try:

1153 if d.ellipsoid.geodesic:

1154 pass

1155 m = _MODS.ellipsoidalKarney

1156 except ImportError:

1157 m = _MODS.ellipsoidalExact

1158 _i = m._intersections2 if s else m._intersection3 # ellispoidalBaseDI

1159 else:

1160 raise _TypeError(datum=datum)

1161 yield _i, d, lat2, lon2, m, n

1162

1163 except (TypeError, ValueError) as x:

1164 raise _xError(x, lat1=lat1, lon1=lon1, datum=datum,

1165 lat2=lat2, lon2=lon2, small=small, **kwds)

1166

1167_idllmn6 = _idllmn6() # PYCHOK singleton

1168

1169

1170 def intersection2(lat1, lon1, bearing1,

1171 lat2, lon2, bearing2, datum=None, wrap=False, small=_100km): # was=True

1172 '''I{Conveniently} compute the intersection of two lines each defined

1173 by a (geodetic) point and a bearing from North, using either ...

1174

1175 1) L{vector3d.intersection3d3} for B{C{small}} distances (below 100 Km

1176 or about 0.88 degrees) or if I{no} B{C{datum}} is specified, or ...

1177

1178 2) L{sphericalTrigonometry.intersection} for a spherical B{C{datum}}

1179 or a C{scalar B{datum}} representing the earth radius, conventionally

1180 in C{meter} or ...

1181

1182 3) L{sphericalNvector.intersection} if B{C{datum}} is a I{negative}

1183 C{scalar}, (negative) earth radius, conventionally in C{meter} or ...

1184

1185 4) L{ellipsoidalKarney.intersection3} for an ellipsoidal B{C{datum}}

1186 and if I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>}

1187 is installed, otherwise ...

1188

1189 5) L{ellipsoidalExact.intersection3}, provided B{C{datum}} is ellipsoidal.

1190

1191 @arg lat1: Latitude of the first point (C{degrees}).

1192 @arg lon1: Longitude of the first point (C{degrees}).

1193 @arg bearing1: Bearing at the first point (compass C{degrees}).

1194 @arg lat2: Latitude of the second point (C{degrees}).

1195 @arg lon2: Longitude of the second point (C{degrees}).

1196 @arg bearing2: Bearing at the second point (compass C{degrees}).

1197 @kwarg datum: Optional datum (L{Datum}) or ellipsoid (L{Ellipsoid},

1198 L{Ellipsoid2} or L{a_f2Tuple}) or C{scalar} earth

1199 radius (C{meter}, same units as B{C{radius1}} and

1200 B{C{radius2}}) or C{None}.

1201 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}}

1202 and B{C{lon2}} (C{bool}).

1203 @kwarg small: Upper limit for small distances (C{meter}).

1204

1205 @return: A L{LatLon2Tuple}C{(lat, lon)} with the lat- and

1206 longitude of the intersection point.

1207

1208 @raise IntersectionError: Ambiguous or infinite intersection

1209 or colinear, parallel or otherwise

1210 non-intersecting lines.

1211

1212 @raise TypeError: Invalid B{C{datum}}.

1213

1214 @raise UnitError: Invalid B{C{lat1}}, B{C{lon1}}, B{C{bearing1}},

1215 B{C{lat2}}, B{C{lon2}} or B{C{bearing2}}.

1216

1217 @see: Method L{RhumbLine.intersection2}.

1218

1219 @note: The returned intersections may be near-antipodal.

1220 '''

1221 b1 = Bearing(bearing1=bearing1)

1222 b2 = Bearing(bearing2=bearing2)

1223 with _idllmn6(datum, lat1, lon1, lat2, lon2,

1224 small, wrap, False, bearing1=b1, bearing2=b2) as t:

1225 _i, d, lat2, lon2, m, n = t

1226 if d is None:

1227 t, _, _ = _i(m.Vector3d(lon1, lat1, 0), b1,

1228 m.Vector3d(lon2, lat2, 0), b2, useZ=False)

1229 t = LatLon2Tuple(t.y, t.x, name=n)

1230

1231 else:

1232 t = _i(m.LatLon(lat1, lon1, datum=d), b1,

1233 m.LatLon(lat2, lon2, datum=d), b2, height=0, wrap=False)

1234 if isinstance(t, Intersection3Tuple): # ellipsoidal

1235 t, _, _ = t

1236 t = LatLon2Tuple(t.lat, t.lon, name=n)

1237 return t

1238

1239

1240 def intersections2(lat1, lon1, radius1,

1241 lat2, lon2, radius2, datum=None, wrap=False, small=_100km): # was=True

1242 '''I{Conveniently} compute the intersections of two circles each defined

1243 by a (geodetic) center point and a radius, using either ...

1244

1245 1) L{vector3d.intersections2} for B{C{small}} distances (below 100 Km

1246 or about 0.88 degrees) or if I{no} B{C{datum}} is specified, or ...

1247

1248 2) L{sphericalTrigonometry.intersections2} for a spherical B{C{datum}}

1249 or a C{scalar B{datum}} representing the earth radius, conventionally

1250 in C{meter} or ...

1251

1252 3) L{ellipsoidalKarney.intersections2} for an ellipsoidal B{C{datum}}

1253 and if I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>}

1254 is installed, otherwise ...

1255

1256 4) L{ellipsoidalExact.intersections2}, provided B{C{datum}} is ellipsoidal.

1257

1258 @arg lat1: Latitude of the first circle center (C{degrees}).

1259 @arg lon1: Longitude of the first circle center (C{degrees}).

1260 @arg radius1: Radius of the first circle (C{meter}, conventionally).

1261 @arg lat2: Latitude of the second circle center (C{degrees}).

1262 @arg lon2: Longitude of the second circle center (C{degrees}).

1263 @arg radius2: Radius of the second circle (C{meter}, same units as B{C{radius1}}).

1264 @kwarg datum: Optional datum (L{Datum}) or ellipsoid (L{Ellipsoid},

1265 L{Ellipsoid2} or L{a_f2Tuple}) or C{scalar} earth

1266 radius (C{meter}, same units as B{C{radius1}} and

1267 B{C{radius2}}) or C{None}.

1268 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}}

1269 and B{C{lon2}} (C{bool}).

1270 @kwarg small: Upper limit for small distances (C{meter}).

1271

1272 @return: 2-Tuple of the intersection points, each a

1273 L{LatLon2Tuple}C{(lat, lon)}. For abutting circles, the

1274 points are the same instance, aka the I{radical center}.

1275

1276 @raise IntersectionError: Concentric, antipodal, invalid or

1277 non-intersecting circles or no

1278 convergence.

1279

1280 @raise TypeError: Invalid B{C{datum}}.

1281

1282 @raise UnitError: Invalid B{C{lat1}}, B{C{lon1}}, B{C{radius1}},

1283 B{C{lat2}}, B{C{lon2}} or B{C{radius2}}.

1284 '''

1285 r1 = Radius_(radius1=radius1)

1286 r2 = Radius_(radius2=radius2)

1287 with _idllmn6(datum, lat1, lon1, lat2, lon2,

1288 small, wrap, True, radius1=r1, radius2=r2) as t:

1289 _i, d, lat2, lon2, m, n = t

1290 if d is None:

1291 r1 = m2degrees(r1, radius=R_M, lat=lat1)

1292 r2 = m2degrees(r2, radius=R_M, lat=lat2)

1293

1294 def _V2T(x, y, _, **unused): # _ == z unused

1295 return LatLon2Tuple(y, x, name=n)

1296

1297 t = _i(m.Vector3d(lon1, lat1, 0), r1,

1298 m.Vector3d(lon2, lat2, 0), r2, sphere=False,

1299 Vector=_V2T)

1300 else:

1301 def _LL2T(lat, lon, **unused):

1302 return LatLon2Tuple(lat, lon, name=n)

1303

1304 t = _i(m.LatLon(lat1, lon1, datum=d), r1,

1305 m.LatLon(lat2, lon2, datum=d), r2,

1306 LatLon=_LL2T, height=0, wrap=False)

1307 return t

1308

1309

1310 def isantipode(lat1, lon1, lat2, lon2, eps=EPS):

1311 '''Check whether two points are I{antipodal}, on diametrically

1312 opposite sides of the earth.

1313

1314 @arg lat1: Latitude of one point (C{degrees}).

1315 @arg lon1: Longitude of one point (C{degrees}).

1316 @arg lat2: Latitude of the other point (C{degrees}).

1317 @arg lon2: Longitude of the other point (C{degrees}).

1318 @kwarg eps: Tolerance for near-equality (C{degrees}).

1319

1320 @return: C{True} if points are antipodal within the

1321 B{C{eps}} tolerance, C{False} otherwise.

1322

1323 @see: Functions L{isantipode_} and L{antipode}.

1324 '''

1325 return (fabs(lat1 + lat2) <= eps and

1326 fabs(lon1 + lon2) <= eps) or _isequalTo(

1327 normal(lat1, lon1), antipode(lat2, lon2), eps)

1328

1329

1330 def isantipode_(phi1, lam1, phi2, lam2, eps=EPS):

1331 '''Check whether two points are I{antipodal}, on diametrically

1332 opposite sides of the earth.

1333

1334 @arg phi1: Latitude of one point (C{radians}).

1335 @arg lam1: Longitude of one point (C{radians}).

1336 @arg phi2: Latitude of the other point (C{radians}).

1337 @arg lam2: Longitude of the other point (C{radians}).

1338 @kwarg eps: Tolerance for near-equality (C{radians}).

1339

1340 @return: C{True} if points are antipodal within the

1341 B{C{eps}} tolerance, C{False} otherwise.

1342

1343 @see: Functions L{isantipode} and L{antipode_}.

1344 '''

1345 return (fabs(phi1 + phi2) <= eps and

1346 fabs(lam1 + lam2) <= eps) or _isequalTo_(

1347 normal_(phi1, lam1), antipode_(phi2, lam2), eps)

1348

1349

1350 def _isequalTo(p1, p2, eps=EPS):

1351 '''Compare 2 point lat-/lons ignoring C{class}.

1352 '''

1353 return (fabs(p1.lat - p2.lat) <= eps and

1354 fabs(p1.lon - p2.lon) <= eps) if eps else (p1.latlon == p2.latlon)

1355

1356

1357 def _isequalTo_(p1, p2, eps=EPS):

1358 '''(INTERNAL) Compare 2 point phi-/lams ignoring C{class}.

1359 '''

1360 return (fabs(p1.phi - p2.phi) <= eps and

1361 fabs(p1.lam - p2.lam) <= eps) if eps else (p1.philam == p2.philam)

1362

1363

1364 def isnormal(lat, lon, eps=0):

1365 '''Check whether B{C{lat}} I{and} B{C{lon}} are within their

1366 respective I{normal} range in C{degrees}.

1367

1368 @arg lat: Latitude (C{degrees}).

1369 @arg lon: Longitude (C{degrees}).

1370 @kwarg eps: Optional tolerance C{degrees}).

1371

1372 @return: C{True} if C{(abs(B{lat}) + B{eps}) <= 90} and

1373 C{(abs(B{lon}) + B{eps}) <= 180}, C{False} othwerwise.

1374

1375 @see: Functions L{isnormal_} and L{normal}.

1376 '''

1377 return (_90_0 - fabs(lat)) >= eps and _loneg(fabs(lon)) >= eps

1378

1379

1380 def isnormal_(phi, lam, eps=0):

1381 '''Check whether B{C{phi}} I{and} B{C{lam}} are within their

1382 respective I{normal} range in C{radians}.

1383

1384 @arg phi: Latitude (C{radians}).

1385 @arg lam: Longitude (C{radians}).

1386 @kwarg eps: Optional tolerance C{radians}).

1387

1388 @return: C{True} if C{(abs(B{phi}) + B{eps}) <= PI/2} and

1389 C{(abs(B{lam}) + B{eps}) <= PI}, C{False} othwerwise.

1390

1391 @see: Functions L{isnormal} and L{normal_}.

1392 '''

1393 return (PI_2 - fabs(phi)) >= eps and (PI - fabs(lam)) >= eps

1394

1395

1396 def latlon2n_xyz(lat, lon, **name):

1397 '''Convert lat-, longitude to C{n-vector} (I{normal} to the

1398 earth's surface) X, Y and Z components.

1399

1400 @arg lat: Latitude (C{degrees}).

1401 @arg lon: Longitude (C{degrees}).

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

1403

1404 @return: A L{Vector3Tuple}C{(x, y, z)}.

1405

1406 @see: Function L{philam2n_xyz}.

1407

1408 @note: These are C{n-vector} x, y and z components,

1409 I{NOT} geocentric ECEF x, y and z coordinates!

1410 '''

1411 return _2n_xyz(name, *sincos2d_(lat, lon))

1412

1413

1414 def _normal2(a, b, n_2, n, n2):

1415 '''(INTERNAL) Helper for C{normal} and C{normal_}.

1416 '''

1417 if fabs(b) > n:

1418 b = remainder(b, n2)

1419 if fabs(a) > n_2:

1420 r = remainder(a, n)

1421 if r != a:

1422 a = -r

1423 b -= n if b > 0 else -n

1424 return float0_(a, b)

1425

1426

1427 def normal(lat, lon, **name):

1428 '''Normalize a lat- I{and} longitude pair in C{degrees}.

1429

1430 @arg lat: Latitude (C{degrees}).

1431 @arg lon: Longitude (C{degrees}).

1432 @kwarg name: Optional, overriding C{B{name}="normal"} (C{str}).

1433

1434 @return: L{LatLon2Tuple}C{(lat, lon)} with C{abs(lat) <= 90}

1435 and C{abs(lon) <= 180}.

1436

1437 @see: Functions L{normal_} and L{isnormal}.

1438 '''

1439 return LatLon2Tuple(*_normal2(lat, lon, _90_0, _180_0, _360_0),

1440 name=_name__(name, name__=normal))

1441

1442

1443 def normal_(phi, lam, **name):

1444 '''Normalize a lat- I{and} longitude pair in C{radians}.

1445

1446 @arg phi: Latitude (C{radians}).

1447 @arg lam: Longitude (C{radians}).

1448 @kwarg name: Optional, overriding C{B{name}="normal_"} (C{str}).

1449

1450 @return: L{PhiLam2Tuple}C{(phi, lam)} with C{abs(phi) <= PI/2}

1451 and C{abs(lam) <= PI}.

1452

1453 @see: Functions L{normal} and L{isnormal_}.

1454 '''

1455 return PhiLam2Tuple(*_normal2(phi, lam, PI_2, PI, PI2),

1456 name=_name__(name, name__=normal_))

1457

1458

1459 def _2n_xyz(name, sa, ca, sb, cb): # name always **name

1460 '''(INTERNAL) Helper for C{latlon2n_xyz} and C{philam2n_xyz}.

1461 '''

1462 # Kenneth Gade eqn 3, but using right-handed

1463 # vector x -> 0°E,0°N, y -> 90°E,0°N, z -> 90°N

1464 return Vector3Tuple(ca * cb, ca * sb, sa, **name)

1465

1466

1467 def n_xyz2latlon(x, y, z, **name):

1468 '''Convert C{n-vector} components to lat- and longitude in C{degrees}.

1469

1470 @arg x: X component (C{scalar}).

1471 @arg y: Y component (C{scalar}).

1472 @arg z: Z component (C{scalar}).

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

1474

1475 @return: A L{LatLon2Tuple}C{(lat, lon)}.

1476

1477 @see: Function L{n_xyz2philam}.

1478 '''

1479 return LatLon2Tuple(atan2d(z, hypot(x, y)), atan2d(y, x), **name)

1480

1481

1482 def n_xyz2philam(x, y, z, **name):

1483 '''Convert C{n-vector} components to lat- and longitude in C{radians}.

1484

1485 @arg x: X component (C{scalar}).

1486 @arg y: Y component (C{scalar}).

1487 @arg z: Z component (C{scalar}).

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

1489

1490 @return: A L{PhiLam2Tuple}C{(phi, lam)}.

1491

1492 @see: Function L{n_xyz2latlon}.

1493 '''

1494 return PhiLam2Tuple(atan2(z, hypot(x, y)), atan2(y, x), **name)

1495

1496

1497 def _opposes(d, m, n, n2):

1498 '''(INTERNAL) Helper for C{opposing} and C{opposing_}.

1499 '''

1500 d = d % n2 # -20 % 360 == 340, -1 % PI2 == PI2 - 1

1501 return False if d < m or d > (n2 - m) else (

1502 True if (n - m) < d < (n + m) else None)

1503

1504

1505 def opposing(bearing1, bearing2, margin=_90_0):

1506 '''Compare the direction of two bearings given in C{degrees}.

1507

1508 @arg bearing1: First bearing (compass C{degrees}).

1509 @arg bearing2: Second bearing (compass C{degrees}).

1510 @kwarg margin: Optional, interior angle bracket (C{degrees}).

1511

1512 @return: C{True} if both bearings point in opposite, C{False} if

1513 in similar or C{None} if in I{perpendicular} directions.

1514

1515 @see: Function L{opposing_}.

1516 '''

1517 m = Degrees_(margin=margin, low=EPS0, high=_90_0)

1518 return _opposes(bearing2 - bearing1, m, _180_0, _360_0)

1519

1520

1521 def opposing_(radians1, radians2, margin=PI_2):

1522 '''Compare the direction of two bearings given in C{radians}.

1523

1524 @arg radians1: First bearing (C{radians}).

1525 @arg radians2: Second bearing (C{radians}).

1526 @kwarg margin: Optional, interior angle bracket (C{radians}).

1527

1528 @return: C{True} if both bearings point in opposite, C{False} if

1529 in similar or C{None} if in perpendicular directions.

1530

1531 @see: Function L{opposing}.

1532 '''

1533 m = Radians_(margin=margin, low=EPS0, high=PI_2)

1534 return _opposes(radians2 - radians1, m, PI, PI2)

1535

1536

1537 def philam2n_xyz(phi, lam, **name):

1538 '''Convert lat-, longitude to C{n-vector} (I{normal} to the

1539 earth's surface) X, Y and Z components.

1540

1541 @arg phi: Latitude (C{radians}).

1542 @arg lam: Longitude (C{radians}).

1543 @kwarg name: Optional name (C{str}).

1544

1545 @return: A L{Vector3Tuple}C{(x, y, z)}.

1546

1547 @see: Function L{latlon2n_xyz}.

1548

1549 @note: These are C{n-vector} x, y and z components,

1550 I{NOT} geocentric ECEF x, y and z coordinates!

1551 '''

1552 return _2n_xyz(name, *sincos2_(phi, lam))

1553

1554

1555 def _Propy(func, nargs, kwds):

1556 '''(INTERNAL) Helper for the C{frechet.[-]Frechet**} and

1557 C{hausdorff.[-]Hausdorff*} classes.

1558 '''

1559 try:

1560 _xcallable(distance=func)

1561 # assert _args_kwds_count2(func)[0] == nargs + int(ismethod(func))

1562 _ = func(*_0_0s(nargs), **kwds)

1563 except Exception as x:

1564 t = unstr(func, **kwds)

1565 raise _TypeError(t, cause=x)

1566 return func

1567

1568

1569 def _radical2(d, r1, r2, **name): # in .ellipsoidalBaseDI, .sphericalTrigonometry, .vector3d

1570 # (INTERNAL) See C{radical2} below

1571 # assert d > EPS0

1572 r = fsumf_(_1_0, (r1 / d)**2, -(r2 / d)**2) * _0_5

1573 n = _name__(name, name__=radical2)

1574 return Radical2Tuple(max(_0_0, min(_1_0, r)), r * d, name=n)

1575

1576

1577 def radical2(distance, radius1, radius2, **name):

1578 '''Compute the I{radical ratio} and I{radical line} of two

1579 U{intersecting circles<https://MathWorld.Wolfram.com/

1580 Circle-CircleIntersection.html>}.

1581

1582 The I{radical line} is perpendicular to the axis thru the

1583 centers of the circles at C{(0, 0)} and C{(B{distance}, 0)}.

1584

1585 @arg distance: Distance between the circle centers (C{scalar}).

1586 @arg radius1: Radius of the first circle (C{scalar}).

1587 @arg radius2: Radius of the second circle (C{scalar}).

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

1589

1590 @return: A L{Radical2Tuple}C{(ratio, xline)} where C{0.0 <=

1591 ratio <= 1.0} and C{xline} is along the B{C{distance}}.

1592

1593 @raise IntersectionError: The B{C{distance}} exceeds the sum

1594 of B{C{radius1}} and B{C{radius2}}.

1595

1596 @raise UnitError: Invalid B{C{distance}}, B{C{radius1}} or

1597 B{C{radius2}}.

1598

1599 @see: U{Circle-Circle Intersection

1600 <https://MathWorld.Wolfram.com/Circle-CircleIntersection.html>}.

1601 '''

1602 d = Distance_(distance, low=_0_0)

1603 r1 = Radius_(radius1=radius1)

1604 r2 = Radius_(radius2=radius2)

1605 if d > (r1 + r2):

1606 raise IntersectionError(distance=d, radius1=r1, radius2=r2,

1607 txt=_too_(_distant_))

1608 return _radical2(d, r1, r2, **name) if d > EPS0 else \

1609 Radical2Tuple(_0_5, _0_0, **name)

1610

1611

1612 class Radical2Tuple(_NamedTuple):

1613 '''2-Tuple C{(ratio, xline)} of the I{radical} C{ratio} and

1614 I{radical} C{xline}, both C{scalar} and C{0.0 <= ratio <= 1.0}

1615 '''

1616 _Names_ = (_ratio_, _xline_)

1617 _Units_ = ( Scalar, Scalar)

1618

1619

1620 def _radistance(inst):

1621 '''(INTERNAL) Helper for the L{frechet._FrechetMeterRadians}

1622 and L{hausdorff._HausdorffMeterRedians} classes.

1623 '''

1624 wrap_, kwds_ = _xkwds_pop2(inst._kwds, wrap=False)

1625 func_ = inst._func_

1626 try: # calling lower-overhead C{func_}

1627 func_(0, _0_25, _0_5, **kwds_)

1628 wrap_ = _Wrap._philamop(wrap_)

1629 except TypeError:

1630 return inst.distance

1631

1632 def _philam(p):

1633 try:

1634 return p.phi, p.lam

1635 except AttributeError: # no .phi or .lam

1636 return radians(p.lat), radians(p.lon)

1637

1638 def _func_wrap(point1, point2):

1639 phi1, lam1 = wrap_(*_philam(point1))

1640 phi2, lam2 = wrap_(*_philam(point2))

1641 return func_(phi2, phi1, lam2 - lam1, **kwds_)

1642

1643 inst._units = inst._units_

1644 return _func_wrap

1645

1646

1647 def _scale_deg(lat1, lat2): # degrees

1648 # scale factor cos(mean of lats) for delta lon

1649 m = fabs(lat1 + lat2) * _0_5

1650 return cos(radians(m)) if m < 90 else _0_0

1651

1652

1653 def _scale_rad(phi1, phi2): # radians, by .frechet, .hausdorff, .heights

1654 # scale factor cos(mean of phis) for delta lam

1655 m = fabs(phi1 + phi2) * _0_5

1656 return cos(m) if m < PI_2 else _0_0

1657

1658

1659 def _sincosa6(phi2, phi1, lam21): # [4] in cosineLaw

1660 '''(INTERNAL) C{sin}es, C{cos}ines and C{acos}ine.

1661 '''

1662 s2, c2, s1, c1, _, c21 = sincos2_(phi2, phi1, lam21)

1663 return s2, c2, s1, c1, acos1(s1 * s2 + c1 * c2 * c21), c21

1664

1665

1666 def thomas(lat1, lon1, lat2, lon2, datum=_WGS84, wrap=False):

1667 '''Compute the distance between two (ellipsoidal) points using

1668 U{Thomas'<https://apps.DTIC.mil/dtic/tr/fulltext/u2/703541.pdf>}

1669 formula.

1670

1671 @arg lat1: Start latitude (C{degrees}).

1672 @arg lon1: Start longitude (C{degrees}).

1673 @arg lat2: End latitude (C{degrees}).

1674 @arg lon2: End longitude (C{degrees}).

1675 @kwarg datum: Datum (L{Datum}) or ellipsoid (L{Ellipsoid},

1676 L{Ellipsoid2} or L{a_f2Tuple}) to use.

1677 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll

1678 B{C{lat2}} and B{C{lon2}} (C{bool}).

1679

1680 @return: Distance (C{meter}, same units as the B{C{datum}}'s

1681 ellipsoid axes).

1682

1683 @raise TypeError: Invalid B{C{datum}}.

1684

1685 @see: Functions L{thomas_}, L{cosineAndoyerLambert}, L{cosineForsytheAndoyerLambert},

1686 L{cosineLaw}, L{equirectangular}, L{euclidean}, L{flatLocal}/L{hubeny},

1687 L{flatPolar}, L{haversine}, L{vincentys} and method L{Ellipsoid.distance2}.

1688 '''

1689 return _dE(thomas_, datum, wrap, lat1, lon1, lat2, lon2)

1690

1691

1692 def thomas_(phi2, phi1, lam21, datum=_WGS84):

1693 '''Compute the I{angular} distance between two (ellipsoidal) points using

1694 U{Thomas'<https://apps.DTIC.mil/dtic/tr/fulltext/u2/703541.pdf>}

1695 formula.

1696

1697 @arg phi2: End latitude (C{radians}).

1698 @arg phi1: Start latitude (C{radians}).

1699 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).

1700 @kwarg datum: Datum or ellipsoid to use (L{Datum}, L{Ellipsoid},

1701 L{Ellipsoid2} or L{a_f2Tuple}).

1702

1703 @return: Angular distance (C{radians}).

1704

1705 @raise TypeError: Invalid B{C{datum}}.

1706

1707 @see: Functions L{thomas}, L{cosineAndoyerLambert_},

1708 L{cosineForsytheAndoyerLambert_}, L{cosineLaw_},

1709 L{euclidean_}, L{flatLocal_}/L{hubeny_}, L{flatPolar_},

1710 L{haversine_} and L{vincentys_} and U{Geodesy-PHP

1711 <https://GitHub.com/jtejido/geodesy-php/blob/master/src/Geodesy/

1712 Distance/ThomasFormula.php>}.

1713 '''

1714 s2, c2, s1, c1, r, _ = _sincosa6(phi2, phi1, lam21)

1715 if r and isnon0(c1) and isnon0(c2):

1716 E = _ellipsoidal(datum, thomas_)

1717 if E.f:

1718 r1 = atan2(E.b_a * s1, c1)

1719 r2 = atan2(E.b_a * s2, c2)

1720

1721 j = (r2 + r1) * _0_5

1722 k = (r2 - r1) * _0_5

1723 sj, cj, sk, ck, h, _ = sincos2_(j, k, lam21 * _0_5)

1724

1725 h = fsumf_(sk**2, (ck * h)**2, -(sj * h)**2)

1726 u = _1_0 - h

1727 if isnon0(u) and isnon0(h):

1728 r = atan(sqrt0(h / u)) * 2 # == acos(1 - 2 * h)

1729 sr, cr = sincos2(r)

1730 if isnon0(sr):

1731 u = 2 * (sj * ck)**2 / u

1732 h = 2 * (sk * cj)**2 / h

1733 x = u + h

1734 y = u - h

1735

1736 s = r / sr

1737 e = 4 * s**2

1738 d = 2 * cr

1739 a = e * d

1740 b = 2 * r

1741 c = s - (a - d) * _0_5

1742 f = E.f * _0_25

1743

1744 t = fsumf_(a * x, -b * y, c * x**2, -d * y**2, e * x * y)

1745 r -= fsumf_(s * x, -y, -t * f * _0_25) * f * sr

1746 return r

1747

1748

1749 def vincentys(lat1, lon1, lat2, lon2, radius=R_M, wrap=False):

1750 '''Compute the distance between two (spherical) points using

1751 U{Vincenty's<https://WikiPedia.org/wiki/Great-circle_distance>}

1752 spherical formula.

1753

1754 @arg lat1: Start latitude (C{degrees}).

1755 @arg lon1: Start longitude (C{degrees}).

1756 @arg lat2: End latitude (C{degrees}).

1757 @arg lon2: End longitude (C{degrees}).

1758 @kwarg radius: Mean earth radius (C{meter}), datum (L{Datum})

1759 or ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or

1760 L{a_f2Tuple}) to use.

1761 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll

1762 B{C{lat2}} and B{C{lon2}} (C{bool}).

1763

1764 @return: Distance (C{meter}, same units as B{C{radius}}).

1765

1766 @raise UnitError: Invalid B{C{radius}}.

1767

1768 @see: Functions L{vincentys_}, L{cosineAndoyerLambert},

1769 L{cosineForsytheAndoyerLambert},L{cosineLaw}, L{equirectangular},

1770 L{euclidean}, L{flatLocal}/L{hubeny}, L{flatPolar},

1771 L{haversine} and L{thomas} and methods L{Ellipsoid.distance2},

1772 C{LatLon.distanceTo*} and C{LatLon.equirectangularTo}.

1773

1774 @note: See note at function L{vincentys_}.

1775 '''

1776 return _dS(vincentys_, radius, wrap, lat1, lon1, lat2, lon2)

1777

1778

1779 def vincentys_(phi2, phi1, lam21):

1780 '''Compute the I{angular} distance between two (spherical) points using

1781 U{Vincenty's<https://WikiPedia.org/wiki/Great-circle_distance>}

1782 spherical formula.

1783

1784 @arg phi2: End latitude (C{radians}).

1785 @arg phi1: Start latitude (C{radians}).

1786 @arg lam21: Longitudinal delta, M{end-start} (C{radians}).

1787

1788 @return: Angular distance (C{radians}).

1789

1790 @see: Functions L{vincentys}, L{cosineAndoyerLambert_},

1791 L{cosineForsytheAndoyerLambert_}, L{cosineLaw_},

1792 L{euclidean_}, L{flatLocal_}/L{hubeny_}, L{flatPolar_},

1793 L{haversine_} and L{thomas_}.

1794

1795 @note: Functions L{vincentys_}, L{haversine_} and L{cosineLaw_}

1796 produce equivalent results, but L{vincentys_} is suitable

1797 for antipodal points and slightly more expensive (M{3 cos,

1798 3 sin, 1 hypot, 1 atan2, 6 mul, 2 add}) than L{haversine_}

1799 (M{2 cos, 2 sin, 2 sqrt, 1 atan2, 5 mul, 1 add}) and

1800 L{cosineLaw_} (M{3 cos, 3 sin, 1 acos, 3 mul, 1 add}).

1801 '''

1802 s1, c1, s2, c2, s21, c21 = sincos2_(phi1, phi2, lam21)

1803

1804 c = c2 * c21

1805 x = s1 * s2 + c1 * c

1806 y = c1 * s2 - s1 * c

1807 return atan2(hypot(c2 * s21, y), x)

1808

1809# **) MIT License

1810#

1811# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

1812#

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

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

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

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

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

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

1819#

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

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

1822#

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

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

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

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

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

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

1829# OTHER DEALINGS IN THE SOFTWARE.