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

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

12 _EPSqrt as _TOL, _0_0, _0_01, _1_0, _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 from pygeodesy.namedTuples import Bearing2Tuple, Destination2Tuple, \

25 Intersection3Tuple, NearestOn2Tuple, \

26 NearestOn8Tuple, _LL4Tuple

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

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

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

30 Radius_, Scalar

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

32 _Wrap, wrap360

33

34 from math import degrees, radians

35

36__all__ = _ALL_LAZY.ellipsoidalBaseDI

37__version__ = '24.11.04'

38

39_polar__ = 'polar?'

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

41

42

43 class LatLonEllipsoidalBaseDI(LatLonEllipsoidalBase):

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

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

46 '''

47

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

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

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

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

52 for more details.

53

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

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

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

57

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

59

60 @raise TypeError: If B{C{other}} not this C{LatLon} class.

61

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

63 ellipsoids are not compatible.

64 '''

65 r = self._Inverse(other, wrap)

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

67

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

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

70 the given distance from this point along a geodesic given

71 by an initial bearing. See method L{destination2} for

72 further details.

73

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

75 '''

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

77

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

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

80 azimuth) after having travelled for the given distance from

81 this point along a geodesic (line) given by an initial bearing

82 at this point.

83

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

85 ellipsoid's axes, conventionally C{meter}. The distance is

86 measured on the surface of the ellipsoid, ignoring this point's

87 height.

88

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

90 are in compass C{degrees360}, clockwise from North.

91

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

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

94

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

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

97 @kwarg height: Optional height, overriding the default

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

99

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

101 '''

102 return self._Direct(distance, bearing, self.classof, height)

103

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

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

106

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

108 L{Destination3Tuple}C{(lat, lon, final)} if C{B{LL} is None}.

109 '''

110 g = self.geodesic

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

112 if LL:

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

114 return r

115

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

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

118 '''

119 h = self._heigHt(height)

120 d = _xkwds_not(None, datum=self.datum, name=self.name,

121 epoch=self.epoch, reframe=self.reframe)

122 d = LL(*_Wrap.latlon(r.lat, r.lon), height=h, **d)

123 return Destination2Tuple(d, wrap360(r.final), name=self.name)

124

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

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

127 a geodesic. See method L{distanceTo3} for more details.

128

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

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

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

132

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

134

135 @raise TypeError: If B{C{other}} not this C{LatLon} class.

136

137 @raise ValueError: This and the B{C{other}} point's L{Datum}

138 ellipsoids are incompatible.

139 '''

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

141

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

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

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

145 C{Inverse} method.

146

147 The distance is in the same units as this point's datum's

148 ellipsoid's axes, conventionally meter. The distance is

149 measured on the surface of the ellipsoid, ignoring this

150 point's height.

151

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

153 are in compass C{degrees360}, clockwise from North.

154

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

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

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

158

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

160

161 @raise TypeError: If B{C{other}} not this C{LatLon} class.

162

163 @raise ValueError: This and the B{C{other}} point's L{Datum}

164 ellipsoids are not compatible.

165 '''

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

167

168 def finalBearingOn(self, distance, bearing):

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

170 travelled for the given distance along a geodesic given

171 by an initial bearing from this point. See method

172 L{destination2} for more details.

173

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

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

176

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

178 '''

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

180

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

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

183 travelled along a geodesic from this point to an other

184 point. See method L{distanceTo3} for more details.

185

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

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

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

189

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

191

192 @raise TypeError: If B{C{other}} not this C{LatLon} class.

193

194 @raise ValueError: This and the B{C{other}} point's L{Datum}

195 ellipsoids are incompatible.

196 '''

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

198

199 @property_RO

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

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

202 '''

203 return None # PYCHOK no cover

204

205 def _g_gl_p3(self, start, end, exact, wrap):

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

207 '''

208 p = _unrollon(self, self.others(start=start), wrap=wrap)

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

210 gl = g._DirectLine( p, end) if _isDegrees(end) else \

211 g._InverseLine(p, self.others(end=end), wrap)

212 return g, gl, p

213

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

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

216 along a geodesic from this point to an other point. See

217 method L{distanceTo3} for more details.

218

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

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

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

222

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

224

225 @raise TypeError: If B{C{other}} not this C{LatLon} class.

226

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

228 ellipsoids are incompatible.

229 '''

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

231

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

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

234 this and an other point.

235

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

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

238 at this and 1.0 at the other point.

239 @kwarg height: Optional height, overriding the fractional

240 height (C{meter}).

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

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

243

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

245

246 @raise TypeError: If B{C{other}} not this C{LatLon} class.

247

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

249

250 @raise ValueError: This and the B{C{other}} point's L{Datum}

251 ellipsoids are incompatible.

252

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

254 C{rhumbMidpointTo}.

255 '''

256 f = Scalar(fraction=fraction)

257 if isnear0(f):

258 r = self

259 elif isnear1(f) and not wrap:

260 r = self.others(other)

261 else: # negative fraction OK

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

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

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

265 return r

266

267 def intersecant2(self, circle, start, end, exact=False, height=None, # PYCHOK signature

268 wrap=False, tol=_TOL):

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

270 points or as a point and a bearing from North.

271

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

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

274 @arg start: Start point of the geodesic (line) (this C{LatLon}).

275 @arg end: End point of the geodesic (line) (this C{LatLon}) or the initial

276 bearing at the B{C{start}} point (compass C{degrees360}).

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

278 method L{geodesic_<Ellipsoid.geodesic_>}.

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

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

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

282 B{C{start}} and/or B{C{end}} (C{bool}).

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

284

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

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

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

288

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

290

291 @raise TypeError: Invalid B{C{circle}}, B{C{start}} or B{C{end}}.

292

293 @raise UnitError: Invalid B{C{circle}}, B{C{end}}, B{C{exact}} or B{C{height}}.

294

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

296 '''

297 try:

298 g, gl, p = self._g_gl_p3(start, end, exact, wrap)

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

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

301

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

303 form=_MODS.dms.F_DMS)

304 return self._intersecend2(p, end, wrap, height, g, P, Q,

305 self.intersecant2)

306 except (TypeError, ValueError) as x:

307 raise _xError(x, center=self, circle=circle, start=start, end=end,

308 exact=exact, wrap=wrap)

309

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

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

312

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

314 '''

315 _ = self.ellipsoids(other)

316 g = self.geodesic

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

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

319

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

321 equidistant=None, tol=_TOL_M):

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

323 to this point.

324

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

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

327 @kwarg height: Optional height, overriding the height of this and all

328 other points (C{meter}, conventionally). If C{B{height}

329 is None}, each point's height is taken into account to

330 compute distances.

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

332 (C{bool}).

333 @kwarg equidistant: An azimuthal equidistant projection (I{class} or function

334 L{pygeodesy.equidistant}) or C{None} for the preferred

335 L{Equidistant<pygeodesy.ellipsoidalBase.Equidistant>}.

336 @kwarg tol: Convergence tolerance (C{meter}, conventionally).

337

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

339 initial, final)} with C{distance} in C{meter}, conventionally

340 and with the C{closest}, the C{start} the C{end} point each

341 an instance of this C{LatLon} class.

342

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

344

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

346

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

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

349

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

351 '''

352 _d3 = self.distanceTo3 # Distance3Tuple

353 _n3 = _nearestOn3

354 try:

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

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

357 _ = self.ellipsoids(p1)

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

359

360 except (TypeError, ValueError) as x:

361 raise _xError(x, Fmt.INDEX(points=0), p1, this=self, tol=tol,

362 closed=closed, height=height, wrap=wrap)

363

364 # get the azimuthal equidistant projection, once

365 A = _Equidistant00(equidistant, c)

366 b = _Box(c, c3.distance)

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

368

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

370 LatLon=self.classof, # this LatLon

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

372 try:

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

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

375 p2 = _unrollon(p1, p2)

376 # skip edge if no overlap with box around closest

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

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

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

380 if d3.distance < c3.distance:

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

382 b = _Box(c, c3.distance)

383 m = max(m, c.iteration)

384 p1, i = p2, j

385

386 except (TypeError, ValueError) as x:

387 raise _xError(x, Fmt.INDEX(points=i), p1,

388 Fmt.INDEX(points=j), p2, this=self, tol=tol,

389 closed=closed, height=height, wrap=wrap)

390

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

392

393 n = self.nearestOn8.__name__ # _DUNDER_nameof

394 c.rename(n)

395 if s is not c:

396 s = s.copy(name=n)

397 if e is not c:

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

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

400 iteration=m) # ._iteration for tests

401

402 def plumbTo(self, start, end, exact=False, height=None, # PYCHOK signature

403 wrap=False, tol=_TOL):

404 '''Compute the intersection of a geodesic from this point I{perpendicular} to

405 a geodesic (line) given as two points or as a point and a bearing from North.

406

407 @arg start: Start point of the geodesic (line) (this C{LatLon}).

408 @arg end: End point of the geodesic (line) (this C{LatLon}) or the initial

409 bearing at the B{C{start}} point (compass C{degrees360}).

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

411 see method L{Ellipsoid.geodesic_}.

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

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

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

415 and/or B{C{end}} point (C{bool}).

416 @kwarg tol: Convergence tolerance (C{meter}).

417

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

419

420 @raise TypeError: If B{C{start}} or B{C{end}} not this C{LatLon} class.

421

422 @raise UnitError: Invalid B{C{end}}, B{C{exact}} or B{C{height}}.

423 '''

424 try:

425 g, gl, p = self._g_gl_p3(start, end, exact, wrap)

426

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

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

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

430 p._iteration = P.iteration

431 except (TypeError, ValueError) as x:

432 raise _xError(x, plumb=self, start=start, end=end,

433 exact=exact, wrap=wrap)

434 return p

435

436

437 class _Box(object):

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

439

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

441 '''

442 _1_01 = _1_0 + _0_01 # ~1% margin

443

444 def __init__(self, center, distance):

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

446

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

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

449 (C{meter}, conventionally)

450 '''

451 m = Radius_(distance=distance)

452 E = center.ellipsoid()

453 d = E.m2degrees(m) * self._1_01

454 self._N = center.lat + d

455 self._S = center.lat - d

456 self._E = center.lon + d

457 self._W = center.lon - d

458

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

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

461

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

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

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

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

466

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

468 otherwise (C{bool}).

469 '''

470 if lat1 > lat2:

471 lat1, lat2 = lat2, lat1

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

473 return False

474 if lon1 > lon2:

475 lon1, lon2 = lon2, lon1

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

477 return False

478 return True

479

480

481 class _Tol(object):

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

483 '''

484 _deg = 0 # tol in degrees

485 _lat = 0

486 _m = 0 # tol in meter

487 _min = MAX # degrees

488 _prev = None

489 _r = 0

490

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

492 '''New L{_Tol}.

493

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

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

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

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

498 '''

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

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

501 self._m = max(EPS, tol_m)

502 self._deg = max(EPS, degrees(self._m / self._r)) # NOT E.m2degrees!

503

504 @property_RO

505 def degrees(self):

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

507 '''

508 return self._deg

509

510 def degrees2m(self, deg):

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

512 '''

513 return self.radius * radians(deg) / PI2 # NOT E.degrees2m!

514

515 def degError(self, Error=_ValueError):

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

517 '''

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

519

520 def done(self, deg):

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

522 '''

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

524 return True

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

526 self._prev = deg

527 return False

528

529 @property_RO

530 def lat(self):

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

532 '''

533 return self._lat

534

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

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

537 '''

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

539 if m2km(m) > self.meter:

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

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

542

543 @property_RO

544 def meter(self):

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

546 '''

547 return self._m

548

549 @property_RO

550 def radius(self):

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

552 '''

553 return self._r

554

555 def reset(self):

556 '''Reset tolerances.

557 '''

558 self._min = MAX # delattrof()

559 self._prev = None # delattrof()

560

561

562 def _Equidistant00(equidistant, p1):

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

564 '''

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

566 equidistant = p1.Equidistant

567 else:

568 _xsubclassof(*_MODS.azimuthal._Equidistants,

569 equidistant=equidistant)

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

571

572

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

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

575 or as a point and (forward) bearing.

576

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

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

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

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

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

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

583 @kwarg exact_height_wrap_tol: Optional keyword arguments C{B{exact}=False},

584 C{B{height}=None}, C{B{wrap}=False} and C{B{tol}}, see method

585 L{intersecant2<LatLonEllipsoidalBaseDI.intersecant2>}.

586

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

588

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

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

591

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

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

594 '''

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

596 raise _IsnotError(_ellipsoidal_, center=center)

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

598

599

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

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

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

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

604 '''

605 _LLS = _MODS.sphericalTrigonometry.LatLon

606 _si = _MODS.sphericalTrigonometry._intersect

607 _vi3 = _MODS.vector3d._intersect3d3

608

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

610 # compute opposing and distance

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

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

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

614 return b, t.distance

615

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

617 # compute an end point along the initial bearing about

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

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

620 # radians in .sphericalTrigonometry._int3d2 and bearing

621 # comparison in .sphericaltrigonometry._intb

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

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

624 d = min(max(d * _1_5, m), m * _3_0)

625 e = s.destination(d, e)

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

627

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

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

630 ll = not _isDegrees(e)

631 if ll:

632 e = s.others(**end)

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

634 e = _unrollon(s, e)

635 return e, ll

636

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

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

639 if not o: # intersection may be on start

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

641 return o

642 return -n if b else n

643

644 E = s1.ellipsoids(s2)

645

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

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

648

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

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

651

652 # get the azimuthal equidistant projection

653 A = _Equidistant00(equidistant, s1)

654

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

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

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

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

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

660

661 if not ll1:

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

663 if not ll2:

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

665

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

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

668 for i in range(1, _TRIPS):

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

670 # convert start and end points to projection

671 # space and compute an intersection there

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

673 eps=e.meter, useZ=False)

674 # convert intersection back to geodetic

675 t, d = A._reverse2(v)

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

677 break

678 else:

679 raise e.degError(Error=IntersectionError)

680

681 # like .sphericalTrigonometry._intersect, if this intersection

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

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

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

685 if b1:

686 t = t.antipodal()

687 b1 = not b1

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

689

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

691 iteration=i, name=n)

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

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

694

695

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

697 **equidistant_tol_LatLon_and_kwds):

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

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

700 point and an initial bearing from North.

701 '''

702 s1 = _xellipsoidal(start1=start1)

703 s2 = s1.others(start2=start2)

704 try:

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

706 **equidistant_tol_LatLon_and_kwds)

707 except (TypeError, ValueError) as x:

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

709

710

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

712 **equidistant_tol_LatLon_and_kwds):

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

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

715 '''

716 c1 = _xellipsoidal(center1=center1)

717 c2 = c1.others(center2=center2)

718 try:

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

720 **equidistant_tol_LatLon_and_kwds)

721 except (TypeError, ValueError) as x:

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

723 center2=center2, radius2=radius2)

724

725

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

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

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

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

730 '''

731 _LLS = _MODS.sphericalTrigonometry.LatLon

732 _si2 = _MODS.sphericalTrigonometry._intersects2

733 _vi2 = _MODS.vector3d._intersects2

734

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

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

737 iteration=t.iteration, name=n)

738

739 r1 = Radius_(radius1=radius1)

740 r2 = Radius_(radius2=radius2)

741

742 E = c1.ellipsoids(c2)

743 # get the azimuthal equidistant projection

744 A = _Equidistant00(equidistant, c1)

745

746 if r1 < r2:

747 c1, c2 = c2, c1

748 r1, r2 = r2, r1

749

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

751 raise _ValueError(_exceed_PI_radians_)

752

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

754 c2 = _unrollon(c1, c2)

755

756 # distance between centers and radii are

757 # measured along the ellipsoid's surface

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

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

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

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

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

763

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

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

766

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

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

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

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

771

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

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

774 ts, ta = [], None

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

776 for i in range(1, _TRIPS):

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

778 # convert centers to projection space and

779 # compute the intersections there

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

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

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

783 sphere=False, too_d=m)

784 # convert intersections back to geodetic

785 if v1 is v2: # abutting

786 t, d = A._reverse2(v1)

787 else: # consider the closer intersection

788 t1, d1 = A._reverse2(v1)

789 t2, d2 = A._reverse2(v2)

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

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

792 t._iteration = i # _NamedTuple._iteration

793 ts.append(t)

794 if v1 is v2: # abutting

795 ta = t

796 break

797 else:

798 raise e.degError(Error=IntersectionError)

799 e.reset()

800

801 if ta: # abutting circles

802 pass # PYCHOK no cover

803 elif len(ts) == 2:

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

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

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

807 ta = ts[0] # assume abutting

808 else: # PYCHOK no cover

809 raise _AssertionError(ts=ts)

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

811 return r, r

812

813

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

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

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

817 separated to allow callers to embellish any exceptions.

818 '''

819 p1 = p.others(point1=point1)

820 p2 = p.others(point2=point2)

821

822 _ = p.ellipsoids(p1)

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

824

825 # get the azimuthal equidistant projection

826 A = _Equidistant00(equidistant, p)

827

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

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

830 tol=tol, **LatLon_and_kwds)

831 return NearestOn2Tuple(r, f)

832

833

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

835 LatLon=None, **LatLon_kwds):

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

837 _LLS = _MODS.sphericalNvector.LatLon

838 _V3d = _MODS.vector3d.Vector3d

839 _vn2 = _MODS.vector3d._nearestOn2

840

841 E = p.ellipsoids(p2)

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

843

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

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

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

847 height=height)

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

849 if height is None:

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

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

852 h0 = favg(h1, h2)

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

854 h0 = h1 = h2 = _0_0

855

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

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

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

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

860 for i in range(1, _TRIPS):

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

862 # convert points to projection space and compute

863 # the nearest one (and its height) there

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

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

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

867 # convert nearest one back to geodetic

868 t, d = A._reverse2(v)

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

870 break

871 else:

872 raise e.degError()

873

874 if height is None:

875 h = v.z # nearest

876 elif _isHeight(height):

877 h = height

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

879 iteration=i, name=n)

880 return r, f, e # fraction or None

881

882

883__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI, intersecant2)

884 del _1_0, _0_01

885

886# **) MIT License

887#

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

889#

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

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

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

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

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

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

896#

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

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

899#

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

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

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

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

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

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

906# OTHER DEALINGS IN THE SOFTWARE.