Working with geographic coordinates in shapely?

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.

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.