Skip to content

Commit a22eeae

Browse files
scivisionscivision
committed
karney: valueerror if too few inputs, use math.
1 parent c96f551 commit a22eeae

2 files changed

Lines changed: 63 additions & 63 deletions

File tree

‎src/pymap3d/karney.py‎

Lines changed: 61 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@
1515

1616
from __future__ import annotations
1717

18-
import math as _math
18+
import math
1919

2020
from .ellipsoid import Ellipsoid, get_cached_geodesic
2121
from .mathfun import copysign, fmod, hypot, isnan, sqrt
@@ -84,8 +84,8 @@ def _sincosd(x):
8484
r = fmod(x, 360.0)
8585
q = 0 if isnan(r) else int(round(r / 90.0))
8686
r -= 90.0 * q
87-
r = _math.radians(r)
88-
s, c = _math.sin(r), _math.cos(r)
87+
r = math.radians(r)
88+
s, c = math.sin(r), math.cos(r)
8989
q = q % 4
9090
if q == 1:
9191
s, c = c, -s
@@ -108,7 +108,7 @@ def _atan2d(y, x):
108108
if x < 0:
109109
x = -x
110110
q += 1
111-
ang = _math.degrees(_math.atan2(y, x))
111+
ang = math.degrees(math.atan2(y, x))
112112
if q == 1:
113113
ang = copysign(180.0, y) - ang
114114
elif q == 2:
@@ -289,8 +289,8 @@ def _Astroid(x, y):
289289
T = T3 ** (1.0 / 3.0) if T3 >= 0 else -((-T3) ** (1.0 / 3.0))
290290
u += T + (r2 / T if T != 0 else 0.0)
291291
else:
292-
ang = _math.atan2(sqrt(-disc), -(S + r3))
293-
u += 2.0 * r * _math.cos(ang / 3.0)
292+
ang = math.atan2(sqrt(-disc), -(S + r3))
293+
u += 2.0 * r * math.cos(ang / 3.0)
294294
v = sqrt(u * u + q)
295295
uv = u + v if u >= 0 else q / (v - u)
296296
w = (uv - q) / (2.0 * v)
@@ -319,10 +319,10 @@ def __init__(self, ell: Ellipsoid | None = None):
319319
self._c2 = self.a ** 2
320320
elif self._e2 > 0:
321321
e = sqrt(self._e2)
322-
self._c2 = (self.a ** 2 + self._b ** 2 * _math.atanh(e) / e) / 2.0
322+
self._c2 = (self.a ** 2 + self._b ** 2 * math.atanh(e) / e) / 2.0
323323
else:
324324
e = sqrt(-self._e2)
325-
self._c2 = (self.a ** 2 + self._b ** 2 * _math.atan(e) / e) / 2.0
325+
self._c2 = (self.a ** 2 + self._b ** 2 * math.atan(e) / e) / 2.0
326326
self._etol2 = 0.1 * _tol2 / sqrt(
327327
max(0.001, abs(self.f)) * min(1.0, 1.0 - self.f / 2.0) / 2.0)
328328

@@ -389,7 +389,7 @@ def _C4f(self, eps, c):
389389

390390
def _Lengths(self, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
391391
cbet1, cbet2, outmask, C1a, C2a):
392-
s12b = m12b = m0 = M12 = M21 = _math.nan
392+
s12b = m12b = m0 = M12 = M21 = math.nan
393393
A1 = _A1m1f(eps)
394394
_C1f(eps, C1a)
395395
if outmask & (_REDUCEDLENGTH | _GEODESICSCALE):
@@ -430,7 +430,7 @@ def _InverseStart(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
430430
lam12, slam12, clam12, C1a, C2a):
431431
sig12 = -1.0
432432
salp1 = 2.0
433-
calp1 = salp2 = calp2 = dnm = _math.nan
433+
calp1 = salp2 = calp2 = dnm = math.nan
434434

435435
sbet12 = sbet2 * cbet1 - cbet2 * sbet1
436436
cbet12 = cbet2 * cbet1 + sbet2 * sbet1
@@ -442,8 +442,8 @@ def _InverseStart(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
442442
sbetm2 /= sbetm2 + (cbet1 + cbet2) ** 2
443443
dnm = sqrt(1.0 + self._ep2 * sbetm2)
444444
omg12 = lam12 / (self._f1 * dnm)
445-
somg12 = _math.sin(omg12)
446-
comg12 = _math.cos(omg12)
445+
somg12 = math.sin(omg12)
446+
comg12 = math.cos(omg12)
447447
else:
448448
somg12 = slam12
449449
comg12 = clam12
@@ -461,26 +461,26 @@ def _InverseStart(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
461461
# Short line case
462462
pass
463463

464-
if not (csig12 >= 0 or ssig12 >= 6 * abs(self.f) * _math.pi * cbet1 ** 2):
464+
if not (csig12 >= 0 or ssig12 >= 6 * abs(self.f) * math.pi * cbet1 ** 2):
465465
if self.f >= 0:
466466
k2 = sbet1 ** 2 * self._ep2
467467
eps = k2 / (2.0 * (1.0 + sqrt(1.0 + k2)) + k2)
468-
lamscale = self.f * cbet1 * self._A3f(eps) * _math.pi
468+
lamscale = self.f * cbet1 * self._A3f(eps) * math.pi
469469
betscale = lamscale * cbet1
470-
x = (lam12 - _math.pi) / lamscale
470+
x = (lam12 - math.pi) / lamscale
471471
y = sbet12a / betscale
472472
else:
473473
cbet12a = cbet2 * cbet1 - sbet2 * sbet1
474-
bet12a = _math.atan2(sbet12a, cbet12a)
474+
bet12a = math.atan2(sbet12a, cbet12a)
475475
_, m12b, m0, _, _ = self._Lengths(
476-
self._n, _math.pi + bet12a,
476+
self._n, math.pi + bet12a,
477477
sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
478478
cbet1, cbet2, _REDUCEDLENGTH, C1a, C2a)
479-
x = -1.0 + m12b / (cbet1 * cbet2 * m0 * _math.pi)
479+
x = -1.0 + m12b / (cbet1 * cbet2 * m0 * math.pi)
480480
betscale = (sbet12a / x if x < -0.01
481-
else -self.f * cbet1 ** 2 * _math.pi)
481+
else -self.f * cbet1 ** 2 * math.pi)
482482
lamscale = betscale / cbet1
483-
y = (lam12 - _math.pi) / lamscale
483+
y = (lam12 - math.pi) / lamscale
484484

485485
if y < -_tol1 and x > -1.0 - _xthresh:
486486
if self.f >= 0:
@@ -494,8 +494,8 @@ def _InverseStart(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
494494
omg12a = lamscale * (
495495
-x * k / (1.0 + k) if self.f >= 0
496496
else -y * (1.0 + k) / k)
497-
somg12 = _math.sin(omg12a)
498-
comg12 = -_math.cos(omg12a)
497+
somg12 = math.sin(omg12a)
498+
comg12 = -math.cos(omg12a)
499499
salp1 = cbet2 * somg12
500500
calp1 = sbet12a - cbet2 * sbet1 * somg12 ** 2 / (1.0 - comg12)
501501

@@ -532,13 +532,13 @@ def _Lambda12(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
532532
csig2 = comg2 = calp2 * cbet2
533533
ssig2, csig2 = _norm2(ssig2, csig2)
534534

535-
sig12 = _math.atan2(
535+
sig12 = math.atan2(
536536
max(0.0, csig1 * ssig2 - ssig1 * csig2),
537537
csig1 * csig2 + ssig1 * ssig2)
538538

539539
somg12 = max(0.0, comg1 * somg2 - somg1 * comg2)
540540
comg12 = comg1 * comg2 + somg1 * somg2
541-
eta = _math.atan2(
541+
eta = math.atan2(
542542
somg12 * clam120 - comg12 * slam120,
543543
comg12 * clam120 + somg12 * slam120)
544544

@@ -560,7 +560,7 @@ def _Lambda12(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
560560
cbet1, cbet2, _REDUCEDLENGTH, C1a, C2a)
561561
dlam12 *= self._f1 / (calp2 * cbet2)
562562
else:
563-
dlam12 = _math.nan
563+
dlam12 = math.nan
564564

565565
return (lam12, salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
566566
eps, somg12, comg12, dlam12)
@@ -570,7 +570,7 @@ def _Inverse(self, lat1, lon1, lat2, lon2):
570570
lonsign = 1.0 if lon12 >= 0 else -1.0
571571
lon12 = lonsign * _AngRound(lon12)
572572
lon12s = _AngRound((180.0 - lon12) - lonsign * lon12s)
573-
lam12 = _math.radians(lon12)
573+
lam12 = math.radians(lon12)
574574
if lon12 > 90:
575575
slam12, clam12 = _sincosd(lon12s)
576576
clam12 = -clam12
@@ -612,8 +612,8 @@ def _Inverse(self, lat1, lon1, lat2, lon2):
612612
C2a = [0.0] * (_nC2 + 1)
613613
C3a = [0.0] * _nC3
614614

615-
s12x = m12x = a12 = _math.nan
616-
M12 = M21 = _math.nan
615+
s12x = m12x = a12 = math.nan
616+
M12 = M21 = math.nan
617617

618618
if meridian:
619619
calp1 = clam12
@@ -626,7 +626,7 @@ def _Inverse(self, lat1, lon1, lat2, lon2):
626626
ssig2 = sbet2
627627
csig2 = calp2 * cbet2
628628

629-
sig12 = _math.atan2(
629+
sig12 = math.atan2(
630630
max(0.0, csig1 * ssig2 - ssig1 * csig2),
631631
csig1 * csig2 + ssig1 * ssig2)
632632

@@ -640,7 +640,7 @@ def _Inverse(self, lat1, lon1, lat2, lon2):
640640
sig12 = m12x = 0.0
641641
m12x *= self._b
642642
s12x *= self._b
643-
a12 = _math.degrees(sig12)
643+
a12 = math.degrees(sig12)
644644
else:
645645
meridian = False
646646

@@ -652,8 +652,8 @@ def _Inverse(self, lat1, lon1, lat2, lon2):
652652
salp1 = salp2 = 1.0
653653
s12x = self.a * lam12
654654
sig12 = lam12 / self._f1
655-
m12x = self._b * _math.sin(sig12)
656-
M12 = M21 = _math.cos(sig12)
655+
m12x = self._b * math.sin(sig12)
656+
M12 = M21 = math.cos(sig12)
657657
a12 = lon12 / self._f1
658658

659659
elif not meridian:
@@ -662,11 +662,11 @@ def _Inverse(self, lat1, lon1, lat2, lon2):
662662
lam12, slam12, clam12, C1a, C2a)
663663

664664
if sig12 >= 0:
665-
B1 = _SinCosSeries(True, _math.sin(sig12), _math.cos(sig12), C1a, _nC1)
665+
B1 = _SinCosSeries(True, math.sin(sig12), math.cos(sig12), C1a, _nC1)
666666
s12x = self._b * (1 + _A1m1f(self._n)) * (sig12 + B1)
667667
m12x = 0.0
668668
M12 = M21 = 1.0
669-
a12 = _math.degrees(sig12)
669+
a12 = math.degrees(sig12)
670670
else:
671671
tripb = False
672672
salp1a = _tiny
@@ -694,8 +694,8 @@ def _Inverse(self, lat1, lon1, lat2, lon2):
694694

695695
if numit < _maxit1 and dlam12 != 0:
696696
dalp1 = -v / dlam12
697-
sdalp1 = _math.sin(dalp1)
698-
cdalp1 = _math.cos(dalp1)
697+
sdalp1 = math.sin(dalp1)
698+
cdalp1 = math.cos(dalp1)
699699
nsalp1 = salp1 * cdalp1 + calp1 * sdalp1
700700
if nsalp1 > 0:
701701
calp1 = calp1 * cdalp1 - salp1 * sdalp1
@@ -716,7 +716,7 @@ def _Inverse(self, lat1, lon1, lat2, lon2):
716716
cbet1, cbet2, _OUT, C1a, C2a)
717717
m12x *= self._b
718718
s12x *= self._b
719-
a12 = _math.degrees(sig12)
719+
a12 = math.degrees(sig12)
720720

721721
if swapp < 0:
722722
salp1, salp2 = salp2, salp1
@@ -764,31 +764,31 @@ def _Direct(self, lat1, lon1, azi1, s12):
764764

765765
A1p1 = 1.0 + A1
766766
B11 = _SinCosSeries(True, ssig1, csig1, C1a, _nC1)
767-
s_B11 = _math.sin(B11)
768-
c_B11 = _math.cos(B11)
767+
s_B11 = math.sin(B11)
768+
c_B11 = math.cos(B11)
769769
stau1 = ssig1 * c_B11 + csig1 * s_B11
770770
ctau1 = csig1 * c_B11 - ssig1 * s_B11
771771

772772
tau12 = s12 / (self._b * A1p1)
773-
s_tau12 = _math.sin(tau12)
774-
c_tau12 = _math.cos(tau12)
773+
s_tau12 = math.sin(tau12)
774+
c_tau12 = math.cos(tau12)
775775

776776
stau2 = stau1 * c_tau12 + ctau1 * s_tau12
777777
ctau2 = ctau1 * c_tau12 - stau1 * s_tau12
778778

779779
B12 = -_SinCosSeries(True, stau2, ctau2, C1pa, _nC1p)
780780
sig12 = tau12 - (B12 - B11)
781-
ssig12 = _math.sin(sig12)
782-
csig12 = _math.cos(sig12)
781+
ssig12 = math.sin(sig12)
782+
csig12 = math.cos(sig12)
783783

784784
if abs(self.f) > 0.01:
785785
ssig2 = ssig1 * csig12 + csig1 * ssig12
786786
csig2 = csig1 * csig12 - ssig1 * ssig12
787787
B12 = _SinCosSeries(True, ssig2, csig2, C1a, _nC1)
788788
serr = (1.0 + A1) * (sig12 + (B12 - B11)) - s12 / self._b
789789
sig12 -= serr / sqrt(1.0 + k2 * ssig2 ** 2)
790-
ssig12 = _math.sin(sig12)
791-
csig12 = _math.cos(sig12)
790+
ssig12 = math.sin(sig12)
791+
csig12 = math.cos(sig12)
792792

793793
ssig2 = ssig1 * csig12 + csig1 * ssig12
794794
csig2 = csig1 * csig12 - ssig1 * ssig12
@@ -814,11 +814,11 @@ def _Direct(self, lat1, lon1, azi1, s12):
814814
# Use LONG_UNROLL method for robustness (handles equatorial case)
815815
E = copysign(1.0, salp0)
816816
omg12 = E * (sig12
817-
- (_math.atan2(ssig2, csig2) - _math.atan2(ssig1, csig1))
818-
+ (_math.atan2(E * somg2, comg2) - _math.atan2(E * somg1, comg1)))
817+
- (math.atan2(ssig2, csig2) - math.atan2(ssig1, csig1))
818+
+ (math.atan2(E * somg2, comg2) - math.atan2(E * somg1, comg1)))
819819
lam12 = omg12 - self.f * A3c * salp0 * (sig12 + (B32 - B31))
820820

821-
lon12 = _math.degrees(lam12)
821+
lon12 = math.degrees(lam12)
822822
lat2 = _atan2d(sbet2, self._f1 * cbet2)
823823
lon2 = _AngNormalize(lon1 + lon12)
824824
azi2 = _atan2d(salp2, calp2)
@@ -871,11 +871,11 @@ def geodesic_inverse(
871871
"""
872872
g = _get_geodesic(ell)
873873
if not deg:
874-
lat1, lon1 = _math.degrees(lat1), _math.degrees(lon1)
875-
lat2, lon2 = _math.degrees(lat2), _math.degrees(lon2)
874+
lat1, lon1 = math.degrees(lat1), math.degrees(lon1)
875+
lat2, lon2 = math.degrees(lat2), math.degrees(lon2)
876876
dist, azi1, azi2 = g._Inverse(lat1, lon1, lat2, lon2)
877877
if not deg:
878-
azi1, azi2 = _math.radians(azi1), _math.radians(azi2)
878+
azi1, azi2 = math.radians(azi1), math.radians(azi2)
879879
return dist, azi1, azi2
880880

881881

@@ -915,11 +915,11 @@ def geodesic_direct(
915915
"""
916916
g = _get_geodesic(ell)
917917
if not deg:
918-
lat1, lon1 = _math.degrees(lat1), _math.degrees(lon1)
919-
azi1 = _math.degrees(azi1)
918+
lat1, lon1 = math.degrees(lat1), math.degrees(lon1)
919+
azi1 = math.degrees(azi1)
920920
lat2, lon2, azi2 = g._Direct(lat1, lon1, azi1, dist)
921921
if not deg:
922-
lat2, lon2, azi2 = _math.radians(lat2), _math.radians(lon2), _math.radians(azi2)
922+
lat2, lon2, azi2 = math.radians(lat2), math.radians(lon2), math.radians(azi2)
923923
return lat2, lon2, azi2
924924

925925

@@ -966,15 +966,16 @@ def geodesic_area(
966966
) -> tuple:
967967
"""
968968
Compute the area and perimeter of a geodesic polygon.
969+
Polygon must have at least 3 vertices.
969970
970971
Counterclockwise vertex ordering gives positive area.
971972
972973
Parameters
973974
----------
974975
lats : array-like
975-
Latitudes of polygon vertices.
976+
Latitudes of polygon vertices: 1-D array or list.
976977
lons : array-like
977-
Longitudes of polygon vertices.
978+
Longitudes of polygon vertices: 1-D array or list.
978979
ell : Ellipsoid, optional
979980
Reference ellipsoid (default WGS84).
980981
deg : bool, optional
@@ -987,9 +988,10 @@ def geodesic_area(
987988
perimeter : float
988989
Perimeter in meters.
989990
"""
991+
990992
n = len(lats)
991993
if n < 3:
992-
return 0.0, 0.0
994+
raise ValueError("Polygon must have at least 3 vertices")
993995

994996
g = _get_geodesic(ell)
995997
perimeter = 0.0
@@ -1012,8 +1014,8 @@ def geodesic_area(
10121014
diff, _ = _AngDiff(azi2_prev, azi1_first)
10131015
area_sum += diff
10141016

1015-
area = g._c2 * _math.radians(area_sum)
1016-
area0 = 4.0 * _math.pi * g._c2
1017+
area = g._c2 * math.radians(area_sum)
1018+
area0 = 4.0 * math.pi * g._c2
10171019
if abs(area) > area0 / 2:
10181020
area -= copysign(area0, area)
10191021
return area, perimeter

‎src/pymap3d/tests/test_karney.py‎

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -115,10 +115,8 @@ def test_triangle(self):
115115
assert perim > 0
116116

117117
def test_too_few_points(self):
118-
area, perim = karney.geodesic_area([0, 1], [0, 1])
119-
assert area == 0.0
120-
assert perim == 0.0
121-
118+
with pytest.raises(ValueError):
119+
area, perim = karney.geodesic_area([0, 1], [0, 1])
122120

123121
class TestGeodesicCache:
124122
def test_default_ellipsoid_cache_reuse(self):

0 commit comments

Comments
 (0)