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 isLatLon, issubclassof, map2

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

12 _EPSqrt as _TOL, _0_0, _0_5, _1_5, _3_0

13# from pygeodesy.dms import F_DMS # _MODS

14 from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase, _TOL_M, property_RO

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 _isequalTo, opposing, _radical2

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

21 _exceed_PI_radians_, _low_, _near_, \

22 _SPACE_, _too_

23 from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS, \

24 _dunder_nameof

25 from pygeodesy.namedTuples import Bearing2Tuple, Destination2Tuple, \

26 Intersection3Tuple, NearestOn2Tuple, \

27 NearestOn8Tuple, _LL4Tuple

28# from pygeodesy.props import property_RO # from .ellipsoidalBase

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

30 from pygeodesy.units import _fi_j2, _isDegrees, _isHeight, _isRadius, \

31 Radius_, Scalar

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

33 _Wrap, wrap360

34

35 from math import degrees, radians

36

37__all__ = _ALL_LAZY.ellipsoidalBaseDI

38__version__ = '24.05.19'

39

40_polar__ = 'polar?'

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

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

43

44

45 class LatLonEllipsoidalBaseDI(LatLonEllipsoidalBase):

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

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

48 '''

49

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

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

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

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

54 for more details.

55

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

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

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

59

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

61

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

63

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

65 ellipsoids are not compatible.

66 '''

67 r = self._Inverse(other, wrap)

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

69

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

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

72 the given distance from this point along a geodesic given

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

74 method L{destination2} for more details.

75

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

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

78 @kwarg height: Optional height, overriding the default

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

80

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

82 '''

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

84

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

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

87 azimuth) after having travelled for the given distance from

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

89 using this C{Direct} method.

90

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

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

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

94

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

96 are in compass C{degrees360}.

97

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

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

100

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

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

103 @kwarg height: Optional height, overriding the default

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

105

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

107 '''

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

109 return self._xnamed(r)

110

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

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

113

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

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

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

117 '''

118 g = self.geodesic

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

120 if LL:

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

122 return r

123

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

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

126 '''

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

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

129 datum=self.datum, name=self.name,

130 epoch=self.epoch, reframe=self.reframe))

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

132

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

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

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

136 L{distanceTo3} for more details.

137

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

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

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

141

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

143

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

145

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

147 ellipsoids are not compatible.

148 '''

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

150

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

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

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

154 C{Inverse} method.

155

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

157 conventionally meter. The distance is measured on the surface

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

159

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

161 are in compass C{degrees360} from North.

162

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

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

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

166

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

168

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

170

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

172 ellipsoids are not compatible.

173 '''

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

175

176 def finalBearingOn(self, distance, bearing):

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

178 travelled for the given distance along a geodesic given by

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

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

181

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

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

184

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

186 '''

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

188

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

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

191 travelled along a geodesic from this point to an other

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

193 L{distanceTo3} for more details.

194

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

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

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

198

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

200

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

202

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

204 ellipsoids are not compatible.

205 '''

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

207

208 @property_RO

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

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

211 '''

212 return None # PYCHOK no cover

213

214 def _g_gl_p3(self, point, other, exact, wrap):

215 '''(INTERNAL) Helper for methods C{.intersecant2} and C{.plumbTo}.

216 '''

217 p = _unrollon(self, self.others(point=point), wrap=wrap)

218 g = self.datum.ellipsoid.geodesic_(exact=exact)

219 gl = g._DirectLine( p, other) if _isDegrees(other) else \

220 g._InverseLine(p, self.others(other), wrap)

221 return g, gl, p

222

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

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

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

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

227 for more details.

228

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

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

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

232

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

234

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

236

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

238 ellipsoids are not compatible.

239 '''

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

241

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

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

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

245 methods.

246

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

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

249 at this and 1.0 at the other point.

250 @kwarg height: Optional height, overriding the fractional

251 height (C{meter}).

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

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

254

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

256

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

258

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

260

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

262 ellipsoids are not compatible.

263

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

265 C{rhumbMidpointTo}.

266 '''

267 f = Scalar(fraction=fraction)

268 if isnear0(f):

269 r = self

270 elif isnear1(f) and not wrap:

271 r = self.others(other)

272 else: # negative fraction OK

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

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

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

276 return r

277

278 def intersecant2(self, circle, point, other, exact=False, height=None, # PYCHOK signature

279 wrap=False, tol=_TOL):

280 '''Compute the intersections of a circle and a geodesic (line) given as

281 two points or as a point and bearing.

282

283 @arg circle: Radius of the circle centered at this location (C{meter},

284 conventionally) or a point on the circle (this C{LatLon}).

285 @arg point: A location on the geodesic (this C{LatLon}).

286 @arg other: An other point on the geodesic (this C{LatLon}) or the

287 (forward) bearing at the B{C{point}} (compass C{degrees}).

288 @kwarg exact: Exact C{geodesic...} to use (C{bool} or C{Geodesic...}), see

289 method L{Ellipsoid.geodesic_}.

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

291 conventionally) or C{None} for interpolated heights.

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

293 B{C{circle}} and/or B{C{other}} (C{bool}).

294 @kwarg tol: Convergence tolerance (C{scalar}).

295

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

297 each an instance of this C{LatLon}. Both points are the same

298 instance if the geodesic (line) is tangential to the circle.

299

300 @raise IntersectionError: The circle and geodesic do not intersect.

301

302 @raise TypeError: If B{C{point}} is not this C{LatLon} or B{C{circle}}

303 or B{C{other}} invalid.

304

305 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{exact}} or

306 B{C{height}}.

307

308 @see: Method L{rhumbIntersecant2<LatLonBase.rhumbIntersecant2>}.

309 '''

310 try:

311 g, gl, p = self._g_gl_p3(point, other, exact, wrap)

312 r = Radius_(circle=circle) if _isRadius(circle) else \

313 g._Inverse(self, self.others(circle=circle), wrap).s12

314

315 P, Q = _MODS.geodesicw._Intersecant2(gl, self.lat, self.lon, r, tol=tol,

316 form=_MODS.dms.F_DMS)

317 return self._intersecend2(p, other, wrap, height, g, P, Q,

318 self.intersecant2)

319 except (TypeError, ValueError) as x:

320 raise _xError(x, center=self, circle=circle, point=point, other=other,

321 exact=exact, wrap=wrap)

322

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

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

325

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

327 '''

328 _ = self.ellipsoids(other)

329 g = self.geodesic

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

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

332

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

334 equidistant=None, tol=_TOL_M):

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

336 to this point.

337

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

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

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

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

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

343 is taken into account for distances.

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

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

346

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

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

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

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

351

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

353

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

355

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

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

358

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

360 '''

361 _d3 = self.distanceTo3 # Distance3Tuple

362 _n3 = _nearestOn3

363 try:

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

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

366 _ = self.ellipsoids(p1)

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

368

369 except (TypeError, ValueError) as x:

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

371 closed=closed, height=height, wrap=wrap)

372

373 # get the azimuthal equidistant projection, once

374 A = _Equidistant00(equidistant, c)

375 b = _Box(c, c3.distance)

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

377

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

379 LatLon=self.classof, # this LatLon

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

381 try:

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

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

384 p2 = _unrollon(p1, p2)

385 # skip edge if no overlap with box around closest

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

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

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

389 if d3.distance < c3.distance:

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

391 b = _Box(c, c3.distance)

392 m = max(m, c.iteration)

393 p1, i = p2, j

394

395 except (TypeError, ValueError) as x:

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

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

398 closed=closed, height=height, wrap=wrap)

399

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

401

402 n = self.nearestOn8.__name__ # _dunder_nameof

403 c.rename(n)

404 if s is not c:

405 s = s.copy(name=n)

406 if e is not c:

407 e = e.copy(name=n) # name__=self.nearestOn8

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

409 iteration=m) # ._iteration for tests

410

411 def plumbTo(self, point, other, exact=False, height=None, # PYCHOK signature

412 wrap=False, tol=_TOL):

413 '''Compute the I{perpendicular} intersection of a geodesic (line) with

414 a geodesic from this point.

415

416 @arg point: A location (C{LatLon}) on the geodesic (line).

417 @arg other: An other point (C{LatLon}) on the geodesic (line) or the

418 (forward) bearing at the B{C{point}} (compass C{degrees}).

419 @kwarg exact: Exact C{geodesic...} to use (C{bool} or C{Geodesic...}),

420 see method L{Ellipsoid.geodesic_}.

421 @kwarg height: Optional height for the intersection point (C{meter},

422 conventionally) or C{None} for an interpolated height.

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

424 and/or B{C{other}} (C{bool}).

425 @kwarg tol: Convergence tolerance (C{scalar}).

426

427 @return: The intersection point, an instance of this C{LatLon}.

428

429 @raise TypeError: If B{C{point}} or B{C{other}} not this C{LatLon}.

430

431 @raise UnitError: Invalid B{C{other}}, B{C{exact}} or B{C{height}}.

432 '''

433 try:

434 g, gl, p = self._g_gl_p3(point, other, exact, wrap)

435

436 P = _MODS.geodesicw._PlumbTo(gl, self.lat, self.lon, tol=tol)

437 h = self._havg(p, h=height)

438 p = self.classof(P.lat2, P.lon2, datum=self.datum, height=h) # name=n

439 p._iteration = P.iteration

440 except (TypeError, ValueError) as x:

441 raise _xError(x, plumb=self, point=point, other=other,

442 exact=exact, wrap=wrap)

443 return p

444

445

446 class _Box(object):

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

448

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

450 '''

451 def __init__(self, center, distance):

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

453

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

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

456 (C{meter}, conventionally)

457 '''

458 m = Radius_(distance=distance)

459 E = center.ellipsoid() # XXX see above

460 d = E.m2degrees(m) + _0_5 # some margin

461 self._N = center.lat + d

462 self._S = center.lat - d

463 self._E = center.lon + d

464 self._W = center.lon - d

465

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

467 '''Check whether this box overlaps an other box.

468

469 @arg lat1: Latitude of a box corner (C{degrees}).

470 @arg lon1: Longitude of a box corner (C{degrees}).

471 @arg lat2: Latitude of the opposing corner (C{degrees}).

472 @arg lon2: Longitude of the opposing corner (C{degrees}).

473

474 @return: C{True} if there is some overlap, C{False}

475 otherwise (C{bool}).

476 '''

477 if lat1 > lat2:

478 lat1, lat2 = lat2, lat1

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

480 return False

481 if lon1 > lon2:

482 lon1, lon2 = lon2, lon1

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

484 return False

485 return True

486

487

488 class _Tol(object):

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

490 '''

491 _deg = 0 # tol in degrees

492 _lat = 0

493 _m = 0 # tol in meter

494 _min = MAX # degrees

495 _prev = None

496 _r = 0

497

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

499 '''New L{_Tol}.

500

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

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

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

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

505 '''

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

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

508 self._m = max(EPS, tol_m)

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

510

511 @property_RO

512 def degrees(self):

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

514 '''

515 return self._deg

516

517 def degrees2m(self, deg):

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

519 '''

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

521

522 def degError(self, Error=_ValueError):

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

524 '''

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

526

527 def done(self, deg):

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

529 '''

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

531 return True

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

533 self._prev = deg

534 return False

535

536 @property_RO

537 def lat(self):

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

539 '''

540 return self._lat

541

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

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

544 '''

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

546 if m2km(m) > self.meter:

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

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

549

550 @property_RO

551 def meter(self):

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

553 '''

554 return self._m

555

556 @property_RO

557 def radius(self):

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

559 '''

560 return self._r

561

562 def reset(self):

563 '''Reset tolerances.

564 '''

565 self._min = MAX

566 self._prev = None

567

568

569 def _Equidistant00(equidistant, p1):

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

571 '''

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

573 equidistant = p1.Equidistant

574 else:

575 Eqs = _MODS.azimuthal._Equidistants

576 if not issubclassof(equidistant, *Eqs): # PYCHOK no cover

577 t = map2(_dunder_nameof, Eqs)

578 raise _IsnotError(*t, equidistant=equidistant)

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

580

581

582 def intersecant2(center, circle, point, other, **exact_height_wrap_tol):

583 '''Compute the intersections of a circle and a geodesic given as two points

584 or as a point and (forward) bearing.

585

586 @arg center: Center of the circle (C{LatLon}).

587 @arg circle: Radius of the circle (C{meter}, conventionally) or a point on

588 the circle (C{LatLon}, as B{C{center}}).

589 @arg point: A point of the geodesic (C{LatLon}, as B{C{center}}).

590 @arg other: An other point of the geodesic (C{LatLon}, as B{C{center}}) or

591 the (forward) bearing at the B{C{point}} (compass C{degrees}).

592 @kwarg exact_height_wrap_tol: Optional keyword arguments, see below.

593

594 @raise NotImplementedError: Method C{intersecant2} not available.

595

596 @raise TypeError: If B{C{center}}, B{C{point}} or B{C{circle}} or B{C{other}}

597 points not ellipsoidal or not compatible with B{C{center}}.

598

599 @see: Method C{LatLon.intersecant2} of class L{ellipsoidalExact.LatLon},

600 L{ellipsoidalKarney.LatLon} or L{ellipsoidalVincenty.LatLon}.

601 '''

602 if not isLatLon(center, ellipsoidal=True): # isinstance(center, LatLonEllipsoidalBase)

603 raise _IsnotError(_ellipsoidal_, center=center)

604 return center.intersecant2(circle, point, other, **exact_height_wrap_tol)

605

606

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

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

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

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

611 '''

612 _LLS = _MODS.sphericalTrigonometry.LatLon

613 _si = _MODS.sphericalTrigonometry._intersect

614 _vi3 = _MODS.vector3d._intersect3d3

615

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

617 # compute opposing and distance

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

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

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

621 return b, t.distance

622

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

624 # compute an end point along the initial bearing about

625 # 1.5 times the distance to the gu-/estimate, at least

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

627 # radians in .sphericalTrigonometry._int3d2 and bearing

628 # comparison in .sphericaltrigonometry._intb

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

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

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

632 e = s.destination(d, e)

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

634

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

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

637 ll = not _isDegrees(e)

638 if ll:

639 e = s.others(**end)

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

641 e = _unrollon(s, e)

642 return e, ll

643

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

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

646 if not o: # intersection may be on start

647 if _isequalTo(s, t, eps=e.degrees):

648 return o

649 return -n if b else n

650

651 E = s1.ellipsoids(s2)

652

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

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

655

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

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

658

659 # get the azimuthal equidistant projection

660 A = _Equidistant00(equidistant, s1)

661

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

663 t = _si(_LLS(s1), (_LLS(e1) if ll1 else e1),

664 _LLS(s2), (_LLS(e2) if ll2 else e2),

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

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

667

668 if not ll1:

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

670 if not ll2:

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

672

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

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

675 for i in range(1, _TRIPS):

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

677 # convert start and end points to projection

678 # space and compute an intersection there

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

680 eps=e.meter, useZ=False)

681 # convert intersection back to geodetic

682 t, d = A._reverse2(v)

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

684 break

685 else:

686 raise e.degError(Error=IntersectionError)

687

688 # like .sphericalTrigonometry._intersect, if this intersection

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

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

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

692 if b1:

693 t = t.antipodal()

694 b1 = not b1

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

696

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

698 iteration=i, name=n)

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

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

701

702

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

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

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

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

707 point and an initial bearing from North.

708 '''

709 s1 = _xellipsoidal(start1=start1)

710 s2 = s1.others(start2=start2)

711 try:

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

713 equidistant=equidistant, tol=tol,

714 LatLon=LatLon, **LatLon_kwds)

715 except (TypeError, ValueError) as x:

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

717

718

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

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

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

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

723 '''

724 c1 = _xellipsoidal(center1=center1)

725 c2 = c1.others(center2=center2)

726 try:

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

728 equidistant=equidistant, tol=tol,

729 LatLon=LatLon, **LatLon_kwds)

730 except (TypeError, ValueError) as x:

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

732 center2=center2, radius2=radius2)

733

734

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

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

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

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

739 '''

740 _LLS = _MODS.sphericalTrigonometry.LatLon

741 _si2 = _MODS.sphericalTrigonometry._intersects2

742 _vi2 = _MODS.vector3d._intersects2

743

744 def _ll4(t, h, n, c):

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

746 iteration=t.iteration, name=n)

747

748 r1 = Radius_(radius1=radius1)

749 r2 = Radius_(radius2=radius2)

750

751 E = c1.ellipsoids(c2)

752 # get the azimuthal equidistant projection

753 A = _Equidistant00(equidistant, c1)

754

755 if r1 < r2:

756 c1, c2 = c2, c1

757 r1, r2 = r2, r1

758

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

760 raise _ValueError(_exceed_PI_radians_)

761

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

763 c2 = _unrollon(c1, c2)

764

765 # distance between centers and radii are

766 # measured along the ellipsoid's surface

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

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

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

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

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

772

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

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

775

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

777 t1, t2 = _si2(_LLS(c1), r1, _LLS(c2), r2, radius=e.radius,

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

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

780

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

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

783 ts, ta = [], None

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

785 for i in range(1, _TRIPS):

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

787 # convert centers to projection space and

788 # compute the intersections there

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

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

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

792 sphere=False, too_d=m)

793 # convert intersections back to geodetic

794 if v1 is v2: # abutting

795 t, d = A._reverse2(v1)

796 else: # consider the closer intersection

797 t1, d1 = A._reverse2(v1)

798 t2, d2 = A._reverse2(v2)

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

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

801 t._iteration = i # _NamedTuple._iteration

802 ts.append(t)

803 if v1 is v2: # abutting

804 ta = t

805 break

806 else:

807 raise e.degError(Error=IntersectionError)

808 e.reset()

809

810 if ta: # abutting circles

811 pass # PYCHOK no cover

812 elif len(ts) == 2:

813 return (_ll4(ts[0], h, n, c1),

814 _ll4(ts[1], h, n, c2))

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

816 ta = ts[0] # assume abutting

817 else: # PYCHOK no cover

818 raise _AssertionError(ts=ts)

819 r = _ll4(ta, h, n, c1)

820 return r, r

821

822

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

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

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

826 separated to allow callers to embellish any exceptions.

827 '''

828 p1 = p.others(point1=point1)

829 p2 = p.others(point2=point2)

830

831 _ = p.ellipsoids(p1)

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

833

834 # get the azimuthal equidistant projection

835 A = _Equidistant00(equidistant, p)

836

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

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

839 tol=tol, **LatLon_and_kwds)

840 return NearestOn2Tuple(r, f)

841

842

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

844 LatLon=None, **LatLon_kwds):

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

846 _LLS = _MODS.sphericalNvector.LatLon

847 _V3d = _MODS.vector3d.Vector3d

848 _vn2 = _MODS.vector3d._nearestOn2

849

850 E = p.ellipsoids(p2)

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

852

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

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

855 t = _LLS(p).nearestOn(_LLS(p1), _LLS(p2), within=within,

856 height=height)

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

858 if height is None:

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

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

861 h0 = favg(h1, h2)

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

863 h0 = h1 = h2 = _0_0

864

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

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

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

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

869 for i in range(1, _TRIPS):

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

871 # convert points to projection space and compute

872 # the nearest one (and its height) there

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

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

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

876 # convert nearest one back to geodetic

877 t, d = A._reverse2(v)

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

879 break

880 else:

881 raise e.degError()

882

883 if height is None:

884 h = v.z # nearest

885 elif _isHeight(height):

886 h = height

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

888 iteration=i, name=n)

889 return r, f, e # fraction or None

890

891

892__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI, intersecant2)

893

894# **) MIT License

895#

896# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

897#

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

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

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

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

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

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

904#

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

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

907#

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

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

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

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

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

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

914# OTHER DEALINGS IN THE SOFTWARE.