1

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

3

4 u'''(INTERNAL) Private spherical base classes C{CartesianSphericalBase} and

5 C{LatLonSphericalBase} for L{sphericalNvector} and L{sphericalTrigonometry}.

6

7 A pure Python implementation of geodetic (lat-/longitude) functions,

8 transcoded in part from JavaScript originals by I{(C) Chris Veness 2011-2016}

9 and published under the same MIT Licence**, see

10 U{Latitude/Longitude<https://www.Movable-Type.co.UK/scripts/latlong.html>}.

11'''

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

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

14

15 from pygeodesy.basics import isbool, isinstanceof, map1

16 from pygeodesy.cartesianBase import Bearing2Tuple, CartesianBase

17 from pygeodesy.constants import EPS, PI, PI2, PI_2, R_M, R_MA, \

18 _umod_360, isnear0, isnon0, _0_0, \

19 _0_5, _1_0, _180_0

20 from pygeodesy.datums import Datums, _spherical_datum

21 from pygeodesy.errors import IntersectionError, _ValueError, _xError

22 from pygeodesy.fmath import favg, fdot, hypot, sqrt_a

23 from pygeodesy.interns import NN, _COMMA_, _concentric_, _datum_, \

24 _distant_, _exceed_PI_radians_, _name_, \

25 _near_, _radius_, _too_

26 from pygeodesy.latlonBase import LatLonBase, _trilaterate5 # PYCHOK passed

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

28# from pygeodesy.namedTuples import Bearing2Tuple # from .cartesianBase

29 from pygeodesy.nvectorBase import NvectorBase, Fmt, _xattrs

30 from pygeodesy.props import deprecated_method, property_doc_, \

31 property_RO, _update_all

32# from pygeodesy.streprs import Fmt, _xattrs # from .nvectorBase

33 from pygeodesy.units import Bearing, Bearing_, Radians_, Radius, \

34 Radius_, Scalar_

35 from pygeodesy.utily import acos1, atan2b, atan2d, degrees90, \

36 degrees180, sincos2, sincos2d, tanPI_2_2, \

37 _unrollon, wrap360, wrapPI

38

39 from math import cos, fabs, log, sin, sqrt

40

41__all__ = _ALL_LAZY.sphericalBase

42__version__ = '23.05.06'

43

44

45 def _angular(distance, radius, low=EPS): # PYCHOK in .spherical*

46 '''(INTERNAL) Return the angular distance in C{radians}.

47

48 @raise UnitError: Invalid B{C{distance}} or B{C{radius}}.

49 '''

50 r = _1_0 if radius is None else Radius_(radius=radius)

51 return Radians_(distance / r, low=low)

52

53

54 def _rads3(rad1, rad2, radius): # in .sphericalTrigonometry

55 '''(INTERNAL) Convert radii to radians.

56 '''

57 r1 = Radius_(rad1=rad1)

58 r2 = Radius_(rad2=rad2)

59 if radius is not None: # convert radii to radians

60 r1 = _angular(r1, radius)

61 r2 = _angular(r2, radius)

62

63 x = r1 < r2

64 if x:

65 r1, r2 = r2, r1

66 if r1 > PI:

67 raise IntersectionError(rad1=rad1, rad2=rad2,

68 txt=_exceed_PI_radians_)

69 return r1, r2, x

70

71

72 class CartesianSphericalBase(CartesianBase):

73 '''(INTERNAL) Base class for spherical C{Cartesian}s.

74 '''

75 _datum = Datums.Sphere # L{Datum}

76

77 def intersections2(self, rad1, other, rad2, radius=R_M):

78 '''Compute the intersection points of two circles each defined

79 by a center point and a radius.

80

81 @arg rad1: Radius of the this circle (C{meter} or C{radians},

82 see B{C{radius}}).

83 @arg other: Center of the other circle (C{Cartesian}).

84 @arg rad2: Radius of the other circle (C{meter} or C{radians},

85 see B{C{radius}}).

86 @kwarg radius: Mean earth radius (C{meter} or C{None} if both

87 B{C{rad1}} and B{C{rad2}} are given in C{radians}).

88

89 @return: 2-Tuple of the intersection points, each C{Cartesian}.

90 For abutting circles, the intersection points are the

91 same C{Cartesian} instance, aka the I{radical center}.

92

93 @raise IntersectionError: Concentric, antipodal, invalid or

94 non-intersecting circles.

95

96 @raise TypeError: If B{C{other}} is not C{Cartesian}.

97

98 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}} or B{C{radius}}.

99

100 @see: U{Calculating intersection of two Circles

101 <https://GIS.StackExchange.com/questions/48937/

102 calculating-intersection-of-two-circles>} and method

103 or function C{trilaterate3d2}.

104 '''

105 x1, x2 = self, self.others(other)

106 r1, r2, x = _rads3(rad1, rad2, radius)

107 if x:

108 x1, x2 = x2, x1

109 try:

110 n, q = x1.cross(x2), x1.dot(x2)

111 n2, q1 = n.length2, (_1_0 - q**2)

112 if n2 < EPS or isnear0(q1):

113 raise ValueError(_near_(_concentric_))

114 c1, c2 = cos(r1), cos(r2)

115 x0 = x1.times((c1 - q * c2) / q1).plus(

116 x2.times((c2 - q * c1) / q1))

117 n1 = _1_0 - x0.length2

118 if n1 < EPS:

119 raise ValueError(_too_(_distant_))

120 except ValueError as x:

121 raise IntersectionError(center=self, rad1=rad1,

122 other=other, rad2=rad2, cause=x)

123 n = n.times(sqrt(n1 / n2))

124 if n.length > EPS:

125 x1 = x0.plus(n)

126 x2 = x0.minus(n)

127 else: # abutting circles

128 x1 = x2 = x0

129

130 return (_xattrs(x1, self, _datum_, _name_),

131 _xattrs(x2, self, _datum_, _name_))

132

133

134 class LatLonSphericalBase(LatLonBase):

135 '''(INTERNAL) Base class for spherical C{LatLon}s.

136 '''

137 _datum = Datums.Sphere # spherical L{Datum}

138

139 def __init__(self, latlonh, lon=None, height=0, datum=None, wrap=False, name=NN):

140 '''Create a spherical C{LatLon} point frome the given lat-, longitude and

141 height on the given datum.

142

143 @arg latlonh: Latitude (C{degrees} or DMS C{str} with N or S suffix) or

144 a previous C{LatLon} instance provided C{B{lon}=None}.

145 @kwarg lon: Longitude (C{degrees} or DMS C{str} with E or W suffix) or

146 C(None), indicating B{C{latlonh}} is a C{LatLon}.

147 @kwarg height: Optional height above (or below) the earth surface (C{meter},

148 same units as the datum's ellipsoid axes or radius).

149 @kwarg datum: Optional, spherical datum to use (L{Datum}, L{Ellipsoid},

150 L{Ellipsoid2}, L{a_f2Tuple}) or earth radius in C{meter},

151 conventionally).

152 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{lat}} and B{C{lon}}

153 (C{bool}).

154 @kwarg name: Optional name (string).

155

156 @raise TypeError: If B{C{latlonh}} is not a C{LatLon} or B{C{datum}} not

157 spherical.

158 '''

159 LatLonBase.__init__(self, latlonh, lon=lon, height=height, wrap=wrap, name=name)

160 if datum not in (None, self.datum):

161 self.datum = datum

162

163 def bearingTo2(self, other, wrap=False, raiser=False):

164 '''Return the initial and final bearing (forward and reverse

165 azimuth) from this to an other point.

166

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

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

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

170

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

172

173 @raise TypeError: The B{C{other}} point is not spherical.

174

175 @see: Methods C{initialBearingTo} and C{finalBearingTo}.

176 '''

177 # .initialBearingTo is inside .-Nvector and .-Trigonometry

178 i = self.initialBearingTo(other, wrap=wrap, raiser=raiser) # PYCHOK .initialBearingTo

179 f = self.finalBearingTo( other, wrap=wrap, raiser=raiser)

180 return Bearing2Tuple(i, f, name=self.name)

181

182 @property_doc_(''' this point's datum (L{Datum}).''')

183 def datum(self):

184 '''Get this point's datum (L{Datum}).

185 '''

186 return self._datum

187

188 @datum.setter # PYCHOK setter!

189 def datum(self, datum):

190 '''Set this point's datum I{without conversion} (L{Datum}, L{Ellipsoid},

191 L{Ellipsoid2}, L{a_f2Tuple}) or C{scalar} spherical earth radius).

192

193 @raise TypeError: If B{C{datum}} invalid or not not spherical.

194 '''

195 d = _spherical_datum(datum, name=self.name, raiser=_datum_)

196 if self._datum != d:

197 _update_all(self)

198 self._datum = d

199

200 def finalBearingTo(self, other, wrap=False, raiser=False):

201 '''Return the final bearing (reverse azimuth) from this to

202 an other point.

203

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

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

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

207

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

209

210 @raise TypeError: The B{C{other}} point is not spherical.

211

212 @example:

213

214 >>> p = LatLon(52.205, 0.119)

215 >>> q = LatLon(48.857, 2.351)

216 >>> b = p.finalBearingTo(q) # 157.9

217 '''

218 p = self.others(other)

219 if wrap:

220 p = _unrollon(self, p, wrap=wrap)

221 # final bearing is the reverse of the other, initial one;

222 # .initialBearingTo is inside .-Nvector and .-Trigonometry

223 b = p.initialBearingTo(self, wrap=False, raiser=raiser)

224 return _umod_360(b + _180_0)

225

226 def intersecant2(self, circle, point, bearing, radius=R_M, exact=False,

227 height=None, wrap=False): # was=True

228 '''Compute the intersections of a circle and a line.

229

230 @arg circle: Radius of the circle centered at this location

231 (C{meter}, same units as B{C{radius}}) or a point

232 on the circle (this C{LatLon}).

233 @arg point: An other point in- or outside the circle (this C{LatLon}).

234 @arg bearing: Bearing at the B{C{point}} (compass C{degrees360})

235 or a second point on the line (this C{LatLon}).

236 @kwarg radius: Mean earth radius (C{meter}, conventionally).

237 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth,

238 destination and distance, if C{False} use the basic

239 rhumb methods (C{bool}) or if C{None} use the I{great

240 circle} methods.

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

242 conventionally) or C{None}.

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

244 B{C{point}}, B{C{circle}} and/or B{C{bearing}} (C{bool}).

245

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

247 each an instance of this class. For a tangent line, each

248 point C{is} this very instance.

249

250 @raise IntersectionError: The circle and line do not intersect.

251

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

253 or B{C{bearing}} invalid.

254

255 @raise ValueError: Invalid B{C{circle}}, B{C{bearing}}, B{C{radius}},

256 B{C{exact}} or B{C{height}}.

257 '''

258 p = self.others(point=point)

259 try:

260 return _intersecant2(self, circle, p, bearing, radius=radius, exact=exact,

261 height=height, wrap=wrap)

262 except (TypeError, ValueError) as x:

263 raise _xError(x, center=self, circle=circle, point=point, bearing=bearing,

264 exact=exact, wrap=wrap)

265

266 def maxLat(self, bearing):

267 '''Return the maximum latitude reached when travelling

268 on a great circle on given bearing from this point

269 based on Clairaut's formula.

270

271 The maximum latitude is independent of longitude

272 and the same for all points on a given latitude.

273

274 Negate the result for the minimum latitude (on the

275 Southern hemisphere).

276

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

278

279 @return: Maximum latitude (C{degrees90}).

280

281 @raise ValueError: Invalid B{C{bearing}}.

282 '''

283 m = acos1(fabs(sin(Bearing_(bearing)) * cos(self.phi)))

284 return degrees90(m)

285

286 def minLat(self, bearing):

287 '''Return the minimum latitude reached when travelling

288 on a great circle on given bearing from this point.

289

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

291

292 @return: Minimum latitude (C{degrees90}).

293

294 @see: Method L{maxLat} for more details.

295

296 @raise ValueError: Invalid B{C{bearing}}.

297 '''

298 return -self.maxLat(bearing)

299

300 def parse(self, strllh, height=0, sep=_COMMA_, name=NN):

301 '''Parse a string representing a similar, spherical C{LatLon}

302 point, consisting of C{"lat, lon[, height]"}.

303

304 @arg strllh: Lat, lon and optional height (C{str}),

305 see function L{pygeodesy.parse3llh}.

306 @kwarg height: Optional, default height (C{meter}).

307 @kwarg sep: Optional separator (C{str}).

308 @kwarg name: Optional instance name (C{str}),

309 overriding this name.

310

311 @return: The similar point (spherical C{LatLon}).

312

313 @raise ParseError: Invalid B{C{strllh}}.

314 '''

315 t = _MODS.dms.parse3llh(strllh, height=height, sep=sep)

316 r = self.classof(*t)

317 if name:

318 r.rename(name)

319 return r

320

321 @property_RO

322 def _radius(self):

323 '''(INTERNAL) Get this sphere's radius.

324 '''

325 return self.datum.ellipsoid.equatoradius

326

327 def _rhumbs3(self, other, wrap, r=False): # != .latlonBase._rhumbx3

328 '''(INTERNAL) Rhumb_ helper function.

329

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

331 '''

332 p = self.others(other, up=2)

333 if wrap:

334 p = _unrollon(self, p, wrap=wrap)

335 a2, b2 = p.philam

336 a1, b1 = self.philam

337 # if |db| > 180 take shorter rhumb

338 # line across the anti-meridian

339 db = wrapPI(b2 - b1)

340 dp = log(tanPI_2_2(a2) / tanPI_2_2(a1))

341 da = a2 - a1

342 if r:

343 # on Mercator projection, longitude distances shrink

344 # by latitude; the 'stretch factor' q becomes ill-

345 # conditioned along E-W line (0/0); use an empirical

346 # tolerance to avoid it

347 q = (da / dp) if fabs(dp) > EPS else cos(self.phi)

348 da = hypot(da, q * db) # angular distance radians

349 return da, db, dp

350

351 def rhumbAzimuthTo(self, other, radius=R_M, exact=False, wrap=False):

352 '''Return the azimuth (bearing) of a rhumb line (loxodrome)

353 between this and an other (spherical) point.

354

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

356 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},

357 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).

358 @kwarg exact: If C{True}, use class L{Rhumb} (C{bool}), default

359 C{False} for backward compatibility.

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

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

362

363 @return: Rhumb line azimuth (compass C{degrees180}).

364

365 @raise TypeError: The B{C{other}} point is incompatible or

366 B{C{radius}} is invalid.

367

368 @example:

369

370 >>> p = LatLon(51.127, 1.338)

371 >>> q = LatLon(50.964, 1.853)

372 >>> b = p.rhumbBearingTo(q) # 116.7

373 '''

374 if exact: # use series, always

375 z = LatLonBase.rhumbAzimuthTo(self, other, exact=False,

376 radius=radius, wrap=wrap)

377 else:

378 _, db, dp = self._rhumbs3(other, wrap)

379 z = atan2d(db, dp) # see .rhumbx.Rhumb.Inverse

380 return z

381

382 @deprecated_method

383 def rhumbBearingTo(self, other): # unwrapped

384 '''DEPRECATED, use method C{.rhumbAzimuthTo}.'''

385 return wrap360(self.rhumbAzimuthTo(other)) # [0..360)

386

387 def rhumbDestination(self, distance, bearing, radius=R_M, height=None, exact=False):

388 '''Return the destination point having travelled the given distance

389 from this point along a rhumb line (loxodrome) at the given bearing.

390

391 @arg distance: Distance travelled (C{meter}, same units as B{C{radius}}),

392 may be negative if C{B{exact}=True}.

393 @arg bearing: Bearing (azimuth) at this point (compass C{degrees360}).

394 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},

395 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}) if

396 C{B{exact}=True}.

397 @kwarg height: Optional height, overriding the default height

398 (C{meter}, same unit as B{C{radius}}).

399 @kwarg exact: If C{True}, use class L{Rhumb} (C{bool}), default

400 C{False} for backward compatibility.

401

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

403

404 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}},

405 B{C{radius}} or B{C{height}}.

406

407 @example:

408

409 >>> p = LatLon(51.127, 1.338)

410 >>> q = p.rhumbDestination(40300, 116.7) # 50.9642°N, 001.8530°E

411 '''

412 if exact: # use series, always

413 r = LatLonBase.rhumbDestination(self, distance, bearing, exact=False,

414 radius=radius, height=height)

415 else: # radius=None from .rhumbMidpointTo

416 if radius in (None, self._radius):

417 d, r = self.datum, radius

418 else:

419 d = _spherical_datum(radius, raiser=_radius_) # spherical only

420 r = d.ellipsoid.equatoradius

421 r = _angular(distance, r, low=_0_0) # distance=0 from .rhumbMidpointTo

422

423 a1, b1 = self.philam

424 sb, cb = sincos2(Bearing_(bearing))

425

426 da = r * cb

427 a2 = a1 + da

428 # normalize latitude if past pole

429 if a2 > PI_2:

430 a2 = PI - a2

431 elif a2 < -PI_2:

432 a2 = -PI - a2

433

434 dp = log(tanPI_2_2(a2) / tanPI_2_2(a1))

435 # q becomes ill-conditioned on E-W course 0/0

436 q = (da / dp) if fabs(dp) > EPS else cos(a1)

437 b2 = (b1 + r * sb / q) if fabs(q) > EPS else b1

438

439 h = self._heigHt(height)

440 r = self.classof(degrees90(a2), degrees180(b2), datum=d, height=h)

441 return r

442

443 def rhumbDistanceTo(self, other, radius=R_M, exact=False, wrap=False):

444 '''Return the distance from this to an other point along

445 a rhumb line (loxodrome).

446

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

448 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},

449 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}) if

450 C{B{exact}=True}.

451 @kwarg exact: If C{True}, use class L{Rhumb} (C{bool}), default

452 C{False} for backward compatibility.

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

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

455

456 @return: Distance (C{meter}, the same units as B{C{radius}}

457 or C{radians} if B{C{radius}} is C{None}).

458

459 @raise TypeError: The B{C{other}} point is incompatible.

460

461 @raise ValueError: Invalid B{C{radius}}.

462

463 @example:

464

465 >>> p = LatLon(51.127, 1.338)

466 >>> q = LatLon(50.964, 1.853)

467 >>> d = p.rhumbDistanceTo(q) # 403100

468 '''

469 if exact: # use series, always

470 r = LatLonBase.rhumbDistanceTo(self, other, exact=False,

471 radius=radius, wrap=wrap)

472 if radius is None: # angular distance in radians

473 r = r / self._radius # /= chokes PyChecker

474 else:

475 # see <https://www.EdWilliams.org/avform.htm#Rhumb>

476 r, _, _ = self._rhumbs3(other, wrap, r=True)

477 if radius is not None:

478 r *= Radius(radius)

479 return r

480

481 def rhumbMidpointTo(self, other, height=None, radius=R_M, exact=False,

482 fraction=_0_5, wrap=False):

483 '''Return the (loxodromic) midpoint on the rhumb line between

484 this and an other point.

485

486 @arg other: The other point (spherical LatLon).

487 @kwarg height: Optional height, overriding the mean height

488 (C{meter}).

489 @kwarg radius: Optional mean earth radius (C{meter}),

490 overriding the default C{R_M}.

491 @kwarg radius: Earth radius (C{meter}) or earth model (L{Datum},

492 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).

493 @kwarg exact: If C{True}, use class L{Rhumb} (C{bool}), default

494 C{False} for backward compatibility.

495 @kwarg fraction: Midpoint location from this point (C{scalar}),

496 may negative if C{B{exact}=True}.

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

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

499

500 @return: The (mid)point at the given B{C{fraction}} along

501 the rhumb line (spherical C{LatLon}).

502

503 @raise TypeError: The B{C{other}} point is incompatible.

504

505 @raise ValueError: Invalid B{C{height}} or B{C{fraction}}

506

507 @example:

508

509 >>> p = LatLon(51.127, 1.338)

510 >>> q = LatLon(50.964, 1.853)

511 >>> m = p.rhumb_midpointTo(q)

512 >>> m.toStr() # '51.0455°N, 001.5957°E'

513 '''

514 if exact: # use series, always

515 r = LatLonBase.rhumbMidpointTo(self, other, exact=False,

516 radius=radius, height=height,

517 fraction=fraction, wrap=wrap)

518 elif fraction is not _0_5:

519 f = Scalar_(fraction=fraction) # low=_0_0

520 r, db, dp = self._rhumbs3(other, wrap, r=True) # radians

521 z = atan2b(db, dp)

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

523 r = self.rhumbDestination(r * f, z, radius=None, height=h)

524

525 else: # for backward compatibility, unwrapped

526 # see <https://MathForum.org/library/drmath/view/51822.html>

527 a1, b1 = self.philam

528 a2, b2 = self.others(other).philam

529

530 if fabs(b2 - b1) > PI:

531 b1 += PI2 # crossing anti-meridian

532

533 a3 = favg(a1, a2)

534 b3 = favg(b1, b2)

535

536 f1 = tanPI_2_2(a1)

537 if isnon0(f1):

538 f2 = tanPI_2_2(a2)

539 f = f2 / f1

540 if isnon0(f):

541 f = log(f)

542 if isnon0(f):

543 f3 = tanPI_2_2(a3)

544 b3 = fdot(map1(log, f1, f2, f3),

545 -b2, b1, b2 - b1) / f

546

547 d = self.datum if radius in (None, self._radius) else \

548 _spherical_datum(radius, name=self.name, raiser=_radius_)

549 h = self._havg(other, h=height)

550 r = self.classof(degrees90(a3), degrees180(b3), datum=d, height=h)

551 return r

552

553 def toNvector(self, Nvector=NvectorBase, **Nvector_kwds): # PYCHOK signature

554 '''Convert this point to C{Nvector} components, I{including

555 height}.

556

557 @kwarg Nvector_kwds: Optional, additional B{C{Nvector}}

558 keyword arguments, ignored if

559 C{B{Nvector} is None}.

560

561 @return: An B{C{Nvector}} or a L{Vector4Tuple}C{(x, y, z, h)}

562 if B{C{Nvector}} is C{None}.

563

564 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}.

565 '''

566 return LatLonBase.toNvector(self, Nvector=Nvector, **Nvector_kwds)

567

568 def toWm(self, radius=R_MA):

569 '''Convert this point to a I{WM} coordinate.

570

571 @kwarg radius: Optional earth radius (C{meter}).

572

573 @return: The WM coordinate (L{Wm}).

574

575 @see: Function L{pygeodesy.toWm} in module L{webmercator} for details.

576 '''

577 return _MODS.webmercator.toWm(self, radius=radius)

578

579

580 def _intersecant2(c, r, p, b, radius=R_M, exact=False,

581 height=None, wrap=False):

582 # (INTERNAL) Intersect a circle and bearing, see L{intersecant2}

583 # above, separated to allow callers to embellish any exceptions

584

585 if wrap:

586 p = _unrollon(c, p, wrap=wrap)

587 nonexact = exact is None

588

589 if not isinstanceof(r, c.__class__, p.__class__):

590 r = Radius_(circle=r)

591 elif nonexact:

592 r = c.distanceTo(r, radius=radius, wrap=wrap)

593 elif isbool(exact):

594 r = c.rhumbDistanceTo(r, radius=radius, exact=exact, wrap=wrap)

595 else:

596 raise _ValueError(exact=exact)

597

598 if not isinstanceof(b, c.__class__, p.__class__):

599 b = Bearing(b)

600 elif nonexact:

601 b = p.initialBearingTo(b, wrap=wrap)

602 else:

603 b = p.rhumbAzimuthTo(b, radius=radius, exact=exact, wrap=wrap)

604

605 d = p.distanceTo(c, radius=radius) if nonexact else \

606 p.rhumbDistanceTo(c, radius=radius, exact=exact)

607 if d > EPS:

608 a = p.initialBearingTo(c) if nonexact else wrap360(

609 p.rhumbAzimuthTo(c, radius=radius, exact=exact))

610 s, c = sincos2d(b - a)

611 s = sqrt_a(r, fabs(s * d))

612 if s > r:

613 raise IntersectionError(_too_(Fmt.distant(s)))

614 elif (r - s) < EPS:

615 return p, p # tangent

616 c *= d

617 else: # coincindent

618 s, c = r, 0

619

620 a = b + _180_0

621 if nonexact:

622 b = p.destination(s + c, b, radius=radius, height=height)

623 a = p.destination(s - c, a, radius=radius, height=height)

624 else:

625 b = p.rhumbDestination(s + c, b, radius=radius, height=height, exact=exact)

626 a = p.rhumbDestination(s - c, a, radius=radius, height=height, exact=exact)

627 return b, a # in bearing direction first

628

629

630__all__ += _ALL_DOCS(CartesianSphericalBase, LatLonSphericalBase)

631

632# **) MIT License

633#

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

635#

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

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

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

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

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

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

642#

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

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

645#

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

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

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

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

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

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

652# OTHER DEALINGS IN THE SOFTWARE.