mrjean1/pygeodesy

Pure Python geodesy tools

pygeodesy's Issues

example ellipsoidalVincenty.py does not work

the following example can be found on the top of the file:

Here's an example usage of Vincenty:
   
      >>> from geodesy.ellipsoidalVincenty import LatLon
      >>> Newport_RI = LatLon(41.49008, -71.312796)
      >>> Cleveland_OH = LatLon(41.499498, -81.695391)
      >>> print(Newport_RI.distance(Cleveland_OH))
      866,457.166175  # meter == 538.3904451566326 miles

LatLon object does not have the distance() function. Probably this should be distanceTo()

Euclidean distance.

There is/was a bug in function euclidean_ affecting the HeightIDW interpolation. Fixed in release PyGeodesy 19.4.2.

Check for convex polygon seems to fail in isenclosedBy of sphericalTrigonometry

Hi,

sphericalTrigonometry.LatLon.isenclosedBy seems to detect polygons with clockwise ordered vertices as concave. Using CCW order or any alternative functions seem to work fine.

Example:

from pygeodesy.sphericalNvector import LatLon as LatLonN
from pygeodesy.sphericalTrigonometry import LatLon as LatLonT
from pygeodesy.points import isenclosedBy

boundsT_CCW = LatLonT(0, 0), LatLonT(0, 1), LatLonT(1, 1), LatLonT(1, 0)
boundsT_CW  = LatLonT(0, 0), LatLonT(1, 0), LatLonT(1, 1), LatLonT(0, 1)
boundsN_CCW = LatLonN(0, 0), LatLonN(0, 1), LatLonN(1, 1), LatLonN(1, 0)
boundsN_CW  = LatLonN(0, 0), LatLonN(1, 0), LatLonN(1, 1), LatLonN(0, 0)

isenclosedBy(LatLonN(0.5, 0.5), boundsN_CCW)  # OK.
isenclosedBy(LatLonN(0.5, 0.5), boundsN_CW)  # OK.
isenclosedBy(LatLonT(0.5, 0.5), boundsT_CCW)  # OK.
isenclosedBy(LatLonT(0.5, 0.5), boundsT_CW)  # OK.
LatLonN(0.5, 0.5).isenclosedBy(boundsN_CCW)  # OK.
LatLonN(0.5, 0.5).isenclosedBy(boundsN_CW)  # OK.
LatLonT(0.5, 0.5).isenclosedBy(boundsT_CCW)  # OK.
LatLonT(0.5, 0.5).isenclosedBy(boundsT_CW)  # "not convex"

Maybe the check in sphericalTrigonometry.py, line 557 should be opposite, when the points are CW?

            # check for convex polygon, the angle between
            # gc vectors, signed by direction of n0
            # (otherwise the test above is not reliable)
            if gc1.angleTo(gc, vSign=n0) < 0:  # > 0 for CW?
                t = _Fmt.SQUARE(points=i)
                raise _ValueError(t, p, txt=_not_(_convex_))

PyGeodesy==21.5.31, Python 3.9.4 64-bit, Windows

Is it possible to calculate elevation above sea level?

Hi mrJean1!
I am trying to retrieve the elevation of many <lat, lon> pairs. I had no problem using PyGeodesy to do the EGM2008 lookups as your great documentation describes. For example Calgary, AB is at an elevation of ~1045 m. The lookup gives a different result.

>>> from pygeodesy import *
>>> predictor = GeoidPGM("/home/caleb/Desktop/egm2008-5.pgm")
>>> pt = points.LatLon_(51.05011, -114.08529)
>>> predictor(pt)
-16.804670290052883

I realize the lookup isn't returning elevation above sea level but rather the height of the geoid above the WGS84 ellipsoid. I also realize the values returned are consistent with other tools.

I was wondering if there was a way to use PyGeodesy to convert that result into elevation above sea level? I apologize if this is out of scope or a dumb question; geoinformatics is not my field by trade.

Sincerely,
Caleb

a.nearestOn(b, c) returns incorrect result when b == c

Example:

In [1]: from pygeodesy.sphericalNvector import LatLon

In [2]: LatLon(1, 1).nearestOn(LatLon(2, 2), LatLon(2, 2))
Out[2]: LatLon(00°00′00.0″N, 000°00′00.0″E)

Expected output is

Out[2]: LatLon(02°00′00.0″N, 002°00′00.0″E)

Tried with versions 17.06.04 and 17.08.24.

Can i use this library on with asyncio

Does this library support async methods with asyncio?

Vectorized numpy operations

Hello, I have large vectors of lat, lon, alt that I would like to convert to a local coordinate system some "enu" frame such that I can do calculations the ecef.Cartesian and the ecef.EcefKarney are doing the job but it is insanely slow to iterate over such a large number ~400,000 points. I believe that this would be much much faster if it could be vectorized with numpy and accept numpy vectors as input. Is this in anyway planned ?

what kind of effort would you estimate it would be to make such a change ? is there an alternative lib that I missed ?

I was also curious to know why the internal data types and not depend on numpy under the hood :), thanks in advance and for the great development :)

GeoidKarney leaves resources open

The GeoidKarney class fails to close the underlying EGM file when the class goes out of scope.

Below is a minimal example to describe the problem.

from pygeodesy.geoids import GeoidKarney
import unittest

class TestGeoidDestructor(unittest.TestCase):
    def test_karney(self):
        geoid = GeoidKarney("/usr/local/share/GeographicLib/geoids/egm96-5.pgm")

if __name__ == "__main__":
    unittest.main()

Which generates the following output and warning...

/usr/lib/python3.8/unittest/case.py:633: ResourceWarning: unclosed file <_io.BufferedReader name='/usr/local/share/GeographicLib/geoids/egm96-5.pgm'>
  method()
ResourceWarning: Enable tracemalloc to get the object allocation traceback
.
----------------------------------------------------------------------
Ran 1 test in 0.003s
OK

A quick look through the source code also showed that there was nowhere that the self._egm was being closed.

Automatic closure of the resource or the addition of an explicit close method to the geoid interface would be a great addition to the project.

---

Technical Details
Python 3.8.5
pygeodesy 21.03.03

intersection

Hi,
Is there a way using PyGeodesy to compute the intersection of two lines defined each by latlon and bearings but for one point, I want to find the intersection point going in Bearing direction and bearing + 180° ?
I know A and bearing_A and C and bearing_C. I want to find G but I don't know if it on the right or left of C.

Thanks

sphericalNvector.LatLon.intermediateTo doesn't interpolate height

There is a bug in pygeodesy.sphericalNvector module LatLon class.
The .intermediateTo method doesn't interpolate the height of points. Behavior differs from the LatLon class in pygeodesy.sphericalTrigonometry module.

The error can be seen from the code snippet below:

>>> import pygeodesy.sphericalTrigonometry as spt
>>> import pygeodesy.sphericalNvector as snv
>>> point_a_trig = spt.LatLon(1, 1, height=100)
>>> point_b_trig = spt.LatLon(2, 2, height=200)
>>> point_a_nvec = snv.LatLon(1, 1, height=100)
>>> point_b_nvec = snv.LatLon(2, 2, height=200)
>>> point_a_trig.intermediateTo(point_b_trig, 0.25)
LatLon(01°15′00.15″N, 001°14′59.71″E, +125.00m)
>>> point_a_nvec.intermediateTo(point_b_nvec, 0.25)
LatLon(01°15′00.15″N, 001°14′59.71″E)

Utm.toLatLon() goes into infinite loop

The function call:

utm.Utm(zone, hemisphere, easting, northing).toLatLon(None)

runs into infinite loop, when called for the following parameters:

zone = 55
hemisphere = 'S'
easting = 321441.0425108216
northing = 5810117.133231169

It cannot reach desired epsilon

Why intern floats?

I hope it's okay to just ask a question and not a bug.

I was reading through the code and I noticed the logic in the interns module relating to maintaining a dict of the floating point constants defined in it. Why is this being done?

I see how len(inters._floats) is used in reporting out the number of constants defined when you python -m pygeodesy, is that the only purpose for maintaining this dict, or is it doing some kind of useful interning of the constants?

Not sync'd

Directories dist, doc and testresults are obsolete and not sync'd: GitHub push blocked due to sudden error "personal email address".

TypeError with sphericalNvector.triangulate

Hello,

Thanks for putting this package together.

I have a TypeError using triangulate:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from geodesy import sphericalNvector
from geodesy.sphericalNvector import LatLon

basse_castouillet = LatLon("47°18.228'N","002°34.326'W")
basse_hergo = LatLon("47°18.664'N","002°31.717'W")

p0 = sphericalNvector.triangulate(basse_castouillet, 7, basse_hergo, 295)
print(str(p0))

I am getting:

Traceback (most recent call last):
  File "./t1.py", line 10, in <module>
    p0 = sphericalNvector.triangulate(basse_castouillet, 7, basse_hergo, 295)
  File "/Users/user/PyGeodesy/geodesy/sphericalNvector.py", line 712, in triangulate
    return Nvector(i.x, i.y, i.z).toLatLon(height=h, datum=point1.datum)
TypeError: toLatLon() got an unexpected keyword argument 'datum'

The package tests are running fine, am I using the package the right way ?

isue with toUtm given LatLon()

This is probably something dumb I am doing, but the toUtm method documentation shows zone as an optional arg; however when trying to use this method, it seems to be asking me for it.

from pygeodesy import Utm
from pygeodesy.sphericalNvector import LatLon
p1Lat = 47.3444439424917
p1Long = -92.6537372154412
p1 = LatLon(p1Lat,p1Long)
p1
LatLon(47°20′40.0″N, 092°39′13.45″W)
u1 = Utm.toUtm(p1)
Traceback (most recent call last):
File "", line 1, in
TypeError: toUtm() missing 1 required positional argument: 'zone'

when I try to give it the correct zone (although I'm not sure how to get that w/o converting to UTM), it still flakes out:

u = Utm.toUtm(p,zone=15)
Traceback (most recent call last):
File "", line 1, in
File "/home/elliot/.local/lib/python3.8/site-packages/pygeodesy/utm.py", line 562, in toUtm
if zone == self.zone and falsed == self.falsed:
AttributeError: 'LatLon' object has no attribute 'zone'
u = Utm.toUtm(p,zone='15')
Traceback (most recent call last):
File "", line 1, in
File "/home/elliot/.local/lib/python3.8/site-packages/pygeodesy/utm.py", line 562, in toUtm
if zone == self.zone and falsed == self.falsed:
AttributeError: 'LatLon' object has no attribute 'zone'

Find the shortest distance between a point and line segments (not line)

Hi,
Is there to find the shortest distance between a point and line segments (not line) like in this image from https://stackoverflow.com/questions/27161533/find-the-shortest-distance-between-a-point-and-line-segments-not-line

The algorithm works fine in regular geometry but I am not sure with LatLon...
Thanks

pip?

Great job that you are converting the JS scripts to Python!
Are you planning to make this library accessible from pip?

EcefKarney.forward method error

Hi there,

I am getting an error when I provide two floats in the EcefKarney.forward method shown below. The two floats are labeled as Ytd (longitude) and Xtd (latitude) as decimal WGS84 points. Please advise. Thanks!

Ytd = degrees(asin(sin(Y0r)*cos(d/r)+cos(Y0r)*sin(d/r)*cos(yaw)))
Xtd = degrees(X0r+atan2(sin(yaw)*sin(d/r)*cos(Y0r), cos(d/r)-sin(Y0r)*sin(radians(Ytd))))
ecef1 = EcefKarney.forward(Ytd, Xtd, 0, 0)

Error message:

Traceback (most recent call last):
File "d:\anaconda\envs\azel\azel\georec.py", line 38, in
ecef1 = EcefKarney.forward(Ytd, Xtd, 0, 0)
File "D:\anaconda\envs\azel\lib\site-packages\pygeodesy\ecef.py", line 341, in forward
return _EcefBase._forward(self, *llhn, M=M)
File "D:\anaconda\envs\azel\lib\site-packages\pygeodesy\ecef.py", line 228, in _forward
E = self.ellipsoid
AttributeError: 'float' object has no attribute 'ellipsoid'

Python 3+ tests on Windows

All PyGeodesy 17.9.22 tests pass with 32-bit Python 2.7.14 on Windows 10 Pro, see the log file in the tests results folder. The tests have also been run with 32-bit Python 3.6.3 now Windows 7 Pro and Windows 10 Pro, but there is still a unicode code to be resolved. In addition, testing with 64-bit Python on Windows is still to be done.

Pygeodesy 20.10.14 broken wheel.

The PyGeodesy-20.10.14-py2.py3-none-any.whl distribution file is still broken due a bad setup.py script. Although other distribution files are OK, release PyGeodesy 20.10.14 has been yanked from PyPI. Release PyGeodesy 20.10.15 is forthcoming.

FAIL: test_WebMercator (test.unitTestSuite.TestSuite)

• Python 2.7 on Linux https://travis-ci.com/cclauss/PyGeodesy/jobs/132401482#L502
• Python 3.7 on Linux https://travis-ci.com/cclauss/PyGeodesy/jobs/132401483#L505

======================================================================
FAIL: test_WebMercator (test.unitTestSuite.TestSuite)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/travis/build/cclauss/PyGeodesy/test/unitTestSuite.py", line 97, in test_WebMercator
    self._run('testWebMercator')
  File "/home/travis/build/cclauss/PyGeodesy/test/unitTestSuite.py", line 31, in _run
    self.assertEqual(x, 0)
AssertionError: 1 != 0
----------------------------------------------------------------------

Pip does not find older versions

It seems that only recent versions of this library are available on Pypi. This breaks my build every week or so, as I pin down the exact version dependencies in my requirements.txt file, like:

PyGeodesy==17.5.25

I would guess this is a reasonably standard practice, so maybe the older versions of the library could be kept on Pypi?

No online docs.

Until further notice, the PyGeodesy docs are no longer available on PythonHosted. The complete documentation is included in the zip distribution file at PyPI. The wheel distribution does not include any documentation, tests and test results.

Issue with simplify functions resolved

There was a serious bug in all simplify functions except one, resulting in incorrect results. The problem has been fixed in PyGeodesy-17.7.25 and later. Additional simplify tests have been added to testSimplify.py, including comparisons with 2 other RDP implementations.

wrong .sizes() output

first_geohash = geohash_freqs[‘geohash’].iloc[0]
sizes = pygeodesy_geohash.sizes(first_geohash)
print(first_geohash,sizes)
u336xv (4890.0, 4890.0)

for precision=6 it’s returning ( 4890, 4890, 2758.887), # 5
but it should return ( 610, 1220, 486.710), # 6

user documentation?

Is there any user documentation for this package? (I have been unable to find any).

Trilaterate by area

Is there any way to trilaterate when some of the circles are not intersecting with the perimeter but with the area?

For example, I would like to trilaterate this circles, but it only works when they are intersecting with the perimeter:

from pygeodesy.sphericalNvector import LatLon

distance = 5110

p1 = LatLon(42.688839, 2.438857)
p2 = LatLon(42.635421, 2.522570)
p3 = LatLon(42.64540, 2.504811)
t = p1.trilaterate(distance, p2, distance, p3, distance)
print('trilaterate:', t.lat, t.lon) # (INCORRECT: 43.051858748236725, 2.944630261790253)


p1 = LatLon(42.688839,2.438857)
p2 = LatLon(42.635421,2.522570)
p3 = LatLon(42.630788,2.500713)
t = p1.trilaterate(distance, p2, distance, p3, distance)
print('trilaterate:', t.lat, t.lon) # (CORRECT: 42.67456065072384, 2.4953950233084456)

Incorrect:

Correct:

Along-track distance

It would be convenient to have a method for computing along-track distance: "distance from the start point to the closest point on the path (to the end point) to the third point". The formula (as provided here) is:

dAt = acos(cos(d13 / R) / cos(dXt / R)) * R

Where dXt is cross-track distance and d13 is the distance from start to the point of interest.

It is also possible to modify it to support negative values, which makes sense in certain cases.

I would be happy to add this method if you're considering accepting pull requests.

LatLon.nearestPointOnSegment() not working properly?

Hi mrJean1,

I think I've found an issue with the nearestPointOnSegment method in sphericalNvector.LatLon.

>>> from PyGeodesy.geodesy.sphericalNvector import LatLon
>>> a = LatLon(0, 0)
>>> b = LatLon(0, 1)
>>> x = LatLon(1, 0)
>>> x.nearestPointOnSegment(a, b)
LatLon(00°00′00.0″N, 180°00′00.0″E)

Cheers,
Eilam

Height not interpolated in nearestOn methods

There is a bug in pygeodesy.sphericalNvector module related to heights.

pygeodesy.sphericalNvector.LatLon.nearestOn and pygeodesy.sphericalNvector.LatLon.nearestOn2 do not interpolate the height on the line, only in case if endpoint is nearest the height is returned. Code to recreate:

>>> import pygeodesy.sphericalNvector as snv
>>> point_a = snv.LatLon(1, 1, height=100)
>>> point_b = snv.LatLon(2, 2, height=200)
>>> # Nearest is close to middle.
>>> point_c = snv.LatLon(1, 2)
>>> point_c.nearestOn(point_a, point_b)
LatLon(01°29′59.93″N, 001°29′59.31″E)
>>> point_c.nearestOn2([point_a, point_b])[0]
LatLon(01°29′59.93″N, 001°29′59.31″E)
>>> # Nearest is endpoint.
>>> point_d = snv.LatLon(0, 0)
>>> point_d.nearestOn(point_a, point_b)
LatLon(01°00′00.0″N, 001°00′00.0″E, +100.00m)
>>> point_d.nearestOn2([point_a, point_b])[0]
LatLon(01°00′00.0″N, 001°00′00.0″E, +100.00m)

trying to get coordinates in the center of a mgrs cell

Thanks for this great library, and sorry for the doom question.

This is what I am doing to convert mgrs reference to both utm and lat/lon coordinates.

mgrs_string = "31TDF3182"
from pygeodesy import mgrs
utm = mgrs.parseMGRS(mgrs_string).toUtm()
print(utm)
latlon = utm.toLatLon()
print(latlon)

output:

31 N 431000 4582000
(41.38657, 2.174726, Datum(name='WGS84', ellipsoid=Ellipsoids.WGS84, transform=Transforms.WGS84), -0.54564, 0.999659)

This works fine, but I am getting the coordinates of the bottom left corner of my mgrs cell.
I want to get the cell's center point.
As in this case the cell size is 1 square Km, it would be just a matter of adding 500 to each utm coordinate, but I don't know what is the best/quickest approach to do it.
So I have some questions related to this:

• Is there any built in parameter to ask mgrs.toUtm() to return the center of mgrs cell?
• Perhaps a built in method to sum a certain number of meters to x and/or y coordinates of the utm object?
• Or at least, could I get Northing and Easting from utm object in a direct way? latlon.lat gives me latitude, but I couldn't get any values from utm.E , utm.Easting , utm['E'] ...

I have not great python skills: what I am doing for now is parsing utm.toStr(), and then splitting the string to get utm Northing and Easting in meters. Then, calculating how big the cell is I know how many meters I should add to N & E ... which looks a bit too much work to me.

Thanks again !!

Availability of GeodesicExact

I need to calculate geodesics on ellipsoids with high flattening, and from my understanding this is possible with the GeodesicExact class of GeographicLib. My problem is that this class isn't available in the Python version of GeographicLib, so a search for alternate options led me here. Is the GeodesicExact class available in PyGeodesy? I saw a reference to it in the PyGeodesy documentation for the elliptic module ("This provides the elliptic functions and integrals needed for Ellipsoid, GeodesicExact, and TransverseMercatorExact."). But I searched throughout the documentation and I couldn't find any other references to it. I apologize if I'm missing something obvious.

If it's not available here, I would appreciate any other suggestions for calculating geodesics on ellipsoids with high flattening using Python, since I'm having a lot of trouble with this.

Intersection of 2 polygons

Hello,

Not an issue, but a question (or a feature request).
Is there in the package a function to know if two polygons intersect ?

Thank you very much.

import problem on LazyImportError

I encountered some problems when using pygeodesy in my own code. The problem comes from these lines:

https://github.com/mrJean1/PyGeodesy/blob/master/pygeodesy/__init__.py#L241-L246

When it cannot import the LazyImportError, it can also not use it the except branch.

Strangely enough i only encountered this issue during compilation with pyinstaller. will see if i can create a reproducable script, but i think it is already clear from the code that this won't work.

environment:

• Mac OS
• python3.6
• only manifests after compilation with pyinstaller

using pygeodesy in standalone application (pyinstaller)

i am having some problems with getting pygeodesy to work when compiled using pyinstaller. it seems that the problems originate from the modifications to sys.path and the checks which are taking place in pygeodesy/__init__.py.

setup:

• Mac OS 10.14.4
• script.py with following content:

from pygeodesy.ellipsoidalVincenty import LatLon
_ = input('Start dist calc?')
print('distance=', LatLon(50, 4).distanceTo(LatLon(51, 4)))

• using system python 3.6 (not virtualenv, as this causes problems with pyinstaller)
  • pip3.6 install pyinstaller
  • pip3.6 install pygeodesy

Running normally

normally running the script with python3.6 script.py works correctly as expected.

Compile

now compiling the script with the following command pyinstaller script.py.
this produces a dist folder in which the script executable can be found. running this executable leads to the following error:

Traceback (most recent call last):
  File "script.py", line 1, in <module>
  File "<frozen importlib._bootstrap>", line 971, in _find_and_load
  File "<frozen importlib._bootstrap>", line 955, in _find_and_load_unlocked
  File "<frozen importlib._bootstrap>", line 665, in _load_unlocked
  File "/usr/local/lib/python3.6/site-packages/PyInstaller/loader/pyimod03_importers.py", line 631, in exec_module
    exec(bytecode, module.__dict__)
  File "pygeodesy/__init__.py", line 245, in <module>
NameError: name 'LazyImportError' is not defined
[73850] Failed to execute script script

This problem originates from this line in the init file, which adds the pygeodesy folder to sys.path. This dynamic behavior is not caught by pyinstaller.

Compile with added path

Explicitly adding the pygeodesy path to pyinstaller fixes the above problem:
PYGEODESY_PATH=$(python3.6 -c "import pygeodesy; print(pygeodesy.__path__[0])")
pyinstaller --paths=$PYGEODEST_PATH script.py

now running the new executable produces the following error:

Traceback (most recent call last):
  File "script.py", line 1, in <module>
  File "/usr/local/lib/python3.6/site-packages/PyInstaller/loader/pyimod03_importers.py", line 631, in exec_module
    exec(bytecode, module.__dict__)
  File "pygeodesy/__init__.py", line 376, in <module>
  File "pygeodesy/__init__.py", line 369, in _ismodule
ImportError: foreign module 'ellipsoidalKarney' from '/Users/matthijs/pygeodesy_bug/dist/script/ellipsoidalKarney.pyc'

It fails on the _ismodule method in the pygeodesy/__init__.py. Adding some prints leads to the following values for dirname(p) and pygeodesy_abspath:

• dirname(p): '.../dist/script'
• pygeodesy_abspath: '.../dist/script/pygeodesy'

These are indeed not equal and thus raising the exception

Conclusion

Checking whether the imported module is indeed part of the pygeodesy distribution seems too restrictive to work with pyinstaller, which layouts the files differently. Would it be an idea to remove this restriction or add an environment flag to control this behaviour?

if os.environ.get('PYGEODESY_SKIP_ISMODULE', 'false') == 'true:
    pass
else:
   # perform _ismodule(...)

Support for Python 3.7 in setup.py?

https://github.com/mrJean1/PyGeodesy/blob/master/setup.py#L98-L100

Extend intermediateTo to support any number as a fraction

Currently it only supports range 0..1 but it is not too difficult to extend it to support an arbitrary distance along great circle defined by self and other including negative values.

PyGeodesy 20.10.9 and 20.10.12 broken.

Both PyGeodesy 20.10.9 and 20.10.12 are bad due to the missing pygeodesy.deprecated package in the MANIFEST.in and all distribution files. Do not use those releases! The issues have been fixed in PyGeodesy 20.10.14.

Conversion from LatLon to OSGR and back may cause systematic offset

Conversion between Lat/Lon and OSGR is not always reversible. See the following examples:
osgr = toOsgr(eV.LatLon(52, -0.12))
osgr.toStr2()
'[G:TQ, E:32014, N:23971]'
ll = osgr.toLatLon(eV.LatLon); (ll.lat, ll.lon)
(51.0, -0.12) # consistent to input
(osgr.easting, osgr.northing)
(532014.2352530259, 123971.60890206805)
Osgr(osgr.easting, osgr.northing)
[G:TQ, E:32014, N:23971]
o = Osgr(532014, 123971).toLatLon(eV.LatLon); (o.lat, o.lon)
(50.99918767451807, -0.12151217235660071) # not consistent to input
https://www.movable-type.co.uk/scripts/latlong-os-gridref.html provides
(50.999995, 0.120004) for 'TQ 32014 23971'.
toOsgr(eV.LatLon(o.lat, o.lon))
[G:TQ, E:31910, N:23878] # consistent to result from https://www.movable-type.co.uk/scripts/latlong-os-gridref.html for (50.99918767451807, -0.12151217235660071)

Test failure on Travis

The testCss failure with Python 2.6.9 and Ubuntu 14.04 on Travis is a minor issue:

test 19 Css.latlon: (50.899999999999999, 1.8) FAILED, expected (50.9, 1.8)

`initialBearingTo` returns wrong values for points on same meridian

I didn't go into too much research reproducing this, but it seems that the method pygeodesy.sphericalNvector.LatLon.initialBearingTo always returns 0.0 if both points are on the same meridian. It should return 180.0 if the point_from is to the north from the point_to. Results seem to be correct as long as points are shifted at least a bit on the Lon axis.

Code snippet to reproduce the issue:

>>> import pygeodesy
>>> pygeodesy.__version__
'21.06.29'
>>> from pygeodesy.sphericalNvector import LatLon
>>> point_a = LatLon(5, 5)
>>> point_b = LatLon(10, 5)
>>> point_a.initialBearingTo(point_b)
0.0
>>> point_b.initialBearingTo(point_a)
0.0
>>> point_c = LatLon(10, 5.0000001)
>>> point_a.initialBearingTo(point_c)
1.1299401739210769e-06
>>> point_c.initialBearingTo(point_a)
180.00000114300522

Trilateration function returns an incorrect result

I ran into some odd results while attempting to locate a point using the trilaterate function (calculated points were almost thousands of meters far from their actual locales), so I decided to check these results against another Python implementation based on the algorithm found in Movable Type website, the JavaScript equivalent to this repository. It turned out that the "same" algorithm was returning points with centimeters precision.

Comparing the equations in the link above to the ones found in the sphericalNvector.py file, the following lines seems a little off to me:

PyGeodesy/pygeodesy/sphericalNvector.py

Lines 1086 to 1087 in 723f27f

	x = fsum_(d12, -d22, d**2) / d # n1->intersection x- and ...
	y = fsum_(d12, -d32, i**2, j**2, -2 * x * i) / j # ... y-component

According to the website, it should be something like this instead:

x = fsum_(d12, -d22, d**2) / (2 * d)  # n1->intersection x- and ...
y = fsum_(d12, -d32, i**2, j**2) / (2 * j) - (x * i / j) # ... y-component

After these changes, both algorithms (the one from the link above and the implementation in this project) returned approximately the same points with centimeters precision, as I said earlier.

I tried to open a pull request with this possible fix, but I have no permission to push to this repository. Nonetheless, I willing to do so if you want me to.

Intersection of circles on ellipsoids.

Pygeodesy 20.7.29 includes 2 functions/methods to intersect circles on ellipsoidal datums, ellipsoidalKarney.intersections2 and ellipsoidalVincenty.intersections2. Both use internal function ellipsoidalBase._intersect2 to do the work based on an iterative procedure involving azimuthal equidistant projections as suggested by Karney. See the documentation and comments for details and references.

Also, Pygeodesy 20.7.29 offers an other, generic function pygeodesy.intersections2 in module pygeodesy.formy to compute circle intersections for 2 sets of plain lat, lon and radius values.

The spherical intersections2 tests in testSpherical.py have been enhanced and tests have been added to testEllipsoidal.py for the new ellipsoidal and to testFormy.py the the generic intersections2 functions.

However, there are still some questions about the test results, in particular about the random tests in testEllipsoidal.py using the azimuthal.Equidistant projection. The results of those differ significantly from same tests using the azimuthal.EquidistantKarney projection and from the random tests in testSpherical.py. That discrepancy is still under investigation.

Do you consider putting docs online?

You probably could just change name of doc.html to index.html and enable GitHub pages from master branch.

Please add some information about how to use and test this library

Please add some information about how to:

1. install
2. use and
3. test this library,

in the README file.

Example runs showing the output would be greatly appreciated.

Paths intersection.

The results from the spericalNvector and -Trigonometry LatLon.intersection methods are inconsistent for paths both given as a point and bearing: the intersection returned by one is occasionally the antipode of the result from the other.

LatLon _Nv is not updating with lat/lon changes

First, thank you so much for converting this to python.

Been using sphericalNvector.LatLon and noticed if I update the lat/lon after I call distanceTo(), the distance is calculated using the old values (cached _Nv vector with old values) and never updates. This probably results in errors for other vector calculations as well since _Nv does not update.

Example:
p1 = LatLon(52.205, 0.119)
p2 = LatLon(48.857, 2.351)
d1 = self.distance(p2) # 404.3 km
p1.lat = 20.0
p1.lon = -60.0
d2 = self.distance(p2) # Remains 404.3 km...should be the distance between 20.0, -60.0 and p2

In the above example, d2 remains 404.3 km because it is still using the cached _Nv from the original p1 location.

Intersection of circles on surface

Is there any way to get the intersections of two circles given by the center LatLon and radius? There is the trilaterate functions which is similar, could it be (ab)used in some way to get the result?
Thank you.

Typo in pygeodesy/sphericalTrigonometry.py

(With thanks for the work:)

"heighth" on line 985 makes intersection() crash with:

 __init__() got an unexpected keyword argument 'heighth'

Trilateration issue

Hi, I am facing to a trilateration exercise, i'm doing this:

from pygeodesy.vector3d import trilaterate3d2, Vector3d
p1 = Vector3d(2, 2, 0)
p2 = Vector3d(3, 3, 0)
p3 = Vector3d(1, 4, 0)
r1 = 1
r2 = 1
r3 = 1.4142
trilaterate3d2(p1, r1, p2, r2, p3, r3)

When i call to trilaterate3d2 i get: pygeodesy.errors.IntersectionError: center1 (Vector3d(2, 2, 0)), center2 (Vector3d(3, 3, 0)), center3 (Vector3d(1, 4, 0)), radius1 (1), radius2 (1) or radius3 (1.4142): no intersection
But it should return P=(2,3,0)
What i'm doing wrong?