Working with geographic coordinates in shapely?

Question

I have a very basic question about working with the shapely module in Python. I have some geographic coordinates received from a GPS. I would like to use shapely to calculate the great cicle distance in meters between two points. I started with:

>>> from shapely.geometry import Point
>>> p1 = Point(43.374880, -78.119956)
>>> p2 = Point(43.374868, -78.119666)

I believe that this gives me two points in the cartesian coordinate system, which isn't going to be very useful. These coordinates come from a GPS which presumably is using the WGS84 CRS. I have seen some example using pyproj and shapely.ops.transform to transform between coordinate systems, but these all seem to involve data points that are already in some valid cartographic coordinate system with an EPSG identifier. I'm not sure what do with naive points like the above that don't have an associated geographic CRS.

How do associate these points with the WGS84 CRS such that p1.distance(p2) will yield the distance in meters? Am I using the wrong tools?

Answer 1

The problem you're trying to solve is the Inverse Geodesic Problem. Happily the Python geographiclib can do this for you:

from geographiclib.geodesic import Geodesic

p1_lat, p1_lon = 43.374880, -78.119956
p2_lat, p2_lon = 43.374868, -78.119666
geod = Geodesic.WGS84

# note the Inverse method expects latitude then longitude
g = geod.Inverse(p1_lat, p1_lon, p2_lat, p2_lon)

print("Distance is {:.2f}m".format(g['s12']))

Giving the result:

Distance is 23.54m

That said, given your points are likely to be only 10s of metres apart, then working in projected units will likely be much easier.

Answer 2

Actually, I think I figured it out. The above two points happen to be on the shores of Lake Ontario, so I can turn the (lat, lon) into a grid reference using EPSG:32117 (NAD83/New York West). I ended up with:

>>> from shapely.geometry import Point
>>> from pyproj import Proj
>>> nys = Proj(init='EPSG:32117')
>>> p1 = Point(43.374880, -78.119956)
>>> p2 = Point(43.374868, -78.119666)
>>> p1_proj = nys(p1.y, p1.x)
>>> p2_proj = nys(p2.y, p2.x)
>>> d = Point(p1_proj).distance(Point(p2_proj))
>>> d
23.539335485921658

Which seems reasonable for my needs (the points I'm using are never going to be more than tens of meters apart).

This is a linear distance rather than a great circle distance, and requires knowing the appropriate CRS for the given coordinates. If there is a more general solution, that would be great.