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' ?

How to reproduce the error?

1. conda create -n myenv python=3.8
2. conda activate myenv
3. pip3 install pymap3d
4. python
5. import pymap3d
6. pymap3d.eci2aer()

OUTPUT:
AttributeError: module 'pymap3d' has no attribute 'eci2aer'

What happened?

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.

Relevant log output

No response

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)

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

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.

What happened?

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

• OS Debian GNU/Linux 32bit ARM processor, little endian
• Python 3.11

See the detailed log at:
https://ci.debian.net/data/autopkgtest/testing/armel/p/pymap3d/34582756/log.gz

Relevant log output

============================= 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 =================

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.

What happened?

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)

Relevant log output

No response

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,

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.

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.

https://github.com/gberrante/map_3d

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

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.

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

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.

What happened?

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?

Relevant log output

No response

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.

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,
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?

The current homepage https://pypi.org/project/pymap3d/ is down

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.

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.

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.

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!

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.

What happened?

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.

Relevant log output

No response

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

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

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 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?

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.

Part of openjournals/joss-reviews#580

The tag for the latest version is v.1.6.0 instead of v1.6.0.

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:

https://answers.opencv.org/question/219040/fastest-way-to-convert-bgr-rgb-aka-do-not-use-numpy-magic-tricks/

The same is almost certainly true for matplotlib too.

With best regards,
Johnny

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.

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.

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?

What happened?

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

Relevant log output

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

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 floats, not three ndarrays.

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.

What happened?

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!

Relevant log output

No response

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

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.

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 ?

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.

Part of openjournals/joss-reviews#580

It would be good to have more detailed install instructions. Things I found missing:

• How to download the source (command to clone the repository or link to download a zip archive from Github)
• How to install the Matlab version

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:

1. Make a release and upload to PyPI. Then the instructions should be: pip install pymap3d
2. Download the zip/clone the repo and install using pip: cd pymap3d; pip install . (no -e)
3. Install directly from master using pip: 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.

Action to publish release v3.0.1 to pypi failed. I think it's because this version string wasn't updated

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.

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.

What happened?

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.

Relevant log output

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 @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:

• Our ReferenceEllipsoid class and handling functions: https://github.com/fatiando/harmonica/blob/master/harmonica/ellipsoid.py
• Conversion from geodetic to geocentric spherical: https://github.com/fatiando/harmonica/blob/master/harmonica/coordinates.py
• Using ReferenceEllipsoid to compute Normal gravity: https://github.com/fatiando/harmonica/blob/master/harmonica/gravity_corrections.py#L10

I 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?