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