geospace-code / pymap3d

pure-Python (Numpy optional) 3D coordinate conversions for geospace ecef enu eci
https://geospace-code.github.io/pymap3d
BSD 2-Clause "Simplified" License
391 stars 85 forks source link

[Bug]: #86

Open AlexAtCCRI opened 1 year ago

AlexAtCCRI commented 1 year ago

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()
Dobatymo commented 1 year ago

Encoutering the same issue.

Dobatymo commented 1 year ago

Here https://github.com/geospace-code/pymap3d/blob/81172e221a4c9884450a38ec8a7ee382198cb7e3/src/pymap3d/vincenty.py#L172 and https://github.com/geospace-code/pymap3d/blob/81172e221a4c9884450a38ec8a7ee382198cb7e3/src/pymap3d/vincenty.py#L232 lamb is set to a float, even when the inputs are numpy arrays. The cause an exception here https://github.com/geospace-code/pymap3d/blob/81172e221a4c9884450a38ec8a7ee382198cb7e3/src/pymap3d/vincenty.py#L264 which is then handled incorrectly.

The fix should be to assign arrays to lamb if the input is arrays.

The whole approach to handle numpy arrays and normal floats with the same code seems a bit brittle to me however. I would create two separate implementations, even that means some code duplication.