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, isscalar, issubclassof

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, Property_RO, \

15 property_RO, _TOL_M

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

17 _or, _ValueError, _xellipsoidal, _xError, _xkwds_not

18 from pygeodesy.fmath import favg, fmean_

19 from pygeodesy.fsums import Fmt, fsumf_

20 from pygeodesy.formy import opposing, _radical2

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

22 _exceed_PI_radians_, _low_, _near_, _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, property_RO # from .ellipsoidalBase

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

29 from pygeodesy.units import _fi_j2, Radius_, Scalar

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

31 _Wrap, wrap360

32

33 from math import degrees, radians

34

35__all__ = _ALL_LAZY.ellipsoidalBaseDI

36__version__ = '23.11.11'

37

38_polar__ = 'polar?'

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

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 (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: The B{C{other}} point is not C{LatLon}.

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, using this C{Direct} method. See

72 method L{destination2} for more details.

73

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

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

76 @kwarg height: Optional height, overriding the default

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

78

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

80 '''

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

82

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

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

85 azimuth) after having travelled for the given distance from

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

87 using this C{Direct} method.

88

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

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

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

92

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

94 are in compass C{degrees360}.

95

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

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

98

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

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

101 @kwarg height: Optional height, overriding the default

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

103

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

105 '''

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

107 return self._xnamed(r)

108

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

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

111

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

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

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

115 '''

116 g = self.geodesic

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

118 if LL:

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

120 return r

121

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

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

124 '''

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

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

127 datum=self.datum, name=self.name,

128 epoch=self.epoch, reframe=self.reframe))

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

130

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

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

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

134 L{distanceTo3} for more details.

135

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

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

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

139

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

141

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

143

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

145 ellipsoids are not compatible.

146 '''

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

148

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

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

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

152 C{Inverse} method.

153

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

155 conventionally meter. The distance is measured on the surface

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

157

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

159 are in compass C{degrees360} from North.

160

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

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

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

164

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

166

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

168

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

170 ellipsoids are not compatible.

171 '''

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

173

174 def finalBearingOn(self, distance, bearing):

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

176 travelled for the given distance along a geodesic given by

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

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

179

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

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

182

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

184 '''

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

186

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

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

189 travelled along a geodesic from this point to an other

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

191 L{distanceTo3} for more details.

192

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

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

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

196

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

198

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

200

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

202 ellipsoids are not compatible.

203 '''

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

205

206 @Property_RO

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

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

209 '''

210 return None # PYCHOK no cover

211

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

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

214 '''

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

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

217 gl = g._DirectLine( p, other) if isscalar(other) else \

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

219 return g, gl, p

220

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

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

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

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

225 for more details.

226

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

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

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

230

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

232

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

234

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

236 ellipsoids are not compatible.

237 '''

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

239

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

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

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

243 methods.

244

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

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

247 at this and 1.0 at the other point.

248 @kwarg height: Optional height, overriding the fractional

249 height (C{meter}).

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

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

252

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

254

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

256

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

258

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

260 ellipsoids are not compatible.

261

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

263 C{rhumbMidpointTo}.

264 '''

265 f = Scalar(fraction=fraction)

266 if isnear0(f):

267 r = self

268 elif isnear1(f) and not wrap:

269 r = self.others(other)

270 else: # negative fraction OK

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

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

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

274 return r

275

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

277 wrap=False, tol=_TOL):

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

279 two points or as a point and bearing.

280

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

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

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

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

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

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

287 method L{Ellipsoid.geodesic_}.

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

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

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

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

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

293

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

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

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

297

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

299

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

301 or B{C{other}} invalid.

302

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

304 B{C{height}}.

305

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

307 '''

308 try:

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

310 r = Radius_(circle=circle) if isscalar(circle) else \

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

312

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

314 form=_MODS.dms.F_DMS)

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

316 self.intersecant2)

317 except (TypeError, ValueError) as x:

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

319 exact=exact, wrap=wrap)

320

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

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

323

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

325 '''

326 _ = self.ellipsoids(other)

327 g = self.geodesic

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

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

330

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

332 equidistant=None, tol=_TOL_M):

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

334 to this point.

335

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

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

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

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

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

341 is taken into account for distances.

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

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

344

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

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

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

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

349

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

351

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

353

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

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

356

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

358 '''

359 _d3 = self.distanceTo3 # Distance3Tuple

360 _n3 = _nearestOn3

361 try:

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

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

364 _ = self.ellipsoids(p1)

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

366

367 except (TypeError, ValueError) as x:

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

369 closed=closed, height=height, wrap=wrap)

370

371 # get the azimuthal equidistant projection, once

372 A = _Equidistant00(equidistant, c)

373 b = _Box(c, c3.distance)

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

375

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

377 LatLon=self.classof, # this LatLon

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

379 try:

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

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

382 p2 = _unrollon(p1, p2)

383 # skip edge if no overlap with box around closest

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

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

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

387 if d3.distance < c3.distance:

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

389 b = _Box(c, c3.distance)

390 m = max(m, c.iteration)

391 p1, i = p2, j

392

393 except (TypeError, ValueError) as x:

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

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

396 closed=closed, height=height, wrap=wrap)

397

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

399

400 n = self.nearestOn8.__name__

401 c.rename(n)

402 if s is not c:

403 s = s.copy(name=n)

404 if e is not c:

405 e = e.copy(name=n)

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

407 iteration=m) # ._iteration for tests

408

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

410 wrap=False, tol=_TOL):

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

412 a geodesic from this point.

413

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

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

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

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

418 nsee method L{Ellipsoid.geodesic_}.

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

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

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

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

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

424

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

426

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

428

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

430 '''

431 try:

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

433

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

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

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

437 p._iteration = P.iteration

438 except (TypeError, ValueError) as x:

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

440 exact=exact, wrap=wrap)

441 return p

442

443

444 class _Box(object):

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

446

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

448 '''

449 def __init__(self, center, distance):

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

451

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

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

454 (C{meter}, conventionally)

455 '''

456 m = Radius_(distance=distance)

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

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

459 self._N = center.lat + d

460 self._S = center.lat - d

461 self._E = center.lon + d

462 self._W = center.lon - d

463

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

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

466

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

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

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

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

471

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

473 otherwise (C{bool}).

474 '''

475 if lat1 > lat2:

476 lat1, lat2 = lat2, lat1

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

478 return False

479 if lon1 > lon2:

480 lon1, lon2 = lon2, lon1

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

482 return False

483 return True

484

485

486 class _Tol(object):

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

488 '''

489 _deg = 0 # tol in degrees

490 _lat = 0

491 _m = 0 # tol in meter

492 _min = MAX # degrees

493 _prev = None

494 _r = 0

495

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

497 '''New L{_Tol}.

498

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

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

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

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

503 '''

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

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

506 self._m = max(EPS, tol_m)

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

508

509 @property_RO

510 def degrees(self):

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

512 '''

513 return self._deg

514

515 def degrees2m(self, deg):

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

517 '''

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

519

520 def degError(self, Error=_ValueError):

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

522 '''

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

524

525 def done(self, deg):

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

527 '''

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

529 return True

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

531 self._prev = deg

532 return False

533

534 @property_RO

535 def lat(self):

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

537 '''

538 return self._lat

539

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

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

542 '''

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

544 if m2km(m) > self.meter:

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

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

547

548 @property_RO

549 def meter(self):

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

551 '''

552 return self._m

553

554 @property_RO

555 def radius(self):

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

557 '''

558 return self._r

559

560 def reset(self):

561 '''Reset tolerances.

562 '''

563 self._min = MAX

564 self._prev = None

565

566

567 def _Equidistant00(equidistant, p1):

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

569 '''

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

571 equidistant = p1.Equidistant

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

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

574 raise _IsnotError(*t, equidistant=equidistant)

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

576

577

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

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

580 or as a point and (forward) bearing.

581

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

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

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

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

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

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

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

589

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

591

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

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

594

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

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

597 '''

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

599 raise _IsnotError(_ellipsoidal_, center=center)

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

601

602

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

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

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

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

607 '''

608 _LLS = _MODS.sphericalTrigonometry.LatLon

609 _si = _MODS.sphericalTrigonometry._intersect

610 _vi3 = _MODS.vector3d._intersect3d3

611

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

613 # compute opposing and distance

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

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

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

617 return b, t.distance

618

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

620 # compute an end point along the initial bearing

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

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

623 # radians in .sphericalTrigonometry._int3d2 and bearing

624 # comparison in .sphericalTrigonometry._intb

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

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

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

628 e = s.destination(d, e)

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

630

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

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

633 ll = not isscalar(e)

634 if ll:

635 e = s.others(**end)

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

637 e = _unrollon(s, e)

638 return e, ll

639

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

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

642

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

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

645 if not o: # intersection may be on start

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

647 return o

648 return -n if b else n

649

650 E = s1.ellipsoids(s2)

651

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

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

654

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

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

657

658 # get the azimuthal equidistant projection

659 A = _Equidistant00(equidistant, s1)

660

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

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

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

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

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

666

667 if not ll1:

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

669 if not ll2:

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

671

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

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

674 for i in range(1, _TRIPS):

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

676 # convert start and end points to projection

677 # space and compute an intersection there

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

679 eps=e.meter, useZ=False)

680 # convert intersection back to geodetic

681 t, d = A._reverse2(v)

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

683 t._iteration = i

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=t._iteration, 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 _latlon4(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 def _llS(ll): # return an C{_LLS} instance

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

750

751 r1 = Radius_(radius1=radius1)

752 r2 = Radius_(radius2=radius2)

753

754 E = c1.ellipsoids(c2)

755 # get the azimuthal equidistant projection

756 A = _Equidistant00(equidistant, c1)

757

758 if r1 < r2:

759 c1, c2 = c2, c1

760 r1, r2 = r2, r1

761

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

763 raise _ValueError(_exceed_PI_radians_)

764

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

766 c2 = _unrollon(c1, c2)

767

768 # distance between centers and radii are

769 # measured along the ellipsoid's surface

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

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

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

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

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

775

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

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

778

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

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

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

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

783

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

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

786 ts, ta = [], None

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

788 for i in range(1, _TRIPS):

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

790 # convert centers to projection space and

791 # compute the intersections there

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

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

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

795 sphere=False, too_d=m)

796 # convert intersections back to geodetic

797 if v1 is v2: # abutting

798 t, d = A._reverse2(v1)

799 else: # consider the closer intersection

800 t1, d1 = A._reverse2(v1)

801 t2, d2 = A._reverse2(v2)

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

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

804 t._iteration = i # _NamedTuple._iteration

805 ts.append(t)

806 if v1 is v2: # abutting

807 ta = t

808 break

809 else:

810 raise e.degError(Error=IntersectionError)

811 e.reset()

812

813 if ta: # abutting circles

814 pass # PYCHOK no cover

815 elif len(ts) == 2:

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

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

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

819 ta = ts[0] # assume abutting

820 else: # PYCHOK no cover

821 raise _AssertionError(ts=ts)

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

823 return r, r

824

825

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

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

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

829 separated to allow callers to embellish any exceptions.

830 '''

831 p1 = p.others(point1=point1)

832 p2 = p.others(point2=point2)

833

834 _ = p.ellipsoids(p1)

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

836

837 # get the azimuthal equidistant projection

838 A = _Equidistant00(equidistant, p)

839

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

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

842 tol=tol, **LatLon_and_kwds)

843 return NearestOn2Tuple(r, f)

844

845

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

847 LatLon=None, **LatLon_kwds):

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

849 _LLS = _MODS.sphericalNvector.LatLon

850 _V3d = _MODS.vector3d.Vector3d

851 _vn2 = _MODS.vector3d._nearestOn2

852

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

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

855

856 E = p.ellipsoids(p2)

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

858

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

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

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

862 height=height)

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

864 if height is None:

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

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

867 h0 = favg(h1, h2)

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

869 h0 = h1 = h2 = _0_0

870

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

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

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

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

875 for i in range(1, _TRIPS):

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

877 # convert points to projection space and compute

878 # the nearest one (and its height) there

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

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

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

882 # convert nearest one back to geodetic

883 t, d = A._reverse2(v)

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

885 t._iteration = i # _NamedTuple._iteration

886 break

887 else:

888 raise e.degError()

889

890 if height is None:

891 h = v.z # nearest

892 elif isscalar(height):

893 h = height

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

895 iteration=t.iteration, name=n)

896 return r, f, e # fraction or None

897

898

899__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI, intersecant2)

900

901# **) MIT License

902#

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

904#

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

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

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

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

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

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

911#

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

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

914#

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

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

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

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

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

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

921# OTHER DEALINGS IN THE SOFTWARE.