1 2# -*- coding: utf-8 -*- 3 4u'''(INTERNAL) Ellipsoidal geodesy base classes C{CartesianEllipsoidalBase} 5and C{LatLonEllipsoidalBase}. 6 7Pure Python implementation of geodesy tools for ellipsoidal earth models, 8transcoded in part from JavaScript originals by I{(C) Chris Veness 2005-2016} 9and published under the same MIT Licence**, see for example U{latlon-ellipsoidal 10<https://www.Movable-Type.co.UK/scripts/geodesy/docs/latlon-ellipsoidal.js.html>}. 11''' 12# make sure int/int division yields float quotient, see .basics 13from __future__ import division 14 15from pygeodesy.basics import _xinstanceof 16from pygeodesy.cartesianBase import CartesianBase 17from pygeodesy.datums import Datum, Datums, _ellipsoidal_datum, _WGS84 18from pygeodesy.errors import _incompatible, _IsnotError, RangeError, \ 19 TRFError, _ValueError, _xError, _xkwds_not 20# from pygeodesy.errors import _xkwds # from .named 21from pygeodesy.interns import _ellipsoidal_ # PYCHOK used! 22from pygeodesy.interns import EPS, EPS0, EPS1, MISSING, NN, _COMMA_, \ 23 _conversion_, _datum_, _DOT_, _N_, _no_, \ 24 _reframe_, _SPACE_, _0_0 25from pygeodesy.latlonBase import LatLonBase, _trilaterate5 26from pygeodesy.lazily import _ALL_DOCS 27from pygeodesy.named import _xnamed, _xkwds 28from pygeodesy.namedTuples import Vector3Tuple 29from pygeodesy.props import deprecated_method, Property_RO, \ 30 property_doc_, property_RO 31from pygeodesy.units import Epoch, _1mm as _TOL_M, Radius_ 32 33__all__ = () 34__version__ = '21.09.14' 35 36 37class CartesianEllipsoidalBase(CartesianBase): 38 '''(INTERNAL) Base class for ellipsoidal C{Cartesian}s. 39 ''' 40 _datum = _WGS84 # L{Datum} 41 42 def _applyHelmerts(self, *transforms): 43 '''(INTERNAL) Apply one I{or more} Helmert transforms. 44 ''' 45 xyz = self.xyz 46 for t in transforms: 47 xyz = t.transform(*xyz) 48 return self.classof(xyz, datum=self.datum) 49 50 @deprecated_method 51 def convertRefFrame(self, reframe2, reframe, epoch=None): 52 '''DEPRECATED, use method L{toRefFrame}.''' 53 return self.toRefFrame(reframe2, reframe, epoch=epoch) 54 55 def intersections2(self, radius, center2, radius2, sphere=True, 56 Vector=None, **Vector_kwds): 57 '''Compute the intersection of two spheres or circles, each defined by a 58 cartesian center point and a radius. 59 60 @arg radius: Radius of this sphere or circle (same units as this point's 61 coordinates). 62 @arg center2: Center of the second sphere or circle (C{Cartesian}, L{Vector3d}, 63 C{Vector3Tuple} or C{Vector4Tuple}). 64 @arg radius2: Radius of the second sphere or circle (same units as this and 65 the B{C{other}} point's coordinates). 66 @kwarg sphere: If C{True} compute the center and radius of the intersection 67 of two I{spheres}. If C{False}, ignore the C{z}-component and 68 compute the intersection of two I{circles} (C{bool}). 69 @kwarg Vector: Class to return intersections (C{Cartesian}, L{Vector3d} or 70 C{Vector3Tuple}) or C{None} for an instance of this (sub-)class. 71 @kwarg Vector_kwds: Optional, additional B{C{Vector}} keyword arguments, 72 ignored if C{B{Vector} is None}. 73 74 @return: If B{C{sphere}} is C{True}, a 2-tuple of the C{center} and C{radius} 75 of the intersection of the I{spheres}. The C{radius} is C{0.0} for 76 abutting spheres (and the C{center} is aka I{radical center}). 77 78 If B{C{sphere}} is C{False}, a 2-tuple with the two intersection 79 points of the I{circles}. For abutting circles, both points are 80 the same instance, aka I{radical center}. 81 82 @raise IntersectionError: Concentric, invalid or non-intersecting spheres or circles. 83 84 @raise TypeError: Invalid B{C{center2}}. 85 86 @raise UnitError: Invalid B{C{radius}} or B{C{radius2}}. 87 88 @see: U{Sphere-Sphere<https://MathWorld.Wolfram.com/Sphere-SphereIntersection.html>}, 89 U{Circle-Circle<https://MathWorld.Wolfram.com/Circle-CircleIntersection.html>} 90 Intersection and function L{radical2}. 91 ''' 92 from pygeodesy.vector3d import _intersects2, _otherV3d 93 try: 94 return _intersects2(_otherV3d(useZ=sphere, center=self), 95 Radius_(radius=radius), 96 center2, Radius_(radius2=radius2), 97 sphere=sphere, clas=self.classof, 98 Vector=Vector, **Vector_kwds) 99 except (TypeError, ValueError) as x: 100 raise _xError(x, center=self, radius=radius, center2=center2, radius2=radius2) 101 102 def toRefFrame(self, reframe2, reframe, epoch=None): 103 '''Convert this cartesian point from one to an other reference frame. 104 105 @arg reframe2: Reference frame to convert I{to} (L{RefFrame}). 106 @arg reframe: Reference frame to convert I{from} (L{RefFrame}). 107 @kwarg epoch: Optional epoch to observe (C{scalar}, fractional 108 calendar year), overriding B{C{reframe}}'s epoch. 109 110 @return: The converted point (C{Cartesian}) or this point if 111 conversion is C{nil}. 112 113 @raise TRFError: No conversion available from B{C{reframe}} 114 to B{C{reframe2}} or invalid B{C{epoch}}. 115 116 @raise TypeError: B{C{reframe2}} or B{C{reframe}} not a 117 L{RefFrame}. 118 ''' 119 from pygeodesy.trf import RefFrame, _reframeTransforms2 120 _xinstanceof(RefFrame, reframe2=reframe2, reframe=reframe) 121 122 _, xs = _reframeTransforms2(reframe2, reframe, epoch) 123 return self._applyHelmerts(*xs) if xs else self 124 125 126class LatLonEllipsoidalBase(LatLonBase): 127 '''(INTERNAL) Base class for ellipsoidal C{LatLon}s. 128 ''' 129 _convergence = None # UTM/UPS meridian convergence (C{degrees}) 130 _datum = _WGS84 # L{Datum} 131 _elevation2to = None # _elevation2 timeout (C{secs}) 132 _epoch = None # overriding .reframe.epoch (C{float}) 133 _geoidHeight2to = None # _geoidHeight2 timeout (C{secs}) 134 _iteration = None # iteration number (C{int} or C{None}) 135 _reframe = None # reference frame (L{RefFrame}) 136 _scale = None # UTM/UPS scale factor (C{float}) 137 138 def __init__(self, lat, lon, height=0, datum=None, reframe=None, 139 epoch=None, name=NN): 140 '''Create an ellipsoidal C{LatLon} point frome the given 141 lat-, longitude and height on the given datum and with 142 the given reference frame and epoch. 143 144 @arg lat: Latitude (C{degrees} or DMS C{[N|S]}). 145 @arg lon: Longitude (C{degrees} or DMS C{str[E|W]}). 146 @kwarg height: Optional elevation (C{meter}, the same units 147 as the datum's half-axes). 148 @kwarg datum: Optional, ellipsoidal datum to use (L{Datum}, 149 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}). 150 @kwarg reframe: Optional reference frame (L{RefFrame}). 151 @kwarg epoch: Optional epoch to observe for B{C{reframe}} 152 (C{scalar}), a non-zero, fractional calendar 153 year; silently ignored if C{B{reframe} is None}. 154 @kwarg name: Optional name (string). 155 156 @raise RangeError: Value of B{C{lat}} or B{C{lon}} outside the valid 157 range and C{rangerrors} set to C{True}. 158 159 @raise TypeError: B{C{datum}} is not a L{datum}, B{C{reframe}} 160 is not a L{RefFrame} or B{C{epoch}} is not 161 C{scalar} non-zero. 162 163 @raise UnitError: Invalid B{C{lat}}, B{C{lon}} or B{C{height}}. 164 165 @example: 166 167 >>> p = LatLon(51.4778, -0.0016) # height=0, datum=Datums.WGS84 168 ''' 169 LatLonBase.__init__(self, lat, lon, height=height, name=name) 170 if datum not in (None, self._datum): 171 self.datum = _ellipsoidal_datum(datum, name=name) 172 if reframe: 173 self.reframe = reframe 174 self.epoch = epoch 175 176 def antipode(self, height=None): 177 '''Return the antipode, the point diametrically opposite 178 to this point. 179 180 @kwarg height: Optional height of the antipode, height 181 of this point otherwise (C{meter}). 182 183 @return: The antipodal point (C{LatLon}). 184 ''' 185 lla = LatLonBase.antipode(self, height=height) 186 if lla.datum != self.datum: 187 lla.datum = self.datum 188 return lla 189 190 @property_RO 191 def convergence(self): 192 '''Get this point's UTM or UPS meridian convergence (C{degrees}) or 193 C{None} if not available or not converted from L{Utm} or L{Ups}. 194 ''' 195 return self._convergence 196 197 @deprecated_method 198 def convertDatum(self, datum2): 199 '''DEPRECATED, use method L{toDatum}.''' 200 return self.toDatum(datum2) 201 202 @deprecated_method 203 def convertRefFrame(self, reframe2): 204 '''DEPRECATED, use method L{toRefFrame}.''' 205 return self.toRefFrame(reframe2) 206 207 @Property_RO 208 def _css(self): 209 '''(INTERNAL) Get this C{LatLon} point as a Cassini-Soldner location (L{Css}). 210 ''' 211 from pygeodesy.css import Css, toCss 212 return toCss(self, height=self.height, Css=Css, name=self.name) 213 214 @property_doc_(''' this points's datum (L{Datum}).''') 215 def datum(self): 216 '''Get this point's datum (L{Datum}). 217 ''' 218 return self._datum 219 220 @datum.setter # PYCHOK setter! 221 def datum(self, datum): 222 '''Set this point's datum I{without conversion}. 223 224 @arg datum: New datum (L{Datum}). 225 226 @raise TypeError: The B{C{datum}} is not a L{Datum} 227 or not ellipsoidal. 228 ''' 229 _xinstanceof(Datum, datum=datum) 230 if not datum.isEllipsoidal: 231 raise _IsnotError(_ellipsoidal_, datum=datum) 232 self._update(datum != self._datum) 233 self._datum = datum 234 235 def distanceTo2(self, other): 236 '''I{Approximate} the distance and (initial) bearing between this 237 and an other (ellipsoidal) point based on the radii of curvature. 238 239 I{Suitable only for short distances up to a few hundred Km 240 or Miles and only between points not near-polar}. 241 242 @arg other: The other point (C{LatLon}). 243 244 @return: An L{Distance2Tuple}C{(distance, initial)}. 245 246 @raise TypeError: The B{C{other}} point is not C{LatLon}. 247 248 @raise ValueError: Incompatible datum ellipsoids. 249 250 @see: Method L{Ellipsoid.distance2} and U{Local, flat earth 251 approximation<https://www.EdWilliams.org/avform.htm#flat>} 252 aka U{Hubeny<https://www.OVG.AT/de/vgi/files/pdf/3781/>} 253 formula. 254 ''' 255 return self.ellipsoids(other).distance2(self.lat, self.lon, 256 other.lat, other.lon) 257 258 @Property_RO 259 def _elevation2(self): 260 '''(INTERNAL) Get elevation and data source. 261 ''' 262 from pygeodesy.elevations import elevation2 263 return elevation2(self.lat, self.lon, timeout=self._elevation2to) 264 265 def elevation2(self, adjust=True, datum=_WGS84, timeout=2): 266 '''Return elevation of this point for its or the given datum. 267 268 @kwarg adjust: Adjust the elevation for a B{C{datum}} other 269 than C{NAD83} (C{bool}). 270 @kwarg datum: Optional datum (L{Datum}). 271 @kwarg timeout: Optional query timeout (C{seconds}). 272 273 @return: An L{Elevation2Tuple}C{(elevation, data_source)} 274 or C{(None, error)} in case of errors. 275 276 @note: The adjustment applied is the difference in geocentric 277 earth radius between the B{C{datum}} and C{NAV83} 278 upon which the L{elevations.elevation2} is based. 279 280 @note: NED elevation is only available for locations within 281 the U{Conterminous US (CONUS) 282 <https://WikiPedia.org/wiki/Contiguous_United_States>}. 283 284 @see: Function L{elevations.elevation2} and method 285 L{Ellipsoid.Rgeocentric} for further details and 286 possible C{error}s. 287 ''' 288 if self._elevation2to != timeout: 289 self._elevation2to = timeout 290 LatLonEllipsoidalBase._elevation2._update(self) 291 return self._Radjust2(adjust, datum, self._elevation2) 292 293 def ellipsoid(self, datum=_WGS84): 294 '''Return the ellipsoid of this point's datum or the given datum. 295 296 @kwarg datum: Default datum (L{Datum}). 297 298 @return: The ellipsoid (L{Ellipsoid} or L{Ellipsoid2}). 299 ''' 300 return getattr(self, _datum_, datum).ellipsoid 301 302 def ellipsoids(self, other): 303 '''Check the type and ellipsoid of this and an other point's datum. 304 305 @arg other: The other point (C{LatLon}). 306 307 @return: This point's datum ellipsoid (L{Ellipsoid} or L{Ellipsoid2}). 308 309 @raise TypeError: The B{C{other}} point is not C{LatLon}. 310 311 @raise ValueError: Incompatible datum ellipsoids. 312 ''' 313 self.others(other, up=2) # ellipsoids' caller 314 315 E = self.ellipsoid() 316 try: # other may be Sphere, etc. 317 e = other.ellipsoid() 318 except AttributeError: 319 try: # no ellipsoid method, try datum 320 e = other.datum.ellipsoid 321 except AttributeError: 322 e = E # no datum, XXX assume equivalent? 323 if e != E: 324 raise _ValueError(e.named2, txt=_incompatible(E.named2)) 325 return E 326 327 @property_doc_(''' this point's observed or C{reframe} epoch (C{float}).''') 328 def epoch(self): 329 '''Get this point's observed or C{reframe} epoch (C{float}) or C{None}. 330 ''' 331 return self._epoch or (self.reframe.epoch if self.reframe else None) 332 333 @epoch.setter # PYCHOK setter! 334 def epoch(self, epoch): 335 '''Set or clear this point's observed epoch. 336 337 @arg epoch: Observed epoch, a fractional calendar year 338 (L{Epoch}, C{scalar}) or C{None}. 339 340 @raise TRFError: Invalid B{C{epoch}}. 341 ''' 342 self._epoch = None if epoch is None else Epoch(epoch) 343 344 @Property_RO 345 def Equidistant(self): 346 '''Get the prefered azimuthal equidistant projection I{class} (L{EquidistantKarney} or L{EquidistantExact}). 347 ''' 348 from pygeodesy.azimuthal import EquidistantKarney, EquidistantExact 349 try: 350 _ = self.datum.ellipsoid.geodesic 351 return EquidistantKarney 352 except ImportError: # no geographiclib 353 return EquidistantExact # XXX no longer L{azimuthal.Equidistant} 354 355 @Property_RO 356 def _etm(self): 357 '''(INTERNAL) Get this C{LatLon} point as an ETM coordinate (L{toEtm8}). 358 ''' 359 from pygeodesy.etm import toEtm8, Etm 360 return toEtm8(self, datum=self.datum, Etm=Etm) 361 362 @Property_RO 363 def _geoidHeight2(self): 364 '''(INTERNAL) Get geoid height and model. 365 ''' 366 from pygeodesy.elevations import geoidHeight2 367 return geoidHeight2(self.lat, self.lon, model=0, timeout=self._geoidHeight2to) 368 369 def geoidHeight2(self, adjust=False, datum=_WGS84, timeout=2): 370 '''Return geoid height of this point for its or the given datum. 371 372 @kwarg adjust: Adjust the geoid height for a B{C{datum}} 373 other than C{NAD83/NADV88} (C{bool}). 374 @kwarg datum: Optional datum (L{Datum}). 375 @kwarg timeout: Optional query timeout (C{seconds}). 376 377 @return: A L{GeoidHeight2Tuple}C{(height, model_name)} or 378 C{(None, error)} in case of errors. 379 380 @note: The adjustment applied is the difference in geocentric 381 earth radius between the B{C{datum}} and C{NAV83/NADV88} 382 upon which the L{elevations.geoidHeight2} is based. 383 384 @note: The geoid height is only available for locations within 385 the U{Conterminous US (CONUS) 386 <https://WikiPedia.org/wiki/Contiguous_United_States>}. 387 388 @see: Function L{elevations.geoidHeight2} and method 389 L{Ellipsoid.Rgeocentric} for further details and 390 possible C{error}s. 391 ''' 392 if self._geoidHeight2to != timeout: 393 self._geoidHeight2to = timeout 394 LatLonEllipsoidalBase._geoidHeight2._update(self) 395 return self._Radjust2(adjust, datum, self._geoidHeight2) 396 397 def intersection3(self, end1, other, end2, height=None, wrap=True, 398 equidistant=None, tol=_TOL_M): 399 '''Interatively compute the intersection point of two paths, each 400 defined by two points or a start point and bearing from North. 401 402 @arg end1: End point of this path (C{LatLon}) or the initial 403 bearing at this point (compass C{degrees360}). 404 @arg other: Start point of the other path (C{LatLon}). 405 @arg end2: End point of the other path (C{LatLon}) or the 406 initial bearing at the other point (compass 407 C{degrees360}). 408 @kwarg height: Optional height at the intersection (C{meter}, 409 conventionally) or C{None} for the mean height. 410 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 411 @kwarg equidistant: An azimuthal equidistant projection (I{class} 412 or function L{equidistant}), or C{None} for 413 this point's preferred C{.Equidistant}. 414 @kwarg tol: Tolerance for skew line distance and length and for 415 convergence (C{meter}, conventionally). 416 417 @return: An L{Intersection3Tuple}C{(point, outside1, outside2)} 418 with C{point} a C{LatLon} instance. 419 420 @raise ImportError: Package U{geographiclib 421 <https://PyPI.org/project/geographiclib>} 422 not installed or not found, but only if 423 C{B{equidistant}=}L{EquidistantKarney}. 424 425 @raise IntersectionError: Skew, colinear, parallel or otherwise 426 non-intersecting paths or no convergence 427 for the B{C{tol}}. 428 429 @raise TypeError: If B{C{end1}}, B{C{other}} or B{C{end2}} point 430 is not C{LatLon}. 431 432 @note: For each path specified with an initial bearing, a pseudo-end 433 point is computed as the C{destination} along that bearing at 434 about 1.5 times the distance from the start point to an initial 435 gu-/estimate of the intersection point (and between 1/8 and 3/8 436 of the authalic earth perimeter). 437 ''' 438 from pygeodesy.ellipsoidalBaseDI import _intersect3 439 try: 440 s2 = self.others(other) 441 return _intersect3(self, end1, s2, end2, height=height, wrap=wrap, 442 equidistant=equidistant, tol=tol, 443 LatLon=self.classof, datum=self.datum) 444 except (TypeError, ValueError) as x: 445 raise _xError(x, start1=self, end1=end1, other=other, end2=end2) 446 447 def intersections2(self, radius1, other, radius2, height=None, wrap=True, 448 equidistant=None, tol=_TOL_M): 449 '''Interatively compute the intersection points of two circles, 450 each defined by a center point and a radius. 451 452 @arg radius1: Radius of this circle (C{meter}, conventionally). 453 @arg other: Center of the other circle (C{LatLon}). 454 @arg radius2: Radius of the other circle (C{meter}, same units as 455 B{C{radius1}}). 456 @kwarg height: Optional height for the intersection points (C{meter}, 457 conventionally) or C{None} for the I{"radical height"} 458 at the I{radical line} between both centers. 459 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 460 @kwarg equidistant: An azimuthal equidistant projection (I{class} 461 or function L{equidistant}), or C{None} for 462 this point's preferred C{.Equidistant}. 463 @kwarg tol: Convergence tolerance (C{meter}, same units as 464 B{C{radius1}} and B{C{radius2}}). 465 466 @return: 2-Tuple of the intersection points, each a C{LatLon} 467 instance. For abutting circles, both intersection 468 points are the same instance, aka I{radical center}. 469 470 @raise ImportError: Package U{geographiclib 471 <https://PyPI.org/project/geographiclib>} 472 not installed or not found, but only if 473 C{B{equidistant}=}L{EquidistantKarney}. 474 475 @raise IntersectionError: Concentric, antipodal, invalid or 476 non-intersecting circles or no 477 convergence for B{C{tol}}. 478 479 @raise TypeError: Invalid B{C{other}} or B{C{equidistant}}. 480 481 @raise UnitError: Invalid B{C{radius1}}, B{C{radius2}} or B{C{height}}. 482 483 @see: U{The B{ellipsoidal} case<https://GIS.StackExchange.com/questions/48937/ 484 calculating-intersection-of-two-circles>}, U{Karney's paper 485 <https://ArXiv.org/pdf/1102.1215.pdf>}, pp 20-21, section B{14. MARITIME BOUNDARIES}, 486 U{circle-circle<https://MathWorld.Wolfram.com/Circle-CircleIntersection.html>} and 487 U{sphere-sphere<https://MathWorld.Wolfram.com/Sphere-SphereIntersection.html>} 488 intersections. 489 ''' 490 from pygeodesy.ellipsoidalBaseDI import _intersections2 491 try: 492 c2 = self.others(other) 493 return _intersections2(self, radius1, c2, radius2, height=height, wrap=wrap, 494 equidistant=equidistant, tol=tol, 495 LatLon=self.classof, datum=self.datum) 496 except (AssertionError, TypeError, ValueError) as x: 497 raise _xError(x, center=self, radius1=radius1, other=other, radius2=radius2) 498 499 @property_RO 500 def iteration(self): 501 '''Get the most recent C{intersections2} or C{nearestOn} iteration 502 number (C{int}) or C{None} if not available/applicable. 503 ''' 504 return self._iteration 505 506 @Property_RO 507 def _lcc(self): 508 '''(INTERNAL) Get this C{LatLon} point as a Lambert location (L{Lcc}). 509 ''' 510 from pygeodesy.lcc import Lcc, toLcc 511 return toLcc(self, height=self.height, Lcc=Lcc, name=self.name) 512 513 def nearestOn(self, point1, point2, within=True, height=None, wrap=True, 514 equidistant=None, tol=_TOL_M): 515 '''Interatively locate the closest point between two other points. 516 517 @arg point1: Start point (C{LatLon}). 518 @arg point2: End point (C{LatLon}). 519 @kwarg within: If C{True} return the closest point I{between} 520 B{C{point1}} and B{C{point2}}, otherwise the 521 closest point elsewhere on the arc (C{bool}). 522 @kwarg height: Optional height for the closest point (C{meter}, 523 conventionally) or C{None} or C{False} to 524 interpolate the height. 525 @kwarg height: Optional height for the closest point (C{meter}, 526 conventionally) or C{None} or C{False} for the 527 interpolated height. If C{False}, the distance 528 between points takes height into account. 529 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 530 @kwarg equidistant: An azimuthal equidistant projection (I{class} 531 or function L{equidistant}), or C{None} for 532 this point's preferred C{.Equidistant}. 533 @kwarg tol: Convergence tolerance (C{meter}, conventionally). 534 535 @return: Closest point (C{LatLon}). 536 537 @raise ImportError: Package U{geographiclib 538 <https://PyPI.org/project/geographiclib>} 539 not installed or not found, but only if 540 C{B{equidistant}=}L{EquidistantKarney}. 541 542 @raise TypeError: Invalid B{C{point1}}, B{C{point2}} or 543 B{C{equidistant}}. 544 545 @raise ValueError: No convergence for the B{C{tol}}. 546 ''' 547 from pygeodesy.ellipsoidalBaseDI import _nearestOne 548 try: 549 return _nearestOne(self, point1, point2, within=within, 550 height=height, wrap=wrap, 551 equidistant=equidistant, tol=tol, 552 LatLon=self.classof, datum=self.datum) 553 except (TypeError, ValueError) as x: 554 raise _xError(x, point=self, point1=point1, point2=point2) 555 556 @Property_RO 557 def _osgr(self): 558 '''(INTERNAL) Get this C{LatLon} point to an OSGR coordinate (L{Osgr}). 559 ''' 560 from pygeodesy.osgr import Osgr, toOsgr 561 return toOsgr(self, datum=self.datum, Osgr=Osgr, name=self.name) 562 563 def parse(self, strllh, height=0, datum=None, epoch=None, reframe=None, 564 sep=_COMMA_, name=NN): 565 '''Parse a string representing a similar, ellipsoidal C{LatLon} 566 point, consisting of C{"lat, lon[, height]"}. 567 568 @arg strllh: Lat, lon and optional height (C{str}), 569 see function L{parse3llh}. 570 @kwarg height: Optional, default height (C{meter} or 571 C{None}). 572 @kwarg datum: Optional datum (L{Datum}), overriding this 573 datum I{without conversion}. 574 @kwarg epoch: Optional datum (L{Epoch}), overriding this 575 epoch I{without conversion}. 576 @kwarg reframe: Optional datum (L{RefFrame}), overriding 577 this reframe I{without conversion}. 578 @kwarg sep: Optional separator (C{str}). 579 @kwarg name: Optional instance name (C{str}), overriding 580 this name. 581 582 @return: The similar point (ellipsoidal C{LatLon}). 583 584 @raise ParseError: Invalid B{C{strllh}}. 585 ''' 586 from pygeodesy.dms import parse3llh 587 a, b, h = parse3llh(strllh, height=height, sep=sep) 588 r = self.classof(a, b, height=h, datum=self.datum) 589 if datum not in (None, self.datum): 590 r.datum = datum 591 if epoch not in (None, self.epoch): 592 r.epoch = epoch 593 if reframe not in (None, self.reframe): 594 r.reframe = reframe 595 if name: 596 r = _xnamed(r, name, force=True) 597 return r 598 599 def _Radjust2(self, adjust, datum, meter_text2): 600 '''(INTERNAL) Adjust elevation or geoidHeight with difference 601 in Gaussian radii of curvature of given datum and NAD83. 602 603 @note: This is an arbitrary, possibly incorrect adjustment. 604 ''' 605 if adjust: # Elevation2Tuple or GeoidHeight2Tuple 606 m, t = meter_text2 607 if isinstance(m, float) and abs(m) > EPS: 608 n = Datums.NAD83.ellipsoid.rocGauss(self.lat) 609 if n > EPS0: 610 # use ratio, datum and NAD83 units may differ 611 e = self.ellipsoid(datum).rocGauss(self.lat) 612 if e > EPS0 and abs(e - n) > EPS: # EPS1 613 m *= e / n 614 meter_text2 = meter_text2.classof(m, t) 615 return self._xnamed(meter_text2) 616 617 @property_doc_(''' this point's reference frame (L{RefFrame}).''') 618 def reframe(self): 619 '''Get this point's reference frame (L{RefFrame}) or C{None}. 620 ''' 621 return self._reframe 622 623 @reframe.setter # PYCHOK setter! 624 def reframe(self, reframe): 625 '''Set or clear this point's reference frame. 626 627 @arg reframe: Reference frame (L{RefFrame}) or C{None}. 628 629 @raise TypeError: The B{C{reframe}} is not a L{RefFrame}. 630 ''' 631 if reframe is not None: 632 from pygeodesy.trf import RefFrame 633 _xinstanceof(RefFrame, reframe=reframe) 634 self._reframe = reframe 635 elif self.reframe is not None: 636 self._reframe = None 637 638 @Property_RO 639 def scale(self): 640 '''Get this point's UTM grid or UPS point scale factor (C{float}) 641 or C{None} if not converted from L{Utm} or L{Ups}. 642 ''' 643 return self._scale 644 645 def toCss(self, **toCss_kwds): 646 '''Convert this C{LatLon} point to a Cassini-Soldner location. 647 648 @kwarg toCss_kwds: Optional L{toCss} keyword arguments. 649 650 @return: The Cassini-Soldner location (L{Css}). 651 652 @see: Function L{toCss}. 653 ''' 654 if toCss_kwds: 655 from pygeodesy.css import toCss 656 r = toCss(self, **_xkwds(toCss_kwds, name=self.name)) 657 else: 658 r = self._css 659 return r 660 661 def toDatum(self, datum2, height=None, name=NN): 662 '''Convert this point to an other datum. 663 664 @arg datum2: Datum to convert I{to} (L{Datum}). 665 @kwarg height: Optional height, overriding the 666 converted height (C{meter}). 667 @kwarg name: Optional name (C{str}), iff converted. 668 669 @return: The converted point (ellipsoidal C{LatLon}) 670 or a copy of this point if B{C{datum2}} 671 matches this point's C{datum}. 672 673 @raise TypeError: Invalid B{C{datum2}}. 674 675 @example: 676 677 >>> p = LatLon(51.4778, -0.0016) # default Datums.WGS84 678 >>> p.toDatum(Datums.OSGB36) # 51.477284°N, 000.00002°E 679 ''' 680 n = name or self.name 681 d2 = _ellipsoidal_datum(datum2, name=n) 682 if self.datum == d2: 683 return self.copy(name=name) 684 r = _xkwds_not(None, epoch=self.epoch, reframe=self.reframe) 685 686 c = self.toCartesian().toDatum(d2) 687 return c.toLatLon(datum=d2, height=height, 688 LatLon=self.classof, name=n, **r) 689 690 def toEtm(self, **toEtm8_kwds): 691 '''Convert this C{LatLon} point to an ETM coordinate. 692 693 @kwarg toEtm8_kwds: Optional L{toEtm8} keyword arguments. 694 695 @return: The ETM coordinate (L{Etm}). 696 697 @see: Function L{toEtm8}. 698 ''' 699 if toEtm8_kwds: 700 from pygeodesy.etm import toEtm8 701 r = toEtm8(self, **_xkwds(toEtm8_kwds, name=self.name)) 702 else: 703 r = self._etm 704 return r 705 706 def toLcc(self, **toLcc_kwds): 707 '''Convert this C{LatLon} point to a Lambert location. 708 709 @kwarg toLcc_kwds: Optional L{toLcc} keyword arguments. 710 711 @return: The Lambert location (L{Lcc}). 712 713 @see: Function L{toLcc}. 714 ''' 715 if toLcc_kwds: 716 from pygeodesy.lcc import toLcc 717 r = toLcc(self, **_xkwds(toLcc_kwds, name=self.name)) 718 else: 719 r = self._lcc 720 return r 721 722 def toMgrs(self, center=False): 723 '''Convert this C{LatLon} point to an MGRS coordinate. 724 725 @kwarg center: If C{True}, try to I{un}-center MGRS 726 to its C{lowerleft} (C{bool}) or by 727 C{B{center} meter} (C{scalar}). 728 729 @return: The MGRS coordinate (L{Mgrs}). 730 731 @see: Method L{Mgrs.toLatLon} and L{toUtm}. 732 ''' 733 return self.toUtm(center=center).toMgrs(center=False) 734 735 def toOsgr(self, **toOsgr_kwds): 736 '''Convert this C{LatLon} point to an OSGR coordinate. 737 738 @kwarg toOsgr_kwds: Optional L{toOsgr} keyword arguments. 739 740 @return: The OSGR coordinate (L{Osgr}). 741 742 @see: Function L{toOsgr}. 743 ''' 744 if toOsgr_kwds: 745 from pygeodesy.osgr import toOsgr 746 r = toOsgr(self, **_xkwds(toOsgr_kwds, name=self.name)) 747 else: 748 r = self._osgr 749 return r 750 751 def toRefFrame(self, reframe2, height=None, name=NN): 752 '''Convert this point to an other reference frame. 753 754 @arg reframe2: Reference frame to convert I{to} (L{RefFrame}). 755 @kwarg height: Optional height, overriding the converted 756 height (C{meter}). 757 @kwarg name: Optional name (C{str}), iff converted. 758 759 @return: The converted point (ellipsoidal C{LatLon}) or this 760 point if conversion is C{nil}, or a copy of this 761 point if the B{C{name}} is non-empty. 762 763 @raise TRFError: This point's C{reframe} is not defined or 764 conversion from this point's C{reframe} to 765 B{C{reframe2}} is not available. 766 767 @raise TypeError: Invalid B{C{reframe2}}, not a L{RefFrame}. 768 769 @example: 770 771 >>> p = LatLon(51.4778, -0.0016, reframe=RefFrames.ETRF2000) # default Datums.WGS84 772 >>> p.toRefFrame(RefFrames.ITRF2014) # 51.477803°N, 000.001597°W, +0.01m 773 >>> p.toRefFrame(RefFrames.ITRF2014, height=0) # 51.477803°N, 000.001597°W 774 ''' 775 if not self.reframe: 776 t = _SPACE_(_DOT_(repr(self), _reframe_), MISSING) 777 raise TRFError(_no_(_conversion_), txt=t) 778 779 from pygeodesy.trf import RefFrame, _reframeTransforms2 780 _xinstanceof(RefFrame, reframe2=reframe2) 781 782 e, xs = _reframeTransforms2(reframe2, self.reframe, self.epoch) 783 if xs: 784 c = self.toCartesian()._applyHelmerts(*xs) 785 n = name or self.name 786 ll = c.toLatLon(datum=self.datum, epoch=e, height=height, 787 LatLon=self.classof, name=n, reframe=reframe2) 788 else: 789 ll = self.copy(name=name) if name else self 790 return ll 791 792 def toUps(self, pole=_N_, falsed=True): 793 '''Convert this C{LatLon} point to a UPS coordinate. 794 795 @kwarg pole: Optional top/center of (stereographic) 796 projection (C{str}, 'N[orth]' or 'S[outh]'). 797 @kwarg falsed: False easting and northing (C{bool}). 798 799 @return: The UPS coordinate (L{Ups}). 800 801 @see: Function L{toUps8}. 802 ''' 803 if self._upsOK(pole, falsed): 804 u = self._ups 805 else: 806 from pygeodesy.ups import toUps8, Ups 807 u = toUps8(self, datum=self.datum, Ups=Ups, 808 pole=pole, falsed=falsed) 809 return u 810 811 def toUtm(self, center=False): 812 '''Convert this C{LatLon} point to a UTM coordinate. 813 814 @kwarg center: If C{True}, I{un}-center the UTM 815 to its C{lowerleft} (C{bool}) or 816 by C{B{center} meter} (C{scalar}). 817 818 @return: The UTM coordinate (L{Utm}). 819 820 @see: Method L{Mgrs.toUtm} and function L{toUtm8}. 821 ''' 822 if center in (False, 0, _0_0): 823 u = self._utm 824 elif center in (True,): 825 u = self._utm._lowerleft 826 else: 827 from pygeodesy.utm import _lowerleft 828 u = _lowerleft(self._utm, center) 829 return u 830 831 def toUtmUps(self, pole=NN): 832 '''Convert this C{LatLon} point to a UTM or UPS coordinate. 833 834 @kwarg pole: Optional top/center of UPS (stereographic) 835 projection (C{str}, 'N[orth]' or 'S[outh]'). 836 837 @return: The UTM or UPS coordinate (L{Utm} or L{Ups}). 838 839 @raise TypeError: Result in L{Utm} or L{Ups}. 840 841 @see: Function L{toUtmUps}. 842 ''' 843 if self._utmOK(): 844 u = self._utm 845 elif self._upsOK(pole): 846 u = self._ups 847 else: # no cover 848 from pygeodesy.utmups import toUtmUps8, Utm, Ups 849 u = toUtmUps8(self, datum=self.datum, Utm=Utm, Ups=Ups, 850 pole=pole, name=self.name) 851 if isinstance(u, Utm): 852 self._overwrite(_utm=u) 853 elif isinstance(u, Ups): 854 self._overwrite(_ups=u) 855 else: 856 _xinstanceof(Utm, Ups, toUtmUps8=u) 857 return u 858 859 def toWm(self, **toWm_kwds): 860 '''Convert this C{LatLon} point to a WM coordinate. 861 862 @kwarg toWm_kwds: Optional L{toWm} keyword arguments. 863 864 @return: The WM coordinate (L{Wm}). 865 866 @see: Function L{toWm}. 867 ''' 868 if toWm_kwds: 869 from pygeodesy.webmercator import toWm 870 r = toWm(self, **_xkwds(toWm_kwds, name=self.name)) 871 else: 872 r = self._wm 873 return r 874 875 @deprecated_method 876 def to3xyz(self): # PYCHOK no cover 877 '''DEPRECATED, use method C{toEcef}. 878 879 @return: A L{Vector3Tuple}C{(x, y, z)}. 880 881 @note: Overloads C{LatLonBase.to3xyz} 882 ''' 883 r = self.toEcef() 884 return Vector3Tuple(r.x, r.y, r.z, name=self.name) 885 886 def trilaterate5(self, distance1, point2, distance2, point3, distance3, 887 area=True, eps=EPS1, wrap=False): 888 '''Trilaterate three points by area overlap or perimeter intersection 889 of three intersecting circles. 890 891 @arg distance1: Distance to this point (C{meter}), same units 892 as B{C{eps}}). 893 @arg point2: Second center point (C{LatLon}). 894 @arg distance2: Distance to point2 (C{meter}, same units as 895 B{C{eps}}). 896 @arg point3: Third center point (C{LatLon}). 897 @arg distance3: Distance to point3 (C{meter}, same units as 898 B{C{eps}}). 899 @kwarg area: If C{True} compute the area overlap, otherwise the 900 perimeter intersection of the circles (C{bool}). 901 @kwarg eps: The required I{minimal overlap} for C{B{area}=True} 902 or the I{intersection margin} for C{B{area}=False} 903 (C{meter}, conventionally). 904 @kwarg wrap: Wrap/unroll angular distances (C{bool}). 905 906 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)} 907 with C{min} and C{max} in C{meter}, same units as B{C{eps}}, 908 the corresponding trilaterated points C{minPoint} and 909 C{maxPoint} as I{ellipsoidal} C{LatLon} and C{n}, the number 910 of trilatered points found for the given B{C{eps}}. 911 912 If only a single trilaterated point is found, C{min I{is} 913 max}, C{minPoint I{is} maxPoint} and C{n = 1}. 914 915 For C{B{area}=True}, C{min} and C{max} are the smallest 916 respectively largest I{radial} overlap found. 917 918 For C{B{area}=False}, C{min} and C{max} represent the 919 nearest respectively farthest intersection margin. 920 921 If C{B{area}=True} and all 3 circles are concentric, C{n=0} 922 and C{minPoint} and C{maxPoint} are the B{C{point#}} with 923 the smallest B{C{distance#}} C{min} respectively C{max} the 924 largest B{C{distance#}}. 925 926 @raise IntersectionError: Trilateration failed for the given B{C{eps}}, 927 insufficient overlap for C{B{area}=True} or 928 no intersection or all (near-)concentric for 929 C{B{area}=False}. 930 931 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}. 932 933 @raise ValueError: Coincident B{C{points}} or invalid B{C{distance1}}, 934 B{C{distance2}} or B{C{distance3}}. 935 936 @note: Ellipsoidal trilateration invokes methods C{LatLon.intersections2} 937 and C{LatLon.nearestOn} based on I{Karney}'s Python U{geographiclib 938 <https://PyPI.org/project/geographiclib>} if installed, otherwise 939 using the accurate (but slower) C{ellipsoidalExact.LatLon} methods. 940 ''' 941 return _trilaterate5(self, distance1, 942 self.others(point2=point2), distance2, 943 self.others(point3=point3), distance3, 944 area=area, eps=eps, wrap=wrap) 945 946 @Property_RO 947 def _ups(self): # __dict__ value overwritten by method C{toUtmUps} 948 '''(INTERNAL) Get this C{LatLon} point as UPS coordinate (L{Ups}), see L{toUps8}. 949 ''' 950 from pygeodesy.ups import toUps8, Ups 951 return toUps8(self, datum=self.datum, Ups=Ups, 952 pole=NN, falsed=True, name=self.name) 953 954 def _upsOK(self, pole=NN, falsed=True): 955 '''(INTERNAL) Check matching C{Ups}. 956 ''' 957 try: 958 u = self._ups 959 except RangeError: 960 return False 961 return falsed and (u.pole == pole[:1].upper() or not pole) 962 963 @Property_RO 964 def _utm(self): # __dict__ value overwritten by method C{toUtmUps} 965 '''(INTERNAL) Get this C{LatLon} point as UTM coordinate (L{Utm}), see L{toUtm8}. 966 ''' 967 from pygeodesy.utm import toUtm8, Utm 968 return toUtm8(self, datum=self.datum, Utm=Utm, name=self.name) 969 970 def _utmOK(self): 971 '''(INTERNAL) Check C{Utm}. 972 ''' 973 try: 974 _ = self._utm 975 except RangeError: 976 return False 977 return True 978 979 @Property_RO 980 def _wm(self): 981 '''(INTERNAL) Get this C{LatLon} point as webmercator (L{Wm}). 982 ''' 983 from pygeodesy.webmercator import toWm 984 return toWm(self) 985 986 987__all__ += _ALL_DOCS(CartesianEllipsoidalBase, LatLonEllipsoidalBase) 988 989# **) MIT License 990# 991# Copyright (C) 2016-2021 -- mrJean1 at Gmail -- All Rights Reserved. 992# 993# Permission is hereby granted, free of charge, to any person obtaining a 994# copy of this software and associated documentation files (the "Software"), 995# to deal in the Software without restriction, including without limitation 996# the rights to use, copy, modify, merge, publish, distribute, sublicense, 997# and/or sell copies of the Software, and to permit persons to whom the 998# Software is furnished to do so, subject to the following conditions: 999# 1000# The above copyright notice and this permission notice shall be included 1001# in all copies or substantial portions of the Software. 1002# 1003# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 1004# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 1005# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 1006# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 1007# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 1008# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 1009# OTHER DEALINGS IN THE SOFTWARE. 1010