geospace-code / pymap3d Goto Github PK
View Code? Open in Web Editor NEWpure-Python (Numpy optional) 3D coordinate conversions for geospace ecef enu eci
Home Page: https://geospace-code.github.io/pymap3d
License: BSD 2-Clause "Simplified" License
pure-Python (Numpy optional) 3D coordinate conversions for geospace ecef enu eci
Home Page: https://geospace-code.github.io/pymap3d
License: BSD 2-Clause "Simplified" License
Consider the following example:
Assuming I'm starting at the point at the bottom left, I want to work out the point in the top right knowing the distance travelled and the azimuth.
To do, I use the following code:
import pymap3d.vincenty
lat2, lon2 = pymap3d.vincenty.vreckon(52.22610277777778, -1.2696583333333333, 839.63, 63.02)
print(lat2, lon2)
The above code returns:
52.22952562862266 358.7412927899441
Note the latitude is correct but the longitude is totally invalid. This appears to be caused by using negative longitudes. Negative latitudes seem to be correctly handled.
No response
Hello.
Can I get a reference to the rsphere.euler function?
I want to understand formulas and principles.
I hope my question wasn't rude.
Thank you.
Hi Team,
This might not be the place for this question, but I'm using pymap3d for this and I'm hoping someone can provide advice - I'm running a simulation that begins at t=0
where the assumption is that at this point ECEF == ECI. As time in the simulation progresses, these two frames naturally deviate. The problem that I'm having is that when using eci2ecef
, I'm not sure what datetime
value to pass as an argument to replicate this assumption.
The docs indicate that the default assumption is the J2000 frame, which presumably implies t=0
at 2000-01-01
when the earths equatorial plane is coincident with the plane defined by its elliptical orbit around the sun - the vernal equinox. This point is obviously not when ECEF == ECI.
So far I've just attempted to search the solution space iterating over a 24 hour period to identify when norm((R,0,0)_ECI) == norm((R,0,0)_ECEF), which has give me a rough-order-of-magnitude response, but it seems as though eci2ecef
only provides a maximum accuracy of 1 second (though it appears it forwards this to astropy, which I assume is capable of providing ~nanosecond accuracy).
I'm wondering if anybody can provide any advice on better ways to identify the time where ECI==ECEF? Is there anything built-in that I'm missing here?
The latest version of pip, and possible several versions prior, raises an error when installing this package due to the discrepancy between the release version in the filename, and the version in the metadata pyproject.toml
.
Is there any release planned to resolve this issue? For now, having to pin to v2.6.1 which isn't ideal.
Thanks
Thanks for porting over these MatLab codes! I was cross checking the eci2aer function with AstroPy before incorperating it into another project I am working on, but I must be missing something or these numbers are very far off...
from pymap3d.aer import eci2aer as eci2aer
from astropy.coordinates import SkyCoord, GCRS, EarthLocation, CartesianRepresentation, AltAz
from astropy import units as u
from datetime import datetime
import numpy as np
eci = np.asarray([-34320604.87943786,2010725.129018341, 2977059.1586334566]) # units in meters
lla = np.asarray([10.0,10.0, 100]) # units in degrees, degrees, meters
t = datetime(2020, 3, 15, 0, 0, 0) # time as a datetime object
sat = SkyCoord(x=eci[0]*u.m, y=eci[1]*u.m, z=eci[2]*u.m, frame=GCRS,
representation_type=CartesianRepresentation)
origin = EarthLocation(lat=lla[0]*u.deg, lon=lla[1]*u.deg, height=lla[2]*u.m)
sataltaz = sat.transform_to(AltAz(obstime=t,location=origin))
az, el, srange = eci2aer(eci[0], eci[1], eci[2], lla[0], lla[1], lla[2],t)
np.array(sataltaz.az)-az
#array([70.90364332])
np.array(sataltaz.alt)-el
#array([-120.90947987])
np.array(sataltaz.distance)-srange
#array([1.77086727e+11])
Is ECI ~ GCRS ? Any feedback is appreciated. If this is my error vice the projects, I will gladly delete this.
When trying to install pymap3d on a fresh Ubuntu 16.04 installation using pip3, the following error is returned:
$ pip3 install pymap3d
Collecting pymap3d
Installing collected packages: pymap3d
Exception:
Traceback (most recent call last):
File "/usr/lib/python3/dist-packages/pip/basecommand.py", line 209, in main
status = self.run(options, args)
File "/usr/lib/python3/dist-packages/pip/commands/install.py", line 335, in run
prefix=options.prefix_path,
File "/usr/lib/python3/dist-packages/pip/req/req_set.py", line 732, in install
**kwargs
File "/usr/lib/python3/dist-packages/pip/req/req_install.py", line 837, in install
self.move_wheel_files(self.source_dir, root=root, prefix=prefix)
File "/usr/lib/python3/dist-packages/pip/req/req_install.py", line 1039, in move_wheel_files
isolated=self.isolated,
File "/usr/lib/python3/dist-packages/pip/wheel.py", line 346, in move_wheel_files
assert info_dir, "%s .dist-info directory not found" % req
AssertionError: pymap3d .dist-info directory not found
I also got the same error on another system using Mint 17.3.
Any ideas ?
Cheers,
Hey @scivision
We are using pymap3d in a GPS-IMU sensor fusion project and noticed that the latest version crashes with Python 2.7
In [1]: import pymap3d
File "/usr/local/lib/python2.7/dist-packages/pymap3d/__init__.py", line 122
ecef[i, :] = _rottrip(gst[i]) @ eci[i, :]
^
SyntaxError: invalid syntax
Best wishes and thanks much for open sourcing the code,
auerl
I followed the link to the API docs from the github project page. The documentation is out of date and still applies to a pre-3.0.0 version. For example, the new code and interface for Ellipsoid is here and has parameters like this:
Parameters
----------
semimajor_axis : float
semimajor axis in meters
semiminor_axis : float
semiminor axis in meters
name: str, optional
Human-friendly name for the ellipsoid
model: str, optional
Short name for the ellipsoid
"""
but the official docs on https://geospace-code.github.io/pymap3d/ellipsoid.html have the old interface, with only a string parameter:
class Ellipsoid(model='wgs84')
and
Parameters
----------
model : str
name of ellipsoid
Can you update the official docs for the new interface (3.0.0) please?
No response
Invoking both geodetic2enu and ecef2enu with the same location and observable produce different output
import pymap3d as pm
lat1 = 45.977
lon1 = 7.658
h1 = 4531
lat0 = 46.017
lon0 = 7.750
h0 = 1673
x,y,z = pm.geodetic2ecef(lat1,lon1,h1)
e,n,u = pm.ecef2enu(x,y,z,lat0,lon0,h0)
e1,n1,u1 = pm.geodetic2enu(lat1,lon1,h1,lat0,lon0,h0)
print ("E=%f N=%f U=%f" % (e,n,u))
print ("E1=%f N1=%f U1=%f" % (e1,n1,u1))
Produce different output. I am expecting same result for enu the U output differs.
If I use the ecef2geodetic function with x >= 521860, y = 0, and z = 0, I get correct outputs of lat = 0, lon =0, and the appropriate alt.
However, if I use x <= 521850, y = 0, z = 0, I get a lat of 90 (incorrect), lon of 0, and an incorrect altitude.
The difference seems to occur in the calculation of u on line 136 in ecef.py, which is:
u = sqrt(0.5 * (r**2 - E**2) + 0.5 * sqrt((r**2 - E**2) ** 2 + 4 * E**2 * z**2))
If r < E, then u = 0.
This can be reproduced by running the following:
import pymap3d
from numpy import linspace
Xs = linspace(521000, 522000, 101)
for x in Xs:
lat, lon, alt = pymap3d.ecef2geodetic(x, 0, 0)
print(x)
print(lat, lon, alt)
print(" ")
You will see that the discontinuity occurs from x=521850 and x=521860
I don't think this should happen, but I'm not sure what the problem is. The equations seem to match the referenced paper.
Any ideas? Thanks!
No response
My following tests both passed in v2.7.0, but now the first one fails in v2.7.1.
import numpy as np
import pymap3d.lox as lox
def test_lox():
lat_ref = 36.0
lon_ref = 140.0
lat = np.arange(30, 40).reshape(2, 5)
lon = np.arange(135, 145).reshape(2, 5)
len_, dir_ = lox.loxodrome_inverse(lat_ref, lon_ref, lat, lon)
def test_lox2():
lat_ref = 36.0
lon_ref = 140.0
lat = np.arange(30, 40).reshape(2, 5)
lon = np.arange(135, 145).reshape(2, 5)
len_, dir_ = lox.loxodrome_inverse(lat, lon, lat_ref, lon_ref)
The failure is like this:
====================================================================== FAILURES =======================================================================
______________________________________________________________________ test_lox _______________________________________________________________________
def test_lox():
lat_ref = 36.0
lon_ref = 140.0
lat = np.arange(30, 40).reshape(2, 5)
lon = np.arange(135, 145).reshape(2, 5)
> len_, dir_ = lox.loxodrome_inverse(lat_ref, lon_ref, lat, lon)
tests/test_pymap3d.py:10:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
../../../../../Library/Caches/pypoetry/virtualenvs/foo-TIPLNhpH-py3.9/lib/python3.9/site-packages/pymap3d/lox.py:173: in loxodrome_inverse
lat2 = broadcast_to(lat2, lat1.shape)
<__array_function__ internals>:5: in broadcast_to
???
../../../../../Library/Caches/pypoetry/virtualenvs/foo-TIPLNhpH-py3.9/lib/python3.9/site-packages/numpy/lib/stride_tricks.py:411: in broadcast_to
return _broadcast_to(array, shape, subok=subok, readonly=True)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
array = array([[0.52359878, 0.54105207, 0.55850536, 0.57595865, 0.59341195],
[0.61086524, 0.62831853, 0.64577182, 0.66322512, 0.68067841]])
shape = (1,), subok = False, readonly = True
def _broadcast_to(array, shape, subok, readonly):
shape = tuple(shape) if np.iterable(shape) else (shape,)
array = np.array(array, copy=False, subok=subok)
if not shape and array.shape:
raise ValueError('cannot broadcast a non-scalar to a scalar array')
if any(size < 0 for size in shape):
raise ValueError('all elements of broadcast shape must be non-'
'negative')
extras = []
> it = np.nditer(
(array,), flags=['multi_index', 'refs_ok', 'zerosize_ok'] + extras,
op_flags=['readonly'], itershape=shape, order='C')
E ValueError: input operand has more dimensions than allowed by the axis remapping
../../../../../Library/Caches/pypoetry/virtualenvs/foo-TIPLNhpH-py3.9/lib/python3.9/site-packages/numpy/lib/stride_tricks.py:348: ValueError
I have a fix, so I will post a PR.
pymap3d.geodetic2enu
has the following type signature:
def geodetic2enu(
lat: ndarray,
lon: ndarray,
h: ndarray,
lat0: ndarray,
lon0: ndarray,
h0: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
) -> tuple[ndarray, ndarray, ndarray]:
...
I'm pretty sure lat0
, lon0
, and h0
must always be scalars. If I provide scalar float
or int
to these args I get type checking errors. Furthermore, lat
, lon
, and h
may be either scalar or n-dimensional. If all these args are scalars, the function returns a tuple of three float
s, not three ndarray
s.
Admittedly numpy and type checking are a tricky business, so it may be better to only annotate non-numeric stuff. Alternatively, consider looking into numpy.typing
and ArrayLike
.
I think this probably applies to quite a number of other conversion functions, but I haven't looked into them.
When azimuth is equal or close to 90 or 270, loxodrome_direct
will return an incorrect longitude value.
The following code will reproduce the problem:
from pymap3d.lox import loxodrome_direct
clat, clon = 35.0, 140.0
az = [
89,
89.9,
89.99,
89.999,
89.9999,
89.99999,
89.999999,
90,
90.000001,
90.00001,
90.0001,
90.001,
90.01,
90.1,
91,
]
tmp = [az_ + 180 for az_ in az]
az += tmp
for az_ in az:
lat, lon = loxodrome_direct(clat, clon, 50_000, az_)
print(f"{az_:.6f}, {lat:.14f}, {lon:.14f}")
This code prints lines like this:
89.000000, 35.00786565022930, 140.54765888297607
89.900000, 35.00078660501723, 140.54771788308574
89.990000, 35.00007866054496, 140.54771634402164
89.999000, 35.00000786605369, 140.54771605442008
89.999900, 35.00000078660448, 140.54771541554422
89.999990, 35.00000007865957, 140.54770927042696
89.999999, 35.00000000786506, 140.54764724167964
90.000000, 34.99999999999901, -19494.22711642675131
90.000001, 34.99999999213296, 140.54778537068364
90.000010, 34.99999992133847, 140.54772297468475
90.000100, 34.99999921339354, 140.54771678232532
90.001000, 34.99999213394431, 140.54771614079507
90.010000, 34.99992133945204, 140.54771583369782
90.100000, 34.99921339487867, 140.54771264295587
91.000000, 34.99213433955717, 140.54760647858132
269.000000, 34.99213433955717, 139.45239352141868
269.900000, 34.99921339487867, 139.45228735704418
269.990000, 34.99992133945204, 139.45228416630223
269.999000, 34.99999213394431, 139.45228385919867
269.999900, 34.99999921339354, 139.45228321768184
269.999990, 34.99999992133847, 139.45227702538708
269.999999, 34.99999999213296, 139.45221462306554
270.000000, 34.99999999999901, -6404.74237214225013
270.000001, 35.00000000786506, 139.45235274366772
270.000010, 35.00000007865956, 139.45229076455408
270.000100, 35.00000078660448, 139.45228458430924
270.001000, 35.00000786605369, 139.45228394556526
270.010000, 35.00007866054496, 139.45228365597760
270.100000, 35.00078660501723, 139.45228211691446
271.000000, 35.00786565022930, 139.45234111702391
At the very least, the values of longitude when the azimuth approaching 90 or 270 to less than 1e-5 look obviously wrong. On the other hand, the latitude values look fine, including when the azimuth is close to 0 or 180.
pymap3d version: 2.7.0
Many thanks.
Hi There,
May I know which ECI frame is currently supported now, is it J2000 ? Is there a way (or any plans) to do this for other ECI frames as mentioned on the ECI wikipedia page ?
Thanks
Shashank
Hi, I have recently released a version of this library written in rust. The attribution link to this GitHub repository is available on the readme
If you are interested in linking also the rust implementation after the matlab and fortran ones, feel free.
Action to publish release v3.0.1 to pypi failed. I think it's because this version string wasn't updated
pymap3d/src/pymap3d/__init__.py
Line 32 in 5802ca3
loxodrome_inverse()
returns (circumference / 2 - expected distance) when lat1 == lat2 and lon1< lon2.
It is clear from the following calculation results:
>>> from pymap3d.lox import loxodrome_inverse
>>> loxodrome_inverse(0,0, 0, 1.0, deg=True)
(19926188.85199597, 90.0)
>>> loxodrome_inverse(0,0, 0, 10.0, deg=True)
(18924313.434856508, 90.0)
>>> loxodrome_inverse(0,0, 0, 170.0, deg=True)
(1113194.9079327343, 90.0)
>>> loxodrome_inverse(0,0, 0, 180.0, deg=True)
(0.0, 90.0)
>>> loxodrome_inverse(0, 1.0, 0, 0, deg=True)
(111319.49079327357, 270.0)
>>> loxodrome_inverse(0, 170.0, 0, 0, deg=True)
(18924313.434856508, 270.0)
At first I thought it was an instability issue around the singularity, but it appears to be a logic error in the singularity-specific processing.
I have made a correction and will send you a PR after this.
When checking a distance along the equator, vdist will return an incorrect result. By adjusting the latitudes so they are close to zero, vdist will correct this behavior.
The following code illustrates the issue:
from pymap3d.vincenty import vdist
print(0, '\t', vdist(0, 0, 0, 1)[0])
for i in range(16, -1, -1):
try:
print(1*10**-i, '\t', vdist(1*10**-i, 0, 1*10**-i, 1)[0])
except ValueError as e:
print(1*10**-i, '\t', f'Failed: {e.args}')
This code prints the following:
0 110946.29557670657
1e-16 110946.29557670657
1e-15 110946.29557670657
1e-14 110946.29557670657
1e-13 110946.29557670657
1e-12 110946.29557670657
1e-11 110946.29557670657
1e-10 110946.29557670657
1e-09 110946.29557670657
1e-08 110946.29557670657
1e-07 110946.29557670657
1e-06 Failed: ('math domain error',)
1e-05 111319.4907932247
0.0001 111319.49079305801
0.001 111319.49077638453
0.01 111319.48910904085
0.1 111319.32237470953
1 111302.64933938583
The values for 1e-05 and greater appear to be correct, and agree with the output from the Matlab distance function. The values for 0 and from 1e-16 through 1e-07 are all off by about the same amount, and 1e-06 produces a math domain error.
Part of openjournals/joss-reviews#580
The tag for the latest version is v.1.6.0
instead of v1.6.0
.
Recently e have updated the Debian package for pymap3d for v2.9.1 to v3.0 and we have started having some test failues on one of our ARM architectures
See the detailed log at:
https://ci.debian.net/data/autopkgtest/testing/armel/p/pymap3d/34582756/log.gz
============================= test session starts ==============================
platform linux -- Python 3.11.4, pytest-7.2.1, pluggy-1.0.0+repack
rootdir: /tmp/autopkgtest-lxc.5m0tpwuk/downtmp
plugins: astropy-header-0.2.2, mock-3.8.2, openfiles-0.5.0, arraydiff-0.5.0, cov-4.0.0, remotedata-0.4.0, astropy-0.10.0, doctestplus-0.12.1, filter-subpackage-0.1.2, hypothesis-6.67.1
collected 321 items
../build.4yE/src/src/pymap3d/tests/test_aer.py .......... [ 3%]
../build.4yE/src/src/pymap3d/tests/test_eci.py ..... [ 4%]
../build.4yE/src/src/pymap3d/tests/test_ellipsoid.py ................... [ 10%]
................. [ 15%]
../build.4yE/src/src/pymap3d/tests/test_enu.py ..... [ 17%]
../build.4yE/src/src/pymap3d/tests/test_geodetic.py .................... [ 23%]
F....F....................F....F................. [ 38%]
../build.4yE/src/src/pymap3d/tests/test_latitude.py .................... [ 45%]
......................... [ 52%]
../build.4yE/src/src/pymap3d/tests/test_look_spheroid.py ......... [ 55%]
../build.4yE/src/src/pymap3d/tests/test_ned.py ... [ 56%]
../build.4yE/src/src/pymap3d/tests/test_pyproj.py .. [ 57%]
../build.4yE/src/src/pymap3d/tests/test_rcurve.py ............. [ 61%]
../build.4yE/src/src/pymap3d/tests/test_rhumb.py ....................... [ 68%]
.......... [ 71%]
../build.4yE/src/src/pymap3d/tests/test_rsphere.py .......... [ 74%]
../build.4yE/src/src/pymap3d/tests/test_sidereal.py .... [ 76%]
../build.4yE/src/src/pymap3d/tests/test_sky.py .... [ 77%]
../build.4yE/src/src/pymap3d/tests/test_spherical.py ................... [ 83%]
....................... [ 90%]
../build.4yE/src/src/pymap3d/tests/test_time.py ..... [ 91%]
../build.4yE/src/src/pymap3d/tests/test_vincenty.py . [ 92%]
../build.4yE/src/src/pymap3d/tests/test_vincenty_dist.py ............... [ 96%]
. [ 97%]
../build.4yE/src/src/pymap3d/tests/test_vincenty_vreckon.py ......... [100%]
=================================== FAILURES ===================================
________________________ test_ecef2geodetic[xyz4-lla4] _________________________
xyz = (6378.137, 0, 0), lla = (0, 0, -6371758.863)
@pytest.mark.parametrize("xyz, lla", xyzlla)
def test_ecef2geodetic(xyz, lla):
lat, lon, alt = pm.ecef2geodetic(*xyz)
> assert lat == approx(lla[0])
E assert nan == 0 ± 1.0e-12
E comparison failed
E Obtained: nan
E Expected: 0 ± 1.0e-12
../build.4yE/src/src/pymap3d/tests/test_geodetic.py:192: AssertionError
________________________ test_ecef2geodetic[xyz9-lla9] _________________________
xyz = (0, 6378.137, 0), lla = (0, 90, -6371758.863)
@pytest.mark.parametrize("xyz, lla", xyzlla)
def test_ecef2geodetic(xyz, lla):
lat, lon, alt = pm.ecef2geodetic(*xyz)
> assert lat == approx(lla[0])
E assert nan == 0 ± 1.0e-12
E comparison failed
E Obtained: nan
E Expected: 0 ± 1.0e-12
../build.4yE/src/src/pymap3d/tests/test_geodetic.py:192: AssertionError
_____________________ test_numpy_ecef2geodetic[xyz4-lla4] ______________________
xyz = (6378.137, 0, 0), lla = (0, 0, -6371758.863)
@pytest.mark.parametrize("xyz, lla", xyzlla)
def test_numpy_ecef2geodetic(xyz, lla):
np = pytest.importorskip("numpy")
lat, lon, alt = pm.ecef2geodetic(
*np.array(
[
[xyz],
],
dtype=np.float32,
).T
)
> assert lat[0] == approx(lla[0])
E assert array([nan], dtype=float32) == 0 ± 1.0e-12
E comparison failed
E Obtained: [nan]
E Expected: 0 ± 1.0e-12
../build.4yE/src/src/pymap3d/tests/test_geodetic.py:263: AssertionError
_____________________ test_numpy_ecef2geodetic[xyz9-lla9] ______________________
xyz = (0, 6378.137, 0), lla = (0, 90, -6371758.863)
@pytest.mark.parametrize("xyz, lla", xyzlla)
def test_numpy_ecef2geodetic(xyz, lla):
np = pytest.importorskip("numpy")
lat, lon, alt = pm.ecef2geodetic(
*np.array(
[
[xyz],
],
dtype=np.float32,
).T
)
> assert lat[0] == approx(lla[0])
E assert array([nan], dtype=float32) == 0 ± 1.0e-12
E comparison failed
E Obtained: [nan]
E Expected: 0 ± 1.0e-12
../build.4yE/src/src/pymap3d/tests/test_geodetic.py:263: AssertionError
=============================== warnings summary ===============================
build.4yE/src/src/pymap3d/tests/test_pyproj.py::test_compare_geodetic
/tmp/autopkgtest-lxc.5m0tpwuk/downtmp/build.4yE/src/src/pymap3d/tests/test_pyproj.py:31: FutureWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
assert pyproj.transform(lla, ecef, lla0[1], lla0[0], lla0[2]) == approx(xyz)
build.4yE/src/src/pymap3d/tests/test_pyproj.py::test_compare_geodetic
/tmp/autopkgtest-lxc.5m0tpwuk/downtmp/build.4yE/src/src/pymap3d/tests/test_pyproj.py:32: FutureWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
assert pyproj.transform(ecef, lla, *xyz) == approx((lla0[1], lla0[0], lla0[2]))
-- Docs: https://docs.pytest.org/en/stable/how-to/capture-warnings.html
=========================== short test summary info ============================
FAILED ../build.4yE/src/src/pymap3d/tests/test_geodetic.py::test_ecef2geodetic[xyz4-lla4]
FAILED ../build.4yE/src/src/pymap3d/tests/test_geodetic.py::test_ecef2geodetic[xyz9-lla9]
FAILED ../build.4yE/src/src/pymap3d/tests/test_geodetic.py::test_numpy_ecef2geodetic[xyz4-lla4]
FAILED ../build.4yE/src/src/pymap3d/tests/test_geodetic.py::test_numpy_ecef2geodetic[xyz9-lla9]
================== 4 failed, 317 passed, 2 warnings in 29.81s =================
Thank you for fixing #42, the lox stability issue, and releasing the fix as v2.7.1.
However, loxodrome_direct
still returns incorrect longitude values, although is is improved before fixing.
The code is the same as one used in #42.
from pymap3d.lox import loxodrome_direct
clat, clon = 35.0, 140.0
az = [
89,
89.9,
89.99,
89.999,
89.9999,
89.99999,
89.999999,
90,
90.000001,
90.00001,
90.0001,
90.001,
90.01,
90.1,
91,
]
tmp = [az_ + 180 for az_ in az]
az += tmp
for az_ in az:
lat, lon = loxodrome_direct(clat, clon, 50_000, az_)
print(f"{az_:.6f}, {lat:.14f}, {lon:.14f}")
This code prints lines like this in PyMap 2.7.1:
89.000000, 35.00786565022930, 140.54765888297607
89.900000, 35.00078660501723, 140.54771788308574
89.990000, 35.00007866054496, 140.54771634402164
89.999000, 35.00000786605369, 140.54771605442008
89.999900, 35.00000078660448, 140.54771541554422
89.999990, 35.00000007865957, 140.54770927042696
89.999999, 35.00000000786506, 140.54764724167964
90.000000, 34.99999999999901, 139.55133723928088
90.000001, 34.99999999213296, 140.54778537068364
90.000010, 34.99999992133847, 140.54772297468475
90.000100, 34.99999921339354, 140.54771678232532
90.001000, 34.99999213394431, 140.54771614079507
90.010000, 34.99992133945204, 140.54771583369782
90.100000, 34.99921339487867, 140.54771264295587
91.000000, 34.99213433955717, 140.54760647858132
269.000000, 34.99213433955717, 139.45239352141868
269.900000, 34.99921339487867, 139.45228735704418
269.990000, 34.99992133945204, 139.45228416630223
269.999000, 34.99999213394431, 139.45228385919867
269.999900, 34.99999921339354, 139.45228321768184
269.999990, 34.99999992133847, 139.45227702538708
269.999999, 34.99999999213296, 139.45221462306554
270.000000, 34.99999999999901, 139.55133723928088
270.000001, 35.00000000786506, 139.45235274366772
270.000010, 35.00000007865956, 139.45229076455408
270.000100, 35.00000078660448, 139.45228458430924
270.001000, 35.00000786605369, 139.45228394556526
270.010000, 35.00007866054496, 139.45228365597760
270.100000, 35.00078660501723, 139.45228211691446
271.000000, 35.00786565022930, 139.45234111702391
Now loxodrome_direct
returns the same value 139.55133723928088 for both 90 and 270, and both are incorrect.
pymap3d version: 2.7.1
Many thanks.
Hi Michael,
I was looking through all of your contact information on scivision.dev, and couldn't find any email address.
I'm just here to report that your article at https://www.scivision.dev/numpy-image-bgr-to-rgb/ is very harmful. There are plenty of people taking that advice and linking to that page, which is why I have to reach you. I used to believe that the same technique was good, but it turns out that there's much more to BGR/RGB conversion than just modifying the "Numpy view" of the RAM storage. That's a "fake" conversion and doesn't actually work in real life.
Here's a full breakdown which explains what happens at a code-level in OpenCV when you give it such "faked" Numpy data:
The same is almost certainly true for matplotlib too.
With best regards,
Johnny
There is a function for converting latitudes from geodetic to geocentric but this is restricted to the surface of the ellipsoid. For many geophysical applications, we need to convert the (latitude, height)
of points into (latitude_spherical, radius)
.
This is implemented in Boule as boule.Ellipsoid.spherical_to_geodetic
. Boule is more focused on gravity calculations than coordinate conversions so I'd be happy to move that implementation to pymap3d if it's of interest to you.
I can volunteer to do this given a bit of guidance as to where to put the functions and what they should be called to fit with pymap3d.
When using columns from a pandas dataframe for the lat0, lon0, h0, az, tilt, I get:
ValueError: The truth value of a Series is ambiguous. Use a.empty, a.bool(), a.item(), a.any() or a.all().
The offending error seems to be caused by:
if h0 < 0:
raise ValueError("Intersection calculation requires altitude [0, Infinity)")
The height values are greater than zero, but that's not the problem exactly. It's testing whether the Series is less than zero, causing the error.
I have numpy (of course, since I'm using pandas) but this doesn't seem to vectorize correctly. Oddly enough, if I use np.vectorize (lookAtSpheroid) and pass the exact same columns from the dataframe it works without error. Something is lost in the wrapper function lookAtSpheroid for lookAtSpheroid_point.
When I call the enu2geodetic method I would expect that I get returned a tuple[float, float, float] instance. However, when I inspect the result what I get is a tuple of numpy.ndarrays. Is this intended behavior? In the documentation it specifies that the result should be a tuple of floats not arrays.
The current homepage https://pypi.org/project/pymap3d/
is down
Hi @scivision, I've been developing a package for gravity and magnetic geophysics called Harmonica with @santisoler. We've been having the need for some coordinate conversions not found in proj4 and ended up implementing our own functions and an ellipsoid class. I recently thought of pymap3d and was poking around the code base again and there are some similarities with what we've been doing. I think there is opportunity for collaboration here and maybe porting some of our code into pymap3d. I always felt like these things were a bit out of place in Harmonica. Here are a few links:
ReferenceEllipsoid
class and handling functions: https://github.com/fatiando/harmonica/blob/master/harmonica/ellipsoid.pyReferenceEllipsoid
to compute Normal gravity: https://github.com/fatiando/harmonica/blob/master/harmonica/gravity_corrections.py#L10I really like our way of handling ellipsoids across the library. We have a global default ellipsoid and ways of setting different (or custom) ellipsoids using context managers. For example, you can change ellipsoids like this:
import harmonica as hm
# Set the ellipsoid globally
hm.set_ellipsoid("GRS80")
# Calculate with GRS80
normal_grav = hm.normal_gravity(latitude, height)
# Set the ellipsoid locally
with hm.set_ellipsoid("WGS84"):
# Calculate with WGS84
normal_grav = hm.normal_gravity(latitude, height)
# Set a custom ellipsoid
ell = ReferenceEllipsoid(name="sphere", semimajor_axis=1, inverse_flattening=1,
geocentric_grav_const=10, angular_velocity=1)
with hm.set_ellipsoid(ell):
normal_grav = hm.normal_gravity(latitude, height)
Functions that need an ellipsoid can get the currently set one by calling get_ellipsoid()
and everything works. See these examples.
I really like this mechanic and I feel like it could work well with pymap3d. I don't what's the best way of moving forward but I wanted to make sure we're all aware of the two projects 🙂 Notice that we need more than a, b, f for our ellipsoids because of the gravity calculations (which could be made optional).
How committed are you to supporting the Matlab toolkit interface? If your primary goal is compatibility then it might not work well but I think there is a lot of room for improvement on the API front (many of these Matlab toolkits are a bit poorly design, IMO).
What do you think?
Hi,
Love pymap3d.
When running my code with mypy without --ignore-missing-imports
, I get the error Skipping analyzing "pymap3d": module is installed, but missing library stubs or py.typed marker
.
Would adding py.typed suffice to add mypy support as described here?
Part of openjournals/joss-reviews#580
I see that there are tests and building the fortran library creates a testmaptran
program that seems to be run fortran tests. But I couldn't find instructions for doing this anywhere in the README, docs or CONTRIBUTING.md.
I found the commands in the .travis.yml
file but this is not accessible to users unfamiliar with continuous integration. Running the Python tests is not trivial because the package needs to be installed with [tests]
apparently and this adds a bunch of extra dependencies.
Even so, I was able to run all the tests on my machine without any problems. It would be nice to have these instructions in the README or CONTRIBUTING files.
The function geodetic2ecef prohibits negative altitude values here.
Can you explain the reasoning behind that choice? What if my altitude values are below the ellipsoid, as they are in many locations?
Thanks for your time!
Part of openjournals/joss-reviews#580
The API docs are complete and probably satisfy most advanced users who are familiar with all the terms and units involved. They can be made more accessible if the function docstrings also include a description of each input and output parameter with units and/or data types and possible accepted values.
For example, it's not clear from the docstring of ecef2geodetic how to specify an ellipsoid (I'm assuming that's what ell
is) or what happens if deg=False
.
This is not a requirement for the review, just a suggestion.
Part of openjournals/joss-reviews#580
It would be good to have more detailed install instructions. Things I found missing:
Regarding the Python instructions, pip install -e
is used for development but not recommended for deploying software because the user can't delete the source directory after install. I would recommend one of the following:
pip install pymap3d
cd pymap3d; pip install .
(no -e
)pip install https://github.com/scivision/pymap3d/archive/master.zip
1 is the best option but takes a bit more work to setup. With 3, you can also include instructions for installing specific versions using the same syntax but replacing master
with the tag name (see #7). The problem with 2 and 3 is that they don't play well with dependency managers. So it's not easy for another package to depend on pymap3d.
There are several syntax errors in Matlab code that cause execution errors. Specifically, the use of 'endif' rather than 'end' in ecef2enu.m and geodetic2ecef.m. The later also has an invalid comment syntax '##' on line 21.
When trying to convert to/from ECI coordinates with dates that don't have whole seconds (i.e. milliseconds are not zero), the sub-second information is ignored in the Julian date conversion: https://github.com/geospace-code/pymap3d/blob/main/src/pymap3d/sidereal.py#L87
When used to for ex. convert satellite trajectories, the resulting trajectory has stair-like shape.
No response
When importing:
import pymap3d as pm
and invoking:
pm.vincenty.vdist(lat1, lon1, lat2, lon2)
raises the exception:
AttributeError: module 'pymap3d' has no attribute 'vincenty'
In fact, the imported module will no longer import the vincenty module, breaking code that needed pymap3d.vincenty.vdist
This behavior was not present in previous installations of the package.
ECEF to ECI velocity conversion
In the tests functions for ECEF/ECI validation values are taken from matlab examples. The functions use position data, but not velocity data.
I found that the resulting velocity was different from that in the matlab documentation after conversion from ECEF to ECI.
I convert (3832 -4024 4837) from ECEF to ECI at datetime (2019, 1, 4, 12)
I expect: (-3384 -4887 4843)
I get: (-3004. -4669. 4842)
Python : 3.9.0
Pymap3d : 2.9.1
pymap3d is awesome, but seems to be missing a one basic coordinate system that would be useful: spherical ECEF. Currently one can do ecef2geodetic and vice versa, but not sph2geodetic and vice versa.
For example, if a user wants to convert geodetic lat/lon/alt to geocentric lat/lon/alt, they need to convert to XYZ ECEF with geodetic2ecef and then manually calculate spherical geocentric coordinates themselves. Is this by design? If not, would it be acceptable for me to submit a PR adding this to the library?
Describe the bug
What is your expected output value(s)?
What program/function are you comparing with?
Hi, I have been actively used an attribute eci2aer
of pymap3d
. However, today I suddenly faces an error when I use the attribute. Am I the only person who see this error: AttributeError: module 'pymap3d' has no attribute 'eci2aer'
?
conda create -n myenv python=3.8
conda activate myenv
pip3 install pymap3d
python
import pymap3d
pymap3d.eci2aer()
OUTPUT:
AttributeError: module 'pymap3d' has no attribute 'eci2aer'
Describe the bug
What is your expected output value(s)?
What program/function are you comparing with?
When I use geodetic2enu
import pymap3d as pm
lat0, lon0, alt0 = 4029.25745, 11146.61129, 1165.2
lat, lon, alt = 4028.25746, 11147.61125, 1165.1
e, n, u = pm.geodetic2enu(lat, lon, alt, lat0, lon0, alt0)
print(e, n, u)
I got an error like
ValueError: -pi/2 <= latitude <= pi/2
Can you help me ?
In this commit you added the line np.seterr(all="raise")
, which is always called when importing pymap3d.
That is breaking many other packages running numpy (e.g. everyone computing the tale of a Gaussian distribution will find a nice surprise). On top of that, it is unnecessary, since you can catch warnings (see https://docs.python.org/3/library/warnings.html#testing-warnings)
No response
for some points in my dataframe the vdist
azimuth solving part has a bug. however, if i look through the rows one at a time ( and pass scalar values instead of arrays to vdist) it works.
as a work-around it would be nice to have a `dist_only=True option which skips all this code.
from pymap3d.vincenty import vdist
n = 100_000
vdist(Lat1 = x.latitude.values[:n],
Lon1 = x.longitude.values[:n],
Lat2 = y.latitude.values[:n],
Lon2 = y.longitude.values[:n])
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
~/anaconda3/lib/python3.9/site-packages/pymap3d/vincenty.py in vdist(Lat1, Lon1, Lat2, Lon2, ell, dist_only)
247 i = sign(sin(lon2 - lon1)) * sign(sin(lamb)) < 0
--> 248 lamb[i] = -lamb[i]
249 except TypeError:
TypeError: 'float' object is not subscriptable
During handling of the above exception, another exception occurred:
ValueError Traceback (most recent call last)
/tmp/ipykernel_1431275/3970588599.py in <module>
1 from pymap3d.vincenty import vdist
2 n = 100_000
----> 3 vdist(Lat1=ship.latitude.values[:n],
4 Lon1 =ship.longitude.values[:n],
5 Lat2=sats.latitude.values[:n],
~/anaconda3/lib/python3.9/site-packages/pymap3d/vincenty.py in vdist(Lat1, Lon1, Lat2, Lon2, ell, dist_only)
248 lamb[i] = -lamb[i]
249 except TypeError:
--> 250 if sign(sin(lon2 - lon1)) * sign(sin(lamb)) < 0:
251 lamb = -lamb
252
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
Hi everyone!
I'd like to ask about how ned2geodetic is suposed to work, since I'm trying to convert points to global coordinates, but it seems that I'm missing something about the API.
I'm given an array of points with the following format:
x_coordinates=[3, 8,-3, 36, 12]
y_coordinates=[4, 45,67, -22, -16]
z_coordinates=[0, 0, 0, 0, 0]
Which represents a sequence of points by their coordinates X,Y,Z (they are in NED format) with respect an inmovile local frame located at the global geodetic coordinates (40.544289, -4.012101):
Geodetic_Origin=[40.544289,-4.012101]
If I undestood correctly, Geodetic_Origin is my observer,
Now I'm trying to get the global geodetic coordinates (i.e Latitude and Longitude) of my points, for which I tried using the ned2geodetic, implemented this way:
import numpy as np
import pymap3d as pm
for i in range (0,length):
geodetic_coordinates=ned2geodetic(x_coordinates[i], y_coordinates[i], z_coordinates[i], Geodetic_Origin=[0], Geodetic_Origin=[1], 0, ell=None, deg=True)
By my understanding, ned2geodetic converts North, East, Down to target latitude, longitude, altitude given the geodetic location of the observer. However, what I'm getting is a circle centered around my local origin, as shown here (used folium for the representation):
[![enter image description here][1]][1]
While I'm expecting something like this:
Any idea of what the problem might be?
EDIT: there was a mistake in the code, it's solved. I proceed to close the issue. Apologies for any inconvenience.
ecef2geodetic returns negative altitudes as positive.
See example code for 1 meter above/below the elliposide.
pymap3d.ecef2geodetic(6378137 - 1, 0, 0, pymap3d.ecef.Ellipsoid('wgs84'), deg=False)
(0.0, 0.0, 1.0)
pymap3d.ecef2geodetic(6378137 + 1, 0, 0, pymap3d.ecef.Ellipsoid('wgs84'), deg=False)
(0.0, 0.0, 1.0)
Invoking geodetic2aer using the following code:
`import pymap3d as pm
lat1 = 45.977
lon1 = 7.658
h1 = 4531
lat0 = 46.017
lon0 = 7.750
h0 = 1673
a,e,r = pm.geodetic2aer(lat1,lon1,h1,lat0,lon0,h0)
print ("A=%f E=%f R=%f" % (a,e,r))
`
Output
A=238.075833 E=-7134.628771 R=8876.843346
It should be (please don't take int account different floating point precision)
A=238.08 E=18.744 R= 8876.8
as you can check here https://it.mathworks.com/help/map/ref/geodetic2aer.html
Part of openjournals/joss-reviews#580
It's great that the project has a contributing guide. Well done!
Currently, the file a bit difficult to find and not very visible. It took me a while to find it and I only looked in the docs folder because I saw the commit message that added it. I would suggest moving it to the repository root and/or including a reference to it in the README so that it can be easily found.
This is not a requirement for the review, just a suggestion.
When attempting to import pymap3d, I get the following TypeError:
TypeError: expected string or bytes-like object
I have attempted to reinstall pymap3d, but the error persists.
Computer: MacOS
Python: 3.9.15
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/Users/e30737/miniconda3/envs/py39/lib/python3.9/site-packages/pymap3d/__init__.py", line 34, in <module>
from .aer import aer2ecef, aer2geodetic, ecef2aer, geodetic2aer
File "/Users/e30737/miniconda3/envs/py39/lib/python3.9/site-packages/pymap3d/aer.py", line 7, in <module>
from .ecef import ecef2enu, ecef2geodetic, enu2uvw, geodetic2ecef
File "/Users/e30737/miniconda3/envs/py39/lib/python3.9/site-packages/pymap3d/ecef.py", line 7, in <module>
from .eci import ecef2eci, eci2ecef
File "/Users/e30737/miniconda3/envs/py39/lib/python3.9/site-packages/pymap3d/eci.py", line 10, in <module>
import astropy.units as u
File "/Users/e30737/miniconda3/envs/py39/lib/python3.9/site-packages/astropy/units/__init__.py", line 17, in <module>
from .quantity import *
File "/Users/e30737/miniconda3/envs/py39/lib/python3.9/site-packages/astropy/units/quantity.py", line 25, in <module>
from astropy.utils.compat import NUMPY_LT_1_22
File "/Users/e30737/miniconda3/envs/py39/lib/python3.9/site-packages/astropy/utils/compat/__init__.py", line 19, in <module>
from .numpycompat import * # noqa
File "/Users/e30737/miniconda3/envs/py39/lib/python3.9/site-packages/astropy/utils/compat/numpycompat.py", line 14, in <module>
NUMPY_LT_1_19 = not minversion('numpy', '1.19')
File "/Users/e30737/miniconda3/envs/py39/lib/python3.9/site-packages/astropy/utils/decorators.py", line 546, in wrapper
return function(*args, **kwargs)
File "/Users/e30737/miniconda3/envs/py39/lib/python3.9/site-packages/astropy/utils/introspection.py", line 167, in minversion
return Version(module_version) >= Version(version)
File "/Users/e30737/miniconda3/envs/py39/lib/python3.9/site-packages/packaging/version.py", line 264, in __init__
match = self._regex.search(version)
TypeError: expected string or bytes-like object
Part of openjournals/joss-reviews#580
It's great that there are examples in the README. They would even more helpful if they included actual values and the expected outcomes so that we can check if we have the correct results.
This is not a requirement for the review, just a suggestion.
Hi folks! I want to install your library, but cannot do it because in the setup.py script it does not have all the information it needs. It creates some UNKNOWN library that I cannot import.
Can you please modify the setup.py script?
D:\Projects\pymap3d>C:/Python36/python setup.py install
running install
running bdist_egg
running egg_info
creating UNKNOWN.egg-info
writing UNKNOWN.egg-info\PKG-INFO
writing dependency_links to UNKNOWN.egg-info\dependency_links.txt
writing top-level names to UNKNOWN.egg-info\top_level.txt
writing manifest file 'UNKNOWN.egg-info\SOURCES.txt'
reading manifest file 'UNKNOWN.egg-info\SOURCES.txt'
reading manifest template 'MANIFEST.in'
writing manifest file 'UNKNOWN.egg-info\SOURCES.txt'
installing library code to build\bdist.win-amd64\egg
running install_lib
warning: install_lib: 'build\lib' does not exist -- no Python modules to install
creating build
creating build\bdist.win-amd64
creating build\bdist.win-amd64\egg
creating build\bdist.win-amd64\egg\EGG-INFO
copying UNKNOWN.egg-info\PKG-INFO -> build\bdist.win-amd64\egg\EGG-INFO
copying UNKNOWN.egg-info\SOURCES.txt -> build\bdist.win-amd64\egg\EGG-INFO
copying UNKNOWN.egg-info\dependency_links.txt -> build\bdist.win-amd64\egg\EGG-INFO
copying UNKNOWN.egg-info\top_level.txt -> build\bdist.win-amd64\egg\EGG-INFO
zip_safe flag not set; analyzing archive contents...
creating dist
creating 'dist\UNKNOWN-0.0.0-py3.6.egg' and adding 'build\bdist.win-amd64\egg' to it
removing 'build\bdist.win-amd64\egg' (and everything under it)
Processing UNKNOWN-0.0.0-py3.6.egg
Copying UNKNOWN-0.0.0-py3.6.egg to c:\python36\lib\site-packages
Adding UNKNOWN 0.0.0 to easy-install.pth file
Installed c:\python36\lib\site-packages\unknown-0.0.0-py3.6.egg
Processing dependencies for UNKNOWN==0.0.0
Finished processing dependencies for UNKNOWN==0.0.0
Some inputs to vdist, noteable equatorial arcs, will return (nan, nan, nan)
Example run using a Python 3 notebook in Google Collab after pip installing pymap3d=1.8.1:
from pymap3d.vincenty import vdist
for i in range(1, 11):
print(f'vdist(0, 0, 0, {i}) -> {vdist(0, 0, 0, i)}')
This returns:
vdist(0, 0, 0, 1) -> (111319.4907932264, 90.0, 270.0)
vdist(0, 0, 0, 2) -> (222638.9815864528, 90.0, 270.0)
vdist(0, 0, 0, 3) -> (333958.4723796792, 90.0, 270.0)
vdist(0, 0, 0, 4) -> (445277.9631729056, 90.0, 270.0)
vdist(0, 0, 0, 5) -> (nan, nan, nan)
vdist(0, 0, 0, 6) -> (nan, nan, nan)
vdist(0, 0, 0, 7) -> (nan, nan, nan)
vdist(0, 0, 0, 8) -> (890555.9263458112, 90.0, 270.0)
vdist(0, 0, 0, 9) -> (1001875.4171390376, 90.0, 270.0)
vdist(0, 0, 0, 10) -> (1113194.907932264, 90.0, 270.0)
/usr/local/lib/python3.6/dist-packages/pymap3d/vincenty.py:175: RuntimeWarning: invalid value encountered in arcsin
sin(lamb[notdone]) / sin(sigma[notdone])))
Altering values by amount less than machine precision will also trigger this behavior:
print(f'vdist( 0, 0, 0, 1) -> {vdist(0, 0, 0, 1)}')
print(f'vdist(1e-16, 1e-16, 1e-16, 1) -> {vdist(1e-16, 1e-16, 1e-16, 1)}')
This returns:
vdist( 0, 0, 0, 1) -> (111319.4907932264, 90.0, 270.0)
vdist(1e-16, 1e-16, 1e-16, 1) -> (nan, nan, nan)
/usr/local/lib/python3.6/dist-packages/pymap3d/vincenty.py:175: RuntimeWarning: invalid value encountered in arcsin
sin(lamb[notdone]) / sin(sigma[notdone])))
It looks like the error correction is causing this and there should be a check that the values going into arcsin are valid.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.