1

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

3

4 u'''(INTERNAL) Private base class L{LatLonBase} for elliposiodal, spherical

5 and N-vectorial C{LatLon}s.

6

7 After I{(C) Chris Veness 2011-2015} published under the same MIT Licence**,

8 see U{https://www.Movable-Type.co.UK/scripts/latlong.html},

9 U{<https://www.Movable-Type.co.UK/scripts/geodesy/docs/latlon-ellipsoidal.js.html>}

10 and U{https://www.Movable-Type.co.UK/scripts/latlong-vectors.html}.

11'''

12

13 from pygeodesy.basics import isscalar, isstr, _xinstanceof

14 from pygeodesy.constants import EPS, EPS0, EPS1, EPS4, INT0, R_M, \

15 _0_0, _0_5, _1_0

16# from pygeodesy.datums import _spherical_datum # from .formy

17 from pygeodesy.dms import F_D, F_DMS, latDMS, lonDMS, parse3llh

18 from pygeodesy.errors import _incompatible, IntersectionError, _TypeError, \

19 _ValueError, _xdatum, _xError, _xkwds, _xkwds_not

20 from pygeodesy.formy import antipode, compassAngle, cosineAndoyerLambert_, \

21 cosineForsytheAndoyerLambert_, cosineLaw, \

22 equirectangular, euclidean, flatLocal_, \

23 flatPolar, hartzell, haversine, isantipode, \

24 isnormal, normal, philam2n_xyz, \

25 _spherical_datum, thomas_, vincentys

26 from pygeodesy.interns import NN, _COMMASPACE_, _concentric_, _height_, \

27 _intersection_, _m_, _no_, _overlap_, _point_

28 from pygeodesy.iters import PointsIter, points2

29 from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS

30 from pygeodesy.named import _NamedBase, notOverloaded

31 from pygeodesy.namedTuples import Bounds2Tuple, LatLon2Tuple, PhiLam2Tuple, \

32 Trilaterate5Tuple, Vector3Tuple

33 from pygeodesy.props import deprecated_method, Property, Property_RO, \

34 property_RO, _update_all

35 from pygeodesy.streprs import Fmt, hstr

36 from pygeodesy.units import Distance_, Lat, Lon, Height, Radius, Radius_, \

37 Scalar, Scalar_

38 from pygeodesy.utily import _unrollon, unrollPI

39 from pygeodesy.vector2d import _circin6, Circin6Tuple, _circum3, Circum3Tuple, \

40 circum4_, Circum4Tuple, _radii11ABC

41 from pygeodesy.vector3d import nearestOn6, Vector3d

42

43 from math import asin, cos, degrees, fabs, radians

44

45__all__ = _ALL_LAZY.latlonBase

46__version__ = '23.03.30'

47

48

49 class LatLonBase(_NamedBase):

50 '''(INTERNAL) Base class for C{LatLon} points on spherical or

51 ellipsoidal earth models.

52 '''

53 _clipid = INT0 # polygonal clip, see .booleans

54 _datum = None # L{Datum}, to be overriden

55 _height = INT0 # height (C{meter}), default

56 _lat = 0 # latitude (C{degrees})

57 _lon = 0 # longitude (C{degrees})

58

59 def __init__(self, lat, lon, height=0, name=NN):

60 '''New C{LatLon}.

61

62 @arg lat: Latitude (C{degrees} or DMS C{str} with N or S suffix).

63 @arg lon: Longitude (C{degrees} or DMS C{str} with E or W suffix).

64 @kwarg height: Optional height (C{meter} above or below the earth surface).

65 @kwarg name: Optional name (C{str}).

66

67 @return: New instance (C{LatLon}).

68

69 @raise RangeError: Value of B{C{lat}} or B{C{lon}} outside the valid

70 range and C{rangerrors} set to C{True}.

71

72 @raise UnitError: Invalid B{C{lat}}, B{C{lon}} or B{C{height}}.

73

74 @example:

75

76 >>> p = LatLon(50.06632, -5.71475)

77 >>> q = LatLon('50°03′59″N', """005°42'53"W""")

78 '''

79 if name:

80 self.name = name

81

82 self._lat = Lat(lat) # parseDMS2(lat, lon)

83 self._lon = Lon(lon) # PYCHOK LatLon2Tuple

84 if height: # elevation

85 self._height = Height(height)

86

87 def __eq__(self, other):

88 return self.isequalTo(other)

89

90 def __ne__(self, other):

91 return not self.isequalTo(other)

92

93 def __str__(self):

94 return self.toStr(form=F_D, prec=6)

95

96 def antipode(self, height=None):

97 '''Return the antipode, the point diametrically opposite

98 to this point.

99

100 @kwarg height: Optional height of the antipode (C{meter}),

101 this point's height otherwise.

102

103 @return: The antipodal point (C{LatLon}).

104 '''

105 a, b = antipode(self.lat, self.lon) # PYCHOK LatLon2Tuple

106 h = self.height if height is None else Height(height)

107 return self.classof(a, b, height=h)

108

109 @deprecated_method

110 def bounds(self, wide, tall, radius=R_M): # PYCHOK no cover

111 '''DEPRECATED, use method C{boundsOf}.'''

112 return self.boundsOf(wide, tall, radius=radius)

113

114 def boundsOf(self, wide, tall, radius=R_M, height=None):

115 '''Return the SW and NE lat-/longitude of a great circle

116 bounding box centered at this location.

117

118 @arg wide: Longitudinal box width (C{meter}, same units as

119 B{C{radius}} or C{degrees} if B{C{radius}} is C{None}).

120 @arg tall: Latitudinal box size (C{meter}, same units as

121 B{C{radius}} or C{degrees} if B{C{radius}} is C{None}).

122 @kwarg radius: Mean earth radius (C{meter}).

123 @kwarg height: Height for C{latlonSW} and C{latlonNE} (C{meter}),

124 overriding the point's height.

125

126 @return: A L{Bounds2Tuple}C{(latlonSW, latlonNE)}, the

127 lower-left and upper-right corner (C{LatLon}).

128

129 @see: U{https://www.Movable-Type.co.UK/scripts/latlong-db.html}

130 '''

131 x = Scalar_(wide=wide) * _0_5

132 y = Scalar_(tall=tall) * _0_5

133 if radius is not None:

134 r = Radius_(radius)

135 c = cos(self.phi)

136 x = degrees(asin(x / r) / c) if fabs(c) > EPS0 else _0_0 # XXX

137 y = degrees(y / r)

138 x, y = fabs(x), fabs(y)

139

140 h = self.height if height is None else Height(height)

141 sw = self.classof(self.lat - y, self.lon - x, height=h)

142 ne = self.classof(self.lat + y, self.lon + x, height=h)

143 return Bounds2Tuple(sw, ne, name=self.name)

144

145 def chordTo(self, other, height=None):

146 '''Compute the length of the chord through the earth between

147 this and an other point.

148

149 @arg other: The other point (C{LatLon}).

150 @kwarg height: Overriding height for both points (C{meter})

151 or C{None} for each point's height.

152

153 @return: The chord length (conventionally C{meter}).

154

155 @raise TypeError: The B{C{other}} point is not C{LatLon}.

156 '''

157 def _v3d(ll):

158 t = ll.toEcef(height=height) # .toVector(Vector=Vector3d)

159 return Vector3d(t.x, t.y, t.z)

160

161 self.others(other)

162 return _v3d(self).minus(_v3d(other)).length

163

164 def circin6(self, point2, point3, eps=EPS4):

165 '''Return the radius and center of the I{inscribed} aka I{In-}circle

166 of the (planar) triangle formed by this and two other points.

167

168 @arg point2: Second point (C{LatLon}).

169 @arg point3: Third point (C{LatLon}).

170 @kwarg eps: Tolerance for function L{pygeodesy.trilaterate3d2}.

171

172 @return: L{Circin6Tuple}C{(radius, center, deltas, cA, cB, cC)}. The

173 C{center} and contact points C{cA}, C{cB} and C{cC}, each an

174 instance of this (sub-)class, are co-planar with this and the

175 two given points, see the B{Note} below.

176

177 @raise ImportError: Package C{numpy} not found, not installed or older

178 than version 1.10.

179

180 @raise IntersectionError: Near-coincident or -colinear points or

181 a trilateration or C{numpy} issue.

182

183 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.

184

185 @note: The C{center} is trilaterated in cartesian (ECEF) space and converted

186 back to geodetic lat-, longitude and height. The latter, conventionally

187 in C{meter} indicates whether the C{center} is above, below or on the

188 surface of the earth model. If C{deltas} is C{None}, the C{center} is

189 I{un}ambigous. Otherwise C{deltas} is a L{LatLon3Tuple}C{(lat, lon,

190 height)} representing the differences between both results from

191 L{pygeodesy.trilaterate3d2} and C{center} is the mean thereof.

192

193 @see: Function L{pygeodesy.circin6}, method L{circum3}, U{Incircle

194 <https://MathWorld.Wolfram.com/Incircle.html>} and U{Contact Triangle

195 <https://MathWorld.Wolfram.com/ContactTriangle.html>}.

196 '''

197 try:

198 cs = self._toCartesian3(point2, point3)

199 r, c, d, cA, cB, cC = _circin6(*cs, eps=eps, useZ=True, dLL3=True,

200 datum=self.datum) # PYCHOK unpack

201 return Circin6Tuple(r, c.toLatLon(), d, cA.toLatLon(), cB.toLatLon(), cC.toLatLon())

202 except (AssertionError, TypeError, ValueError) as x:

203 raise _xError(x, point=self, point2=point2, point3=point3)

204

205 def circum3(self, point2, point3, circum=True, eps=EPS4):

206 '''Return the radius and center of the smallest circle I{through} or I{containing}

207 this and two other points.

208

209 @arg point2: Second point (C{LatLon}).

210 @arg point3: Third point (C{LatLon}).

211 @kwarg circum: If C{True} return the C{circumradius} and C{circumcenter},

212 always, ignoring the I{Meeus}' Type I case (C{bool}).

213 @kwarg eps: Tolerance for function L{pygeodesy.trilaterate3d2}.

214

215 @return: A L{Circum3Tuple}C{(radius, center, deltas)}. The C{center}, an

216 instance of this (sub-)class, is co-planar with this and the two

217 given points. If C{deltas} is C{None}, the C{center} is

218 I{un}ambigous. Otherwise C{deltas} is a L{LatLon3Tuple}C{(lat,

219 lon, height)} representing the difference between both results

220 from L{pygeodesy.trilaterate3d2} and C{center} is the mean thereof.

221

222 @raise ImportError: Package C{numpy} not found, not installed or older than

223 version 1.10.

224

225 @raise IntersectionError: Near-concentric, -coincident or -colinear points,

226 incompatible C{Ecef} classes or a trilateration

227 or C{numpy} issue.

228

229 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.

230

231 @note: The C{center} is trilaterated in cartesian (ECEF) space and converted

232 back to geodetic lat-, longitude and height. The latter, conventionally

233 in C{meter} indicates whether the C{center} is above, below or on the

234 surface of the earth model. If C{deltas} is C{None}, the C{center} is

235 I{un}ambigous. Otherwise C{deltas} is a L{LatLon3Tuple}C{(lat, lon,

236 height)} representing the difference between both results from

237 L{pygeodesy.trilaterate3d2} and C{center} is the mean thereof.

238

239 @see: Function L{pygeodesy.circum3} and methods L{circin6} and L{circum4_}.

240 '''

241 try:

242 cs = self._toCartesian3(point2, point3)

243 r, c, d = _circum3(*cs, circum=circum, eps=eps, useZ=True, dLL3=True, # XXX -3d2

244 clas=cs[0].classof, datum=self.datum) # PYCHOK unpack

245 return Circum3Tuple(r, c.toLatLon(), d)

246 except (AssertionError, TypeError, ValueError) as x:

247 raise _xError(x, point=self, point2=point2, point3=point3, circum=circum)

248

249 def circum4_(self, *points):

250 '''Best-fit a sphere through this and two or more other points.

251

252 @arg points: The other points (each a C{LatLon}).

253

254 @return: L{Circum4Tuple}C{(radius, center, rank, residuals)} with C{center}

255 an instance of this (sub-)class.

256

257 @raise ImportError: Package C{numpy} not found, not installed or older than

258 version 1.10.

259

260 @raise NumPyError: Some C{numpy} issue.

261

262 @raise TypeError: One of the B{C{points}} invalid.

263

264 @raise ValueError: Too few B{C{points}}.

265

266 @see: Function L{pygeodesy.circum4_} and L{circum3}.

267 '''

268 C = self._toCartesianEcef

269 c = C(point=self)

270 t = circum4_(c, Vector=c.classof, *(C(i=i, points=p) for i, p in enumerate(points)))

271 c = t.center.toLatLon(LatLon=self.classof, name=t.name)

272 return Circum4Tuple(t.radius, c, t.rank, t.residuals, name=c.name)

273

274 @property

275 def clipid(self):

276 '''Get the (polygonal) clip (C{int}).

277 '''

278 return self._clipid

279

280 @clipid.setter # PYCHOK setter!

281 def clipid(self, clipid):

282 '''Get the (polygonal) clip (C{int}).

283 '''

284 self._clipid = int(clipid)

285

286 @deprecated_method

287 def compassAngle(self, other, adjust=True, wrap=False): # PYCHOK no cover

288 '''DEPRECATED, use method L{compassAngleTo}.'''

289 return self.compassAngleTo(other, adjust=adjust, wrap=wrap)

290

291 def compassAngleTo(self, other, adjust=True, wrap=False):

292 '''Return the angle from North for the direction vector between

293 this and an other point.

294

295 Suitable only for short, non-near-polar vectors up to a few

296 hundred Km or Miles. Use method C{initialBearingTo} for

297 larger distances.

298

299 @arg other: The other point (C{LatLon}).

300 @kwarg adjust: Adjust the longitudinal delta by the

301 cosine of the mean latitude (C{bool}).

302 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes and

303 longitudinal delta (C{bool}).

304

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

306

307 @raise TypeError: The B{C{other}} point is not C{LatLon}.

308

309 @note: Courtesy of Martin Schultz.

310

311 @see: U{Local, flat earth approximation

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

313 '''

314 self.others(other)

315 return compassAngle(self.lat, self.lon, other.lat, other.lon,

316 adjust=adjust, wrap=wrap)

317

318 def cosineAndoyerLambertTo(self, other, wrap=False):

319 '''Compute the distance between this and an other point using

320 the U{Andoyer-Lambert correction<https://navlib.net/wp-content/uploads/

321 2013/10/admiralty-manual-of-navigation-vol-1-1964-english501c.pdf>} of

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

323 formula.

324

325 @arg other: The other point (C{LatLon}).

326 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).

327

328 @return: Distance (C{meter}, same units as the axes of

329 this point's datum ellipsoid).

330

331 @raise TypeError: The B{C{other}} point is not C{LatLon}.

332

333 @see: Function L{pygeodesy.cosineAndoyerLambert} and methods

334 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo},

335 C{distanceTo*}, L{equirectangularTo}, L{euclideanTo},

336 L{flatLocalTo}/L{hubenyTo}, L{flatPolarTo}, L{haversineTo},

337 L{thomasTo} and L{vincentysTo}.

338 '''

339 return self._distanceTo_(cosineAndoyerLambert_, other, wrap=wrap)

340

341 def cosineForsytheAndoyerLambertTo(self, other, wrap=False):

342 '''Compute the distance between this and an other point using

343 the U{Forsythe-Andoyer-Lambert correction

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

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

346 formula.

347

348 @arg other: The other point (C{LatLon}).

349 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).

350

351 @return: Distance (C{meter}, same units as the axes of

352 this point's datum ellipsoid).

353

354 @raise TypeError: The B{C{other}} point is not C{LatLon}.

355

356 @see: Function L{pygeodesy.cosineForsytheAndoyerLambert} and methods

357 L{cosineAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},

358 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},

359 L{flatPolarTo}, L{haversineTo}, L{thomasTo} and L{vincentysTo}.

360 '''

361 return self._distanceTo_(cosineForsytheAndoyerLambert_, other, wrap=wrap)

362

363 def cosineLawTo(self, other, radius=None, wrap=False):

364 '''Compute the distance between this and an other point using the

365 U{spherical Law of Cosines

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

367 formula.

368

369 @arg other: The other point (C{LatLon}).

370 @kwarg radius: Mean earth radius (C{meter}) or C{None}

371 for the mean radius of this point's datum

372 ellipsoid.

373 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).

374

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

376

377 @raise TypeError: The B{C{other}} point is not C{LatLon}.

378

379 @see: Function L{pygeodesy.cosineLaw} and methods L{cosineAndoyerLambertTo},

380 L{cosineForsytheAndoyerLambertTo}, C{distanceTo*}, L{equirectangularTo},

381 L{euclideanTo}, L{flatLocalTo}/L{hubenyTo}, L{flatPolarTo},

382 L{haversineTo}, L{thomasTo} and L{vincentysTo}.

383 '''

384 return self._distanceTo(cosineLaw, other, radius, wrap=wrap)

385

386 @property_RO

387 def datum(self): # PYCHOK no cover

388 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.

389 '''

390 notOverloaded(self)

391

392 def destinationXyz(self, delta, LatLon=None, **LatLon_kwds):

393 '''Calculate the destination using a I{local} delta from this point.

394

395 @arg delta: Local delta to the destination (L{XyzLocal}, L{Enu},

396 L{Ned} or L{Local9Tuple}).

397 @kwarg LatLon: Optional (geodetic) class to return the destination

398 or C{None}.

399 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword

400 arguments, ignored if C{B{LatLon} is None}.

401

402 @return: Destination as a C{B{LatLon}(lat, lon, **B{LatLon_kwds})}

403 instance or if C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat,

404 lon, height)} respectively L{LatLon4Tuple}C{(lat, lon,

405 height, datum)} depending on whether a C{datum} keyword

406 is un-/specified.

407

408 @raise TypeError: Invalid B{C{delta}}, B{C{LatLon}} or B{C{LatLon_kwds}}.

409 '''

410 t = self._ltp._local2ecef(delta, nine=True)

411 return t.toLatLon(LatLon=LatLon, **_xkwds(LatLon_kwds, name=self.name))

412

413 def _distanceTo(self, func, other, radius, **options):

414 '''(INTERNAL) Helper for methods C{<func>To}.

415 '''

416 self.others(other) # up=2

417 if radius is None:

418 radius = self._datum.ellipsoid.R1 if self._datum else R_M

419 return func(self.lat, self.lon, other.lat, other.lon,

420 radius=radius, **options)

421

422 def _distanceTo_(self, func_, other, wrap=False):

423 '''(INTERNAL) Helper for (ellipsoidal) methods C{<func>To}.

424 '''

425 self.others(other) # up=2

426 r, _ = unrollPI(self.lam, other.lam, wrap=wrap)

427 r = func_(other.phi, self.phi, r, datum=self.datum)

428 return r * self.datum.ellipsoid.a

429

430 @Property_RO

431 def Ecef(self):

432 '''Get the ECEF I{class} (L{EcefKarney}), I{lazily}.

433 '''

434 return _MODS.ecef.EcefKarney # default

435

436 @Property_RO

437 def _Ecef_forward(self):

438 '''(INTERNAL) Helper for L{_ecef9} and L{toEcef} (C{callable}).

439 '''

440 return self.Ecef(self.datum, name=self.name).forward

441

442 @Property_RO

443 def _ecef9(self):

444 '''(INTERNAL) Helper for L{toCartesian}, L{toEcef} and L{toCartesian} (L{Ecef9Tuple}).

445 '''

446 return self._Ecef_forward(self, M=True)

447

448 @deprecated_method

449 def equals(self, other, eps=None): # PYCHOK no cover

450 '''DEPRECATED, use method L{isequalTo}.'''

451 return self.isequalTo(other, eps=eps)

452

453 @deprecated_method

454 def equals3(self, other, eps=None): # PYCHOK no cover

455 '''DEPRECATED, use method L{isequalTo3}.'''

456 return self.isequalTo3(other, eps=eps)

457

458 def equirectangularTo(self, other, radius=None, **options):

459 '''Compute the distance between this and an other point

460 using the U{Equirectangular Approximation / Projection

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

462

463 Suitable only for short, non-near-polar distances up to a

464 few hundred Km or Miles. Use method L{haversineTo} or

465 C{distanceTo*} for more accurate and/or larger distances.

466

467 See function L{pygeodesy.equirectangular_} for more details,

468 the available B{C{options}} and errors raised.

469

470 @arg other: The other point (C{LatLon}).

471 @kwarg radius: Mean earth radius (C{meter}) or C{None} for

472 the mean radius of this point's datum ellipsoid.

473 @kwarg options: Optional keyword arguments for function

474 L{pygeodesy.equirectangular}.

475

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

477

478 @raise TypeError: The B{C{other}} point is not C{LatLon}.

479

480 @see: Function L{pygeodesy.equirectangular} and methods L{cosineAndoyerLambertTo},

481 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},

482 C{euclideanTo}, L{flatLocalTo}/L{hubenyTo}, L{flatPolarTo},

483 L{haversineTo}, L{thomasTo} and L{vincentysTo}.

484 '''

485 return self._distanceTo(equirectangular, other, radius, **options)

486

487 def euclideanTo(self, other, radius=None, **options):

488 '''Approximate the C{Euclidian} distance between this and

489 an other point.

490

491 See function L{pygeodesy.euclidean} for the available B{C{options}}.

492

493 @arg other: The other point (C{LatLon}).

494 @kwarg radius: Mean earth radius (C{meter}) or C{None} for

495 the mean radius of this point's datum ellipsoid.

496 @kwarg options: Optional keyword arguments for function

497 L{pygeodesy.euclidean}.

498

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

500

501 @raise TypeError: The B{C{other}} point is not C{LatLon}.

502

503 @see: Function L{pygeodesy.euclidean} and methods L{cosineAndoyerLambertTo},

504 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},

505 L{equirectangularTo}, L{flatLocalTo}/L{hubenyTo}, L{flatPolarTo},

506 L{haversineTo}, L{thomasTo} and L{vincentysTo}.

507 '''

508 return self._distanceTo(euclidean, other, radius, **options)

509

510 def flatLocalTo(self, other, radius=None, wrap=False):

511 '''Compute the distance between this and an other point using the

512 U{ellipsoidal Earth to plane projection

513 <https://WikiPedia.org/wiki/Geographical_distance#Ellipsoidal_Earth_projected_to_a_plane>}

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

515

516 @arg other: The other point (C{LatLon}).

517 @kwarg radius: Mean earth radius (C{meter}) or C{None} for the

518 equatorial radius of this point's datum/ellipsoid.

519 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).

520

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

522

523 @raise TypeError: The B{C{other}} point is not C{LatLon}.

524

525 @raise ValueError: Invalid B{C{radius}}.

526

527 @see: Function L{pygeodesy.flatLocal}/L{pygeodesy.hubeny}, methods

528 L{cosineAndoyerLambertTo}, L{cosineForsytheAndoyerLambertTo},

529 L{cosineLawTo}, C{distanceTo*}, L{equirectangularTo}, L{euclideanTo},

530 L{flatPolarTo}, L{haversineTo}, L{thomasTo} and L{vincentysTo} and

531 U{local, flat Earth approximation<https://www.edwilliams.org/avform.htm#flat>}.

532 '''

533 E = self.datum.ellipsoid

534 r = self._distanceTo_(flatLocal_, other, wrap=wrap) * E.a2_

535 a = E.a if radius in (None, 1, _1_0) else Radius(radius)

536 return r * a

537

538 hubenyTo = flatLocalTo # for Karl Hubeny

539

540 def flatPolarTo(self, other, radius=None, wrap=False):

541 '''Compute the distance between this and an other point using

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

543 Geographical_distance#Polar_coordinate_flat-Earth_formula>}formula.

544

545 @arg other: The other point (C{LatLon}).

546 @kwarg radius: Mean earth radius (C{meter}) or C{None} for the

547 mean radius of this point's datum ellipsoid.

548 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).

549

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

551

552 @raise TypeError: The B{C{other}} point is not C{LatLon}.

553

554 @see: Function L{pygeodesy.flatPolar} and methods L{cosineAndoyerLambertTo},

555 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},

556 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},

557 L{haversineTo}, L{thomasTo} and L{vincentysTo}.

558 '''

559 return self._distanceTo(flatPolar, other, radius, wrap=wrap)

560

561 def hartzell(self, los=None, earth=None):

562 '''Compute the intersection of a Line-Of-Sight (los) from this Point-Of-View

563 (pov) with this point's ellipsoid surface.

564

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

566 C{None} to point to the ellipsoid's center.

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

568 L{a_f2Tuple} or C{scalar} radius in C{meter}) overriding

569 this point's C{datum} ellipsoid.

570

571 @return: The ellipsoid intersection (C{LatLon}) or this very instance

572 if this C{pov's height} is C{0}.

573

574 @raise IntersectionError: Null C{pov} or B{C{los}} vector, this

575 C{pov's height} is negative or B{C{los}}

576 points outside the ellipsoid or in an

577 opposite direction.

578

579 @raise TypeError: Invalid B{C{los}}.

580

581 @see: Function C{hartzell} for further details.

582 '''

583 h = self.height

584 if not h:

585 r = self

586 elif h < 0:

587 raise IntersectionError(pov=self, los=los, height=h, txt=_no_(_height_))

588 elif los is None:

589 d = self.datum if earth is None else _spherical_datum(earth)

590 r = self.dup(datum=d, height=0, name=self.hartzell.__name__)

591 else:

592 c = self.toCartesian()

593 r = hartzell(c, los=los, earth=earth or self.datum, LatLon=self.classof)

594 return r

595

596 def haversineTo(self, other, radius=None, wrap=False):

597 '''Compute the distance between this and an other point using the

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

599 formula.

600

601 @arg other: The other point (C{LatLon}).

602 @kwarg radius: Mean earth radius (C{meter}) or C{None} for

603 the mean radius of this point's datum ellipsoid.

604 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).

605

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

607

608 @raise TypeError: The B{C{other}} point is not C{LatLon}.

609

610 @see: Function L{pygeodesy.haversine} and methods L{cosineAndoyerLambertTo},

611 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},

612 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},

613 L{flatPolarTo}, L{thomasTo} and L{vincentysTo}.

614 '''

615 return self._distanceTo(haversine, other, radius, wrap=wrap)

616

617 def _havg(self, other, f=_0_5):

618 '''(INTERNAL) Weighted, average height.

619

620 @arg other: An other point (C{LatLon}).

621 @kwarg f: Optional fraction (C{float}).

622

623 @return: Average, fractional height (C{float}).

624 '''

625 return _MODS.fmath.favg(self.height, other.height, f=f)

626

627 @Property

628 def height(self):

629 '''Get the height (C{meter}).

630 '''

631 return self._height

632

633 @height.setter # PYCHOK setter!

634 def height(self, height):

635 '''Set the height (C{meter}).

636

637 @raise TypeError: Invalid B{C{height}} C{type}.

638

639 @raise ValueError: Invalid B{C{height}}.

640 '''

641 h = Height(height)

642 if self._height != h:

643 _update_all(self)

644 self._height = h

645

646 def height4(self, earth=None, normal=True, LatLon=None, **LatLon_kwds):

647 '''Compute the height above or below and the projection on this datum's

648 ellipsoid surface.

649

650 @kwarg earth: A datum, ellipsoid, triaxial ellipsoid or earth radius

651 I{overriding} this datum (L{Datum}, L{Ellipsoid},

652 L{Ellipsoid2}, L{a_f2Tuple}, L{Triaxial}, L{Triaxial_},

653 L{JacobiConformal} or C{meter}, conventionally).

654 @kwarg normal: If C{True} the projection is the nearest point on the

655 ellipsoid's surface, otherwise the intersection of the

656 radial line to the center and the ellipsoid's surface.

657 @kwarg LatLon: Optional class to return the height and projection

658 (C{LatLon}) or C{None}.

659 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword arguments,

660 ignored if C{B{LatLon} is None}.

661

662 @note: Use keyword argument C{height=0} to override C{B{LatLon}.height}

663 to {0} or any other C{scalar}, conventionally in C{meter}.

664

665 @return: An instance of B{C{LatLon}} or if C{B{LatLon} is None}, a

666 L{Vector4Tuple}C{(x, y, z, h)} with the I{projection} C{x}, C{y}

667 and C{z} coordinates and height C{h} in C{meter}, conventionally.

668

669 @raise TriaxialError: No convergence in triaxial root finding.

670

671 @raise TypeError: Invalid B{C{earth}}.

672

673 @see: L{Ellipsoid.height4} and L{Triaxial_.height4} for more information.

674 '''

675 c = self.toCartesian()

676 if LatLon is None:

677 r = c.height4(earth=earth, normal=normal)

678 else:

679 r = c.height4(earth=earth, normal=normal, Cartesian=c.classof, height=0)

680 r = r.toLatLon(LatLon=LatLon, **_xkwds(LatLon_kwds, height=r.height))

681 return r

682

683 def heightStr(self, prec=-2, m=_m_):

684 '''Return this B{C{height}} as C{str}ing.

685

686 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).

687 @kwarg m: Optional unit of the height (C{str}).

688

689 @see: Function L{pygeodesy.hstr}.

690 '''

691 return hstr(self.height, prec=prec, m=m)

692

693 @deprecated_method

694 def isantipode(self, other, eps=EPS): # PYCHOK no cover

695 '''DEPRECATED, use method L{isantipodeTo}.'''

696 return self.isantipodeTo(other, eps=eps)

697

698 def isantipodeTo(self, other, eps=EPS):

699 '''Check whether this and an other point are antipodal,

700 on diametrically opposite sides of the earth.

701

702 @arg other: The other point (C{LatLon}).

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

704

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

706 tolerance, C{False} otherwise.

707 '''

708 return isantipode(self.lat, self.lon,

709 other.lat, other.lon, eps=eps)

710

711 @Property_RO

712 def isEllipsoidal(self):

713 '''Check whether this point is ellipsoidal (C{bool} or C{None} if unknown).

714 '''

715 return self.datum.isEllipsoidal if self._datum else None

716

717 @Property_RO

718 def isEllipsoidalLatLon(self):

719 '''Get C{LatLon} base.

720 '''

721 return False

722

723 def isequalTo(self, other, eps=None):

724 '''Compare this point with an other point, I{ignoring} height.

725

726 @arg other: The other point (C{LatLon}).

727 @kwarg eps: Tolerance for equality (C{degrees}).

728

729 @return: C{True} if both points are identical,

730 I{ignoring} height, C{False} otherwise.

731

732 @raise TypeError: The B{C{other}} point is not C{LatLon}

733 or mismatch of the B{C{other}} and

734 this C{class} or C{type}.

735

736 @raise UnitError: Invalid B{C{eps}}.

737

738 @see: Method L{isequalTo3}.

739 '''

740 self.others(other)

741

742 return _isequalTo(self, other, eps=Scalar_(eps=eps)) if eps else \

743 (self.lat == other.lat and self.lon == other.lon)

744

745 def isequalTo3(self, other, eps=None):

746 '''Compare this point with an other point, I{including} height.

747

748 @arg other: The other point (C{LatLon}).

749 @kwarg eps: Tolerance for equality (C{degrees}).

750

751 @return: C{True} if both points are identical

752 I{including} height, C{False} otherwise.

753

754 @raise TypeError: The B{C{other}} point is not C{LatLon}

755 or mismatch of the B{C{other}} and

756 this C{class} or C{type}.

757

758 @see: Method L{isequalTo}.

759 '''

760 return self.height == other.height and self.isequalTo(other, eps=eps)

761

762 @Property_RO

763 def isnormal(self):

764 '''Return C{True} if this point is normal (C{bool}),

765 meaning C{abs(lat) <= 90} and C{abs(lon) <= 180}.

766

767 @see: Functions L{pygeodesy.isnormal} and

768 L{pygeodesy.normal}.

769 '''

770 return isnormal(self.lat, self.lon, eps=0)

771

772 @Property_RO

773 def isSpherical(self):

774 '''Check whether this point is spherical (C{bool} or C{None} if unknown).

775 '''

776 return self.datum.isSpherical if self._datum else None

777

778 @Property_RO

779 def lam(self):

780 '''Get the longitude (B{C{radians}}).

781 '''

782 return radians(self.lon)

783

784 @Property

785 def lat(self):

786 '''Get the latitude (C{degrees90}).

787 '''

788 return self._lat

789

790 @lat.setter # PYCHOK setter!

791 def lat(self, lat):

792 '''Set the latitude (C{str[N|S]} or C{degrees}).

793

794 @raise ValueError: Invalid B{C{lat}}.

795 '''

796 lat = Lat(lat) # parseDMS(lat, suffix=_NS_, clip=90)

797 if self._lat != lat:

798 _update_all(self)

799 self._lat = lat

800

801 @Property

802 def latlon(self):

803 '''Get the lat- and longitude (L{LatLon2Tuple}C{(lat, lon)}).

804 '''

805 return LatLon2Tuple(self._lat, self._lon, name=self.name)

806

807 @latlon.setter # PYCHOK setter!

808 def latlon(self, latlonh):

809 '''Set the lat- and longitude and optionally the height

810 (2- or 3-tuple or comma- or space-separated C{str}

811 of C{degrees90}, C{degrees180} and C{meter}).

812

813 @raise TypeError: Height of B{C{latlonh}} not C{scalar} or

814 B{C{latlonh}} not C{list} or C{tuple}.

815

816 @raise ValueError: Invalid B{C{latlonh}} or M{len(latlonh)}.

817

818 @see: Function L{pygeodesy.parse3llh} to parse a B{C{latlonh}}

819 string into a 3-tuple (lat, lon, h).

820 '''

821 if isstr(latlonh):

822 latlonh = parse3llh(latlonh, height=self.height)

823 else:

824 _xinstanceof(list, tuple, latlonh=latlonh)

825 if len(latlonh) == 3:

826 h = Height(latlonh[2], name=Fmt.SQUARE(latlonh=2))

827 elif len(latlonh) != 2:

828 raise _ValueError(latlonh=latlonh)

829 else:

830 h = self.height

831

832 llh = Lat(latlonh[0]), Lon(latlonh[1]), h # parseDMS2(latlonh[0], latlonh[1])

833 if (self._lat, self._lon, self._height) != llh:

834 _update_all(self)

835 self._lat, self._lon, self._height = llh

836

837 def latlon2(self, ndigits=0):

838 '''Return this point's lat- and longitude in C{degrees}, rounded.

839

840 @kwarg ndigits: Number of (decimal) digits (C{int}).

841

842 @return: A L{LatLon2Tuple}C{(lat, lon)}, both C{float}

843 and rounded away from zero.

844

845 @note: The C{round}ed values are always C{float}, also

846 if B{C{ndigits}} is omitted.

847 '''

848 return LatLon2Tuple(round(self.lat, ndigits),

849 round(self.lon, ndigits), name=self.name)

850

851 @deprecated_method

852 def latlon_(self, ndigits=0): # PYCHOK no cover

853 '''DEPRECATED, use method L{latlon2}.'''

854 return self.latlon2(ndigits=ndigits)

855

856 latlon2round = latlon_ # PYCHOK no cover

857

858 @Property_RO

859 def latlonheight(self):

860 '''Get the lat-, longitude and height (L{LatLon3Tuple}C{(lat, lon, height)}).

861 '''

862 return self.latlon.to3Tuple(self.height)

863

864 @Property

865 def lon(self):

866 '''Get the longitude (C{degrees180}).

867 '''

868 return self._lon

869

870 @lon.setter # PYCHOK setter!

871 def lon(self, lon):

872 '''Set the longitude (C{str[E|W]} or C{degrees}).

873

874 @raise ValueError: Invalid B{C{lon}}.

875 '''

876 lon = Lon(lon) # parseDMS(lon, suffix=_EW_, clip=180)

877 if self._lon != lon:

878 _update_all(self)

879 self._lon = lon

880

881 @Property_RO

882 def _ltp(self):

883 '''(INTERNAL) Cache for L{toLtp}.

884 '''

885 return _MODS.ltp.Ltp(self, ecef=self.Ecef(self.datum), name=self.name)

886

887 def nearestOn6(self, points, closed=False, height=None, wrap=False):

888 '''Locate the point on a path or polygon closest to this point.

889

890 Points are converted to and distances are computed in

891 I{geocentric}, cartesian space.

892

893 @arg points: The path or polygon points (C{LatLon}[]).

894 @kwarg closed: Optionally, close the polygon (C{bool}).

895 @kwarg height: Optional height, overriding the height of

896 this and all other points (C{meter}). If

897 C{None}, take the height of points into

898 account for distances.

899 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes

900 (C{bool}).

901

902 @return: A L{NearestOn6Tuple}C{(closest, distance, fi, j,

903 start, end)} with the C{closest}, the C{start}

904 and the C{end} point each an instance of this

905 C{LatLon} and C{distance} in C{meter}, same

906 units as the cartesian axes.

907

908 @raise PointsError: Insufficient number of B{C{points}}.

909

910 @raise TypeError: Some B{C{points}} or some B{C{points}}'

911 C{Ecef} invalid.

912

913 @raise ValueError: Some B{C{points}}' C{Ecef} is incompatible.

914

915 @see: Function L{pygeodesy.nearestOn6}.

916 '''

917 def _cs(Ps, h, w, C):

918 p = None # not used

919 for i, q in Ps.enumerate():

920 if w and i != 0:

921 q = _unrollon(p, q)

922 yield C(height=h, i=i, up=3, points=q)

923 p = q

924

925 C = self._toCartesianEcef # to verify datum and Ecef

926 Ps = self.PointsIter(points)

927

928 c = C(height=height, this=self) # this Cartesian

929 t = nearestOn6(c, _cs(Ps, height, wrap, C), closed=closed)

930 c, s, e = t.closest, t.start, t.end

931

932 kwds = _xkwds_not(None, LatLon=self.classof, # this LatLon

933 height=height)

934 r = self.Ecef(self.datum).reverse

935 p = r(c).toLatLon(**kwds)

936 s = r(s).toLatLon(**kwds) if s is not c else p

937 e = r(e).toLatLon(**kwds) if e is not c else p

938 return t.dup(closest=p, start=s, end=e)

939

940 def normal(self):

941 '''Normalize this point to C{abs(lat) <= 90} and C{abs(lon) <= 180}.

942

943 @return: C{True} if this point was I{normal}, C{False}

944 if it wasn't (but is now).

945

946 @see: Property L{isnormal} and function L{pygeodesy.normal}.

947 '''

948 n = self.isnormal

949 if not n:

950 self.latlon = normal(self.lat, self.lon)

951 return n

952

953 @Property_RO

954 def _N_vector(self):

955 '''(INTERNAL) Get the (C{nvectorBase._N_vector_})

956 '''

957 return _MODS.nvectorBase._N_vector_(*self.xyzh)

958

959 @Property_RO

960 def phi(self):

961 '''Get the latitude (B{C{radians}}).

962 '''

963 return radians(self.lat)

964

965 @Property_RO

966 def philam(self):

967 '''Get the lat- and longitude (L{PhiLam2Tuple}C{(phi, lam)}).

968 '''

969 return PhiLam2Tuple(self.phi, self.lam, name=self.name)

970

971 def philam2(self, ndigits=0):

972 '''Return this point's lat- and longitude in C{radians}, rounded.

973

974 @kwarg ndigits: Number of (decimal) digits (C{int}).

975

976 @return: A L{PhiLam2Tuple}C{(phi, lam)}, both C{float}

977 and rounded away from zero.

978

979 @note: The C{round}ed values are always C{float}, also

980 if B{C{ndigits}} is omitted.

981 '''

982 return PhiLam2Tuple(round(self.phi, ndigits),

983 round(self.lam, ndigits), name=self.name)

984

985 @Property_RO

986 def philamheight(self):

987 '''Get the lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}).

988 '''

989 return self.philam.to3Tuple(self.height)

990

991 @deprecated_method

992 def points(self, points, closed=True): # PYCHOK no cover

993 '''DEPRECATED, use method L{points2}.'''

994 return self.points2(points, closed=closed)

995

996 def points2(self, points, closed=True):

997 '''Check a path or polygon represented by points.

998

999 @arg points: The path or polygon points (C{LatLon}[])

1000 @kwarg closed: Optionally, consider the polygon closed,

1001 ignoring any duplicate or closing final

1002 B{C{points}} (C{bool}).

1003

1004 @return: A L{Points2Tuple}C{(number, points)}, an C{int}

1005 and C{list} or C{tuple}.

1006

1007 @raise PointsError: Insufficient number of B{C{points}}.

1008

1009 @raise TypeError: Some B{C{points}} are not C{LatLon}.

1010 '''

1011 return points2(points, closed=closed, base=self)

1012

1013 def PointsIter(self, points, loop=0, dedup=False):

1014 '''Return a C{PointsIter} iterator.

1015

1016 @arg points: The path or polygon points (C{LatLon}[])

1017 @kwarg loop: Number of loop-back points (non-negative C{int}).

1018 @kwarg dedup: Skip duplicate points (C{bool}).

1019

1020 @return: A new C{PointsIter} iterator.

1021

1022 @raise PointsError: Insufficient number of B{C{points}}.

1023 '''

1024 return PointsIter(points, base=self, loop=loop, dedup=dedup)

1025

1026 def radii11(self, point2, point3):

1027 '''Return the radii of the C{Circum-}, C{In-}, I{Soddy} and C{Tangent}

1028 circles of a (planar) triangle formed by this and two other points.

1029

1030 @arg point2: Second point (C{LatLon}).

1031 @arg point3: Third point (C{LatLon}).

1032

1033 @return: L{Radii11Tuple}C{(rA, rB, rC, cR, rIn, riS, roS, a, b, c, s)}.

1034

1035 @raise IntersectionError: Near-coincident or -colinear points.

1036

1037 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.

1038

1039 @see: Function L{pygeodesy.radii11}, U{Incircle

1040 <https://MathWorld.Wolfram.com/Incircle.html>}, U{Soddy Circles

1041 <https://MathWorld.Wolfram.com/SoddyCircles.html>} and U{Tangent

1042 Circles<https://MathWorld.Wolfram.com/TangentCircles.html>}.

1043 '''

1044 try:

1045 cs = self._toCartesian3(point2, point3)

1046 return _radii11ABC(*cs, useZ=True)[0]

1047 except (TypeError, ValueError) as x:

1048 raise _xError(x, point=self, point2=point2, point3=point3)

1049

1050 def _rhumb2(self, exact, radius):

1051 '''(INTERNAL) Get the C{rhumb} for this point's datum or for

1052 the earth model or earth B{C{radius}} if not C{None}.

1053 '''

1054 D = self.datum if radius is None else _spherical_datum(radius) # ellipsoidal OK

1055 r = D.ellipsoid.rhumbx if exact else \

1056 _MODS.rhumbx.Rhumb(D, exact=False, name=D.name)

1057 return r, D

1058

1059 def rhumbAzimuthTo(self, other, exact=False, radius=None):

1060 '''Return the azimuth (bearing) of a rhumb line (loxodrome)

1061 between this and an other (ellipsoidal) point.

1062

1063 @arg other: The other point (C{LatLon}).

1064 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),

1065 default C{False}.

1066 @kwarg radius: Optional earth radius (C{meter}) or earth model

1067 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or

1068 L{a_f2Tuple}), overriding this point's datum.

1069

1070 @return: Rhumb azimuth (compass C{degrees360}).

1071

1072 @raise TypeError: The B{C{other}} point is incompatible or

1073 B{C{radius}} is invalid.

1074 '''

1075 self.others(other)

1076 r, _ = self._rhumb2(exact, radius)

1077 C = _MODS.rhumbx.Caps

1078 return r.Inverse(self.lat, self.lon, other.lat, other.lon,

1079 outmask=C.AZIMUTH).azi12

1080

1081 def rhumbDestination(self, distance, azimuth, exact=False, radius=None, height=None):

1082 '''Return the destination point having travelled the given distance

1083 from this point along a rhumb line (loxodrome) at the given azimuth.

1084

1085 @arg distance: Distance travelled (C{meter}, same units as this

1086 point's datum (ellipsoid) axes or B{C{radius}},

1087 may be negative.

1088 @arg azimuth: Azimuth (bearing) at this point (compass C{degrees}).

1089 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),

1090 default C{False}.

1091 @kwarg radius: Optional earth radius (C{meter}) or earth model

1092 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or

1093 L{a_f2Tuple}), overriding this point's datum.

1094 @kwarg height: Optional height, overriding the default height

1095 (C{meter}).

1096

1097 @return: The destination point (ellipsoidal C{LatLon}).

1098

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

1100

1101 @raise ValueError: Invalid B{C{distance}}, B{C{azimuth}},

1102 B{C{radius}} or B{C{height}}.

1103 '''

1104 r, D = self._rhumb2(exact, radius)

1105 d = r.Direct(self.lat, self.lon, azimuth, distance)

1106 h = self.height if height is None else Height(height)

1107 return self.classof(d.lat2, d.lon2, datum=D, height=h)

1108

1109 def rhumbDistanceTo(self, other, exact=False, radius=None):

1110 '''Return the distance from this to an other point along

1111 a rhumb line (loxodrome).

1112

1113 @arg other: The other point (C{LatLon}).

1114 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),

1115 default C{False}.

1116 @kwarg radius: Optional earth radius (C{meter}) or earth model

1117 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or

1118 L{a_f2Tuple}), overriding this point's datum.

1119

1120 @return: Distance (C{meter}, the same units as this point's

1121 datum (ellipsoid) axes or B{C{radius}}.

1122

1123 @raise TypeError: The B{C{other}} point is incompatible or

1124 B{C{radius}} is invalid.

1125

1126 @raise ValueError: Invalid B{C{radius}}.

1127 '''

1128 self.others(other)

1129 r, _ = self._rhumb2(exact, radius)

1130 C = _MODS.rhumbx.Caps

1131 return r.Inverse(self.lat, self.lon, other.lat, other.lon,

1132 outmask=C.DISTANCE).s12

1133

1134 def rhumbLine(self, azimuth_other, exact=False, radius=None, name=NN, **caps):

1135 '''Get a rhumb line through this point at a given azimuth or

1136 through this and an other point.

1137

1138 @arg azimuth_other: The azimuth of the rhumb line (compass)

1139 C{degrees} or the other point (C{LatLon}).

1140 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),

1141 default C{False}.

1142 @kwarg radius: Optional earth radius (C{meter}) or earth model

1143 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or

1144 L{a_f2Tuple}), overriding this point's datum.

1145 @kwarg name: Optional name (C{str}).

1146 @kwarg caps: Optional C{caps}, see L{RhumbLine} C{B{caps}}.

1147

1148 @return: A L{RhumbLine} instance.

1149

1150 @raise TypeError: Invalid B{C{radius}} or BC{C{azimuth_other}}

1151 not a C{scalar} nor a C{LatLon}.

1152

1153 @see: Classes L{RhumbLine} and L{Rhumb}, property L{Rhumb.exact}

1154 and methods L{Rhumb.DirectLine} and L{Rhumb.InverseLine}.

1155 '''

1156 r, _ = self._rhumb2(exact, radius)

1157 a = azimuth_other

1158 if isscalar(a):

1159 r = r.DirectLine(self.lat, self.lon, a,

1160 name=name or self.name, **caps)

1161 elif isinstance(a, LatLonBase):

1162 self.others(a)

1163 r = r.InverseLine(self.lat, self.lon, a.lat, a.lon,

1164 name=name or self.name, **caps)

1165 else:

1166 raise _TypeError(azimuth_other=a)

1167 return r

1168

1169 def rhumbMidpointTo(self, other, exact=False, radius=None,

1170 height=None, fraction=_0_5):

1171 '''Return the (loxodromic) midpoint on the rhumb line between

1172 this and an other point.

1173

1174 @arg other: The other point (C{LatLon}).

1175 @kwarg exact: If C{True}, use the I{exact} L{Rhumb} (C{bool}),

1176 default C{False}.

1177 @kwarg radius: Optional earth radius (C{meter}) or earth model

1178 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or

1179 L{a_f2Tuple}), overriding this point's datum.

1180 @kwarg height: Optional height, overriding the mean height

1181 (C{meter}).

1182 @kwarg fraction: Midpoint location from this point (C{scalar}),

1183 may be negative or greater than 1.0.

1184

1185 @return: The midpoint at the given B{C{fraction}} along the

1186 rhumb line (C{LatLon}).

1187

1188 @raise TypeError: The B{C{other}} point is incompatible or

1189 B{C{radius}} is invalid.

1190

1191 @raise ValueError: Invalid B{C{height}} or B{C{fraction}}.

1192 '''

1193 self.others(other)

1194 r, D = self._rhumb2(exact, radius)

1195 f = Scalar(fraction=fraction)

1196 d = r.Inverse(self.lat, self.lon, other.lat, other.lon)

1197 d = r.Direct( self.lat, self.lon, d.azi12, d.s12 * f)

1198 h = self._havg(other, f=f) if height is None else Height(height)

1199 return self.classof(d.lat2, d.lon2, datum=D, height=h)

1200

1201 def thomasTo(self, other, wrap=False):

1202 '''Compute the distance between this and an other point using

1203 U{Thomas'<https://apps.DTIC.mil/dtic/tr/fulltext/u2/703541.pdf>}

1204 formula.

1205

1206 @arg other: The other point (C{LatLon}).

1207 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).

1208

1209 @return: Distance (C{meter}, same units as the axes of

1210 this point's datum ellipsoid).

1211

1212 @raise TypeError: The B{C{other}} point is not C{LatLon}.

1213

1214 @see: Function L{pygeodesy.thomas} and methods L{cosineAndoyerLambertTo},

1215 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},

1216 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},

1217 L{flatPolarTo}, L{haversineTo} and L{vincentysTo}.

1218 '''

1219 return self._distanceTo_(thomas_, other, wrap=wrap)

1220

1221 @deprecated_method

1222 def to2ab(self): # PYCHOK no cover

1223 '''DEPRECATED, use property L{philam}.'''

1224 return self.philam

1225

1226 def toCartesian(self, height=None, Cartesian=None, **Cartesian_kwds):

1227 '''Convert this point to cartesian, I{geocentric} coordinates,

1228 also known as I{Earth-Centered, Earth-Fixed} (ECEF).

1229

1230 @kwarg height: Optional height, overriding this point's height

1231 (C{meter}, conventionally).

1232 @kwarg Cartesian: Optional class to return the geocentric

1233 coordinates (C{Cartesian}) or C{None}.

1234 @kwarg Cartesian_kwds: Optional, additional B{C{Cartesian}}

1235 keyword arguments, ignored if

1236 C{B{Cartesian} is None}.

1237

1238 @return: A B{C{Cartesian}} or if B{C{Cartesian}} is C{None},

1239 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M,

1240 datum)} with C{C=0} and C{M} if available.

1241

1242 @raise TypeError: Invalid B{C{Cartesian}} or B{C{Cartesian_kwds}}.

1243 '''

1244 r = self._ecef9 if height is None else self.toEcef(height=height)

1245 if Cartesian is not None: # class or .classof

1246 r = self._xnamed(Cartesian(r, **Cartesian_kwds))

1247 _xdatum(r.datum, self.datum)

1248 return r

1249

1250 def _toCartesian3(self, point2, point3):

1251 '''(INTERNAL) Convert this and 2 other points.

1252 '''

1253 return (self. toCartesian().copy(name=_point_), # copy to rename

1254 self._toCartesianEcef(up=3, point2=point2),

1255 self._toCartesianEcef(up=3, point3=point3))

1256

1257 def _toCartesianEcef(self, height=None, i=None, up=2, **name_point):

1258 '''(INTERNAL) Convert to cartesian and check Ecef's before and after.

1259 '''

1260 p = self.others(up=up, **name_point)

1261 c = p.toCartesian(height=height)

1262 E = self.Ecef

1263 if E:

1264 for p in (p, c):

1265 e = getattr(p, LatLonBase.Ecef.name, None)

1266 if e not in (None, E): # PYCHOK no cover

1267 n, _ = name_point.popitem()

1268 if i is not None:

1269 Fmt.SQUARE(n, i)

1270 raise _ValueError(n, e, txt=_incompatible(E.__name__))

1271 return c

1272

1273 def toEcef(self, height=None, M=False):

1274 '''Convert this point to I{geocentric} coordinates, also known as

1275 I{Earth-Centered, Earth-Fixed} (U{ECEF<https://WikiPedia.org/wiki/ECEF>}).

1276

1277 @kwarg height: Optional height, overriding this point's height

1278 (C{meter}, conventionally).

1279 @kwarg M: Optionally, include the rotation L{EcefMatrix} (C{bool}).

1280

1281 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)}

1282 with C{C=0} and C{M} if available.

1283

1284 @raise EcefError: A C{.datum} or an ECEF issue.

1285 '''

1286 return self._ecef9 if height in (None, self.height) else \

1287 self._Ecef_forward(self.lat, self.lon, height=height, M=M)

1288

1289 @deprecated_method

1290 def to3llh(self, height=None): # PYCHOK no cover

1291 '''DEPRECATED, use property L{latlonheight} or C{latlon.to3Tuple(B{height})}.'''

1292 return self.latlonheight if height in (None, self.height) else \

1293 self.latlon.to3Tuple(height)

1294

1295 def toLocal(self, Xyz=None, ltp=None, **Xyz_kwds):

1296 '''Convert this I{geodetic} point to I{local} C{X}, C{Y} and C{Z}.

1297

1298 @kwarg Xyz: Optional class to return C{X}, C{Y} and C{Z}

1299 (L{XyzLocal}, L{Enu}, L{Ned}) or C{None}.

1300 @kwarg ltp: The I{local tangent plane} (LTP) to use,

1301 overriding this point's LTP (L{Ltp}).

1302 @kwarg Xyz_kwds: Optional, additional B{C{Xyz}} keyword

1303 arguments, ignored if C{B{Xyz} is None}.

1304

1305 @return: An B{C{Xyz}} instance or if C{B{Xyz} is None},

1306 a L{Local9Tuple}C{(x, y, z, lat, lon, height,

1307 ltp, ecef, M)} with C{M=None}, always.

1308

1309 @raise TypeError: Invalid B{C{ltp}}.

1310 '''

1311 p = _MODS.ltp._xLtp(ltp, self._ltp)

1312 return p._ecef2local(self._ecef9, Xyz, Xyz_kwds)

1313

1314 def toLtp(self, Ecef=None):

1315 '''Return the I{local tangent plane} (LTP) for this point.

1316

1317 @kwarg Ecef: Optional ECEF I{class} (L{EcefKarney}, ...

1318 L{EcefYou}), overriding this point's C{Ecef}.

1319 '''

1320 return self._ltp if Ecef in (None, self.Ecef) else _MODS.ltp.Ltp(

1321 self, ecef=Ecef(self.datum), name=self.name)

1322

1323 def toNvector(self, h=None, Nvector=None, **Nvector_kwds):

1324 '''Convert this point to C{n-vector} (normal to the earth's

1325 surface) components, I{including height}.

1326

1327 @kwarg h: Optional height, overriding this point's

1328 height (C{meter}).

1329 @kwarg Nvector: Optional class to return the C{n-vector}

1330 components (C{Nvector}) or C{None}.

1331 @kwarg Nvector_kwds: Optional, additional B{C{Nvector}}

1332 keyword arguments, ignored if

1333 C{B{Nvector} is None}.

1334

1335 @return: A B{C{Nvector}} or a L{Vector4Tuple}C{(x, y, z, h)}

1336 if B{C{Nvector}} is C{None}.

1337

1338 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}.

1339 '''

1340 return self.toVector(Vector=Nvector, h=self.height if h is None else h,

1341 ll=self, **Nvector_kwds)

1342

1343 def toStr(self, form=F_DMS, joined=_COMMASPACE_, m=_m_, **prec_sep_s_D_M_S): # PYCHOK expected

1344 '''Convert this point to a "lat, lon[, +/-height]" string, formatted

1345 in the given C{B{form}at}.

1346

1347 @kwarg form: The lat-/longitude C{B{form}at} to use (C{str}), see

1348 functions L{pygeodesy.latDMS} or L{pygeodesy.lonDMS}.

1349 @kwarg joined: Separator to join the lat-, longitude and heigth

1350 strings (C{str} or C{None} or C{NN} for non-joined).

1351 @kwarg m: Optional unit of the height (C{str}), use C{None} to

1352 exclude height from the returned string.

1353 @kwarg prec_sep_s_D_M_S: Optional C{B{prec}ision}, C{B{sep}arator},

1354 B{C{s_D}}, B{C{s_M}}, B{C{s_S}} and B{C{s_DMS}} keyword

1355 arguments, see function L{pygeodesy.latDMS} or

1356 L{pygeodesy.lonDMS}.

1357

1358 @return: This point in the specified C{B{form}at}, etc. (C{str} or

1359 a 2- or 3-tuple C{(lat_str, lon_str[, height_str])} if

1360 C{B{joined}=NN} or C{B{joined}=None}).

1361

1362 @see: Function L{pygeodesy.latDMS} or L{pygeodesy.lonDMS} for more

1363 details about keyword arguments C{B{form}at}, C{B{prec}ision},

1364 C{B{sep}arator}, B{C{s_D}}, B{C{s_M}}, B{C{s_S}} and B{C{s_DMS}}.

1365

1366 @example:

1367

1368 >>> LatLon(51.4778, -0.0016).toStr() # 51°28′40″N, 000°00′06″W

1369 >>> LatLon(51.4778, -0.0016).toStr(F_D) # 51.4778°N, 000.0016°W

1370 >>> LatLon(51.4778, -0.0016, 42).toStr() # 51°28′40″N, 000°00′06″W, +42.00m

1371 '''

1372 t = (latDMS(self.lat, form=form, **prec_sep_s_D_M_S),

1373 lonDMS(self.lon, form=form, **prec_sep_s_D_M_S))

1374 if self.height and m is not None:

1375 t += (self.heightStr(m=m),)

1376 return joined.join(t) if joined else t

1377

1378 def toVector(self, Vector=None, **Vector_kwds):

1379 '''Convert this point to C{n-vector} (normal to the earth's

1380 surface) components, I{ignoring height}.

1381

1382 @kwarg Vector: Optional class to return the C{n-vector}

1383 components (L{Vector3d}) or C{None}.

1384 @kwarg Vector_kwds: Optional, additional B{C{Vector}}

1385 keyword arguments, ignored if

1386 C{B{Vector} is None}.

1387

1388 @return: A B{C{Vector}} or a L{Vector3Tuple}C{(x, y, z)}

1389 if B{C{Vector}} is C{None}.

1390

1391 @raise TypeError: Invalid B{C{Vector}} or B{C{kwds}}.

1392

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

1394 I{NOT} geocentric (ECEF) x, y and z coordinates!

1395 '''

1396 r = self._vector3tuple

1397 if Vector is not None:

1398 r = Vector(*r, **_xkwds(Vector_kwds, name=self.name))

1399 return r

1400

1401 def toVector3d(self):

1402 '''Convert this point to C{n-vector} (normal to the earth's

1403 surface) components, I{ignoring height}.

1404

1405 @return: Unit vector (L{Vector3d}).

1406

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

1408 I{NOT} geocentric (ECEF) x, y and z coordinates!

1409 '''

1410 return self._vector3d # XXX .unit()

1411

1412 @deprecated_method

1413 def to3xyz(self): # PYCHOK no cover

1414 '''DEPRECATED, use property L{xyz} or method L{toNvector}, L{toVector},

1415 L{toVector3d} or perhaps (geocentric) L{toEcef}.'''

1416 return self.xyz # self.toVector()

1417

1418 @Property_RO

1419 def _vector3d(self):

1420 '''(INTERNAL) Cache for L{toVector3d}.

1421 '''

1422 return self.toVector(Vector=Vector3d) # XXX .unit()

1423

1424 @Property_RO

1425 def _vector3tuple(self):

1426 '''(INTERNAL) Cache for L{toVector}.

1427 '''

1428 return philam2n_xyz(self.phi, self.lam, name=self.name)

1429

1430 def vincentysTo(self, other, radius=None, wrap=False):

1431 '''Compute the distance between this and an other point using

1432 U{Vincenty's<https://WikiPedia.org/wiki/Great-circle_distance>}

1433 spherical formula.

1434

1435 @arg other: The other point (C{LatLon}).

1436 @kwarg radius: Mean earth radius (C{meter}) or C{None}

1437 for the mean radius of this point's datum

1438 ellipsoid.

1439 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).

1440

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

1442

1443 @raise TypeError: The B{C{other}} point is not C{LatLon}.

1444

1445 @see: Function L{pygeodesy.vincentys} and methods L{cosineAndoyerLambertTo},

1446 L{cosineForsytheAndoyerLambertTo}, L{cosineLawTo}, C{distanceTo*},

1447 L{equirectangularTo}, L{euclideanTo}, L{flatLocalTo}/L{hubenyTo},

1448 L{flatPolarTo}, L{haversineTo} and L{thomasTo}.

1449 '''

1450 return self._distanceTo(vincentys, other, radius, wrap=wrap)

1451

1452 @Property_RO

1453 def xyz(self):

1454 '''Get the C{n-vector} X, Y and Z components (L{Vector3Tuple}C{(x, y, z)})

1455

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

1457 geocentric (ECEF) x, y and z coordinates!

1458 '''

1459 return self.toVector(Vector=Vector3Tuple)

1460

1461 @Property_RO

1462 def xyzh(self):

1463 '''Get the C{n-vector} X, Y, Z and H components (L{Vector4Tuple}C{(x, y, z, h)})

1464

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

1466 geocentric (ECEF) x, y and z coordinates!

1467 '''

1468 return self.xyz.to4Tuple(self.height)

1469

1470

1471 def _isequalTo(point1, point2, eps=EPS): # in .ellipsoidalBaseDI._intersect3._on, .formy

1472 '''(INTERNAL) Compare point lat-/lon without type.

1473 '''

1474 return fabs(point1.lat - point2.lat) <= eps and \

1475 fabs(point1.lon - point2.lon) <= eps

1476

1477

1478 def _isequalTo_(point1, point2, eps=EPS): # PYCHOK in .formy

1479 '''(INTERNAL) Compare point phi-/lam without type.

1480 '''

1481 return fabs(point1.phi - point2.phi) <= eps and \

1482 fabs(point1.lam - point2.lam) <= eps

1483

1484

1485 def _trilaterate5(p1, d1, p2, d2, p3, d3, area=True, eps=EPS1,

1486 radius=R_M, wrap=False):

1487 '''(INTERNAL) Trilaterate three points by area overlap or by

1488 perimeter intersection of three circles.

1489

1490 @note: The B{C{radius}} is only needed for both the n-vectorial

1491 and C{sphericalTrigonometry.LatLon.distanceTo} methods and

1492 silently ignored by the C{ellipsoidalExact}, C{-GeodSolve},

1493 C{-Karney} and C{-Vincenty.LatLon.distanceTo} methods.

1494 '''

1495 r1 = Distance_(distance1=d1)

1496 r2 = Distance_(distance2=d2)

1497 r3 = Distance_(distance3=d3)

1498

1499 m = 0 if area else (r1 + r2 + r3)

1500 pc = 0

1501 t = []

1502 for _ in range(3):

1503 try: # intersection of circle (p1, r1) and (p2, r2)

1504 c1, c2 = p1.intersections2(r1, p2, r2, wrap=wrap)

1505

1506 if area: # check overlap

1507 if c1 is c2: # abutting

1508 c = c1

1509 else: # nearest point on radical

1510 c = p3.nearestOn(c1, c2, within=True, wrap=wrap)

1511 d = r3 - p3.distanceTo(c, radius=radius, wrap=wrap)

1512 if d > eps: # sufficient overlap

1513 t.append((d, c))

1514 m = max(m, d)

1515

1516 else: # check intersection

1517 for c in ((c1,) if c1 is c2 else (c1, c2)):

1518 d = fabs(r3 - p3.distanceTo(c, radius=radius, wrap=wrap))

1519 if d < eps: # below margin

1520 t.append((d, c))

1521 m = min(m, d)

1522

1523 except IntersectionError as x:

1524 if _concentric_ in str(x): # XXX ConcentricError?

1525 pc += 1

1526

1527 p1, r1, p2, r2, p3, r3 = p2, r2, p3, r3, p1, r1 # rotate

1528

1529 if t: # get min, max, points and count ...

1530 t = tuple(sorted(t))

1531 n = (len(t),) # as 1-tuple

1532 # ... or for a single trilaterated result,

1533 # min *is* max, min- *is* maxPoint and n=1

1534 return Trilaterate5Tuple(t[0] + t[-1] + n) # *(t[0] + ...)

1535

1536 elif area and pc == 3: # all pairwise concentric ...

1537 r, p = min((r1, p1), (r2, p2), (r3, p3))

1538 m = max(r1, r2, r3)

1539 # ... return "smallest" point twice, the smallest

1540 # and largest distance and n=0 for concentric

1541 return Trilaterate5Tuple(float(r), p, float(m), p, 0)

1542

1543 n, f = (_overlap_, max) if area else (_intersection_, min)

1544 t = '%s (%s %.3f)' % (_no_(n), f.__name__, m)

1545 raise IntersectionError(area=area, eps=eps, wrap=wrap, txt=t)

1546

1547

1548__all__ += _ALL_DOCS(LatLonBase)

1549

1550# **) MIT License

1551#

1552# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.

1553#

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

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

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

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

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

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

1560#

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

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

1563#

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

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

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

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

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

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

1570# OTHER DEALINGS IN THE SOFTWARE.