Coverage for pygeodesy/sphericalBase.py: 92%

211 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-05 15:46 -0400

1 

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

3 

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

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

6 

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

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

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

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

11''' 

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

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

14 

15from pygeodesy.basics import isbool, isinstanceof, map1 

16from pygeodesy.cartesianBase import Bearing2Tuple, CartesianBase 

17from 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 

20from pygeodesy.datums import Datums, _spherical_datum 

21from pygeodesy.errors import IntersectionError, _ValueError, _xError 

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

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

24 _distant_, _exceed_PI_radians_, _name_, \ 

25 _near_, _radius_, _too_ 

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

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

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

29from pygeodesy.nvectorBase import NvectorBase, Fmt, _xattrs 

30from pygeodesy.props import deprecated_method, property_doc_, \ 

31 property_RO, _update_all 

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

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

34 Radius, Radius_, Scalar_ 

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

36 sincos2, sincos2d, tanPI_2_2, wrap360, wrapPI 

37 

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

39 

40__all__ = _ALL_LAZY.sphericalBase 

41__version__ = '23.03.22' 

42 

43 

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

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

46 

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

48 ''' 

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

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

51 

52 

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

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

55 ''' 

56 r1 = Radius_(rad1=rad1) 

57 r2 = Radius_(rad2=rad2) 

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

59 r1 = _angular(r1, radius) 

60 r2 = _angular(r2, radius) 

61 

62 x = r1 < r2 

63 if x: 

64 r1, r2 = r2, r1 

65 if r1 > PI: 

66 raise IntersectionError(rad1=rad1, rad2=rad2, 

67 txt=_exceed_PI_radians_) 

68 return r1, r2, x 

69 

70 

71class CartesianSphericalBase(CartesianBase): 

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

73 ''' 

74 _datum = Datums.Sphere # L{Datum} 

75 

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

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

78 by a center point and a radius. 

79 

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

81 see B{C{radius}}). 

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

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

84 see B{C{radius}}). 

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

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

87 

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

89 The intersection points are the same C{Cartesian} 

90 instance for abutting circles, aka I{radical center}. 

91 

92 @raise IntersectionError: Concentric, antipodal, invalid or 

93 non-intersecting circles. 

94 

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

96 

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

98 

99 @see: U{Calculating intersection of two Circles 

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

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

102 or function C{trilaterate3d2}. 

103 ''' 

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

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

106 if x: 

107 x1, x2 = x2, x1 

108 try: 

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

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

111 if n2 < EPS or isnear0(q1): 

112 raise ValueError(_near_(_concentric_)) 

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

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

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

116 n1 = _1_0 - x0.length2 

117 if n1 < EPS: 

118 raise ValueError(_too_(_distant_)) 

119 except ValueError as x: 

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

121 other=other, rad2=rad2, cause=x) 

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

123 if n.length > EPS: 

124 x1 = x0.plus(n) 

125 x2 = x0.minus(n) 

126 else: # abutting circles 

127 x1 = x2 = x0 

128 

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

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

131 

132 

133class LatLonSphericalBase(LatLonBase): 

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

135 ''' 

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

137 

138 def __init__(self, lat, lon, height=0, datum=None, name=NN): 

139 '''Create a spherical C{LatLon} point frome the given 

140 lat-, longitude and height on the given datum. 

141 

142 @arg lat: Latitude (C{degrees} or DMS C{[N|S]}). 

143 @arg lon: Longitude (C{degrees} or DMS C{str[E|W]}). 

144 @kwarg height: Optional elevation (C{meter}, the same units 

145 as the datum's half-axes). 

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

147 L{Ellipsoid}, L{Ellipsoid2}, L{a_f2Tuple}) 

148 or C{scalar} earth radius). 

149 @kwarg name: Optional name (string). 

150 

151 @raise TypeError: If B{C{datum}} invalid or not 

152 not spherical. 

153 

154 @example: 

155 

156 >>> p = LatLon(51.4778, -0.0016) # height=0, datum=Datums.WGS84 

157 ''' 

158 LatLonBase.__init__(self, lat, lon, height=height, name=name) 

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

160 self.datum = datum 

161 

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

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

164 azimuth) from this to an other point. 

165 

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

167 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

168 @kwarg raiser: Optionally, raise L{CrossError} (C{bool}). 

169 

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

171 

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

173 

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

175 ''' 

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

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

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

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

180 

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

182 def datum(self): 

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

184 ''' 

185 return self._datum 

186 

187 @datum.setter # PYCHOK setter! 

188 def datum(self, datum): 

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

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

191 

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

193 ''' 

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

195 if self._datum != d: 

196 _update_all(self) 

197 self._datum = d 

198 

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

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

201 an other point. 

202 

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

204 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

205 @kwarg raiser: Optionally, raise L{CrossError} (C{bool}). 

206 

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

208 

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

210 

211 @example: 

212 

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

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

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

216 ''' 

217 self.others(other) 

218 

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

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

221 b = other.initialBearingTo(self, wrap=wrap, raiser=raiser) 

222 return _umod_360(b + _180_0) 

223 

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

225 height=None, wrap=True): 

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

227 

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

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

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

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

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

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

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

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

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

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

238 circle} methods. 

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

240 conventionally) or C{None}. 

241 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

242 

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

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

245 point C{is} this very instance. 

246 

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

248 

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

250 or B{C{bearing}} invalid. 

251 

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

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

254 ''' 

255 p = self.others(point=point) 

256 try: 

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

258 height=height, wrap=wrap) 

259 except (TypeError, ValueError) as x: 

260 raise _xError(x, center=self, circle=circle, point=point, bearing=bearing, exact=exact) 

261 

262 def maxLat(self, bearing): 

263 '''Return the maximum latitude reached when travelling 

264 on a great circle on given bearing from this point 

265 based on Clairaut's formula. 

266 

267 The maximum latitude is independent of longitude 

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

269 

270 Negate the result for the minimum latitude (on the 

271 Southern hemisphere). 

272 

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

274 

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

276 

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

278 ''' 

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

280 return degrees90(m) 

281 

282 def minLat(self, bearing): 

283 '''Return the minimum latitude reached when travelling 

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

285 

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

287 

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

289 

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

291 

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

293 ''' 

294 return -self.maxLat(bearing) 

295 

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

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

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

299 

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

301 see function L{pygeodesy.parse3llh}. 

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

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

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

305 overriding this name. 

306 

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

308 

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

310 ''' 

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

312 r = self.classof(*t) 

313 if name: 

314 r.rename(name) 

315 return r 

316 

317 @property_RO 

318 def _radius(self): 

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

320 ''' 

321 return self.datum.ellipsoid.equatoradius 

322 

323 def _rhumb3(self, other, r=False): 

324 '''(INTERNAL) Rhumb_ helper function. 

325 

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

327 ''' 

328 self.others(other) 

329 

330 a1, b1 = self.philam 

331 a2, b2 = other.philam 

332 # if |db| > 180 take shorter rhumb 

333 # line across the anti-meridian 

334 db = wrapPI(b2 - b1) 

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

336 da = a2 - a1 

337 if r: 

338 # on Mercator projection, longitude distances shrink 

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

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

341 # tolerance to avoid it 

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

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

344 return da, db, dp 

345 

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

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

348 between this and an other (spherical) point. 

349 

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

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

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

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

354 C{False} for backward compatibility. 

355 

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

357 

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

359 B{C{radius}} is invalid. 

360 

361 @example: 

362 

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

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

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

366 ''' 

367 if exact: # use series, always 

368 z = LatLonBase.rhumbAzimuthTo(self, other, exact=False, radius=radius) 

369 else: 

370 _, db, dp = self._rhumb3(other) 

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

372 return z 

373 

374 @deprecated_method 

375 def rhumbBearingTo(self, other): 

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

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

378 

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

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

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

382 

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

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

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

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

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

388 C{B{exact}=True}. 

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

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

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

392 C{False} for backward compatibility. 

393 

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

395 

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

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

398 

399 @example: 

400 

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

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

403 ''' 

404 if exact: # use series, always 

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

406 radius=radius, height=height) 

407 else: # radius=None from .rhumbMidpointTo 

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

409 d, r = self.datum, radius 

410 else: 

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

412 r = d.ellipsoid.equatoradius 

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

414 

415 a1, b1 = self.philam 

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

417 

418 da = r * cb 

419 a2 = a1 + da 

420 # normalize latitude if past pole 

421 if a2 > PI_2: 

422 a2 = PI - a2 

423 elif a2 < -PI_2: 

424 a2 = -PI - a2 

425 

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

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

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

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

430 

431 h = self.height if height is None else Height(height) 

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

433 return r 

434 

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

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

437 a rhumb line (loxodrome). 

438 

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

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

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

442 C{B{exact}=True}. 

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

444 C{False} for backward compatibility. 

445 

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

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

448 

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

450 

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

452 

453 @example: 

454 

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

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

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

458 ''' 

459 if exact: # use series, always 

460 r = LatLonBase.rhumbDistanceTo(self, other, exact=False, radius=radius) 

461 if radius is None: # angular distance in radians 

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

463 else: 

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

465 r, _, _ = self._rhumb3(other, r=True) 

466 if radius is not None: 

467 r *= Radius(radius) 

468 return r 

469 

470 def rhumbMidpointTo(self, other, height=None, radius=R_M, 

471 exact=False, fraction=_0_5): 

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

473 this and an other point. 

474 

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

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

477 (C{meter}). 

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

479 overriding the default C{R_M}. 

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

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

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

483 C{False} for backward compatibility. 

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

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

486 

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

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

489 

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

491 

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

493 

494 @example: 

495 

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

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

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

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

500 ''' 

501 if exact: # use series, always 

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

503 radius=radius, height=height, fraction=fraction) 

504 elif fraction is not _0_5: 

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

506 r, db, dp = self._rhumb3(other, r=True) # radians 

507 z = atan2b(db, dp) 

508 h = self._havg(other, f=f) if height is None else height 

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

510 

511 else: # for backward compatibility 

512 self.others(other) 

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

514 a1, b1 = self.philam 

515 a2, b2 = other.philam 

516 if fabs(b2 - b1) > PI: 

517 b1 += PI2 # crossing anti-meridian 

518 

519 a3 = favg(a1, a2) 

520 b3 = favg(b1, b2) 

521 

522 f1 = tanPI_2_2(a1) 

523 if isnon0(f1): 

524 f2 = tanPI_2_2(a2) 

525 f = f2 / f1 

526 if isnon0(f): 

527 f = log(f) 

528 if isnon0(f): 

529 f3 = tanPI_2_2(a3) 

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

531 -b2, b1, b2 - b1) / f 

532 

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

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

535 h = self._havg(other) if height is None else Height(height) 

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

537 return r 

538 

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

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

541 height}. 

542 

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

544 keyword arguments, ignored if 

545 C{B{Nvector} is None}. 

546 

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

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

549 

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

551 ''' 

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

553 

554 def toWm(self, radius=R_MA): 

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

556 

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

558 

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

560 

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

562 ''' 

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

564 

565 

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

567 height=None, wrap=False): 

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

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

570 

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

572 r = Radius_(circle=r) 

573 elif exact is None: 

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

575 elif isbool(exact): 

576 r = c.rhumbDistanceTo(r, radius=radius, exact=exact) 

577 else: 

578 raise _ValueError(exact=exact) 

579 

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

581 b = Bearing(b) 

582 elif exact is None: 

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

584 else: 

585 b = p.rhumbAzimuthTo(b, radius=radius, exact=exact) 

586 

587 if exact is None: 

588 d = p.distanceTo(c, radius=radius, wrap=wrap) 

589 a = p.initialBearingTo(c, wrap=wrap) if d > EPS else b 

590 else: 

591 d = p.rhumbDistanceTo(c, radius=radius, exact=exact) 

592 a = wrap360(p.rhumbAzimuthTo( c, radius=radius, exact=exact)) if d > EPS else b 

593 

594 if d > EPS: 

595 s, c = sincos2d(b - a) 

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

597 if s > r: 

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

599 elif (r - s) < EPS: 

600 return p, p # tangent 

601 c *= d 

602 else: # coincindent 

603 s, c = r, 0 

604 a = b + _180_0 

605 if exact is None: 

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

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

608 else: 

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

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

611 return b, a # in bearing direction first 

612 

613 

614__all__ += _ALL_DOCS(CartesianSphericalBase, LatLonSphericalBase) 

615 

616# **) MIT License 

617# 

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

619# 

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

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

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

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

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

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

626# 

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

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

629# 

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

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

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

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

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

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

636# OTHER DEALINGS IN THE SOFTWARE.