1

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

3

4 u'''Spherical, C{N-vector}-based geodesy.

5

6 N-vector-based classes geodetic (lat-/longitude) L{LatLon}, geocentric

7(ECEF) L{Cartesian} and L{Nvector} and functions L{areaOf}, L{intersection},

8 L{meanOf}, L{nearestOn3}, L{perimeterOf}, L{sumOf}, L{triangulate} and

9 L{trilaterate}, I{all spherical}.

10

11 Pure Python implementation of n-vector-based spherical geodetic (lat-/longitude)

12 methods, transcoded from JavaScript originals by I{(C) Chris Veness 2011-2016},

13 published under the same MIT Licence**. See U{Vector-based geodesy

14<https://www.Movable-Type.co.UK/scripts/latlong-vectors.html>} and

15 U{Module latlon-nvector-spherical

16<https://www.Movable-Type.co.UK/scripts/geodesy/docs/module-latlon-nvector-spherical.html>}.

17

18 Tools for working with points and lines on (a spherical model of) the

19 earth’s surface using using n-vectors rather than the more common

20 spherical trigonometry. N-vectors make many calculations much simpler,

21 and easier to follow, compared with the trigonometric equivalents.

22

23 Based on Kenneth Gade’s U{‘Non-singular Horizontal Position Representation’

24<https://www.NavLab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf>},

25 The Journal of Navigation (2010), vol 63, nr 3, pp 395-417.

26

27 Note that the formulations below take x => 0°N,0°E, y => 0°N,90°E and

28 z => 90°N while Gade uses x => 90°N, y => 0°N,90°E, z => 0°N,0°E.

29

30 Also note that on a spherical earth model, an n-vector is equivalent

31 to a normalised version of an (ECEF) cartesian coordinate.

32'''

33# make sure int/int division yields float quosient, see .basics

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

35

36 from pygeodesy.basics import isscalar, _xinstanceof

37 from pygeodesy.constants import EPS, EPS0, PI, PI2, PI_2, R_M, \

38 _0_0, _0_5, _1_0

39# from pygeodesy.datums import Datums # from .sphericalBase

40 from pygeodesy.errors import PointsError, VectorError, _xError, _xkwds

41 from pygeodesy.fmath import fmean, fsum

42# from pygeodesy.fsums import fsum # from .fmath

43 from pygeodesy.interns import _composite_, _end_, _Nv00_, _other_, \

44 _point_, _pole_

45 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS, _ALL_OTHER

46# from pygeodesy.named import notImplemented # from .nvectorBase

47# from pygeodesy.namedTuples import NearestOn3Tuple # from .points

48 from pygeodesy.nvectorBase import LatLonNvectorBase, NorthPole, \

49 notImplemented, NvectorBase, _nsumOf, \

50 _triangulate, _trilaterate

51 from pygeodesy.points import NearestOn3Tuple, ispolar # PYCHOK exported

52 from pygeodesy.props import deprecated_function, deprecated_method

53 from pygeodesy.sphericalBase import _angular, CartesianSphericalBase, \

54 _intersecant2, LatLonSphericalBase, \

55 Datums, _r2m

56 from pygeodesy.units import Bearing, Bearing_, Radius, Scalar

57 from pygeodesy.utily import atan2, degrees360, fabs, sincos2, sincos2_, \

58 sincos2d, _unrollon, _Wrap

59

60# from math import atan2, fabs # from utily

61

62__all__ = _ALL_LAZY.sphericalNvector

63__version__ = '23.06.30'

64

65_lines_ = 'lines'

66

67

68 class Cartesian(CartesianSphericalBase):

69 '''Extended to convert geocentric, L{Cartesian} points to

70 L{Nvector} and n-vector-based, spherical L{LatLon}.

71 '''

72

73 def toLatLon(self, **LatLon_and_kwds): # PYCHOK LatLon=LatLon

74 '''Convert this cartesian to an C{Nvector}-based geodetic point.

75

76 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword

77 arguments, like C{datum}. Use C{B{LatLon}=...}

78 to override this L{LatLon} class or specify

79 C{B{LatLon}=None}.

80

81 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is set

82 to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,

83 C, M, datum)} with C{C} and C{M} if available.

84

85 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument.

86 '''

87 kwds = _xkwds(LatLon_and_kwds, LatLon=LatLon, datum=self.datum)

88 return CartesianSphericalBase.toLatLon(self, **kwds)

89

90 def toNvector(self, **Nvector_and_kwds): # PYCHOK Datums.WGS84

91 '''Convert this cartesian to L{Nvector} components, I{including height}.

92

93 @kwarg Nvector_and_kwds: Optional L{Nvector} and L{Nvector} keyword

94 arguments, like C{datum}. Use C{B{Nvector}=...}

95 to override this L{Nvector} class or specify

96 C{B{Nvector}=None}.

97

98 @return: The C{n-vector} components (L{Nvector}) or if B{C{Nvector}}

99 is set to C{None}, a L{Vector4Tuple}C{(x, y, z, h)}

100

101 @raise TypeError: Invalid B{C{Nvector_and_kwds}} argument.

102 '''

103 # ll = CartesianBase.toLatLon(self, LatLon=LatLon,

104 # datum=datum or self.datum)

105 # kwds = _xkwds(kwds, Nvector=Nvector)

106 # return ll.toNvector(**kwds)

107 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector, datum=self.datum)

108 return CartesianSphericalBase.toNvector(self, **kwds)

109

110

111 class LatLon(LatLonNvectorBase, LatLonSphericalBase):

112 '''New n-vector based point on a spherical earth model.

113

114 Tools for working with points, lines and paths on (a spherical

115 model of) the earth's surface using vector-based methods.

116

117 @example:

118

119 >>> from sphericalNvector import LatLon

120 >>> p = LatLon(52.205, 0.119)

121 '''

122 _Nv = None # cached_toNvector L{Nvector})

123

124 def _update(self, updated, *attrs, **setters): # PYCHOK args

125 '''(INTERNAL) Zap cached attributes if updated.

126 '''

127 if updated: # reset caches

128 LatLonNvectorBase._update(self, updated, _Nv=self._Nv) # special case

129 LatLonSphericalBase._update(self, updated, *attrs, **setters)

130

131 def alongTrackDistanceTo(self, start, end, radius=R_M, wrap=False):

132 '''Compute the (signed) distance from the start to the closest

133 point on the great circle line defined by a start and an

134 end point.

135

136 That is, if a perpendicular is drawn from this point to the

137 great circle line, the along-track distance is the distance

138 from the start point to the point where the perpendicular

139 crosses the line.

140

141 @arg start: Start point of great circle line (L{LatLon}).

142 @arg end: End point of great circle line (L{LatLon}) or

143 initial bearing from start point (compass

144 C{degrees360}).

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

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

147 the B{C{start}} and B{C{end}} points (C{bool}).

148

149 @return: Distance along the great circle line (C{radians}

150 if C{B{radius} is None} else C{meter}, same units

151 as B{C{radius}}), positive if "after" the start

152 toward the end point of the line or negative if

153 "before" the start point.

154

155 @raise TypeError: If B{C{start}} or B{C{end}} point is not L{LatLon}.

156

157 @raise Valuerror: Some points coincide.

158

159 @example:

160

161 >>> p = LatLon(53.2611, -0.7972)

162

163 >>> s = LatLon(53.3206, -1.7297)

164 >>> e = LatLon(53.1887, 0.1334)

165 >>> d = p.alongTrackDistanceTo(s, e) # 62331.58

166 '''

167 p = self.others(start=start)

168 n = self.toNvector()

169

170 gc, _, _ = self._gc3(p, end, _end_, wrap=wrap)

171 a = gc.cross(n).cross(gc) # along-track point gc × p × gc

172 return _r2m(start.toNvector().angleTo(a, vSign=gc), radius)

173

174 @deprecated_method

175 def bearingTo(self, other, **unused): # PYCHOK no cover

176 '''DEPRECATED, use method L{initialBearingTo}.

177 '''

178 return self.initialBearingTo(other)

179

180 def crossTrackDistanceTo(self, start, end, radius=R_M, wrap=False):

181 '''Compute the (signed) distance from this point to great circle

182 defined by a start and end point.

183

184 @arg start: Start point of great circle line (L{LatLon}).

185 @arg end: End point of great circle line (L{LatLon}) or

186 initial bearing from start point (compass

187 C{degrees360}).

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

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

190 the B{C{start}} and B{C{end}} points (C{bool}).

191

192 @return: Distance to great circle (C{radians} if C{B{radius}

193 is None} else C{meter}, same units as B{C{radius}}),

194 negative if to the left or positive if to the right

195 of the line .

196

197 @raise TypeError: If B{C{start}} or B{C{end}} point is not L{LatLon}.

198

199 @raise Valuerror: Some points coincide.

200

201 @example:

202

203 >>> p = LatLon(53.2611, -0.7972)

204

205 >>> s = LatLon(53.3206, -1.7297)

206 >>> d = p.crossTrackDistanceTo(s, 96) # -305.7

207

208 >>> e = LatLon(53.1887, 0.1334)

209 >>> d = p.crossTrackDistanceTo(s, e) # -307.5

210 '''

211 p = self.others(start=start)

212 n = self.toNvector()

213

214 gc, _, _ = self._gc3(p, end, _end_, wrap=wrap)

215 return _r2m(gc.angleTo(n) - PI_2, radius)

216

217 def destination(self, distance, bearing, radius=R_M, height=None):

218 '''Locate the destination from this point after having

219 travelled the given distance on the given bearing.

220

221 @arg distance: Distance travelled (C{meter}, same units

222 as B{C{radius}}).

223 @arg bearing: Bearing from this point (compass C{degrees360}).

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

225 @kwarg height: Optional height at destination, overriding the

226 default height (C{meter}, same units as B{C{radius}}).

227

228 @return: Destination point (L{LatLon}).

229

230 @raise Valuerror: Polar coincidence or invalid B{C{distance}},

231 B{C{bearing}}, B{C{radius}} or B{C{height}}.

232

233 @example:

234

235 >>> p = LatLon(51.4778, -0.0015)

236 >>> q = p.destination(7794, 300.7)

237 >>> q.toStr() # 51.513546°N, 000.098345°W

238 '''

239 b = Bearing_(bearing)

240 a = _angular(distance, radius, low=None)

241 sa, ca, sb, cb = sincos2_(a, b)

242

243 n = self.toNvector()

244 e = NorthPole.cross(n, raiser=_pole_).unit() # east vector at n

245 x = n.cross(e) # north vector at n

246 d = x.times(cb).plus(e.times(sb)) # direction vector @ n

247 n = n.times(ca).plus(d.times(sa))

248 return n.toLatLon(height=height, LatLon=self.classof) # Nvector(n.x, n.y, n.z).toLatLon(...)

249

250 def distanceTo(self, other, radius=R_M, wrap=False):

251 '''Compute the distance from this to an other point.

252

253 @arg other: The other point (L{LatLon}).

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

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

256 the B{C{other}} point (C{bool}).

257

258 @return: Distance between this and the B{C{other}} point

259 (C{meter}, same units as B{C{radius}} or C{radians}

260 if B{C{radius}} is C{None}).

261

262 @raise TypeError: Invalid B{C{other}} point.

263

264 @example:

265

266 >>> p = LatLon(52.205, 0.119)

267 >>> q = LatLon(48.857, 2.351);

268 >>> d = p.distanceTo(q) # 404.3 Km

269 '''

270 p = self.others(other)

271 if wrap:

272 p = _unrollon(self, p)

273 n = p.toNvector()

274 r = fabs(self.toNvector().angleTo(n, wrap=wrap))

275 return r if radius is None else (Radius(radius) * r)

276

277# @Property_RO

278# def Ecef(self):

279# '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}.

280# '''

281# return _ALL_MODS.ecef.EcefKarney

282

283 def _gc3(self, start, end, namend, raiser=_point_, wrap=False):

284 '''(INTERNAL) Return great circle, start and end Nvectors.

285 '''

286 s = start.toNvector()

287 if isscalar(end): # bearing

288 gc = s.greatCircle(end)

289 e = None

290 else: # point

291 p = self.others(end, name=namend)

292 if wrap:

293 p = _unrollon(start, p, wrap=wrap)

294 e = p.toNvector()

295 gc = s.cross(e, raiser=raiser) # XXX .unit()?

296 return gc, s, e

297

298 def greatCircle(self, bearing):

299 '''Compute the vector normal to great circle obtained by

300 heading on the given bearing from this point.

301

302 Direction of vector is such that initial bearing vector

303 b = c × n, where n is an n-vector representing this point.

304

305 @arg bearing: Bearing from this point (compass C{degrees360}).

306

307 @return: N-vector representing the great circle (L{Nvector}).

308 '''

309 t = Bearing_(bearing)

310 a, b = self.philam

311

312 sa, ca, sb, cb, st, ct = sincos2_(a, b, t)

313 return Nvector(sb * ct - sa * cb * st,

314 -cb * ct - sa * sb * st,

315 ca * st, name=self.name) # XXX .unit()

316

317 def greatCircleTo(self, other, wrap=False):

318 '''Compute the vector normal to great circle obtained by

319 heading from this to an other point or on a given bearing.

320

321 Direction of vector is such that initial bearing vector

322 b = c × n, where n is an n-vector representing this point.

323

324 @arg other: The other point (L{LatLon}) or the bearing from

325 this point (compass C{degrees360}).

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

327 the B{C{other}} point (C{bool}).

328

329 @return: N-vector representing the great circle (L{Nvector}).

330

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

332

333 @raise Valuerror: Points coincide.

334

335 @example:

336

337 >>> p = LatLon(53.3206, -1.7297)

338 >>> gc = p.greatCircle(96.0)

339 >>> gc.toStr() # (-0.79408, 0.12856, 0.59406)

340

341 >>> q = LatLon(53.1887, 0.1334)

342 >>> g = p.greatCircleTo(q)

343 >>> g.toStr() # (-0.79408, 0.12859, 0.59406)

344 '''

345 gc, _, _ = self._gc3(self, other, _other_, wrap=wrap)

346 return gc.unit()

347

348 def initialBearingTo(self, other, wrap=False, **unused): # raiser=...

349 '''Compute the initial bearing (forward azimuth) from this

350 to an other point.

351

352 @arg other: The other point (L{LatLon}).

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

354 the B{C{other}} point (C{bool}).

355

356 @return: Initial bearing (compass C{degrees360}).

357

358 @raise Crosserror: This point coincides with the B{C{other}}

359 point or the C{NorthPole}, provided

360 L{pygeodesy.crosserrors} is C{True}.

361

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

363

364 @example:

365

366 >>> p1 = LatLon(52.205, 0.119)

367 >>> p2 = LatLon(48.857, 2.351)

368 >>> b = p1.initialBearingTo(p2) # 156.2

369 '''

370 n = self.toNvector()

371 p = self.others(other)

372 if wrap:

373 p = _unrollon(self, p, wrap=wrap)

374 p = p.toNvector()

375 # see <https://MathForum.org/library/drmath/view/55417.html>

376# gc1 = self.greatCircleTo(other)

377 gc1 = n.cross(p, raiser=_point_) # .unit()

378# gc2 = self.greatCircleTo(NorthPole)

379 gc2 = n.cross(NorthPole, raiser=_pole_) # .unit()

380 return degrees360(gc1.angleTo(gc2, vSign=n))

381

382 def intermediateChordTo(self, other, fraction, height=None, wrap=False):

383 '''Locate the point projected from the point at given fraction

384 on a straight line (chord) between this and an other point.

385

386 @arg other: The other point (L{LatLon}).

387 @arg fraction: Fraction between both points (float, between

388 0.0 for this and 1.0 for the other point).

389 @kwarg height: Optional height at the intermediate point,

390 overriding the fractional height (C{meter}).

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

392 the B{C{other}} point (C{bool}).

393

394 @return: Intermediate point (L{LatLon}).

395

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

397

398 @example:

399

400 >>> p = LatLon(52.205, 0.119)

401 >>> q = LatLon(48.857, 2.351)

402 >>> i = p.intermediateChordTo(q, 0.25) # 51.3723°N, 000.7072°E

403 '''

404 n = self.toNvector()

405 p = self.others(other)

406 if wrap:

407 p = _unrollon(self, p, wrap=wrap)

408

409 f = Scalar(fraction=fraction)

410 i = p.toNvector().times(f).plus(n.times(1 - f))

411# i = p.toNvector() * f + self.toNvector() * (1 - f))

412

413 h = self._havg(other, f=f, h=height)

414 return i.toLatLon(height=h, LatLon=self.classof) # Nvector(i.x, i.y, i.z).toLatLon(...)

415

416 def intermediateTo(self, other, fraction, height=None, wrap=False):

417 '''Locate the point at a given fraction between this and an

418 other point.

419

420 @arg other: The other point (L{LatLon}).

421 @arg fraction: Fraction between both points (C{float}, between

422 0.0 for this and 1.0 for the other point).

423 @kwarg height: Optional height at the intermediate point,

424 overriding the fractional height (C{meter}).

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

426 the B{C{other}} point (C{bool}).

427

428 @return: Intermediate point (L{LatLon}).

429

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

431

432 @raise Valuerror: Points coincide or invalid B{C{height}}.

433

434 @see: Methods C{midpointTo} and C{rhumbMidpointTo}.

435

436 @example:

437

438 >>> p = LatLon(52.205, 0.119)

439 >>> q = LatLon(48.857, 2.351)

440 >>> i = p.intermediateTo(q, 0.25) # 51.3721°N, 000.7074°E

441 '''

442 n = self.toNvector()

443 p = self.others(other)

444 if wrap:

445 p = _unrollon(self, p, wrap=wrap)

446 p = p.toNvector()

447 f = Scalar(fraction=fraction)

448

449 x = n.cross(p, raiser=_point_)

450 d = x.unit().cross(n) # unit(n × p) × n

451 # angular distance α, tan(α) = |n × p| / n ⋅ p

452 s, c = sincos2(atan2(x.length, n.dot(p)) * f) # interpolated

453 i = n.times(c).plus(d.times(s)) # n * cosα + d * sinα

454

455 h = self._havg(other, f=f, h=height)

456 return i.toLatLon(height=h, LatLon=self.classof) # Nvector(i.x, i.y, i.z).toLatLon(...)

457

458 def intersection(self, end1, start2, end2, height=None, wrap=False):

459 '''Locate the intersection point of two lines each defined

460 by two points or a start point and bearing from North.

461

462 @arg end1: End point of the first line (L{LatLon}) or the

463 initial bearing at this point (compass C{degrees360}).

464 @arg start2: Start point of the second line (L{LatLon}).

465 @arg end2: End point of the second line (L{LatLon}) or the

466 initial bearing at the second point (compass

467 C{degrees}).

468 @kwarg height: Optional height at the intersection point,

469 overriding the mean height (C{meter}).

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

471 start and end points (C{bool}).

472

473 @return: The intersection point (L{LatLon}).

474

475 @raise TypeError: If B{C{start2}}, B{C{end1}} or B{C{end2}}

476 point is not L{LatLon}.

477

478 @raise ValueError: Intersection is ambiguous or infinite or

479 the lines are parallel, coincident or null.

480

481 @see: Function L{sphericalNvector.intersection} and method

482 L{intersection2}.

483

484 @example:

485

486 >>> s = LatLon(51.8853, 0.2545)

487 >>> e = LatLon(49.0034, 2.5735)

488 >>> i = s.intersection(108.55, e, 32.44) # 50.9076°N, 004.5086°E

489 '''

490 return intersection(self, end1, start2, end2, height=height,

491 wrap=wrap, LatLon=self.classof)

492

493 def intersection2(self, end1, start2, end2, height=None, wrap=False):

494 '''Locate the intersections of two (great circle) lines each defined

495 by two points or by a start point and an (initial) bearing.

496

497 @arg end1: End point of the first line (L{LatLon}) or the

498 initial bearing at this point (compass C{degrees360}).

499 @arg start2: Start point of the second line (L{LatLon}).

500 @arg end2: End point of the second line (L{LatLon}) or the

501 initial bearing at the second start point (compass

502 C{degrees360}).

503 @kwarg height: Optional height at the intersection and antipodal

504 point, overriding the mean height (C{meter}).

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

506 B{C{start2}} and both B{C{end*}} points (C{bool}).

507

508 @return: 2-Tuple C{(intersection, antipode)}, each a B{C{LatLon}}.

509

510 @raise TypeError: If B{C{start2}}, B{C{end1}} or B{C{end2}}

511 point is not L{LatLon}.

512

513 @raise ValueError: Intersection is ambiguous or infinite or

514 the lines are parallel, coincident or null.

515

516 @see: Function L{sphericalNvector.intersection2} and method

517 L{intersection}.

518 '''

519 return intersection2(self, end1, start2, end2, height=height,

520 wrap=wrap, LatLon=self.classof)

521

522 def isenclosedBy(self, points, wrap=False):

523 '''Check whether a (convex) polygon or composite encloses this point.

524

525 @arg points: The polygon points or composite (L{LatLon}[],

526 L{BooleanFHP} or L{BooleanGH}).

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

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

529

530 @return: C{True} if this point is inside the polygon or composite,

531 C{False} otherwise.

532

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

534

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

536

537 @see: Functions L{pygeodesy.isconvex}, L{pygeodesy.isenclosedBy}

538 and L{pygeodesy.ispolar} especially if the B{C{points}} may

539 enclose a pole or wrap around the earth I{longitudinally}.

540 '''

541 if _MODS.booleans.isBoolean(points):

542 return points._encloses(self.lat, self.lon, wrap=wrap)

543

544 # sum subtended angles of each edge (using n0, the

545 # normal vector to this point for sign of α)

546 def _subtangles(ps, w):

547 Ps = self.PointsIter(ps, loop=1, wrap=w)

548 n0 = self.toNvector()

549 _m0 = n0.minus

550 p1 = Ps[0]

551 vs1 = _m0(p1.toNvector())

552 for p2 in Ps.iterate(closed=True):

553 if w and not Ps.looped:

554 p2 = _unrollon(p1, p2)

555 p1 = p2

556 vs2 = _m0(p2.toNvector())

557 yield vs1.angleTo(vs2, vSign=n0) # PYCHOK false

558 vs1 = vs2

559

560 # Note, this method uses angle summation test: on a plane,

561 # angles for an enclosed point will sum to 360°, angles for

562 # an exterior point will sum to 0°. On a sphere, enclosed

563 # point angles will sum to less than 360° (due to spherical

564 # excess), exterior point angles will be small but non-zero.

565 s = fsum(_subtangles(points, wrap), floats=True) # normal vector

566 # XXX are winding number optimisations equally applicable to

567 # spherical surface?

568 return fabs(s) > PI

569

570 @deprecated_method

571 def isEnclosedBy(self, points): # PYCHOK no cover

572 '''DEPRECATED, use method C{isenclosedBy}.'''

573 return self.isenclosedBy(points)

574

575 def iswithin(self, point1, point2, wrap=False):

576 '''Check whether this point is between two other points.

577

578 If this point is not on the great circle arc defined by

579 both points, return whether it is within the area bound

580 by perpendiculars to the great circle at each point (in

581 the same hemispere).

582

583 @arg point1: Start point of the arc (L{LatLon}).

584 @arg point2: End point of the arc (L{LatLon}).

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

586 B{C{point1}} and B{C{point2}} (C{bool}).

587

588 @return: C{True} if this point is within the (great circle)

589 arc, C{False} otherwise.

590

591 @raise TypeError: If B{C{point1}} or B{C{point2}} is not

592 L{LatLon}.

593 '''

594 p1 = self.others(point1=point1)

595 p2 = self.others(point2=point2)

596 if wrap:

597 p1 = _Wrap.point(p1)

598 p2 = _unrollon(p1, p2, wrap=wrap)

599 n, n1, n2 = (_.toNvector() for _ in (self, p1, p2))

600

601 # corner case, null arc

602 if n1.isequalTo(n2):

603 return n.isequalTo(n1) or n.isequalTo(n2) # PYCHOK returns

604

605 if n.dot(n1) < 0 or n.dot(n2) < 0: # different hemisphere

606 return False # PYCHOK returns

607

608 # get vectors representing d0=p0->p1 and d2=p2->p1 and the

609 # dot product d0⋅d2 tells us if p0 is on the p2 side of p1 or

610 # on the other side (similarly for d0=p0->p2 and d1=p1->p2

611 # and dot product d0⋅d1 and p0 on the p1 side of p2 or not)

612 return n.minus(n1).dot(n2.minus(n1)) >= 0 and \

613 n.minus(n2).dot(n1.minus(n2)) >= 0

614

615 @deprecated_method

616 def isWithin(self, point1, point2): # PYCHOK no cover

617 '''DEPRECATED, use method C{iswithin}.'''

618 return self.iswithin(point1, point2)

619

620 def midpointTo(self, other, height=None, fraction=_0_5, wrap=False):

621 '''Find the midpoint between this and an other point.

622

623 @arg other: The other point (L{LatLon}).

624 @kwarg height: Optional height at the midpoint, overriding

625 the mean height (C{meter}).

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

627 may be negative or greater than 1.0.

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

629 B{C{other}} point (C{bool}).

630

631 @return: Midpoint (L{LatLon}).

632

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

634

635 @see: Methods C{intermediateTo} and C{rhumbMidpointTo}.

636

637 @example:

638

639 >>> p1 = LatLon(52.205, 0.119)

640 >>> p2 = LatLon(48.857, 2.351)

641 >>> m = p1.midpointTo(p2) # '50.5363°N, 001.2746°E'

642 '''

643 if fraction is _0_5:

644 p = self.others(other)

645 if wrap:

646 p = _unrollon(self, p, wrap=wrap)

647 m = self.toNvector().plus(p.toNvector())

648 h = self._havg(other, f=fraction, h=height)

649 r = m.toLatLon(height=h, LatLon=self.classof)

650 else:

651 r = self.intermediateTo(other, fraction, height=height, wrap=wrap)

652 return r

653

654 def nearestOn(self, point1, point2, height=None, within=True, wrap=False):

655 '''Locate the point on the great circle arc between two points

656 closest to this point.

657

658 @arg point1: Start point of the arc (L{LatLon}).

659 @arg point2: End point of the arc (L{LatLon}).

660 @kwarg height: Optional height, overriding the mean height for

661 the point within the arc (C{meter}), or C{None}

662 to interpolate the height.

663 @kwarg within: If C{True}, return the closest point between both

664 given points, otherwise the closest point

665 elsewhere on the great circle arc (C{bool}).

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

667 B{C{point1}} and B{C{point2}} (C{bool}).

668

669 @return: Closest point on the arc (L{LatLon}).

670

671 @raise NotImplementedError: Keyword argument C{B{wrap}=True}

672 not supported.

673

674 @raise TypeError: Invalid B{C{point1}} or B{C{point2}}.

675

676 @example:

677

678 >>> s1 = LatLon(51.0, 1.0)

679 >>> s2 = LatLon(51.0, 2.0)

680

681 >>> s = LatLon(51.0, 1.9)

682 >>> p = s.nearestOn(s1, s2) # 51.0004°N, 001.9000°E

683

684 >>> d = p.distanceTo(s) # 42.71 m

685

686 >>> s = LatLon(51.0, 2.1)

687 >>> p = s.nearestOn(s1, s2) # 51.0000°N, 002.0000°E

688 '''

689 p1 = self.others(point1=point1)

690 p2 = self.others(point2=point2)

691 if wrap:

692 p1 = _Wrap.point(p1)

693 p2 = _unrollon(p1, p2, wrap=wrap)

694 p0 = self

695

696 if p0.iswithin(p1, p2) and not p1.isequalTo(p2, EPS):

697 # closer to arc than to its endpoints,

698 # find the closest point on the arc

699 gc1 = p1.toNvector().cross(p2.toNvector())

700 gc2 = p0.toNvector().cross(gc1)

701 n = gc1.cross(gc2)

702

703 elif within: # for backward compatibility, XXX unwrapped

704 return point1 if (self.distanceTo(point1) <

705 self.distanceTo(point2)) else point2

706

707 else: # handle beyond arc extent by .vector3d.nearestOn

708 n1 = p1.toNvector()

709 n2 = p2.toNvector()

710 n = p0.toNvector().nearestOn(n1, n2, within=False)

711 if n is n1:

712 return p1 # is point1

713 elif n is n2:

714 return p2 # is point2 if not wrap

715

716 p = n.toLatLon(height=height or 0, LatLon=self.classof)

717 if height in (None, False): # interpolate height within extent

718 d = p1.distanceTo(p2)

719 f = (p1.distanceTo(p) / d) if d > EPS0 else _0_5

720 p.height = p1._havg(p2, f=max(_0_0, min(f, _1_0)))

721 return p

722

723 # @deprecated_method

724 def nearestOn2(self, points, **closed_radius_height): # PYCHOK no cover

725 '''DEPRECATED, use method L{sphericalNvector.LatLon.nearestOn3}.

726

727 @return: ... 2-Tuple C{(closest, distance)} of the C{closest}

728 point (L{LatLon}) on the polygon and the C{distance}

729 to that point from this point ...

730 '''

731 r = self.nearestOn3(points, **closed_radius_height)

732 return r.closest, r.distance

733

734 def nearestOn3(self, points, closed=False, radius=R_M, height=None, wrap=False):

735 '''Locate the point on a path or polygon (with great circle

736 arcs joining consecutive points) closest to this point.

737

738 The closest point is either on within the extent of any great

739 circle arc or the nearest of the arc's end points.

740

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

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

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

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

745 for a point within the arc (C{meter}).

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

747 the B{C{points}} (C{bool}).

748

749 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of

750 the C{closest} point (L{LatLon}), the C{distance}

751 between this and the C{closest} point in C{meter},

752 same units as B{C{radius}} or in C{radians} if

753 B{C{radius}} is C{None} and the C{angle} from this

754 to the C{closest} point in compass C{degrees360}.

755

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

757

758 @raise ValueError: No B{C{points}}.

759 '''

760 Ps = self.PointsIter(points, loop=1, wrap=wrap)

761 _r = self.distanceTo

762 _n = self.nearestOn

763

764 c = p1 = Ps[0]

765 r = _r(c, radius=None) # radians

766 for p2 in Ps.iterate(closed=closed):

767 if wrap and not Ps.looped:

768 p2 = _unrollon(p1, p2)

769 p = _n(p1, p2, height=height)

770 d = _r(p, radius=None) # radians

771 if d < r:

772 c, r = p, d

773 p1 = p2

774 d = r if radius is None else (Radius(radius) * r)

775 return NearestOn3Tuple(c, d, degrees360(r))

776

777 def toCartesian(self, **Cartesian_and_kwds): # PYCHOK Cartesian=Cartesian, datum=None

778 '''Convert this point to C{Nvector}-based cartesian (ECEF) coordinates.

779

780 @kwarg Cartesian_and_kwds: Optional L{Cartesian} and L{Cartesian} keyword

781 arguments, like C{datum}. Use C{B{Cartesian}=...}

782 to override this L{Cartesian} class or specify

783 C{B{Cartesian}=None}.

784

785 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}} is

786 set to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,

787 C, M, datum)} with C{C} and C{M} if available.

788

789 @raise TypeError: Invalid B{C{Cartesian_and_kwds}} argument.

790 '''

791 kwds = _xkwds(Cartesian_and_kwds, Cartesian=Cartesian, datum=self.datum)

792 return LatLonSphericalBase.toCartesian(self, **kwds)

793

794 def toNvector(self, **Nvector_and_kwds): # PYCHOK signature

795 '''Convert this point to L{Nvector} components, I{including height}.

796

797 @kwarg Nvector_and_kwds: Optional L{Nvector} and L{Nvector} keyword

798 arguments. Specify C{B{Nvector}=...} to

799 override this L{Nvector} class or use

800 C{B{Nvector}=None}.

801

802 @return: The C{n-vector} components (L{Nvector}) or if B{C{Nvector}} is

803 set to C{None}, a L{Vector4Tuple}C{(x, y, z, h)}.

804

805 @raise TypeError: Invalid B{C{Nvector_and_kwds}} argument.

806

807 @example:

808

809 >>> p = LatLon(45, 45)

810 >>> n = p.toNvector()

811 >>> n.toStr() # [0.50, 0.50, 0.70710]

812 '''

813 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector)

814 return LatLonNvectorBase.toNvector(self, **kwds)

815

816

817 class Nvector(NvectorBase):

818 '''An n-vector is a position representation using a (unit) vector

819 normal to the earth's surface. Unlike lat-/longitude points,

820 n-vectors have no singularities or discontinuities.

821

822 For many applications, n-vectors are more convenient to work

823 with than other position representations like lat-/longitude,

824 earth-centred earth-fixed (ECEF) vectors, UTM coordinates, etc.

825

826 On a spherical model earth, an n-vector is equivalent to an

827 earth-centred earth-fixed (ECEF) vector.

828

829 Note commonality with L{ellipsoidalNvector.Nvector}.

830 '''

831 _datum = Datums.Sphere # default datum (L{Datum})

832

833 def toCartesian(self, **Cartesian_and_kwds): # PYCHOK Cartesian=Cartesian

834 '''Convert this n-vector to C{Nvector}-based cartesian

835 (ECEF) coordinates.

836

837 @kwarg Cartesian_and_kwds: Optional L{Cartesian} and L{Cartesian} keyword

838 arguments, like C{h}. Use C{B{Cartesian}=...}

839 to override this L{Cartesian} class or specify

840 C{B{Cartesian}=None}.

841

842 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}} is

843 set to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,

844 C, M, datum)} with C{C} and C{M} if available.

845

846 @raise TypeError: Invalid B{C{Cartesian_and_kwds}} argument.

847 '''

848 kwds = _xkwds(Cartesian_and_kwds, h=self.h, Cartesian=Cartesian)

849 return NvectorBase.toCartesian(self, **kwds) # class or .classof

850

851 def toLatLon(self, **LatLon_and_kwds): # PYCHOK height=None, LatLon=LatLon

852 '''Convert this n-vector to an C{Nvector}-based geodetic point.

853

854 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword

855 arguments, like C{height}. Use C{B{LatLon}=...}

856 to override this L{LatLon} class or specify

857 C{B{LatLon}=None}.

858

859 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is set

860 to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,

861 C, M, datum)} with C{C} and C{M} if available.

862

863 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument.

864

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

866 '''

867 kwds = _xkwds(LatLon_and_kwds, height=self.h, LatLon=LatLon)

868 return NvectorBase.toLatLon(self, **kwds) # class or .classof

869

870 def greatCircle(self, bearing):

871 '''Compute the n-vector normal to great circle obtained by

872 heading on given compass bearing from this point as its

873 n-vector.

874

875 Direction of vector is such that initial bearing vector

876 b = c × p.

877

878 @arg bearing: Initial compass bearing (C{degrees}).

879

880 @return: N-vector representing great circle (L{Nvector}).

881

882 @raise Valuerror: Polar coincidence.

883

884 @example:

885

886 >>> n = LatLon(53.3206, -1.7297).toNvector()

887 >>> gc = n.greatCircle(96.0) # [-0.794, 0.129, 0.594]

888 '''

889 s, c = sincos2d(Bearing(bearing))

890

891 e = NorthPole.cross(self, raiser=_pole_) # easting

892 n = self.cross(e, raiser=_point_) # northing

893

894 e = e.times(c / e.length)

895 n = n.times(s / n.length)

896 return n.minus(e)

897

898

899_Nvll = LatLon(_0_0, _0_0, name=_Nv00_) # reference instance (L{LatLon})

900

901

902 def areaOf(points, radius=R_M, wrap=False):

903 '''Calculate the area of a (spherical) polygon or composite (with

904 great circle arcs joining consecutive points).

905

906 @arg points: The polygon points or clips (C{LatLon}[],

907 L{BooleanFHP} or L{BooleanGH}).

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

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

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

911

912 @return: Polygon area (C{meter} I{squared} , same units as

913 B{C{radius}}, or C{radians} if B{C{radius}} is C{None}).

914

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

916

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

918

919 @see: Functions L{pygeodesy.areaOf}, L{sphericalTrigonometry.areaOf}

920 and L{ellipsoidalKarney.areaOf}.

921

922 @example:

923

924 >>> b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1)

925 >>> areaOf(b) # 8666058750.718977

926 '''

927 def _interangles(ps, w): # like .karney._polygon

928 Ps = _Nvll.PointsIter(ps, loop=2, wrap=w)

929 # use vector to 1st point as plane normal for sign of α

930 n0 = Ps[0].toNvector()

931

932 v2 = Ps[0]._N_vector # XXX v2 == no?

933 p1 = Ps[1]

934 v1 = p1._N_vector

935 gc = v2.cross(v1)

936 for p2 in Ps.iterate(closed=True):

937 if w and not Ps.looped:

938 p2 = _unrollon(p1, p2)

939 p1 = p2

940 v2 = p2._N_vector

941 gc1 = v1.cross(v2)

942 v1 = v2

943 yield gc.angleTo(gc1, vSign=n0)

944 gc = gc1

945

946 if _MODS.booleans.isBoolean(points):

947 r = points._sum2(LatLon, areaOf, radius=None, wrap=wrap)

948 else:

949 # sum interior angles: depending on whether polygon is cw or ccw,

950 # angle between edges is π−α or π+α, where α is angle between

951 # great-circle vectors; so sum α, then take n·π − |Σα| (cannot

952 # use Σ(π−|α|) as concave polygons would fail)

953 s = fsum(_interangles(points, wrap), floats=True)

954 # using Girard’s theorem: A = [Σθᵢ − (n−2)·π]·R²

955 # (PI2 - abs(s) == (n*PI - abs(s)) - (n-2)*PI)

956 r = fabs(PI2 - fabs(s))

957 return r if radius is None else (r * Radius(radius)**2)

958

959

960 def intersecant2(center, circle, point, bearing, radius=R_M, exact=False,

961 height=None, wrap=False): # was=True

962 '''Compute the intersections of a circle and a line.

963

964 @arg center: Center of the circle (L{LatLon}).

965 @arg circle: Radius of the circle (C{meter}, same units as B{C{radius}})

966 or a point on the circle (L{LatLon}).

967 @arg point: A point in- or outside the circle (L{LatLon}).

968 @arg bearing: Bearing at the B{C{point}} (compass C{degrees360}) or

969 a second point on the line (L{LatLon}).

970 @kwarg radius: Mean earth radius (C{meter}, conventionally).

971 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth,

972 destination and distance, if C{False} use the basic

973 rhumb methods (C{bool}) or if C{None} use the I{great

974 circle} methods.

975 @kwarg height: Optional height for the intersection points (C{meter},

976 conventionally) or C{None}.

977 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{point}}

978 and the B{C{circle}} and B{C{bearing}} points (C{bool}).

979

980 @return: 2-Tuple of the intersection points (representing a chord), each

981 an instance of this class. For a tangent line, both points are

982 the same instance, the B{C{point}} or wrapped or I{normalized}.

983

984 @raise IntersectionError: The circle and line do not intersect.

985

986 @raise TypeError: If B{C{center}} or B{C{point}} not L{LatLon} or

987 B{C{circle}} or B{C{bearing}} invalid.

988

989 @raise ValueError: Invalid B{C{circle}}, B{C{bearing}}, B{C{radius}},

990 B{C{exact}} or B{C{height}}.

991 '''

992 c = _Nvll.others(center=center)

993 p = _Nvll.others(point=point)

994 try:

995 return _intersecant2(c, circle, p, bearing, radius=radius, exact=exact,

996 height=height, wrap=wrap)

997 except (TypeError, ValueError) as x:

998 raise _xError(x, center=center, circle=circle, point=point, bearing=bearing,

999 exact=exact, wrap=wrap)

1000

1001

1002 def intersection(start1, end1, start2, end2, height=None, wrap=False,

1003 LatLon=LatLon, **LatLon_kwds):

1004 '''Locate the intersections of two (great circle) lines each defined

1005 by two points or by a start point and an (initial) bearing.

1006

1007 @arg start1: Start point of the first line (L{LatLon}).

1008 @arg end1: End point of the first line (L{LatLon}) or the initial

1009 bearing at the first start point (compass C{degrees360}).

1010 @arg start2: Start point of the second line (L{LatLon}).

1011 @arg end2: End point of the second line (L{LatLon}) or the initial

1012 bearing at the second start point (compass C{degrees360}).

1013 @kwarg height: Optional height at the intersection point,

1014 overriding the mean height (C{meter}).

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

1016 and both B{C{end*}} points (C{bool}).

1017 @kwarg LatLon: Optional class to return the intersection point

1018 (L{LatLon}).

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

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

1021

1022 @return: The intersection point (B{C{LatLon}}) or if C{B{LatLon}

1023 is None}, a cartesian L{Ecef9Tuple}C{(x, y, z, lat, lon,

1024 height, C, M, datum)} with C{C} and C{M} if available.

1025

1026 @raise TypeError: If B{C{start*}} or B{C{end*}} is not L{LatLon}.

1027

1028 @raise ValueError: Intersection is ambiguous or infinite or

1029 the lines are parallel, coincident or null.

1030

1031 @see: Function L{sphericalNvector.intersection2}.

1032

1033 @example:

1034

1035 >>> p = LatLon(51.8853, 0.2545)

1036 >>> q = LatLon(49.0034, 2.5735)

1037 >>> i = intersection(p, 108.55, q, 32.44) # 50.9076°N, 004.5086°E

1038 '''

1039 i, _, h = _intersect3(start1, end1, start2, end2, height, wrap)

1040 kwds = _xkwds(LatLon_kwds, height=h, LatLon=LatLon)

1041 return i.toLatLon(**kwds)

1042

1043

1044 def intersection2(start1, end1, start2, end2, height=None, wrap=False,

1045 LatLon=LatLon, **LatLon_kwds):

1046 '''Locate the intersections of two (great circle) lines each defined

1047 by two points or by a start point and an (initial) bearing.

1048

1049 @arg start1: Start point of the first line (L{LatLon}).

1050 @arg end1: End point of the first line (L{LatLon}) or the

1051 initial bearing at the first start point

1052 (compass C{degrees360}).

1053 @arg start2: Start point of the second line (L{LatLon}).

1054 @arg end2: End point of the second line (L{LatLon}) or the

1055 initial bearing at the second start point

1056 (compass C{degrees360}).

1057 @kwarg height: Optional height at the intersection and antipodal

1058 point, overriding the mean height (C{meter}).

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

1060 and both B{C{end*}} points (C{bool}).

1061 @kwarg LatLon: Optional class to return the intersection and

1062 antipodal points (L{LatLon}).

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

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

1065

1066 @return: 2-Tuple C{(intersection, antipode)}, each a (B{C{LatLon}})

1067 or if C{B{LatLon} is None}, a cartesian L{Ecef9Tuple}C{(x,

1068 y, z, lat, lon, height, C, M, datum)} with C{C} and C{M}

1069 if available.

1070

1071 @raise TypeError: If B{C{start*}} or B{C{end*}} is not L{LatLon}.

1072

1073 @raise ValueError: Intersection is ambiguous or infinite or

1074 the lines are parallel, coincident or null.

1075

1076 @see: Function L{sphericalNvector.intersection}.

1077 '''

1078 i, a, h = _intersect3(start1, end1, start2, end2, height, wrap)

1079 kwds = _xkwds(LatLon_kwds, height=h, LatLon=LatLon)

1080 return i.toLatLon(**kwds), a.toLatLon(**kwds)

1081

1082

1083 def _intersect3(start1, end1, start2, end2, height, wrap):

1084 '''(INTERNAL) Return the intersection and antipodal points for

1085 functions C{intersection} and C{intersection2}.

1086 '''

1087 p1 = _Nvll.others(start1=start1)

1088 p2 = _Nvll.others(start2=start2)

1089 if wrap:

1090 p2 = _unrollon(p1, p2, wrap=wrap)

1091 # If gc1 and gc2 are great circles through start and end points

1092 # (or defined by start point and bearing), then the candidate

1093 # intersections are simply gc1 × gc2 and gc2 × gc1. Most of the

1094 # work is deciding the correct intersection point to select! If

1095 # bearing is given, that determines the intersection, but if both

1096 # lines are defined by start/end points, take closer intersection.

1097 gc1, s1, e1 = _Nvll._gc3(p1, end1, 'end1', wrap=wrap)

1098 gc2, s2, e2 = _Nvll._gc3(p2, end2, 'end2', wrap=wrap)

1099

1100 hs = start1.height, start2.height

1101 # there are two (antipodal) candidate intersection

1102 # points ... we have to choose the one to return

1103 i1 = gc1.cross(gc2, raiser=_lines_)

1104 i2 = gc2.cross(gc1, raiser=_lines_)

1105

1106 # selection of intersection point depends on how

1107 # lines are defined (by bearings or endpoints)

1108 if e1 and e2: # endpoint+endpoint

1109 d = sumOf((s1, s2, e1, e2)).dot(i1)

1110 hs += end1.height, end2.height

1111 elif e1 and not e2: # endpoint+bearing

1112 # gc2 x v2 . i1 +ve means v2 bearing points to i1

1113 d = gc2.cross(s2).dot(i1)

1114 hs += end1.height,

1115 elif e2 and not e1: # bearing+endpoint

1116 # gc1 x v1 . i1 +ve means v1 bearing points to i1

1117 d = gc1.cross(s1).dot(i1)

1118 hs += end2.height,

1119 else: # bearing+bearing

1120 # if gc x v . i1 is +ve, initial bearing is

1121 # towards i1, otherwise towards antipodal i2

1122 d1 = gc1.cross(s1).dot(i1) # +ve means p1 bearing points to i1

1123 d2 = gc2.cross(s2).dot(i1) # +ve means p2 bearing points to i1

1124 if d1 > 0 and d2 > 0:

1125 d = 1 # both point to i1

1126 elif d1 < 0 and d2 < 0:

1127 d = -1 # both point to i2

1128 else: # d1, d2 opposite signs

1129 # intersection is at further-away intersection point,

1130 # take opposite intersection from mid- point of v1

1131 # and v2 [is this always true?] XXX changed to always

1132 # get intersection p1 bearing points to, aka being

1133 # located "after" p1 along the bearing at p1, like

1134 # function .sphericalTrigonometry._intersect and

1135 # .ellipsoidalBaseDI._intersect3

1136 d = d1 # neg(s1.plus(s2).dot(i1))

1137

1138 h = fmean(hs) if height is None else height

1139 return (i1, i2, h) if d > 0 else (i2, i1, h)

1140

1141

1142 def meanOf(points, height=None, wrap=False, LatLon=LatLon, **LatLon_kwds):

1143 '''Compute the I{geographic} mean of the supplied points.

1144

1145 @arg points: Array of points to be averaged (L{LatLon}[]).

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

1147 (C{meter}).

1148 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{points}} (C{bool}).

1149 @kwarg LatLon: Optional class to return the mean point (L{LatLon}).

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

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

1152

1153 @return: Point at geographic mean and mean height (B{C{LatLon}}).

1154

1155 @raise PointsError: Insufficient number of B{C{points}} or

1156 some B{C{points}} are not C{LatLon}.

1157 '''

1158 def _N_vs(ps, w):

1159 Ps = _Nvll.PointsIter(ps, wrap=w)

1160 for p in Ps.iterate(closed=False):

1161 yield p._N_vector

1162

1163 try:

1164 # geographic mean

1165 n = _nsumOf(_N_vs(points, wrap), height, Nvector, {})

1166 except (TypeError, ValueError) as x:

1167 raise PointsError(points=points, wrap=wrap, LatLon=LatLon, cause=x)

1168 return n.toLatLon(**_xkwds(LatLon_kwds, LatLon=LatLon, height=n.h,

1169 name=meanOf.__name__))

1170

1171

1172@deprecated_function

1173 def nearestOn2(point, points, **closed_radius_height): # PYCHOK no cover

1174 '''DEPRECATED, use method L{sphericalNvector.nearestOn3}.

1175

1176 @return: ... 2-Tuple C{(closest, distance)} of the C{closest}

1177 point (L{LatLon}) on the polygon and the C{distance}

1178 between the C{closest} and the given B{C{point}} ...

1179 '''

1180 r = nearestOn3(point, points, **closed_radius_height)

1181 return r.closest, r.distance

1182

1183

1184 def nearestOn3(point, points, closed=False, radius=R_M, height=None, wrap=False):

1185 '''Locate the point on a polygon (with great circle arcs joining

1186 consecutive points) closest to an other point.

1187

1188 If the given point is between the end points of a great circle

1189 arc, the closest point is on that arc. Otherwise, the closest

1190 point is the nearest of the arc's end points.

1191

1192 @arg point: The other, reference point (L{LatLon}).

1193 @arg points: The polygon points (L{LatLon}[]).

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

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

1196 @kwarg height: Optional height, overriding the mean height for

1197 a point within the (great circle) arc (C{meter}).

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

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

1200

1201 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of

1202 the C{closest} point (L{LatLon}) on the polygon, the

1203 C{distance} and the C{angle} between the C{closest}

1204 and the given B{C{point}}. The C{distance} is in

1205 C{meter}, same units as B{C{radius}} or in C{radians}

1206 if B{C{radius}} is C{None}, the C{angle} is in compass

1207 C{degrees360}.

1208

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

1210

1211 @raise TypeError: Some B{C{points}} or B{C{point}} not C{LatLon}.

1212 '''

1213 _xinstanceof(LatLon, point=point)

1214

1215 return point.nearestOn3(points, closed=closed, radius=radius,

1216 height=height, wrap=wrap)

1217

1218

1219 def perimeterOf(points, closed=False, radius=R_M, wrap=False):

1220 '''Compute the perimeter of a (spherical) polygon or composite (with

1221 great circle arcs joining consecutive points).

1222

1223 @arg points: The polygon points (L{LatLon}[]).

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

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

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

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

1228

1229 @return: Polygon perimeter (C{meter}, same units as B{C{radius}}

1230 or C{radians} if B{C{radius}} is C{None}).

1231

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

1233

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

1235

1236 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with

1237 C{B{points}} a composite.

1238

1239 @see: Functions L{pygeodesy.perimeterOf}, L{ellipsoidalKarney.perimeterOf}

1240 and L{sphericalTrigonometry.perimeterOf}.

1241 '''

1242 def _rads(ps, c, w): # angular edge lengths in radians

1243 Ps = _Nvll.PointsIter(ps, loop=1, wrap=w)

1244 p1 = Ps[0]

1245 v1 = p1._N_vector

1246 for p2 in Ps.iterate(closed=c):

1247 if w and not (c and Ps.looped):

1248 p2 = _unrollon(p1, p2)

1249 p1 = p2

1250 v2 = p2._N_vector

1251 yield v1.angleTo(v2)

1252 v1 = v2

1253

1254 if _MODS.booleans.isBoolean(points):

1255 if not closed:

1256 notImplemented(None, closed=closed, points=_composite_)

1257 r = points._sum2(LatLon, perimeterOf, closed=True, radius=None, wrap=wrap)

1258 else:

1259 r = fsum(_rads(points, closed, wrap), floats=True)

1260 return r if radius is None else (Radius(radius) * r)

1261

1262

1263 def sumOf(nvectors, Vector=Nvector, h=None, **Vector_kwds):

1264 '''Return the I{vectorial} sum of two or more n-vectors.

1265

1266 @arg nvectors: Vectors to be added (L{Nvector}[]).

1267 @kwarg Vector: Optional class for the vectorial sum (L{Nvector}).

1268 @kwarg h: Optional height, overriding the mean height (C{meter}).

1269 @kwarg Vector_kwds: Optional, additional B{C{Vector}} keyword arguments.

1270

1271 @return: Vectorial sum (B{C{Vector}}).

1272

1273 @raise VectorError: No B{C{nvectors}}.

1274 '''

1275 try:

1276 return _nsumOf(nvectors, h, Vector, Vector_kwds)

1277 except (TypeError, ValueError) as x:

1278 raise VectorError(nvectors=nvectors, Vector=Vector, cause=x)

1279

1280

1281 def triangulate(point1, bearing1, point2, bearing2,

1282 height=None, wrap=False,

1283 LatLon=LatLon, **LatLon_kwds):

1284 '''Locate a point given two known points and the initial bearings

1285 from those points.

1286

1287 @arg point1: First reference point (L{LatLon}).

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

1289 @arg point2: Second reference point (L{LatLon}).

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

1291 @kwarg height: Optional height at the triangulated point, overriding

1292 the mean height (C{meter}).

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

1294 (C{bool}).

1295 @kwarg LatLon: Optional class to return the triangulated point (L{LatLon}).

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

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

1298

1299 @return: Triangulated point (B{C{LatLon}}).

1300

1301 @raise TypeError: If B{C{point1}} or B{C{point2}} is not L{LatLon}.

1302

1303 @raise Valuerror: Points coincide.

1304

1305 @example:

1306

1307 >>> p = LatLon("47°18.228'N","002°34.326'W") # Basse Castouillet

1308 >>> q = LatLon("47°18.664'N","002°31.717'W") # Basse Hergo

1309 >>> t = triangulate(p, 7, q, 295) # 47.323667°N, 002.568501°W'

1310 '''

1311 return _triangulate(_Nvll.others(point1=point1), bearing1,

1312 _Nvll.others(point2=point2), bearing2,

1313 height=height, wrap=wrap,

1314 LatLon=LatLon, **LatLon_kwds)

1315

1316

1317 def trilaterate(point1, distance1, point2, distance2, point3, distance3, # PYCHOK args

1318 radius=R_M, height=None, useZ=False, wrap=False,

1319 LatLon=LatLon, **LatLon_kwds):

1320 '''Locate a point at given distances from three other points.

1321

1322 @arg point1: First point (L{LatLon}).

1323 @arg distance1: Distance to the first point (C{meter}, same units

1324 as B{C{radius}}).

1325 @arg point2: Second point (L{LatLon}).

1326 @arg distance2: Distance to the second point (C{meter}, same units

1327 as B{C{radius}}).

1328 @arg point3: Third point (L{LatLon}).

1329 @arg distance3: Distance to the third point (C{meter}, same units

1330 as B{C{radius}}).

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

1332 @kwarg height: Optional height at the trilaterated point, overriding

1333 the IDW height (C{meter}, same units as B{C{radius}}).

1334 @kwarg useZ: Include Z component iff non-NaN, non-zero (C{bool}).

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

1336 and B{C{point3}} (C{bool}).

1337 @kwarg LatLon: Optional class to return the trilaterated point (L{LatLon}).

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

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

1340

1341 @return: Trilaterated point (B{C{LatLon}}).

1342

1343 @raise IntersectionError: No intersection, trilateration failed.

1344

1345 @raise TypeError: Invalid B{C{point1}}, B{C{point2}} or B{C{point3}}.

1346

1347 @raise ValueError: Coincident B{C{points}} or invalid B{C{distance1}},

1348 B{C{distance2}}, B{C{distance3}} or B{C{radius}}.

1349

1350 @see: U{Trilateration<https://WikiPedia.org/wiki/Trilateration>}.

1351 '''

1352 return _trilaterate(_Nvll.others(point1=point1), distance1,

1353 _Nvll.others(point2=point2), distance2,

1354 _Nvll.others(point3=point3), distance3,

1355 radius=radius, height=height, useZ=useZ,

1356 wrap=wrap, LatLon=LatLon, **LatLon_kwds)

1357

1358

1359__all__ += _ALL_OTHER(Cartesian, LatLon, Nvector, # classes

1360 areaOf, # functions

1361 intersecant2, intersection, intersection2, ispolar,

1362 meanOf,

1363 nearestOn2, nearestOn3,

1364 perimeterOf,

1365 sumOf,

1366 triangulate, trilaterate)

1367

1368# **) MIT License

1369#

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

1371#

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

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

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

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

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

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

1378#

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

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

1381#

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

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

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

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

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

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

1388# OTHER DEALINGS IN THE SOFTWARE.