Skip to content

Commit

Permalink
bugfix: catch Numpy ArithmeticError
Browse files Browse the repository at this point in the history
fixes #78
  • Loading branch information
scivision committed Feb 26, 2023
1 parent 162d5db commit dde1b9c
Show file tree
Hide file tree
Showing 6 changed files with 39 additions and 15 deletions.
6 changes: 5 additions & 1 deletion src/pymap3d/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@
* Fortran: [Maptran3D](https://github.com/geospace-code/maptran3d)
"""

__version__ = "2.10.0"
__version__ = "3.0.0"

from .aer import aer2ecef, aer2geodetic, ecef2aer, geodetic2aer
from .ecef import (
Expand Down Expand Up @@ -94,6 +94,10 @@
]

try:
import numpy as np

np.seterr(all="raise")

This comment has been minimized.

Copy link
@davidrigie

davidrigie Mar 4, 2023

Contributor

This caused me a lot of grief


from .aer import aer2eci, eci2aer
from .azelradec import azel2radec, radec2azel
from .eci import ecef2eci, eci2ecef
Expand Down
14 changes: 8 additions & 6 deletions src/pymap3d/ecef.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
from math import pi

from .ellipsoid import Ellipsoid
from .mathfun import atan, atan2, cos, degrees, hypot, radians, sin, sqrt, tan
from .mathfun import atan, atan2, cos, degrees, hypot, isclose, radians, sin, sqrt, tan
from .utils import sanitize

__all__ = [
Expand Down Expand Up @@ -70,8 +70,8 @@ def geodetic2ecef(
lon = radians(lon)

# radius of curvature of the prime vertical section
N = ell.semimajor_axis**2 / sqrt(
ell.semimajor_axis**2 * cos(lat) ** 2 + ell.semiminor_axis**2 * sin(lat) ** 2
N = ell.semimajor_axis**2 / hypot(
ell.semimajor_axis * cos(lat), ell.semiminor_axis * sin(lat)
)
# Compute cartesian (geocentric) coordinates given (curvilinear) geodetic coordinates.
x = (N + alt) * cos(lat) * cos(lon)
Expand Down Expand Up @@ -133,7 +133,7 @@ def ecef2geodetic(
E = sqrt(ell.semimajor_axis**2 - ell.semiminor_axis**2)

# eqn. 4a
u = sqrt(0.5 * (r**2 - E**2) + 0.5 * sqrt((r**2 - E**2) ** 2 + 4 * E**2 * z**2))
u = sqrt(0.5 * (r**2 - E**2) + 0.5 * hypot(r**2 - E**2, 2 * E * z))

Q = hypot(x, y)

Expand All @@ -142,8 +142,10 @@ def ecef2geodetic(
# eqn. 4b
try:
Beta = atan(huE / u * z / hypot(x, y))
except ZeroDivisionError:
if z >= 0:
except ArithmeticError:
if isclose(z, 0):
Beta = 0
elif z > 0:
Beta = pi / 2
else:
Beta = -pi / 2
Expand Down
2 changes: 1 addition & 1 deletion src/pymap3d/latitude.py
Original file line number Diff line number Diff line change
Expand Up @@ -344,7 +344,7 @@ def geodetic2conformal(geodetic_lat, ell: Ellipsoid = None, deg: bool = True):
# compute conformal latitudes with correction for points at +90
try:
conformal_lat = 2 * atan(sqrt((f4 / f3) * ((f1 / f2) ** e))) - (pi / 2)
except ZeroDivisionError:
except ArithmeticError:
conformal_lat = pi / 2

return degrees(conformal_lat) if deg else conformal_lat
Expand Down
14 changes: 10 additions & 4 deletions src/pymap3d/mathfun.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
exp,
hypot,
inf,
isclose,
isnan,
log,
power,
Expand All @@ -36,6 +37,7 @@
exp,
hypot,
inf,
isclose,
isnan,
log,
radians,
Expand All @@ -48,7 +50,7 @@ def power(x, y): # type: ignore
return pow(x, y)

def sign(x) -> float: # type: ignore
"""signum function"""
"""signum"""
if x < 0:
y = -1.0
elif x > 0:
Expand All @@ -58,9 +60,12 @@ def sign(x) -> float: # type: ignore

return y

def cbrt(x) -> float: # type: ignore
"""math.cbrt was added in Python 3.11"""
return x ** (1 / 3)
try:
import math.cbrt as cbrt # type: ignore
except ImportError:

def cbrt(x) -> float: # type: ignore
return x ** (1 / 3)


__all__ = [
Expand All @@ -75,6 +80,7 @@ def cbrt(x) -> float: # type: ignore
"exp",
"hypot",
"inf",
"isclose",
"isnan",
"log",
"power",
Expand Down
14 changes: 13 additions & 1 deletion src/pymap3d/tests/test_geodetic.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,9 +14,21 @@
B = ELL.semiminor_axis

xyzlla = [
((A, 0, 0), (0, 0, 0)),
((A - 1, 0, 0), (0, 0, -1)),
((A + 1, 0, 0), (0, 0, 1)),
((0.1 * A, 0, 0), (0, 0, -0.9 * A)),
((0.001 * A, 0, 0), (0, 0, -0.999 * A)),
((0, A, 0), (0, 90, 0)),
((0, A - 1, 0), (0, 90, -1)),
((0, A + 1, 0), (0, 90, 1)),
((0, 0.1 * A, 0), (0, 90, -0.9 * A)),
((0, 0.001 * A, 0), (0, 90, -0.999 * A)),
((0, 0, B), (90, 0, 0)),
((0, 0, B + 1), (90, 0, 1)),
((0, 0, B - 1), (90, 0, -1)),
((0, 0, 0.1 * B), (90, 0, -0.9 * B)),
((0, 0, 0.001 * B), (90, 0, -0.999 * B)),
((0, 0, B - 1), (89.999999, 0, -1)),
((0, 0, B - 1), (89.99999, 0, -1)),
((0, 0, -B + 1), (-90, 0, -1)),
Expand Down Expand Up @@ -179,7 +191,7 @@ def test_ecef2geodetic(xyz, lla):
lat, lon, alt = pm.ecef2geodetic(*xyz)
assert lat == approx(lla[0])
assert lon == approx(lla[1])
assert alt == approx(lla[2])
assert alt == approx(lla[2], abs=1e-9)


@pytest.mark.parametrize(
Expand Down
4 changes: 2 additions & 2 deletions src/pymap3d/vincenty.py
Original file line number Diff line number Diff line change
Expand Up @@ -190,10 +190,10 @@ def vdist(
alpha[isnan(sinAlpha)] = 0
alpha[(sinAlpha > 1) | (abs(sinAlpha - 1) < 1e-16)] = pi / 2

except (ZeroDivisionError, TypeError, ValueError):
except (ArithmeticError, TypeError, ValueError):
try:
sinAlpha = cos(U1) * cos(U2) * sin(lamb) / sin(sigma)
except ZeroDivisionError:
except ArithmeticError:
sinAlpha = 0.0

if isnan(sinAlpha):
Expand Down

0 comments on commit dde1b9c

Please sign in to comment.