1 2# -*- coding: utf-8 -*- 3 4u'''(INTERNAL) Spherical geodesy bases. 5 6Base classes C{CartesianSphericalBase} and C{LatLonSphericalBase} 7imported by L{sphericalNvector} or L{sphericalTrigonometry}. 8 9Pure Python implementation of geodetic (lat-/longitude) functions, 10transcoded in part from JavaScript originals by I{(C) Chris Veness 2011-2016} 11and published under the same MIT Licence**, see 12U{Latitude/Longitude<https://www.Movable-Type.co.UK/scripts/latlong.html>}. 13''' 14# make sure int/int division yields float quotient, see .basics 15from __future__ import division 16 17from pygeodesy.basics import isnear0, isnon0, map1 18from pygeodesy.cartesianBase import CartesianBase 19from pygeodesy.datums import Datums, _spherical_datum 20from pygeodesy.ellipsoids import R_M, R_MA 21from pygeodesy.errors import IntersectionError 22from pygeodesy.fmath import favg, fdot, hypot 23from pygeodesy.interns import EPS, NN, PI, PI2, PI_2, _COMMA_, \ 24 _concentric_, _datum_, _distant_, \ 25 _exceed_PI_radians_, _name_, _near_, \ 26 _too_, _1_0, _180_0, _360_0 27from pygeodesy.latlonBase import LatLonBase, _trilaterate5 # PYCHOK passed 28from pygeodesy.lazily import _ALL_DOCS 29from pygeodesy.namedTuples import Bearing2Tuple 30from pygeodesy.nvectorBase import NvectorBase, _xattrs # streprs 31from pygeodesy.props import property_doc_ 32from pygeodesy.units import Bearing_, Height, Radians_, Radius, Radius_ 33from pygeodesy.utily import acos1, atan2b, degrees90, degrees180, \ 34 sincos2, tanPI_2_2, wrapPI 35 36from math import cos, log, sin, sqrt 37 38__all__ = () 39__version__ = '21.08.28' 40 41 42def _angular(distance, radius): # PYCHOK for export 43 '''(INTERNAL) Return the angular distance in C{radians}. 44 45 @raise UnitError: Invalid B{C{distance}} or B{C{radius}}. 46 ''' 47 return Radians_(distance / Radius_(radius=radius), low=EPS) 48 49 50def _rads3(rad1, rad2, radius): # in .sphericalTrigonometry 51 '''(INTERNAL) Convert radii to radians. 52 ''' 53 r1 = Radius_(rad1=rad1) 54 r2 = Radius_(rad2=rad2) 55 if radius is not None: # convert radii to radians 56 r1 = _angular(r1, radius) 57 r2 = _angular(r2, radius) 58 59 x = r1 < r2 60 if x: 61 r1, r2 = r2, r1 62 if r1 > PI: 63 raise IntersectionError(rad1=rad1, rad2=rad2, 64 txt=_exceed_PI_radians_) 65 return r1, r2, x 66 67 68class CartesianSphericalBase(CartesianBase): 69 '''(INTERNAL) Base class for spherical C{Cartesian}s. 70 ''' 71 _datum = Datums.Sphere # L{Datum} 72 73 def intersections2(self, rad1, other, rad2, radius=R_M): 74 '''Compute the intersection points of two circles each defined 75 by a center point and a radius. 76 77 @arg rad1: Radius of the this circle (C{meter} or C{radians}, 78 see B{C{radius}}). 79 @arg other: Center of the other circle (C{Cartesian}). 80 @arg rad2: Radius of the other circle (C{meter} or C{radians}, 81 see B{C{radius}}). 82 @kwarg radius: Mean earth radius (C{meter} or C{None} if both 83 B{C{rad1}} and B{C{rad2}} are given in C{radians}). 84 85 @return: 2-Tuple of the intersection points, each C{Cartesian}. 86 The intersection points are the same C{Cartesian} 87 instance for abutting circles, aka I{radical center}. 88 89 @raise IntersectionError: Concentric, antipodal, invalid or 90 non-intersecting circles. 91 92 @raise TypeError: If B{C{other}} is not C{Cartesian}. 93 94 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}} or B{C{radius}}. 95 96 @see: U{Calculating intersection of two Circles 97 <https://GIS.StackExchange.com/questions/48937/ 98 calculating-intersection-of-two-circles>} and method 99 or function C{trilaterate3d2}. 100 ''' 101 x1, x2 = self, self.others(other) 102 r1, r2, x = _rads3(rad1, rad2, radius) 103 if x: 104 x1, x2 = x2, x1 105 try: 106 n, q = x1.cross(x2), x1.dot(x2) 107 n2, q1 = n.length2, (_1_0 - q**2) 108 if n2 < EPS or isnear0(q1): 109 raise ValueError(_near_(_concentric_)) 110 c1, c2 = cos(r1), cos(r2) 111 x0 = x1.times((c1 - q * c2) / q1).plus( 112 x2.times((c2 - q * c1) / q1)) 113 n1 = _1_0 - x0.length2 114 if n1 < EPS: 115 raise ValueError(_too_(_distant_)) 116 except ValueError as x: 117 raise IntersectionError(center=self, rad1=rad1, 118 other=other, rad2=rad2, txt=str(x)) 119 n = n.times(sqrt(n1 / n2)) 120 if n.length > EPS: 121 x1 = x0.plus(n) 122 x2 = x0.minus(n) 123 else: # abutting circles 124 x1 = x2 = x0 125 126 return (_xattrs(x1, self, _datum_, _name_), 127 _xattrs(x2, self, _datum_, _name_)) 128 129 130class LatLonSphericalBase(LatLonBase): 131 '''(INTERNAL) Base class for spherical C{LatLon}s. 132 ''' 133 _datum = Datums.Sphere # spherical L{Datum} 134 135 def __init__(self, lat, lon, height=0, datum=None, name=NN): 136 '''Create a spherical C{LatLon} point frome the given 137 lat-, longitude and height on the given datum. 138 139 @arg lat: Latitude (C{degrees} or DMS C{[N|S]}). 140 @arg lon: Longitude (C{degrees} or DMS C{str[E|W]}). 141 @kwarg height: Optional elevation (C{meter}, the same units 142 as the datum's half-axes). 143 @kwarg datum: Optional, spherical datum to use (L{Datum}, 144 L{Ellipsoid}, L{Ellipsoid2}, L{a_f2Tuple}) 145 or C{scalar} earth radius). 146 @kwarg name: Optional name (string). 147 148 @raise TypeError: If B{C{datum}} invalid or not 149 not spherical. 150 151 @example: 152 153 >>> p = LatLon(51.4778, -0.0016) # height=0, datum=Datums.WGS84 154 ''' 155 LatLonBase.__init__(self, lat, lon, height=height, name=name) 156 if datum not in (None, self.datum): 157 self._datum = _spherical_datum(datum, name=self.name, raiser=True) 158 159 def bearingTo2(self, other, wrap=False, raiser=False): 160 '''Return the initial and final bearing (forward and reverse 161 azimuth) from this to an other point. 162 163 @arg other: The other point (C{LatLon}). 164 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 165 @kwarg raiser: Optionally, raise L{CrossError} (C{bool}). 166 167 @return: A L{Bearing2Tuple}C{(initial, final)}. 168 169 @raise TypeError: The B{C{other}} point is not spherical. 170 171 @see: Methods C{initialBearingTo} and C{finalBearingTo}. 172 ''' 173 # .initialBearingTo is inside .-Nvector and .-Trigonometry 174 i = self.initialBearingTo(other, wrap=wrap, raiser=raiser) # PYCHOK .initialBearingTo 175 f = self.finalBearingTo( other, wrap=wrap, raiser=raiser) 176 return Bearing2Tuple(i, f, name=self.name) 177 178 @property_doc_(''' this point's datum (L{Datum}).''') 179 def datum(self): 180 '''Get this point's datum (L{Datum}). 181 ''' 182 return self._datum 183 184 @datum.setter # PYCHOK setter! 185 def datum(self, datum): 186 '''Set this point's datum I{without conversion}. 187 188 @arg datum: New spherical datum (L{Datum}, L{Ellipsoid}, 189 L{Ellipsoid2}, L{a_f2Tuple}) or C{scalar} 190 earth radius). 191 192 @raise TypeError: If B{C{datum}} invalid or not 193 not spherical. 194 ''' 195 d = _spherical_datum(datum, name=self.name, raiser=True) 196 self._update(d != self._datum) 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 (b + _180_0) % _360_0 # == wrap360 since b >= 0 223 224 def maxLat(self, bearing): 225 '''Return the maximum latitude reached when travelling 226 on a great circle on given bearing from this point 227 based on Clairaut's formula. 228 229 The maximum latitude is independent of longitude 230 and the same for all points on a given latitude. 231 232 Negate the result for the minimum latitude (on the 233 Southern hemisphere). 234 235 @arg bearing: Initial bearing (compass C{degrees360}). 236 237 @return: Maximum latitude (C{degrees90}). 238 239 @raise ValueError: Invalid B{C{bearing}}. 240 241 @JSname: I{maxLatitude}. 242 ''' 243 m = acos1(abs(sin(Bearing_(bearing)) * cos(self.phi))) 244 return degrees90(m) 245 246 def minLat(self, bearing): 247 '''Return the minimum latitude reached when travelling 248 on a great circle on given bearing from this point. 249 250 @arg bearing: Initial bearing (compass C{degrees360}). 251 252 @return: Minimum latitude (C{degrees90}). 253 254 @see: Method L{maxLat} for more details. 255 256 @raise ValueError: Invalid B{C{bearing}}. 257 258 @JSname: I{minLatitude}. 259 ''' 260 return -self.maxLat(bearing) 261 262 def parse(self, strllh, height=0, sep=_COMMA_, name=NN): 263 '''Parse a string representing a similar, spherical C{LatLon} 264 point, consisting of C{"lat, lon[, height]"}. 265 266 @arg strllh: Lat, lon and optional height (C{str}), 267 see function L{parse3llh}. 268 @kwarg height: Optional, default height (C{meter}). 269 @kwarg sep: Optional separator (C{str}). 270 @kwarg name: Optional instance name (C{str}), 271 overriding this name. 272 273 @return: The similar point (spherical C{LatLon}). 274 275 @raise ParseError: Invalid B{C{strllh}}. 276 ''' 277 from pygeodesy.dms import parse3llh 278 r = self.classof(*parse3llh(strllh, height=height, sep=sep)) 279 if name: 280 r.rename(name) 281 return r 282 283 def _rhumb3(self, other): 284 '''(INTERNAL) Rhumb_ helper function. 285 286 @arg other: The other point (spherical C{LatLon}). 287 ''' 288 self.others(other) 289 290 a1, b1 = self.philam 291 a2, b2 = other.philam 292 # if |db| > 180 take shorter rhumb 293 # line across the anti-meridian 294 db = wrapPI(b2 - b1) 295 dp = log(tanPI_2_2(a2) / tanPI_2_2(a1)) 296 return (a2 - a1), db, dp 297 298 def rhumbBearingTo(self, other): 299 '''Return the initial bearing (forward azimuth) from this to 300 an other point along a rhumb (loxodrome) line. 301 302 @arg other: The other point (spherical C{LatLon}). 303 304 @return: Initial bearing (compass C{degrees360}). 305 306 @raise TypeError: The B{C{other}} point is not spherical. 307 308 @example: 309 310 >>> p = LatLon(51.127, 1.338) 311 >>> q = LatLon(50.964, 1.853) 312 >>> b = p.rhumbBearingTo(q) # 116.7 313 ''' 314 _, db, dp = self._rhumb3(other) 315 return atan2b(db, dp) 316 317 def rhumbDestination(self, distance, bearing, radius=R_M, height=None): 318 '''Return the destination point having travelled along a rhumb 319 (loxodrome) line from this point the given distance on the 320 given bearing. 321 322 @arg distance: Distance travelled (C{meter}, same units as 323 B{C{radius}}). 324 @arg bearing: Bearing from this point (compass C{degrees360}). 325 @kwarg radius: Mean earth radius (C{meter}). 326 @kwarg height: Optional height, overriding the default height 327 (C{meter}, same unit as B{C{radius}}). 328 329 @return: The destination point (spherical C{LatLon}). 330 331 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}}, 332 B{C{radius}} or B{C{height}}. 333 334 @example: 335 336 >>> p = LatLon(51.127, 1.338) 337 >>> q = p.rhumbDestination(40300, 116.7) # 50.9642°N, 001.8530°E 338 339 @JSname: I{rhumbDestinationPoint} 340 ''' 341 r = _angular(distance, radius) 342 343 a1, b1 = self.philam 344 sb, cb = sincos2(Bearing_(bearing)) 345 346 da = r * cb 347 a2 = a1 + da 348 # normalize latitude if past pole 349 if a2 > PI_2: 350 a2 = PI - a2 351 elif a2 < -PI_2: 352 a2 = -PI - a2 353 354 dp = log(tanPI_2_2(a2) / tanPI_2_2(a1)) 355 # E-W course becomes ill-conditioned with 0/0 356 q = (da / dp) if abs(dp) > EPS else cos(a1) 357 b2 = (b1 + r * sb / q) if abs(q) > EPS else b1 358 359 h = self.height if height is None else Height(height) 360 return self.classof(degrees90(a2), degrees180(b2), height=h) 361 362 def rhumbDistanceTo(self, other, radius=R_M): 363 '''Return the distance from this to an other point along a rhumb 364 (loxodrome) line. 365 366 @arg other: The other point (spherical C{LatLon}). 367 @kwarg radius: Mean earth radius (C{meter}) or C{None}. 368 369 @return: Distance (C{meter}, the same units as B{C{radius}} 370 or C{radians} if B{C{radius}} is C{None}). 371 372 @raise TypeError: The B{C{other}} point is not spherical. 373 374 @raise ValueError: Invalid B{C{radius}}. 375 376 @example: 377 378 >>> p = LatLon(51.127, 1.338) 379 >>> q = LatLon(50.964, 1.853) 380 >>> d = p.rhumbDistanceTo(q) # 403100 381 ''' 382 # see <https://www.EdWilliams.org/avform.htm#Rhumb> 383 da, db, dp = self._rhumb3(other) 384 385 # on Mercator projection, longitude distances shrink 386 # by latitude; the 'stretch factor' q becomes ill- 387 # conditioned along E-W line (0/0); use an empirical 388 # tolerance to avoid it 389 q = (da / dp) if abs(dp) > EPS else cos(self.phi) 390 r = hypot(da, q * db) 391 return r if radius is None else (Radius(radius) * r) 392 393 def rhumbMidpointTo(self, other, height=None): 394 '''Return the (loxodromic) midpoint between this and 395 an other point. 396 397 @arg other: The other point (spherical LatLon). 398 @kwarg height: Optional height, overriding the mean height 399 (C{meter}). 400 401 @return: The midpoint (spherical C{LatLon}). 402 403 @raise TypeError: The B{C{other}} point is not spherical. 404 405 @raise ValueError: Invalid B{C{height}}. 406 407 @example: 408 409 >>> p = LatLon(51.127, 1.338) 410 >>> q = LatLon(50.964, 1.853) 411 >>> m = p.rhumb_midpointTo(q) 412 >>> m.toStr() # '51.0455°N, 001.5957°E' 413 ''' 414 self.others(other) 415 416 # see <https://MathForum.org/library/drmath/view/51822.html> 417 a1, b1 = self.philam 418 a2, b2 = other.philam 419 if abs(b2 - b1) > PI: 420 b1 += PI2 # crossing anti-meridian 421 422 a3 = favg(a1, a2) 423 b3 = favg(b1, b2) 424 425 f1 = tanPI_2_2(a1) 426 if isnon0(f1): 427 f2 = tanPI_2_2(a2) 428 f = f2 / f1 429 if isnon0(f): 430 f = log(f) 431 if isnon0(f): 432 f3 = tanPI_2_2(a3) 433 b3 = fdot(map1(log, f1, f2, f3), 434 -b2, b1, b2 - b1) / f 435 436 h = self._havg(other) if height is None else Height(height) 437 return self.classof(degrees90(a3), degrees180(b3), height=h) 438 439 def toNvector(self, Nvector=NvectorBase, **Nvector_kwds): # PYCHOK signature 440 '''Convert this point to C{Nvector} components, I{including 441 height}. 442 443 @kwarg Nvector_kwds: Optional, additional B{C{Nvector}} 444 keyword arguments, ignored if 445 C{B{Nvector} is None}. 446 447 @return: An B{C{Nvector}} or a L{Vector4Tuple}C{(x, y, z, h)} 448 if B{C{Nvector}} is C{None}. 449 450 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}. 451 ''' 452 return LatLonBase.toNvector(self, Nvector=Nvector, **Nvector_kwds) 453 454 def toWm(self, radius=R_MA): 455 '''Convert this point to a I{WM} coordinate. 456 457 @kwarg radius: Optional earth radius (C{meter}). 458 459 @return: The WM coordinate (L{Wm}). 460 461 @see: Function L{toWm} in module L{webmercator} for details. 462 ''' 463 from pygeodesy.webmercator import toWm 464 return toWm(self, radius=radius) 465 466 467__all__ += _ALL_DOCS(CartesianSphericalBase, LatLonSphericalBase) 468 469# **) MIT License 470# 471# Copyright (C) 2016-2021 -- mrJean1 at Gmail -- All Rights Reserved. 472# 473# Permission is hereby granted, free of charge, to any person obtaining a 474# copy of this software and associated documentation files (the "Software"), 475# to deal in the Software without restriction, including without limitation 476# the rights to use, copy, modify, merge, publish, distribute, sublicense, 477# and/or sell copies of the Software, and to permit persons to whom the 478# Software is furnished to do so, subject to the following conditions: 479# 480# The above copyright notice and this permission notice shall be included 481# in all copies or substantial portions of the Software. 482# 483# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 484# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 485# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 486# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 487# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 488# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 489# OTHER DEALINGS IN THE SOFTWARE. 490