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, fsumf_

19 from pygeodesy.formy import opposing, _radical2

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

21 _low_, _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, Radius_, Scalar

29 from pygeodesy.utily import m2km, unroll180, _unrollon, _unrollon3, \

30 _Wrap, wrap360

31

32 from math import degrees, radians

33

34__all__ = _ALL_LAZY.ellipsoidalBaseDI

35__version__ = '23.05.15'

36

37_polar__ = 'polar?'

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: If C{True}, wrap or I{normalize} and unroll the

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

56

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

58

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

60

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

62 ellipsoids are not compatible.

63 '''

64 r = self._Inverse(other, wrap)

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

66

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

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

69 the given distance from this point along a geodesic given

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

71 method L{destination2} for more details.

72

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

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

75 @kwarg height: Optional height, overriding the default

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

77

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

79 '''

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

81

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

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

84 azimuth) after having travelled for the given distance from

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

86 using this C{Direct} method.

87

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

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

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

91

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

93 are in compass C{degrees360}.

94

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

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

97

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

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

100 @kwarg height: Optional height, overriding the default

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

102

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

104 '''

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

106 return self._xnamed(r)

107

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

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

110

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

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

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

114 '''

115 g = self.geodesic

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

117 if LL:

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

119 return r

120

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

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

123 '''

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

125 d = LL(*_Wrap.latlon(r.lat, r.lon), height=h, **_xkwds_not(None,

126 datum=self.datum, name=self.name,

127 epoch=self.epoch, reframe=self.reframe))

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

129

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

131 '''Compute the distance between this and an other point along

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

133 L{distanceTo3} for more details.

134

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

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

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

138

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

140

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

142

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

144 ellipsoids are not compatible.

145 '''

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

147

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

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

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

151 C{Inverse} method.

152

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

154 conventionally meter. The distance is measured on the surface

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

156

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

158 are in compass C{degrees360} from North.

159

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

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

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

163

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

165

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

167

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

169 ellipsoids are not compatible.

170 '''

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

172

173 def finalBearingOn(self, distance, bearing):

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

175 travelled for the given distance along a geodesic given by

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

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

178

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

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

181

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

183 '''

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

185

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

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

188 travelled along a geodesic from this point to an other

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

190 L{distanceTo3} for more details.

191

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

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

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

195

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

197

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

199

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

201 ellipsoids are not compatible.

202 '''

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

204

205 @Property_RO

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

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

208 '''

209 return None # PYCHOK no cover

210

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

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

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

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

215 for more details.

216

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

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

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

220

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

222

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

224

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

226 ellipsoids are not compatible.

227 '''

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

229

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

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

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

233 methods.

234

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

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

237 at this and 1.0 at the other point.

238 @kwarg height: Optional height, overriding the fractional

239 height (C{meter}).

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

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

242

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

244

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

246

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

248

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

250 ellipsoids are not compatible.

251

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

253 C{rhumbMidpointTo}.

254 '''

255 f = Scalar(fraction=fraction)

256 if isnear0(f):

257 r = self

258 elif isnear1(f) and not wrap:

259 r = self.others(other)

260 else: # negative fraction OK

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

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

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

264 return r

265

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

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

268

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

270 '''

271 _ = self.ellipsoids(other)

272 g = self.geodesic

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

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

275

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

277 equidistant=None, tol=_TOL_M):

278 '''I{Iteratively} locate the point on a path or polygon closest

279 to this point.

280

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

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

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

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

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

286 is taken into account for distances.

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

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

289

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

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

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

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

294

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

296

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

298

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

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

301

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

303 '''

304 _d3 = self.distanceTo3 # Distance3Tuple

305 _n3 = _nearestOn3

306 try:

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

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

309 _ = self.ellipsoids(p1)

310 c3 = _d3(c, wrap=wrap) # XXX wrap=False?

311

312 except (TypeError, ValueError) as x:

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

314 closed=closed, height=height, wrap=wrap)

315

316 # get the azimuthal equidistant projection, once

317 A = _Equidistant00(equidistant, c)

318 b = _Box(c, c3.distance)

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

320

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

322 LatLon=self.classof, # this LatLon

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

324 try:

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

326 if wrap and j != 0: # not Ps.looped

327 p2 = _unrollon(p1, p2)

328 # skip edge if no overlap with box around closest

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

330 p, t, _ = _n3(self, p1, p2, A, **kwds)

331 d3 = _d3(p, wrap=False) # already unrolled

332 if d3.distance < c3.distance:

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

334 b = _Box(c, c3.distance)

335 m = max(m, c.iteration)

336 p1, i = p2, j

337

338 except (TypeError, ValueError) as x:

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

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

341 closed=closed, height=height, wrap=wrap)

342

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

344

345 n = self.nearestOn8.__name__

346 c.rename(n)

347 if s is not c:

348 s = s.copy(name=n)

349 if e is not c:

350 e = e.copy(name=n)

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

352 iteration=m) # ._iteration for tests

353

354

355 class _Box(object):

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

357

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

359 '''

360 def __init__(self, center, distance):

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

362

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

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

365 (C{meter}, conventionally)

366 '''

367 E = center.ellipsoid()

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

369 self._N = center.lat + d

370 self._S = center.lat - d

371 self._E = center.lon + d

372 self._W = center.lon - d

373

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

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

376

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

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

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

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

381

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

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

384 '''

385 if lat1 > lat2:

386 lat1, lat2 = lat2, lat1

387 if lat1 > self._N or lat2 < self._S:

388 return False

389 if lon1 > lon2:

390 lon1, lon2 = lon2, lon1

391 if lon1 > self._E or lon2 < self._W:

392 return False

393 return True

394

395

396 class _Tol(object):

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

398 '''

399 _deg = 0 # tol in degrees

400 _lat = 0

401 _m = 0 # tol in meter

402 _min = MAX # degrees

403 _prev = None

404 _r = 0

405

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

407 '''New L{_Tol}.

408

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

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

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

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

413 '''

414 self._lat = fmean_(lat, *lats) if lats else lat

415 self._r = max(EPS, E.rocMean(self._lat))

416 self._m = max(EPS, tol_m)

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

418

419 @property_RO

420 def degrees(self):

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

422 '''

423 return self._deg

424

425 def degrees2m(self, deg):

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

427 '''

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

429

430 def degError(self, Error=_ValueError):

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

432 '''

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

434

435 def done(self, deg):

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

437 '''

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

439 return True

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

441 self._prev = deg

442 return False

443

444 @property_RO

445 def lat(self):

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

447 '''

448 return self._lat

449

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

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

452 '''

453 t = _SPACE_(Fmt.tolerance(self.meter), _too_(_low_))

454 if m2km(m) > self.meter:

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

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

457

458 @property_RO

459 def meter(self):

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

461 '''

462 return self._m

463

464 @property_RO

465 def radius(self):

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

467 '''

468 return self._r

469

470 def reset(self):

471 '''Reset tolerances.

472 '''

473 self._min = MAX

474 self._prev = None

475

476

477 def _Equidistant00(equidistant, p1):

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

479 '''

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

481 equidistant = p1.Equidistant

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

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

484 raise _IsnotError(*t, equidistant=equidistant)

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

486

487

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

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

490 '''(INTERNAL) Intersect two (ellipsoidal) lines, see ellipsoidal method

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

492 '''

493 _LLS = _MODS.sphericalTrigonometry.LatLon

494 _si = _MODS.sphericalTrigonometry._intersect

495 _vi3 = _MODS.vector3d._intersect3d3

496

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

498 # compute opposing and distance

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

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

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

502 return b, t.distance

503

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

505 # compute an end point along the initial bearing

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

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

508 # radians in .sphericalTrigonometry._int3d2 and bearing

509 # comparison in .sphericalTrigonometry._intb

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

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

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

513 e = s.destination(d, e)

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

515

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

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

518 ll = not isscalar(e)

519 if ll:

520 e = s.others(**end)

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

522 e = _unrollon(s, e)

523 return e, ll

524

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

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

527

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

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

530 if not o: # intersection may be on start

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

532 return o

533 return -n if b else n

534

535 E = s1.ellipsoids(s2)

536

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

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

539

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

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

542

543 # get the azimuthal equidistant projection

544 A = _Equidistant00(equidistant, s1)

545

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

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

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

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

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

551

552 if not ll1:

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

554 if not ll2:

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

556

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

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

559 for i in range(1, _TRIPS):

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

561 # convert start and end points to projection

562 # space and compute an intersection there

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

564 eps=e.meter, useZ=False)

565 # convert intersection back to geodetic

566 t, d = A._reverse2(v)

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

568 t._iteration = i

569 break

570 else:

571 raise e.degError(Error=IntersectionError)

572

573 # like .sphericalTrigonometry._intersect, if this intersection

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

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

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

577 if b1:

578 t = t.antipodal()

579 b1 = not b1

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

581

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

583 iteration=t._iteration, name=n)

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

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

586

587

588 def _intersection3(start1, end1, start2, end2, height=None, wrap=False, # was=True

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

590 '''(INTERNAL) Iteratively compute the intersection point of two lines,

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

592 point and an initial bearing from North.

593 '''

594 s1 = _xellipsoidal(start1=start1)

595 s2 = s1.others(start2=start2)

596 try:

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

598 equidistant=equidistant, tol=tol,

599 LatLon=LatLon, **LatLon_kwds)

600 except (TypeError, ValueError) as x:

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

602

603

604 def _intersections2(center1, radius1, center2, radius2, height=None, wrap=False, # was=True

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

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

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

608 '''

609 c1 = _xellipsoidal(center1=center1)

610 c2 = c1.others(center2=center2)

611 try:

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

613 equidistant=equidistant, tol=tol,

614 LatLon=LatLon, **LatLon_kwds)

615 except (TypeError, ValueError) as x:

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

617 center2=center2, radius2=radius2)

618

619

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

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

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

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

624 '''

625 _LLS = _MODS.sphericalTrigonometry.LatLon

626 _si2 = _MODS.sphericalTrigonometry._intersects2

627 _vi2 = _MODS.vector3d._intersects2

628

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

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

631 iteration=t.iteration, name=n)

632

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

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

635

636 r1 = Radius_(radius1=radius1)

637 r2 = Radius_(radius2=radius2)

638

639 E = c1.ellipsoids(c2)

640 # get the azimuthal equidistant projection

641 A = _Equidistant00(equidistant, c1)

642

643 if r1 < r2:

644 c1, c2 = c2, c1

645 r1, r2 = r2, r1

646

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

648 raise _ValueError(_exceed_PI_radians_)

649

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

651 c2 = _unrollon(c1, c2)

652

653 # distance between centers and radii are

654 # measured along the ellipsoid's surface

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

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

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

658 if fsumf_(r1, r2, -m) < 0:

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

660

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

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

663

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

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

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

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

668

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

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

671 ts, ta = [], None

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

673 for i in range(1, _TRIPS):

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

675 # convert centers to projection space and

676 # compute the intersections there

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

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

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

680 sphere=False, too_d=m)

681 # convert intersections back to geodetic

682 if v1 is v2: # abutting

683 t, d = A._reverse2(v1)

684 else: # consider the closer intersection

685 t1, d1 = A._reverse2(v1)

686 t2, d2 = A._reverse2(v2)

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

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

689 t._iteration = i # _NamedTuple._iteration

690 ts.append(t)

691 if v1 is v2: # abutting

692 ta = t

693 break

694 else:

695 raise e.degError(Error=IntersectionError)

696 e.reset()

697

698 if ta: # abutting circles

699 pass # PYCHOK no cover

700 elif len(ts) == 2:

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

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

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

704 ta = ts[0] # assume abutting

705 else: # PYCHOK no cover

706 raise _AssertionError(ts=ts)

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

708 return r, r

709

710

711 def _nearestOn2(p, point1, point2, within=True, height=None, wrap=False, # was=True

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

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

714 separated to allow callers to embellish any exceptions.

715 '''

716 p1 = p.others(point1=point1)

717 p2 = p.others(point2=point2)

718

719 _ = p.ellipsoids(p1)

720# E = p.ellipsoids(p2) # done in _nearestOn3

721

722 # get the azimuthal equidistant projection

723 A = _Equidistant00(equidistant, p)

724

725 p1, p2, _ = _unrollon3(p, p1, p2, wrap) # XXX don't unroll?

726 r, f, _ = _nearestOn3(p, p1, p2, A, within=within, height=height,

727 tol=tol, **LatLon_and_kwds)

728 return NearestOn2Tuple(r, f)

729

730

731 def _nearestOn3(p, p1, p2, A, within=True, height=None, tol=_TOL_M,

732 LatLon=None, **LatLon_kwds):

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

734 _LLS = _MODS.sphericalNvector.LatLon

735 _V3d = _MODS.vector3d.Vector3d

736 _vn2 = _MODS.vector3d._nearestOn2

737

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

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

740

741 E = p.ellipsoids(p2)

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

743

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

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

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

747 height=height)

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

749 if height is None:

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

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

752 h0 = favg(h1, h2)

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

754 h0 = h1 = h2 = _0_0

755

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

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

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

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

760 for i in range(1, _TRIPS):

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

762 # convert points to projection space and compute

763 # the nearest one (and its height) there

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

765 v, f = _vn2(vp, _V3d(s.x, s.y, h1),

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

767 # convert nearest one back to geodetic

768 t, d = A._reverse2(v)

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

770 t._iteration = i # _NamedTuple._iteration

771 break

772 else:

773 raise e.degError()

774

775 if height is None:

776 h = v.z # nearest

777 elif isscalar(height):

778 h = height

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

780 iteration=t.iteration, name=n)

781 return r, f, e # fraction or None

782

783

784__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI)

785

786# **) MIT License

787#

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

789#

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

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

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

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

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

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

796#

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

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

799#

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

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

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

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

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

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

806# OTHER DEALINGS IN THE SOFTWARE.