1

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

3

4 u'''(INTERNAL) Private, ellipsoidal Direct/Inverse geodesy base

5 class C{LatLonEllipsoidalBaseDI} and functions.

6'''

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

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

9

10 from pygeodesy.basics import isscalar, issubclassof

11 from pygeodesy.constants import EPS, MAX, PI, PI2, PI_4, isnear0, isnear1, \

12 _0_0, _0_5, _1_5, _3_0

13 from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase, Property_RO, \

14 property_RO, _TOL_M

15 from pygeodesy.errors import _AssertionError, IntersectionError, _IsnotError, \

16 _or, _ValueError, _xellipsoidal, _xError, _xkwds_not

17 from pygeodesy.fmath import favg, fmean_

18 from pygeodesy.fsums import Fmt, fsum_

19 from pygeodesy.formy import opposing, _radical2

20 from pygeodesy.interns import _antipodal_, _concentric_, _exceed_PI_radians_, \

21 _near_, _SPACE_, _too_

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

23 from pygeodesy.namedTuples import Bearing2Tuple, Destination2Tuple, \

24 Intersection3Tuple, NearestOn2Tuple, \

25 NearestOn8Tuple, _LL4Tuple

26# from pygeodesy.props import Property_RO, property_RO # from .ellipsoidalBase

27# from pygeodesy.streprs import Fmt # from .fsums

28 from pygeodesy.units import _fi_j2, Height, Radius_, Scalar

29 from pygeodesy.utily import m2km, unroll180, _unrollon, wrap90, wrap180, wrap360

30

31 from math import degrees, radians

32

33__all__ = _ALL_LAZY.ellipsoidalBaseDI

34__version__ = '23.03.22'

35

36_polar__ = 'polar?'

37_too_low_ = _too_('low')

38_B2END = _1_5 # _intersect3 bearing to pseudo-end point factor

39_TRIPS = 33 # _intersect3, _intersects2, _nearestOn interations, 6..9 sufficient?

40

41

42 class LatLonEllipsoidalBaseDI(LatLonEllipsoidalBase):

43 '''(INTERNAL) Base class for C{ellipsoidal*.LatLon} classes

44 with I{overloaded} C{Direct} and C{Inverse} methods.

45 '''

46

47 def bearingTo2(self, other, wrap=False):

48 '''Compute the initial and final bearing (forward and reverse

49 azimuth) from this to an other point, using this C{Inverse}

50 method. See methods L{initialBearingTo} and L{finalBearingTo}

51 for more details.

52

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

54 @kwarg wrap: Wrap and unroll longitudes (C{bool}).

55

56 @return: A L{Bearing2Tuple}C{(initial, final)}.

57

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

59

60 @raise ValueError: If this and the B{C{other}} point's L{Datum}

61 ellipsoids are not compatible.

62 '''

63 r = self._Inverse(other, wrap)

64 return Bearing2Tuple(r.initial, r.final, name=self.name)

65

66 def destination(self, distance, bearing, height=None):

67 '''Compute the destination point after having travelled for

68 the given distance from this point along a geodesic given

69 by an initial bearing, using this C{Direct} method. See

70 method L{destination2} for more details.

71

72 @arg distance: Distance (C{meter}).

73 @arg bearing: Initial bearing in (compass C{degrees360}).

74 @kwarg height: Optional height, overriding the default

75 height (C{meter}, same units as C{distance}).

76

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

78 '''

79 return self._Direct(distance, bearing, self.classof, height).destination

80

81 def destination2(self, distance, bearing, height=None):

82 '''Compute the destination point and the final bearing (reverse

83 azimuth) after having travelled for the given distance from

84 this point along a geodesic given by an initial bearing,

85 using this C{Direct} method.

86

87 The distance must be in the same units as this point's datum

88 axes, conventionally C{meter}. The distance is measured on

89 the surface of the ellipsoid, ignoring this point's height.

90

91 The initial and final bearing (forward and reverse azimuth)

92 are in compass C{degrees360}.

93

94 The destination point's height and datum are set to this

95 point's height and datum, unless the former is overridden.

96

97 @arg distance: Distance (C{meter}).

98 @arg bearing: Initial bearing (compass C{degrees360}).

99 @kwarg height: Optional height, overriding the default

100 height (C{meter}, same units as C{distance}).

101

102 @return: A L{Destination2Tuple}C{(destination, final)}.

103 '''

104 r = self._Direct(distance, bearing, self.classof, height)

105 return self._xnamed(r)

106

107 def _Direct(self, distance, bearing, LL, height): # overloaded by I{Vincenty}

108 '''(INTERNAL) I{Karney}'s C{Direct} method.

109

110 @return: A L{Destination2Tuple}C{(destination, final)} or

111 a L{Destination3Tuple}C{(lat, lon, final)} if

112 B{C{LL}} is C{None}.

113 '''

114 g = self.geodesic

115 r = g.Direct3(self.lat, self.lon, bearing, distance)

116 if LL:

117 r = self._Direct2Tuple(LL, height, r)

118 return r

119

120 def _Direct2Tuple(self, LL, height, r):

121 '''(INTERNAL) Helper for C{._Direct} result L{Destination2Tuple}.

122 '''

123 h = self.height if height is None else height

124 d = LL(wrap90(r.lat), wrap180(r.lon), height=h, datum=self.datum, name=self.name,

125 **_xkwds_not(None, epoch=self.epoch, reframe=self.reframe))

126 return Destination2Tuple(d, wrap360(r.final))

127

128 def distanceTo(self, other, wrap=False, **unused): # ignore radius=R_M

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

130 along a geodesic, using this C{Inverse} method. See method

131 L{distanceTo3} for more details.

132

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

134 @kwarg wrap: Wrap and unroll longitudes (C{bool}).

135

136 @return: Distance (C{meter}).

137

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

139

140 @raise ValueError: If this and the B{C{other}} point's L{Datum}

141 ellipsoids are not compatible.

142 '''

143 return self._Inverse(other, wrap, azis=False).distance

144

145 def distanceTo3(self, other, wrap=False):

146 '''Compute the distance, the initial and final bearing along

147 a geodesic between this and an other point, using this

148 C{Inverse} method.

149

150 The distance is in the same units as this point's datum axes,

151 conventionally meter. The distance is measured on the surface

152 of the ellipsoid, ignoring this point's height.

153

154 The initial and final bearing (forward and reverse azimuth)

155 are in compass C{degrees360} from North.

156

157 @arg other: Destination point (C{LatLon}).

158 @kwarg wrap: Wrap and unroll longitudes (C{bool}).

159

160 @return: A L{Distance3Tuple}C{(distance, initial, final)}.

161

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

163

164 @raise ValueError: If this and the B{C{other}} point's L{Datum}

165 ellipsoids are not compatible.

166 '''

167 return self._xnamed(self._Inverse(other, wrap))

168

169 def finalBearingOn(self, distance, bearing):

170 '''Compute the final bearing (reverse azimuth) after having

171 travelled for the given distance along a geodesic given by

172 an initial bearing from this point, using this C{Direct}

173 method. See method L{destination2} for more details.

174

175 @arg distance: Distance (C{meter}).

176 @arg bearing: Initial bearing (compass C{degrees360}).

177

178 @return: Final bearing (compass C{degrees360}).

179 '''

180 return self._Direct(distance, bearing, None, None).final

181

182 def finalBearingTo(self, other, wrap=False):

183 '''Compute the final bearing (reverse azimuth) after having

184 travelled along a geodesic from this point to an other

185 point, using this C{Inverse} method. See method

186 L{distanceTo3} for more details.

187

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

189 @kwarg wrap: Wrap and unroll longitudes (C{bool}).

190

191 @return: Final bearing (compass C{degrees360}).

192

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

194

195 @raise ValueError: If this and the B{C{other}} point's L{Datum}

196 ellipsoids are not compatible.

197 '''

198 return self._Inverse(other, wrap).final

199

200 @Property_RO

201 def geodesic(self): # overloaded by I{Karney}'s, N/A for I{Vincenty}

202 '''N/A, invalid (C{None} I{always}).

203 '''

204 return None # PYCHOK no cover

205

206 def initialBearingTo(self, other, wrap=False):

207 '''Compute the initial bearing (forward azimuth) to travel

208 along a geodesic from this point to an other point,

209 using this C{Inverse} method. See method L{distanceTo3}

210 for more details.

211

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

213 @kwarg wrap: Wrap and unroll longitudes (C{bool}).

214

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

216

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

218

219 @raise ValueError: If this and the B{C{other}} point's L{Datum}

220 ellipsoids are not compatible.

221 '''

222 return self._Inverse(other, wrap).initial

223

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

225 '''Return the point at given fraction along the geodesic between

226 this and an other point, using this C{Direct} and C{Inverse}

227 methods.

228

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

230 @arg fraction: Fraction between both points (C{scalar},

231 0.0 at this and 1.0 at the other point.

232 @kwarg height: Optional height, overriding the fractional

233 height (C{meter}).

234 @kwarg wrap: Wrap and unroll longitudes (C{bool}).

235

236 @return: Intermediate point (C{LatLon}).

237

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

239

240 @raise UnitError: Invalid B{C{fraction}} or B{C{height}}.

241

242 @raise ValueError: If this and the B{C{other}} point's L{Datum}

243 ellipsoids are not compatible.

244

245 @see: Methods L{distanceTo3}, L{destination}, C{midpointTo} and

246 C{rhumbMidpointTo}.

247 '''

248 f = Scalar(fraction=fraction)

249 if isnear0(f):

250 r = self

251 elif isnear1(f) and not wrap:

252 r = self.others(other)

253 else: # negative fraction OK

254 t = self.distanceTo3(other, wrap=wrap)

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

256 r = self.destination(t.distance * f, t.initial, height=h)

257 return r

258

259 def _Inverse(self, other, wrap, **unused): # azis=False, overloaded by I{Vincenty}

260 '''(INTERNAL) I{Karney}'s C{Inverse} method.

261

262 @return: A L{Distance3Tuple}C{(distance, initial, final)}.

263

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

265

266 @raise ValueError: If this and the B{C{other}} point's

267 L{Datum} ellipsoids are not compatible.

268 '''

269 _ = self.ellipsoids(other)

270 g = self.geodesic

271 _, lon = unroll180(self.lon, other.lon, wrap=wrap)

272 return g.Inverse3(self.lat, self.lon, other.lat, lon)

273

274 def nearestOn8(self, points, closed=False, height=None, wrap=False,

275 equidistant=None, tol=_TOL_M):

276 '''Iteratively locate the point on a path or polygon closest

277 to this point.

278

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

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

281 @kwarg height: Optional height, overriding the height of this and

282 all other points (C{meter}, conventionally). If

283 B{C{height}} is C{None}, the height of each point

284 is taken into account for distances.

285

286 @return: A L{NearestOn8Tuple}C{(closest, distance, fi, j, start,

287 end, initial, final)} with C{distance} in C{meter},

288 conventionally and with the C{closest}, the C{start}

289 the C{end} point each an instance of this C{LatLon}.

290

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

292

293 @raise TypeError: Some B{C{points}} or B{C{equidistant}} invalid.

294

295 @raise ValueError: Some B{C{points}}' datum or ellipsoid incompatible

296 or no convergence for the given B{C{tol}}.

297

298 @see: Function L{pygeodesy.nearestOn6} and method C{nearestOn6}.

299 '''

300 D3 = self.distanceTo3 # Distance3Tuple

301

302 try:

303 Ps = self.PointsIter(points, loop=1)

304 p1 = c = s = e = Ps[0]

305 _ = self.ellipsoids(p1)

306 c3 = D3(c, wrap=wrap) # XXX wrap=False?

307

308 except (TypeError, ValueError) as x:

309 raise _xError(x, Fmt.SQUARE(points=0), p1, this=self, tol=tol,

310 closed=closed, height=height, wrap=wrap)

311

312 # get the azimuthal equidistant projection, once

313 A = _Equidistant00(equidistant, c)

314 b = _Box(c, c3.distance)

315 m = f = i = 0 # p1..p2 == points[i]..[j]

316

317 kwds = dict(within=True, height=height, tol=tol,

318 LatLon=self.classof, # this LatLon

319 datum=self.datum, epoch=self.epoch, reframe=self.reframe)

320 try:

321 for j, p2 in Ps.enumerate(closed=closed):

322 if wrap and j != 0:

323 p2 = _unrollon(p1, p2)

324 # skip edge if no overlap with box around closest

325 if j < 4 or b.overlaps(p1.lat, p1.lon, p2.lat, p2.lon):

326 p, t, _ = _nearestOn3_(self, p1, p2, A, **kwds)

327 d3 = D3(p, wrap=False) # already unrolled

328 if d3.distance < c3.distance:

329 c3, c, s, e, f = d3, p, p1, p2, (i + t)

330 b = _Box(c, c3.distance)

331 m = max(m, c.iteration)

332 p1, i = p2, j

333

334 except (TypeError, ValueError) as x:

335 raise _xError(x, Fmt.SQUARE(points=i), p1,

336 Fmt.SQUARE(points=j), p2, this=self, tol=tol,

337 closed=closed, height=height, wrap=wrap)

338

339 f, j = _fi_j2(f, len(Ps)) # like .vector3d.nearestOn6

340

341 n = self.nearestOn8.__name__

342 c.rename(n)

343 if s is not c:

344 s = s.copy(name=n)

345 if e is not c:

346 e = e.copy(name=n)

347 return NearestOn8Tuple(c, c3.distance, f, j, s, e, c3.initial, c3.final,

348 iteration=m) # ._iteration for tests

349

350

351 class _Box(object):

352 '''Bounding box around a C{LatLon} point.

353

354 @see: Function C{_box4} in .clipy.py.

355 '''

356 def __init__(self, center, distance):

357 '''New L{_Box} around point.

358

359 @arg center: The center point (C{LatLon}).

360 @arg distance: Radius, half-size of the box

361 (C{meter}, conventionally)

362 '''

363 E = center.ellipsoid()

364 d = degrees(distance / max(E.a, E.b)) + _0_5 # some margin

365 self._N = center.lat + d

366 self._S = center.lat - d

367 self._E = center.lon + d

368 self._W = center.lon - d

369

370 def overlaps(self, lat1, lon1, lat2, lon2):

371 '''Check whether this box overlaps a line between 2 points.

372

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

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

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

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

377

378 @return: C{False} if there is certainly no overlap,

379 C{True} otherwise (C{bool}).

380 '''

381 non_ = ((lat1 > self._N or lat2 < self._S) if lat1 < lat2 else

382 (lat2 > self._N or lat1 < self._S)) or \

383 ((lon1 > self._E or lon2 < self._W) if lon1 < lon2 else

384 (lon2 > self._E or lon1 < self._W))

385 return not non_

386

387

388 class _Tol(object):

389 '''Handle a tolerance in C{meter} as C{degrees} and C{meter}.

390 '''

391 _deg = 0 # tol in degrees

392 _lat = 0

393 _m = 0 # tol in meter

394 _min = MAX # degrees

395 _prev = None

396 _r = 0

397

398 def __init__(self, tol_m, E, lat, *lats):

399 '''New L{_Tol}.

400

401 @arg tol_m: Tolerance (C{meter}, only).

402 @arg E: Earth ellipsoid (L{Ellipsoid}).

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

404 @arg lats: Additional latitudes (C{degrees}).

405 '''

406 if lats:

407 lat = fmean_(lat, *lats)

408 self._lat = lat

409 self._r = max(EPS, E.rocMean(lat))

410 self._m = max(EPS, tol_m)

411 self._deg = max(EPS, degrees(self._m / self._r)) # avoid m2degrees!

412

413 @property_RO

414 def degrees(self):

415 '''Get this tolerance in C{degrees}.

416 '''

417 return self._deg

418

419 def degrees2m(self, deg):

420 '''Convert B{C{deg}} to meter at the same C{lat} and earth radius.

421 '''

422 return self.radius * radians(deg) / PI2 # avoid degrees2m!

423

424 def degError(self, Error=_ValueError):

425 '''Compose an error with the C{deg}rees minimum.

426 '''

427 return self.mError(self.degrees2m(self._min), Error=Error)

428

429 def done(self, deg):

430 '''Check C{deg} vs. tolerance and previous value.

431 '''

432 if deg < self._deg or deg == self._prev:

433 return True

434 self._min = min(self._min, deg)

435 self._prev = deg

436 return False

437

438 @property_RO

439 def lat(self):

440 '''Get the mean latitude in C{degrees}.

441 '''

442 return self._lat

443

444 def mError(self, m, Error=_ValueError):

445 '''Compose an error with B{C{m}}eter minimum.

446 '''

447 t = _SPACE_(Fmt.tolerance(self.meter), _too_low_)

448 if m2km(m) > self.meter:

449 t = _or(t, _antipodal_, _near_(_polar__))

450 return Error(Fmt.no_convergence(m), txt=t)

451

452 @property_RO

453 def meter(self):

454 '''Get this tolerance in C{meter}.

455 '''

456 return self._m

457

458 @property_RO

459 def radius(self):

460 '''Get the earth radius in C{meter}.

461 '''

462 return self._r

463

464 def reset(self):

465 '''Reset tolerances.

466 '''

467 self._min = MAX

468 self._prev = None

469

470

471 def _Equidistant00(equidistant, p1):

472 '''(INTERNAL) Get an C{Equidistant*(0, 0, ...)} instance.

473 '''

474 if equidistant is None or not callable(equidistant):

475 equidistant = p1.Equidistant

476 elif not issubclassof(equidistant, *_MODS.azimuthal._Equidistants): # PYCHOK no cover

477 t = tuple(_.__name__ for _ in _MODS.azimuthal._Equidistants)

478 raise _IsnotError(*t, equidistant=equidistant)

479 return equidistant(0, 0, p1.datum)

480

481

482 def _intersect3(s1, end1, s2, end2, height=None, wrap=True, # MCCABE 16

483 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):

484 '''(INTERNAL) Intersect two (ellipsoidal) path, see ellipsoidal method

485 L{intersection3}, separated to allow callers to embellish any exceptions.

486 '''

487 _LLS = _MODS.sphericalTrigonometry.LatLon

488 _si = _MODS.sphericalTrigonometry._intersect

489 _vi3 = _MODS.vector3d._intersect3d3

490

491 def _b_d(s, e, w, t, h=_0_0):

492 # compute opposing and distance

493 t = s.classof(t.lat, t.lon, height=h, name=t.name)

494 t = s.distanceTo3(t, wrap=w) # Distance3Tuple

495 b = opposing(e, t.initial) # "before" start

496 return b, t.distance

497

498 def _b_e(s, e, w, t):

499 # compute an end point along the initial bearing

500 # about 1.5 times the distance to the gu-/estimate, at

501 # least 1/8 and at most 3/8 of the earth perimeter like

502 # radians in .sphericalTrigonometry._int3d2 and bearing

503 # comparison in .sphericalTrigonometry._intb

504 b, d = _b_d(s, e, w, t, h=t.height)

505 m = s.ellipsoid().R2x * PI_4 # authalic exact

506 d = min(max(d * _B2END, m), m * _3_0)

507 e = s.destination(d, e)

508 return b, (_unrollon(s, e) if w else e)

509

510 def _e_ll(s, e, w, **end):

511 # return 2-tuple (end, False if bearing else True)

512 ll = not isscalar(e)

513 if ll:

514 e = s.others(**end)

515 if w: # unroll180 == .karney._unroll2

516 e = _unrollon(s, e)

517 return e, ll

518

519 def _llS(ll): # return an C{_LLS} instance

520 return _LLS(ll.lat, ll.lon, height=ll.height)

521

522 def _o(o, b, n, s, t, e):

523 # determine C{o}utside before, on or after start point

524 if not o: # intersection may be on start

525 if _MODS.latlonBase._isequalTo(s, t, eps=e.degrees):

526 return o

527 return -n if b else n

528

529 E = s1.ellipsoids(s2)

530

531 e1, ll1 = _e_ll(s1, end1, wrap, end1=end1)

532 e2, ll2 = _e_ll(s2, end2, wrap, end2=end2)

533

534 e = _Tol(tol, E, s1.lat, (e1.lat if ll1 else s1.lat),

535 s2.lat, (e2.lat if ll2 else s2.lat))

536

537 # get the azimuthal equidistant projection

538 A = _Equidistant00(equidistant, s1)

539

540 # gu-/estimate initial intersection, spherically ...

541 t = _si(_llS(s1), (_llS(e1) if ll1 else e1),

542 _llS(s2), (_llS(e2) if ll2 else e2),

543 height=height, wrap=False, LatLon=_LLS) # unrolled already

544 h, n = t.height, t.name

545

546 if not ll1:

547 b1, e1 = _b_e(s1, e1, wrap, t)

548 if not ll2:

549 b2, e2 = _b_e(s2, e2, wrap, t)

550

551 # ... and iterate as Karney describes, for references

552 # @see: Function L{ellipsoidalKarney.intersection3}.

553 for i in range(1, _TRIPS):

554 A.reset(t.lat, t.lon) # gu-/estimate as origin

555 # convert start and end points to projection

556 # space and compute an intersection there

557 v, o1, o2 = _vi3(*A._forwards(s1, e1, s2, e2),

558 eps=e.meter, useZ=False)

559 # convert intersection back to geodetic

560 t, d = A._reverse2(v)

561 if e.done(d): # below tol or unchanged?

562 t._iteration = i

563 break

564 else:

565 raise e.degError(Error=IntersectionError)

566

567 # like .sphericalTrigonometry._intersect, if this intersection

568 # is "before" the first point, use the antipodal intersection

569 if not (ll1 or ll2): # end1 and end2 are an initial bearing

570 b1, _ = _b_d(s1, end1, wrap, t)

571 if b1:

572 t = t.antipodal()

573 b1 = not b1

574 b2, _ = _b_d(s2, end2, wrap, t)

575

576 r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=s1,

577 iteration=t._iteration, name=n)

578 return Intersection3Tuple(r, (o1 if ll1 else _o(o1, b1, 1, s1, t, e)),

579 (o2 if ll2 else _o(o2, b2, 2, s2, t, e)))

580

581

582 def _intersection3(start1, end1, start2, end2, height=None, wrap=True,

583 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):

584 '''(INTERNAL) Iteratively compute the intersection point of two paths,

585 each defined by two (ellipsoidal) points or an (ellipsoidal) start

586 point and an initial bearing from North.

587 '''

588 s1 = _xellipsoidal(start1=start1)

589 s2 = s1.others(start2=start2)

590 try:

591 return _intersect3(s1, end1, s2, end2, height=height, wrap=wrap,

592 equidistant=equidistant, tol=tol,

593 LatLon=LatLon, **LatLon_kwds)

594 except (TypeError, ValueError) as x:

595 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)

596

597

598 def _intersections2(center1, radius1, center2, radius2, height=None, wrap=True,

599 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):

600 '''(INTERNAL) Iteratively compute the intersection points of two circles,

601 each defined by an (ellipsoidal) center point and a radius.

602 '''

603 c1 = _xellipsoidal(center1=center1)

604 c2 = c1.others(center2=center2)

605 try:

606 return _intersects2(c1, radius1, c2, radius2, height=height, wrap=wrap,

607 equidistant=equidistant, tol=tol,

608 LatLon=LatLon, **LatLon_kwds)

609 except (TypeError, ValueError) as x:

610 raise _xError(x, center1=center1, radius1=radius1,

611 center2=center2, radius2=radius2)

612

613

614 def _intersects2(c1, radius1, c2, radius2, height=None, wrap=True, # MCCABE 16

615 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds):

616 '''(INTERNAL) Intersect two (ellipsoidal) circles, see L{_intersections2}

617 above, separated to allow callers to embellish any exceptions.

618 '''

619 _LLS = _MODS.sphericalTrigonometry.LatLon

620 _si2 = _MODS.sphericalTrigonometry._intersects2

621 _vi2 = _MODS.vector3d._intersects2

622

623 def _latlon4(t, h, n, c):

624 return _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=c,

625 iteration=t.iteration, name=n)

626

627 def _llS(ll): # return an C{_LLS} instance

628 return _LLS(ll.lat, ll.lon, height=ll.height)

629

630 r1 = Radius_(radius1=radius1)

631 r2 = Radius_(radius2=radius2)

632

633 E = c1.ellipsoids(c2)

634 # get the azimuthal equidistant projection

635 A = _Equidistant00(equidistant, c1)

636

637 if r1 < r2:

638 c1, c2 = c2, c1

639 r1, r2 = r2, r1

640

641 if r1 > (min(E.b, E.a) * PI):

642 raise _ValueError(_exceed_PI_radians_)

643

644 if wrap: # unroll180 == .karney._unroll2

645 c2 = _unrollon(c1, c2)

646

647 # distance between centers and radii are

648 # measured along the ellipsoid's surface

649 m = c1.distanceTo(c2, wrap=False) # meter

650 if m < max(r1 - r2, EPS):

651 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?

652 if fsum_(r1, r2, -m) < 0:

653 raise IntersectionError(_too_(Fmt.distant(m)))

654

655 f = _radical2(m, r1, r2).ratio # "radical fraction"

656 e = _Tol(tol, E, favg(c1.lat, c2.lat, f=f))

657

658 # gu-/estimate initial intersections, spherically ...

659 t1, t2 = _si2(_llS(c1), r1, _llS(c2), r2, radius=e.radius,

660 height=height, wrap=False, too_d=m) # unrolled already

661 h, n = t1.height, t1.name

662

663 # ... and iterate as Karney describes, for references

664 # @see: Function L{ellipsoidalKarney.intersections2}.

665 ts, ta = [], None

666 for t in ((t1,) if t1 is t2 else (t1, t2)):

667 for i in range(1, _TRIPS):

668 A.reset(t.lat, t.lon) # gu-/estimate as origin

669 # convert centers to projection space and

670 # compute the intersections there

671 t1, t2 = A._forwards(c1, c2)

672 v1, v2 = _vi2(t1, r1, # XXX * t1.scale?,

673 t2, r2, # XXX * t2.scale?,

674 sphere=False, too_d=m)

675 # convert intersections back to geodetic

676 if v1 is v2: # abutting

677 t, d = A._reverse2(v1)

678 else: # consider the closer intersection

679 t1, d1 = A._reverse2(v1)

680 t2, d2 = A._reverse2(v2)

681 t, d = (t1, d1) if d1 < d2 else (t2, d2)

682 if e.done(d): # below tol or unchanged?

683 t._iteration = i # _NamedTuple._iteration

684 ts.append(t)

685 if v1 is v2: # abutting

686 ta = t

687 break

688 else:

689 raise e.degError(Error=IntersectionError)

690 e.reset()

691

692 if ta: # abutting circles

693 pass # PYCHOK no cover

694 elif len(ts) == 2:

695 return (_latlon4(ts[0], h, n, c1),

696 _latlon4(ts[1], h, n, c2))

697 elif len(ts) == 1: # PYCHOK no cover

698 ta = ts[0] # assume abutting

699 else: # PYCHOK no cover

700 raise _AssertionError(ts=ts)

701 r = _latlon4(ta, h, n, c1)

702 return r, r

703

704

705 def _nearestOn2(p, point1, point2, within=True, height=None, wrap=True,

706 equidistant=None, tol=_TOL_M, **LatLon_and_kwds):

707 '''(INTERNAL) Closest point and fraction, like L{_intersects2} above,

708 separated to allow callers to embellish any exceptions.

709 '''

710 p1 = p.others(point1=point1)

711 p2 = p.others(point2=point2)

712

713 _ = p.ellipsoids(p1)

714# E = p.ellipsoids(p2) # done in __nearestOn2__

715

716 # get the azimuthal equidistant projection

717 A = _Equidistant00(equidistant, p)

718

719 if wrap:

720 p1 = _unrollon(p, p1) # XXX do not unroll?

721 p2 = _unrollon(p, p2) # XXX do not unroll?

722 p2 = _unrollon(p1, p2)

723

724 r, f, e = _nearestOn3_(p, p1, p2, A, within=within, height=height,

725 tol=tol, **LatLon_and_kwds)

726 return NearestOn2Tuple(r, f)

727

728

729 def _nearestOn3_(p, p1, p2, A, within=True, height=None, tol=_TOL_M,

730 LatLon=None, **LatLon_kwds):

731 # Only in function C{_nearestOn2} and method C{nearestOn8} above

732 _LLS = _MODS.sphericalNvector.LatLon

733 _V3d = _MODS.vector3d.Vector3d

734 _vnOn2 = _MODS.vector3d._nearestOn2

735

736 def _llS(ll): # return an C{_LLS} instance

737 return _LLS(ll.lat, ll.lon, height=ll.height)

738

739 E = p.ellipsoids(p2)

740 e = _Tol(tol, E, p.lat, p1.lat, p2.lat)

741

742 # gu-/estimate initial nearestOn, spherically ... wrap=False, only!

743 # using sphericalNvector.LatLon.nearestOn for within=False support

744 t = _llS(p).nearestOn(_llS(p1), _llS(p2), within=within,

745 height=height)

746 n, h = t.name, t.height

747 if height is None:

748 h1 = p1.height # use heights as pseudo-Z in projection space

749 h2 = p2.height # to be included in the closest function

750 h0 = favg(h1, h2)

751 else: # ignore heights in distances, Z=0

752 h0 = h1 = h2 = _0_0

753

754 # ... and iterate to find the closest (to the origin with .z

755 # to interpolate height) as Karney describes, for references

756 # @see: Function L{ellipsoidalKarney.nearestOn}.

757 vp, f = _V3d(_0_0, _0_0, h0), None

758 for i in range(1, _TRIPS):

759 A.reset(t.lat, t.lon) # gu-/estimate as origin

760 # convert points to projection space and compute

761 # the nearest one (and its height) there

762 s, t = A._forwards(p1, p2)

763 v, f = _vnOn2(vp, _V3d(s.x, s.y, h1),

764 _V3d(t.x, t.y, h2), within=within)

765 # convert nearest one back to geodetic

766 t, d = A._reverse2(v)

767 if e.done(d): # below tol or unchanged?

768 t._iteration = i # _NamedTuple._iteration

769 break

770 else:

771 raise e.degError()

772

773 if height is None:

774 h = v.z # nearest

775 elif isscalar(height):

776 h = height

777 r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=p,

778 iteration=t.iteration, name=n)

779 return r, f, e # fraction or None

780

781

782__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI)

783

784# **) MIT License

785#

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

787#

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

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

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

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

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

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

794#

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

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

797#

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

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

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

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

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

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

804# OTHER DEALINGS IN THE SOFTWARE.