-
Notifications
You must be signed in to change notification settings - Fork 57
/
sphericalNvector.py
1304 lines (996 loc) · 52.5 KB
/
sphericalNvector.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# -*- coding: utf-8 -*-
u'''N-vector-based classes geodetic (lat-/longitude) L{LatLon}, geocentric
(ECEF) L{Cartesian} and L{Nvector} and functions L{areaOf}, L{intersection},
L{meanOf}, L{nearestOn2}, L{triangulate} and L{trilaterate}, I{all spherical}.
Pure Python implementation of n-vector-based spherical geodetic (lat-/longitude)
methods, transcribed from JavaScript originals by I{(C) Chris Veness 2011-2016},
published under the same MIT Licence**. See U{Vector-based geodesy
<https://www.Movable-Type.co.UK/scripts/latlong-vectors.html>} and
U{Module latlon-nvector-spherical
<https://www.Movable-Type.co.UK/scripts/geodesy/docs/module-latlon-nvector-spherical.html>}.
Tools for working with points and paths on (a spherical model of) the
earth’s surface using using n-vectors rather than the more common
spherical trigonometry. N-vectors make many calculations much simpler,
and easier to follow, compared with the trigonometric equivalents.
Based on Kenneth Gade’s U{‘Non-singular Horizontal Position Representation’
<https://www.NavLab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf>},
The Journal of Navigation (2010), vol 63, nr 3, pp 395-417.
Note that the formulations below take x => 0°N,0°E, y => 0°N,90°E and
z => 90°N while Gade uses x => 90°N, y => 0°N,90°E, z => 0°N,0°E.
Also note that on a spherical earth model, an n-vector is equivalent
to a normalised version of an (ECEF) cartesian coordinate.
@newfield example: Example, Examples
'''
from pygeodesy.basics import EPS, EPS_2, PI, PI2, PI_2, R_M, \
isscalar, map1, _xinstanceof, _xkwds
from pygeodesy.datum import Datums
from pygeodesy.ecef import EcefKarney
from pygeodesy.errors import _ValueError
from pygeodesy.fmath import fidw, fmean, fsum, fsum_
from pygeodesy.interns import _1_, _2_, _bearing_, _coincident_, \
_distance_, _end_, _fraction_, \
_other_, _point_, _points_, _pole_, \
_start_, _start1_, _start2_
from pygeodesy.lazily import _ALL_LAZY, _ALL_OTHER
from pygeodesy.named import NearestOn3Tuple
from pygeodesy.nvectorBase import NvectorBase, NorthPole, LatLonNvectorBase, \
sumOf as _sumOf
from pygeodesy.points import _imdex2, ispolar # PYCHOK exported
from pygeodesy.sphericalBase import _angular, CartesianSphericalBase, \
LatLonSphericalBase
from pygeodesy.streprs import unstr
from pygeodesy.units import Bearing, Bearing_, Height, Radius, Radius_, Scalar
from pygeodesy.utily import degrees360, iterNumpy2, sincos2, sincos2d
from math import atan2, fabs, sqrt
__all__ = _ALL_LAZY.sphericalNvector
__version__ = '20.07.24'
_paths_ = 'paths'
class Cartesian(CartesianSphericalBase):
'''Extended to convert geocentric, L{Cartesian} points to
L{Nvector} and n-vector-based, spherical L{LatLon}.
'''
def toLatLon(self, **LatLon_datum_kwds): # PYCHOK LatLon=LatLon
'''Convert this cartesian to an C{Nvector}-based geodetic point.
@kwarg LatLon_datum_kwds: Optional L{LatLon}, B{C{datum}} and
other keyword arguments, ignored if B{C{LatLon=None}}.
Use B{C{LatLon=...}} to override the L{LatLon} class
or specify B{C{LatLon=None}}.
@return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is
C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
C, M, datum)} with C{C} and C{M} if available.
@raise TypeError: Invalid B{C{LatLon_datum_kwds}}.
'''
kwds = _xkwds(LatLon_datum_kwds, LatLon=LatLon, datum=self.datum)
return CartesianSphericalBase.toLatLon(self, **kwds)
def toNvector(self, **Nvector_datum_kwds): # PYCHOK Datums.WGS84
'''Convert this cartesian to L{Nvector} components,
I{including height}.
@kwarg Nvector_datum_kwds: Optional L{Nvector}, B{C{datum}} and
other keyword arguments, ignored if B{C{Nvector=None}}.
Use B{C{Nvector=...}} to override the L{Nvector} class
or specify B{C{Nvector=None}}.
@return: The C{n-vector} components (L{Nvector}) or a
L{Vector4Tuple}C{(x, y, z, h)} if B{C{Nvector}}
is C{None}.
@raise TypeError: Invalid B{C{Nvector_datum_kwds}}.
'''
# ll = CartesianBase.toLatLon(self, LatLon=LatLon,
# datum=datum or self.datum)
# kwds = _xkwds(kwds, Nvector=Nvector)
# return ll.toNvector(**kwds)
kwds = _xkwds(Nvector_datum_kwds, Nvector=Nvector, datum=self.datum)
return CartesianSphericalBase.toNvector(self, **kwds)
class LatLon(LatLonNvectorBase, LatLonSphericalBase):
'''New n-vector based point on a spherical earth model.
Tools for working with points and paths on (a spherical
model of) the earth's surface using vector-based methods.
@example:
>>> from sphericalNvector import LatLon
>>> p = LatLon(52.205, 0.119)
'''
_Nv = None #: (INTERNAL) cache _toNvector L{Nvector}).
def _gc3(self, start, end, namend, raiser=_points_):
'''(INTERNAL) Return great circle, start and end Nvectors.
'''
s = start.toNvector()
if isscalar(end): # bearing
gc = s.greatCircle(end)
e = None
else:
self.others(end, name=namend)
e = end.toNvector()
gc = s.cross(e, raiser=raiser) # XXX .unit()?
return gc, s, e
def _update(self, updated, *attrs): # PYCHOK args
'''(INTERNAL) Zap cached attributes if updated.
'''
if updated: # reset caches
LatLonNvectorBase._update(self, updated, _Nv=self._Nv) # special case
LatLonSphericalBase._update(self, updated, *attrs)
def alongTrackDistanceTo(self, start, end, radius=R_M):
'''Compute the (signed) distance from the start to the closest
point on the great circle path defined by a start and an
end point.
That is, if a perpendicular is drawn from this point to the
great circle path, the along-track distance is the distance
from the start point to the point where the perpendicular
crosses the path.
@arg start: Start point of great circle path (L{LatLon}).
@arg end: End point of great circle path (L{LatLon}) or
initial bearing from start point (compass
C{degrees360}).
@kwarg radius: Mean earth radius (C{meter}).
@return: Distance along the great circle path (positive if
after the start toward the end point of the path
or negative if before the start point).
@raise TypeError: If B{C{start}} or B{C{end}} point is not L{LatLon}.
@raise Valuerror: Some points coincide.
@example:
>>> p = LatLon(53.2611, -0.7972)
>>> s = LatLon(53.3206, -1.7297)
>>> e = LatLon(53.1887, 0.1334)
>>> d = p.alongTrackDistanceTo(s, e) # 62331.58
'''
self.others(start, name=_start_)
gc, _, _ = self._gc3(start, end, _end_)
p = self.toNvector()
a = gc.cross(p).cross(gc) # along-track point gc × p × gc
return start.toNvector().angleTo(a, vSign=gc) * radius
def bearingTo(self, other, **unused): # PYCHOK no cover
'''DEPRECATED, use method C{initialBearingTo}.
'''
return self.initialBearingTo(other)
def crossTrackDistanceTo(self, start, end, radius=R_M):
'''Compute the (signed) distance from this point to great circle
defined by a start and end point.
@arg start: Start point of great circle path (L{LatLon}).
@arg end: End point of great circle path (L{LatLon}) or
initial bearing from start point (compass
C{degrees360}).
@kwarg radius: Mean earth radius (C{meter}).
@return: Distance to great circle (negative if to the
left or positive if to the right of the path).
@raise TypeError: If B{C{start}} or B{C{end}} point is not L{LatLon}.
@raise Valuerror: Some points coincide.
@example:
>>> p = LatLon(53.2611, -0.7972)
>>> s = LatLon(53.3206, -1.7297)
>>> d = p.crossTrackDistanceTo(s, 96) # -305.7
>>> e = LatLon(53.1887, 0.1334)
>>> d = p.crossTrackDistanceTo(s, e) # -307.5
'''
self.others(start, name=_start_)
gc, _, _ = self._gc3(start, end, _end_)
p = self.toNvector()
return (gc.angleTo(p) - PI_2) * radius
def destination(self, distance, bearing, radius=R_M, height=None):
'''Locate the destination from this point after having
travelled the given distance on the given bearing.
@arg distance: Distance travelled (C{meter}, same units
as B{C{radius}}).
@arg bearing: Bearing from this point (compass C{degrees360}).
@kwarg radius: Mean earth radius (C{meter}).
@kwarg height: Optional height at destination, overriding the
default height (C{meter}, same units as B{C{radius}}).
@return: Destination point (L{LatLon}).
@raise Valuerror: Polar coincidence ior invalid B{C{distance}},
B{C{bearing}}, B{C{radius}} or B{C{height}}.
@example:
>>> p = LatLon(51.4778, -0.0015)
>>> q = p.destination(7794, 300.7)
>>> q.toStr() # 51.513546°N, 000.098345°W
@JSname: I{destinationPoint}.
'''
a = _angular(distance, radius)
sa, ca, sb, cb = sincos2(a, Bearing_(bearing))
p = self.toNvector()
e = NorthPole.cross(p, raiser=_pole_).unit() # east vector at p
n = p.cross(e) # north vector at p
q = n.times(cb).plus(e.times(sb)) # direction vector @ p
n = p.times(ca).plus(q.times(sa))
return n.toLatLon(height=height, LatLon=self.classof) # Nvector(n.x, n.y, n.z).toLatLon(...)
def distanceTo(self, other, radius=R_M, **unused): # for -DistanceTo
'''Compute the distance from this to an other point.
@arg other: The other point (L{LatLon}).
@kwarg radius: Mean earth radius (C{meter}).
@return: Distance between this and the B{C{other}} point
(C{meter}, same units as B{C{radius}}).
@raise TypeError: The B{C{other}} point is not L{LatLon}.
@example:
>>> p = LatLon(52.205, 0.119)
>>> q = LatLon(48.857, 2.351);
>>> d = p.distanceTo(q) # 404.3 km
'''
self.others(other)
return self.toNvector().angleTo(other.toNvector()) * Radius(radius)
def greatCircle(self, bearing):
'''Compute the vector normal to great circle obtained by
heading on the given bearing from this point.
Direction of vector is such that initial bearing vector
b = c × n, where n is an n-vector representing this point.
@arg bearing: Bearing from this point (compass C{degrees360}).
@return: N-vector representing the great circle (L{Nvector}).
'''
a, b = self.philam
t = Bearing_(bearing)
sa, ca, sb, cb, st, ct = sincos2(a, b, t)
return Nvector(sb * ct - sa * cb * st,
-cb * ct - sa * sb * st,
ca * st, name=self.name) # XXX .unit()
def greatCircleTo(self, other):
'''Compute the vector normal to great circle obtained by
heading from this to an other point or on a given bearing.
Direction of vector is such that initial bearing vector
b = c × n, where n is an n-vector representing this point.
@arg other: The other point (L{LatLon}) or the bearing from
this point (compass C{degrees360}).
@return: N-vector representing the great circle (L{Nvector}).
@raise TypeError: The B{C{other}} point is not L{LatLon}.
@raise Valuerror: Points coincide.
@example:
>>> p = LatLon(53.3206, -1.7297)
>>> gc = p.greatCircle(96.0)
>>> gc.toStr() # (-0.79408, 0.12856, 0.59406)
>>> q = LatLon(53.1887, 0.1334)
>>> g = p.greatCircleTo(q)
>>> g.toStr() # (-0.79408, 0.12859, 0.59406)
'''
gc, _, _ = self._gc3(self, other, _other_)
return gc.unit()
def initialBearingTo(self, other, **unused):
'''Compute the initial bearing (forward azimuth) from this
to an other point.
@arg other: The other point (L{LatLon}).
@arg unused: Optional keyword argument B{C{wrap}} ignored.
@return: Initial bearing (compass C{degrees360}).
@raise Crosserror: This point coincides with the B{C{other}}
point or the C{NorthPole}, provided
L{crosserrors} is C{True}.
@raise TypeError: The B{C{other}} point is not L{LatLon}.
@example:
>>> p1 = LatLon(52.205, 0.119)
>>> p2 = LatLon(48.857, 2.351)
>>> b = p1.initialBearingTo(p2) # 156.2
@JSname: I{bearingTo}.
'''
self.others(other, name=_other_)
# see <https://MathForum.org/library/drmath/view/55417.html>
n = self.toNvector()
# gc1 = self.greatCircleTo(other)
gc1 = n.cross(other.toNvector(), raiser=_points_) # .unit()
# gc2 = self.greatCircleTo(NorthPole)
gc2 = n.cross(NorthPole, raiser=_pole_) # .unit()
return degrees360(gc1.angleTo(gc2, vSign=n))
def intermediateChordTo(self, other, fraction, height=None):
'''Locate the point projected from the point at given fraction
on a straight line (chord) between this and an other point.
@arg other: The other point (L{LatLon}).
@arg fraction: Fraction between both points (float, between
0.0 for this and 1.0 for the other point).
@kwarg height: Optional height at the intermediate point,
overriding the fractional height (C{meter}).
@return: Intermediate point (L{LatLon}).
@raise TypeError: The B{C{other}} point is not L{LatLon}.
@example:
>>> p = LatLon(52.205, 0.119)
>>> q = LatLon(48.857, 2.351)
>>> i = p.intermediateChordTo(q, 0.25) # 51.3723°N, 000.7072°E
@JSname: I{intermediatePointOnChordTo}, I{intermediatePointDirectlyTo}.
'''
self.others(other)
f = Scalar(fraction, name=_fraction_)
i = other.toNvector().times(f).plus(
self.toNvector().times(1 - f))
# i = other.toNvector() * f + \
# self.toNvector() * (1 - f))
h = self._havg(other, f=f) if height is None else Height(height)
return i.toLatLon(height=h, LatLon=self.classof) # Nvector(i.x, i.y, i.z).toLatLon(...)
def intermediateTo(self, other, fraction, height=None):
'''Locate the point at a given fraction between this and an
other point.
@arg other: The other point (L{LatLon}).
@arg fraction: Fraction between both points (float, between
0.0 for this and 1.0 for the other point).
@kwarg height: Optional height at the intermediate point,
overriding the fractional height (C{meter}).
@return: Intermediate point (L{LatLon}).
@raise TypeError: The B{C{other}} point is not L{LatLon}.
@raise Valuerror: Points coincide or invalid B{C{height}}.
@example:
>>> p = LatLon(52.205, 0.119)
>>> q = LatLon(48.857, 2.351)
>>> i = p.intermediateTo(q, 0.25) # 51.3721°N, 000.7074°E
@JSname: I{intermediatePointTo}.
'''
q = self.others(other).toNvector()
p = self.toNvector()
f = Scalar(fraction, name=_fraction_)
x = p.cross(q, raiser=_points_)
d = x.unit().cross(p) # unit(p × q) × p
# angular distance α, tan(α) = |p × q| / p ⋅ q
s, c = sincos2(atan2(x.length, p.dot(q)) * f) # interpolated
i = p.times(c).plus(d.times(s)) # p * cosα + d * sinα
h = self._havg(other, f=f) if height is None else Height(height)
return i.toLatLon(height=h, LatLon=self.classof) # Nvector(i.x, i.y, i.z).toLatLon(...)
def intersection(self, end1, start2, end2, height=None):
'''Locate the intersection point of two paths each defined
by two points or a start point and bearing from North.
@arg end1: End point of the first path (L{LatLon}) or the
initial bearing at this point (compass C{degrees360}).
@arg start2: Start point of the second path (L{LatLon}).
@arg end2: End point of the second path (L{LatLon}) or the
initial bearing at the second point (compass
C{degrees}).
@kwarg height: Optional height at the intersection point,
overriding the mean height (C{meter}).
@return: The intersection point (L{LatLon}) or C{None}
if no unique intersection exists.
@raise TypeError: If B{C{start2}}, B{C{end1}} or B{C{end2}} point
is not L{LatLon}.
@raise ValueError: Intersection is ambiguous or infinite or
the paths are parallel, coincident or null.
@example:
>>> s = LatLon(51.8853, 0.2545)
>>> e = LatLon(49.0034, 2.5735)
>>> i = s.intersection(108.55, e, 32.44) # 50.9076°N, 004.5086°E
'''
return intersection(self, end1, start2, end2,
height=height, LatLon=self.classof)
def isenclosedBy(self, points):
'''Check whether this point is enclosed by a (convex) polygon.
@arg points: The polygon points (L{LatLon}[]).
@return: C{True} if the polygon encloses this point,
C{False} otherwise.
@raise PointsError: Insufficient number of B{C{points}}.
@raise TypeError: Some B{C{points}} are not L{LatLon}.
@example:
>>> b = LatLon(45,1), LatLon(45,2), LatLon(46,2), LatLon(46,1)
>>> p = LatLon(45.1, 1.1)
>>> inside = p.isEnclosedBy(b) # True
@JSname: I{enclosedBy}.
'''
n, points = self.points2(points, closed=True)
# use normal vector to this point for sign of α
n0 = self.toNvector()
if iterNumpy2(points):
def _subtangles(n, points, n0): # iterate
vs = n0.minus(points[n-1].toNvector())
for i in range(n):
vs1 = n0.minus(points[i].toNvector())
yield vs.angleTo(vs1, vSign=n0) # PYCHOK false
vs = vs1
# sum subtended angles
s = fsum(_subtangles(n, points, n0))
else:
# get vectors from this to each point
vs = [n0.minus(points[i].toNvector()) for i in range(n)]
# close the polygon
vs.append(vs[0])
# sum subtended angles of each edge (using n0 to determine sign)
s = fsum(vs[i].angleTo(vs[i+1], vSign=n0) for i in range(n))
# Note, this method uses angle summation test: on a plane,
# angles for an enclosed point will sum to 360°, angles for
# an exterior point will sum to 0°. On a sphere, enclosed
# point angles will sum to less than 360° (due to spherical
# excess), exterior point angles will be small but non-zero.
# XXX are winding number optimisations equally applicable to
# spherical surface?
return abs(s) > PI
def isEnclosedBy(self, points): # PYCHOK no cover
'''DEPRECATED, use method C{isenclosedBy}.
'''
return self.isenclosedBy(points)
def iswithin(self, point1, point2):
'''Check whether this point is between two other points.
If this point is not on the great circle arc defined by
both points, return whether it is within the area bound
by perpendiculars to the great circle at each point (in
the same hemispere).
@arg point1: Start point of the arc (L{LatLon}).
@arg point2: End point of the arc (L{LatLon}).
@return: C{True} if this point is within the arc,
C{False} otherwise.
@raise TypeError: If B{C{point1}} or B{C{point2}} is not L{LatLon}.
@JSname: I{isBetween}.
'''
n0 = self.toNvector()
n1 = self.others(point1, name='point1').toNvector()
n2 = self.others(point2, name='point2').toNvector()
# corner case, null arc
if n1.isequalTo(n2):
return n0.isequalTo(n1) or n0.isequalTo(n2) # PYCHOK returns
if n0.dot(n1) < 0 or n0.dot(n2) < 0: # different hemisphere
return False # PYCHOK returns
# get vectors representing d0=p0->p1 and d2=p2->p1
# and dot product d0⋅d2 tells us if p0 is on the
# p2 side of p1 or on the other side (similarly
# for d0=p0->p2 and d1=p1->p2 and dot product
# d0⋅d1 and p0 on the p1 side of p2 or not)
return n0.minus(n1).dot(n2.minus(n1)) >= 0 and \
n0.minus(n2).dot(n1.minus(n2)) >= 0
def isWithin(self, point1, point2): # PYCHOK no cover
'''DEPRECATED, use method C{iswithin}.
'''
return self.iswithin(point1, point2)
def midpointTo(self, other, height=None):
'''Find the midpoint between this and an other point.
@arg other: The other point (L{LatLon}).
@kwarg height: Optional height at the midpoint, overriding
the mean height (C{meter}).
@return: Midpoint (L{LatLon}).
@raise TypeError: The B{C{other}} point is not L{LatLon}.
@example:
>>> p1 = LatLon(52.205, 0.119)
>>> p2 = LatLon(48.857, 2.351)
>>> m = p1.midpointTo(p2) # '50.5363°N, 001.2746°E'
'''
self.others(other)
m = self.toNvector().plus(other.toNvector())
h = self._havg(other) if height is None else height
return m.toLatLon(height=h, LatLon=self.classof)
def nearestOn(self, point1, point2, height=None):
'''Locate the point on the great circle arc between two
points closest to this point.
If this point is within the extent of the arc between both
end points, return the closest point on the arc. Otherwise,
return the closest of the arc's end points.
@arg point1: Start point of the arc (L{LatLon}).
@arg point2: End point of the arc (L{LatLon}).
@kwarg height: Optional height, overriding the mean height
for the point within the arc (C{meter}).
@return: Closest point on the arc (L{LatLon}).
@raise TypeError: If B{C{point1}} or B{C{point2}} is not L{LatLon}.
@example:
>>> s1 = LatLon(51.0, 1.0)
>>> s2 = LatLon(51.0, 2.0)
>>> s = LatLon(51.0, 1.9)
>>> p = s.nearestOn(s1, s2) # 51.0004°N, 001.9000°E
>>> d = p.distanceTo(s) # 42.71 m
>>> s = LatLon(51.0, 2.1)
>>> p = s.nearestOn(s1, s2) # 51.0000°N, 002.0000°E
@JSname: I{nearestPointOnSegment}.
'''
if self.isWithin(point1, point2) and not point1.isequalTo(point2, EPS):
# closer to arc than to its endpoints,
# find the closest point on the arc
gc1 = point1.toNvector().cross(point2.toNvector())
gc2 = self.toNvector().cross(gc1)
p = gc1.cross(gc2).toLatLon(height=height or 0, LatLon=self.classof)
if height is None: # interpolate height
d = point1.distanceTo(point2)
f = 0.5 if d < EPS else point1.distanceTo(p) / d
p.height = point1._havg(point2, f=f)
# beyond arc extent, take closer endpoint
elif self.distanceTo(point1) < self.distanceTo(point2):
p = point1
else:
p = point2
return p
def nearestOn2(self, points, **closed_radius_height): # PYCHOK no cover
'''DEPRECATED, use method L{sphericalNvector.LatLon.nearestOn3}.
@return: ... 2-Tuple C{(closest, distance)} of the C{closest}
point (L{LatLon}) on the polygon and the C{distance}
to that point from this point ...
'''
r = self.nearestOn3(points, **closed_radius_height)
return r.closest, r.distance
def nearestOn3(self, points, closed=False, radius=R_M, height=None):
'''Locate the point on a polygon (with great circle arcs
joining consecutive points) closest to this point.
If this point is within the extent of any great circle
arc, the closest point is on that arc. Otherwise,
the closest is the nearest of the arc's end points.
@arg points: The polygon points (L{LatLon}[]).
@kwarg closed: Optionally, close the polygon (C{bool}).
@kwarg radius: Mean earth radius (C{meter}).
@kwarg height: Optional height, overriding the mean height
for a point within the arc (C{meter}).
@return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of
the C{closest} point (L{LatLon}), the C{distance}
between this and the C{closest} point in C{meter},
same units as B{C{radius}} and the C{angle} from
this to the C{closest} point in compass C{degrees360}.
@raise TypeError: Some B{C{points}} are not C{LatLon}.
@raise ValueError: No B{C{points}}.
'''
n, points = self.points2(points, closed=closed)
i, m = _imdex2(closed, n)
c = p2 = points[i]
r = self.distanceTo(c, radius=1) # force radians
for i in range(m, n):
p1, p2 = p2, points[i]
p = self.nearestOn(p1, p2, height=height)
t = self.distanceTo(p, radius=1) # force radians
if t < r:
c, r = p, t
return NearestOn3Tuple(c, r * Radius(radius), degrees360(r))
def toCartesian(self, **Cartesian_datum_kwds): # PYCHOK Cartesian=Cartesian, datum=None
'''Convert this point to C{Nvector}-based cartesian (ECEF)
coordinates.
@kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}} or
other keyword arguments, ignored if
B{C{Cartesian=None}}. Use
B{C{Cartesian=...}} to override this
L{Cartesian} class or specify
B{C{Cartesian=None}}.
@return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}}
is C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
C, M, datum)} with C{C} and C{M} if available.
@raise TypeError: Invalid B{C{Cartesian_datum_kwds}}.
'''
kwds = _xkwds(Cartesian_datum_kwds, Cartesian=Cartesian,
datum=self.datum)
return LatLonSphericalBase.toCartesian(self, **kwds)
def toNvector(self, **Nvector_kwds): # PYCHOK signature
'''Convert this point to L{Nvector} components, I{including
height}.
@kwarg Nvector_kwds: Optional L{Nvector} keyword arguments,
ignored if B{C{Nvector=None}}. Use
B{C{Nvector=...}} to override this
L{Nvector} class or specify
B{C{Nvector=None}}
@return: The C{n-vector} components (L{Nvector}) or a
L{Vector4Tuple}C{(x, y, z, h)} if B{C{Nvector}}
is C{None}.
@raise TypeError: Invalid B{C{Nvector_kwds}}.
@example:
>>> p = LatLon(45, 45)
>>> n = p.toNvector()
>>> n.toStr() # [0.50, 0.50, 0.70710]
@JSname: I{toVector}.
'''
kwds = _xkwds(Nvector_kwds, Nvector=Nvector)
return LatLonNvectorBase.toNvector(self, **kwds)
def triangulate(self, bearing1, other, bearing2, height=None):
'''Locate a point given this and an other point and a bearing
at this and the other point.
@arg bearing1: Bearing at this point (compass C{degrees360}).
@arg other: The other point (L{LatLon}).
@arg bearing2: Bearing at the other point (compass C{degrees360}).
@kwarg height: Optional height at the triangulated point,
overriding the mean height (C{meter}).
@return: Triangulated point (L{LatLon}).
@raise TypeError: The B{C{other}} point is not L{LatLon}.
@raise Valuerror: Points coincide.
@example:
>>> p = LatLon("47°18.228'N","002°34.326'W") # Basse Castouillet
>>> q = LatLon("47°18.664'N","002°31.717'W") # Basse Hergo
>>> t = p.triangulate(7, q, 295) # 47.323667°N, 002.568501°W'
'''
return triangulate(self, bearing1, other, bearing2,
height=height, LatLon=self.classof)
def trilaterate(self, distance1, point2, distance2, point3, distance3,
radius=R_M, height=None, useZ=False):
'''Locate a point at given distances from this and two other points.
See also U{Trilateration<https://WikiPedia.org/wiki/Trilateration>}.
@arg distance1: Distance to this point (C{meter}, same units
as B{C{radius}}).
@arg point2: Second reference point (L{LatLon}).
@arg distance2: Distance to point2 (C{meter}, same units as
B{C{radius}}).
@arg point3: Third reference point (L{LatLon}).
@arg distance3: Distance to point3 (C{meter}, same units as
B{C{radius}}).
@kwarg radius: Mean earth radius (C{meter}).
@kwarg height: Optional height at trilaterated point, overriding
the mean height (C{meter}, same units as B{C{radius}}).
@kwarg useZ: Include Z component iff non-NaN, non-zero (C{bool}).
@return: Trilaterated point (L{LatLon}).
@raise TypeError: Some B{C{points}} are not L{LatLon}.
@raise ValueError: Distance(s) exceeds trilateration or
some B{C{points}} coincide.
'''
return trilaterate(self, distance1, point2, distance2,
point3, distance3,
radius=radius, height=height,
LatLon=self.classof, useZ=useZ)
class Nvector(NvectorBase):
'''An n-vector is a position representation using a (unit) vector
normal to the earth's surface. Unlike lat-/longitude points,
n-vectors have no singularities or discontinuities.
For many applications, n-vectors are more convenient to work
with than other position representations like lat-/longitude,
earth-centred earth-fixed (ECEF) vectors, UTM coordinates, etc.
On a spherical model earth, an n-vector is equivalent to an
earth-centred earth-fixed (ECEF) vector.
Note commonality with L{ellipsoidalNvector.Nvector}.
'''
_datum = Datums.Sphere #: (INTERNAL) Default datum (L{Datum}).
_Ecef = EcefKarney #: (INTERNAL) Preferred C{Ecef...} class.
def toCartesian(self, **Cartesian_h_kwds): # PYCHOK Cartesian=Cartesian
'''Convert this n-vector to C{Nvector}-based cartesian
(ECEF) coordinates.
@kwarg Cartesian_h_kwds: Optional L{Cartesian}, B{C{h}} and
other keyword arguments, ignored if
B{C{Cartesian=None}}. Use
B{C{Cartesian=...}} to override this
L{Cartesian} class or specify
B{C{Cartesian=None}}.
@return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}}
is C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
C, M, datum)} with C{C} and C{M} if available.
@raise TypeError: Invalid B{C{Cartesian}}, B{C{h}} or other
B{C{Cartesian_h_kwds}}.
'''
kwds = _xkwds(Cartesian_h_kwds, h=self.h, Cartesian=Cartesian)
return NvectorBase.toCartesian(self, **kwds) # class or .classof
def toLatLon(self, **LatLon_height_kwds): # PYCHOK height=None, LatLon=LatLon
'''Convert this n-vector to an C{Nvector}-based geodetic point.
@kwarg LatLon_height_kwds: Optional L{LatLon}, B{C{height}} and
other keyword arguments, ignored if
B{C{LatLon=None}}. Use
B{C{LatLon=...}} to override this
L{LatLon} class or specify
B{C{LatLon=None}}.
@return: The geodetic point (L{LatLon}) or if B{C{LatLon}}
is C{None}, a L{LatLon3Tuple}C{(lat, lon, height)}.
@raise TypeError: Invalid B{C{LatLon}}, B{C{height}} or
other B{C{LatLon_height_kwds}}.
@raise ValueError: Invalid B{C{height}}.
'''
kwds = _xkwds(LatLon_height_kwds, height=self.h, LatLon=LatLon)
return NvectorBase.toLatLon(self, **kwds) # class or .classof
def greatCircle(self, bearing):
'''Compute the n-vector normal to great circle obtained by
heading on given compass bearing from this point as its
n-vector.
Direction of vector is such that initial bearing vector
b = c × p.
@arg bearing: Initial compass bearing (C{degrees}).
@return: N-vector representing great circle (L{Nvector}).
@raise Valuerror: Polar coincidence.
@example:
>>> n = LatLon(53.3206, -1.7297).toNvector()
>>> gc = n.greatCircle(96.0) # [-0.794, 0.129, 0.594]
'''
s, c = sincos2d(Bearing(bearing))
e = NorthPole.cross(self, raiser=_pole_) # easting
n = self.cross(e, raiser=_point_) # northing
e = e.times(c / e.length)
n = n.times(s / n.length)
return n.minus(e)
_Nvll = LatLon(0, 0, name='Nv00') #: (INTERNAL) Reference instance (L{LatLon}).
def areaOf(points, radius=R_M):
'''Calculate the area of a (spherical) polygon (with great circle
arcs joining consecutive points).
@arg points: The polygon points (L{LatLon}[]).
@kwarg radius: Mean earth radius (C{meter}).
@return: Polygon area (C{meter}, same units as B{C{radius}}, squared).
@raise PointsError: Insufficient number of B{C{points}}.
@raise TypeError: Some B{C{points}} are not L{LatLon}.
@see: L{pygeodesy.areaOf}, L{sphericalTrigonometry.areaOf} and
L{ellipsoidalKarney.areaOf}.
@example:
>>> b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1)
>>> areaOf(b) # 8666058750.718977
'''
n, points = _Nvll.points2(points, closed=True)
# use vector to 1st point as plane normal for sign of α
n0 = points[0].toNvector()
if iterNumpy2(points):
def _interangles(n, points, n0): # iterate
v2 = points[n-2]._N_vector
v1 = points[n-1]._N_vector
gc = v2.cross(v1)
for i in range(n):
v2 = points[i]._N_vector
gc1 = v1.cross(v2)
v1 = v2
yield gc.angleTo(gc1, vSign=n0)
gc = gc1
# sum interior angles
s = fsum(_interangles(n, points, n0))
else:
# get great-circle vector for each edge
gc, v1 = [], points[n-1]._N_vector
for i in range(n):
v2 = points[i]._N_vector
gc.append(v1.cross(v2)) # PYCHOK false, does have .cross
v1 = v2
gc.append(gc[0]) # XXX needed?
# sum interior angles: depending on whether polygon is cw or ccw,
# angle between edges is π−α or π+α, where α is angle between
# great-circle vectors; so sum α, then take n·π − |Σα| (cannot
# use Σ(π−|α|) as concave polygons would fail)
s = fsum(gc[i].angleTo(gc[i+1], vSign=n0) for i in range(n))
# using Girard’s theorem: A = [Σθᵢ − (n−2)·π]·R²
# (PI2 - abs(s) == (n*PI - abs(s)) - (n-2)*PI)
return abs(PI2 - abs(s)) * Radius(radius)**2
def intersection(start1, end1, start2, end2,
height=None, LatLon=LatLon, **LatLon_kwds):
'''Locate the intersection of two paths each defined by two
points or by a start point and an initial bearing.
@arg start1: Start point of the first path (L{LatLon}).
@arg end1: End point of the first path (L{LatLon}) or the
initial bearing at the first start point
(compass C{degrees360}).
@arg start2: Start point of the second path (L{LatLon}).
@arg end2: End point of the second path (L{LatLon}) or the
initial bearing at the second start point
(compass C{degrees360}).
@kwarg height: Optional height at the intersection point,
overriding the mean height (C{meter}).
@kwarg LatLon: Optional class to return the intersection
point (L{LatLon}).
@kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
arguments, ignored if B{C{LatLon=None}}.
@return: The intersection point (B{C{LatLon}}) or 3-tuple
(C{degrees90}, C{degrees180}, height) if B{C{LatLon}}
is C{None} or C{None} if no unique intersection
exists.
@raise TypeError: If B{C{start*}} or B{C{end*}} is not L{LatLon}.
@raise ValueError: Intersection is ambiguous or infinite or
the paths are parallel, coincident or null.
@example:
>>> p = LatLon(51.8853, 0.2545)
>>> q = LatLon(49.0034, 2.5735)
>>> i = intersection(p, 108.55, q, 32.44) # 50.9076°N, 004.5086°E
'''
_Nvll.others(start1, name=_start1_)
_Nvll.others(start2, name=_start2_)
# If gc1 and gc2 are great circles through start and end points
# (or defined by start point and bearing), then the candidate
# intersections are simply gc1 × gc2 and gc2 × gc1. Most of the
# work is deciding the correct intersection point to select! If
# bearing is given, that determines the intersection, but if both
# paths are defined by start/end points, take closer intersection.
gc1, s1, e1 = _Nvll._gc3(start1, end1, 'end1')
gc2, s2, e2 = _Nvll._gc3(start2, end2, 'end2')
hs = start1.height, start2.height
# there are two (antipodal) candidate intersection
# points ... we have to choose the one to return
i1 = gc1.cross(gc2, raiser=_paths_)
# postpone computing i2 until needed
# i2 = gc2.cross(gc1, raiser=_paths_)
# selection of intersection point depends on how
# paths are defined (by bearings or endpoints)
if e1 and e2: # endpoint+endpoint
d = sumOf((s1, s2, e1, e2)).dot(i1)
hs += end1.height, end2.height
elif e1 and not e2: # endpoint+bearing
# gc2 x v2 . i1 +ve means v2 bearing points to i1
d = gc2.cross(s2).dot(i1)
hs += end1.height,
elif e2 and not e1: # bearing+endpoint