# Finding distance between two points with GeoSeries.distance

When I want the distance between two points `[(117.454361,38.8459879),(117.459880 ,38.846255)]` (longitude,latitude) on the earth, I take the `GeoSeries.distance` method, but the method does not give me the right answer.

The real distance is about 479 meters, however the method give the result far from it, why?

``````import geopandas as gpd
from shapely.geometry import Point

geom = [Point(xy) for xy in zip([117.454361,117.459880],[38.8459879,38.846255])]
gdf = gpd.GeoDataFrame(geometry=geom, crs={'init':'epsg:4326'})
dis = gdf.distance(gdf.shift())
print(dis)
``````

``````import geopandas as gpd
from shapely.geometry import Point

geom = [Point(xy) for xy in zip([117.454361,117.459880], [38.8459879,38.846255])]
gdf = gpd.GeoDataFrame(geometry=geom, crs={'init':'epsg:4326'})
gdf.to_crs(epsg=3310, inplace=True)
l = gdf.distance(gdf.shift())
print(l)
``````

The result is 479.450134 meters.

• What is crs 3310? Why are you using it? Commented Jun 3, 2019 at 14:33
• I believe it's a conversion to UTM, a rectilinear projection coordinate system where distance can be calculated with pythagorean theorem. But I don't know why it's EPSG 3310 in particular. Is that a general solution for anywhere on Earth? Commented Jul 31, 2019 at 20:43
• EPSG 32663 will also be a good candidate (WGS 84 / World Equidistant Cylindrical) spatialreference.org/ref/epsg/32663 Commented Aug 11, 2020 at 11:42

This is a distance in degrees, the coordinates of your data. I can get this using Pythagoras' theorem from your coordinates:

``````>>> a = [117.454361, 38.8459879]
>>> b = [117.45988, 38.846255]
>>> math.sqrt((a[0]-b[0])**2 + (a[1]-b[1])**2)
0.005525459565494833
``````

to get a distance in metres, transform to a metric coordinate system, such as the UTM zone for your data.

Possible help here:

Calculate distance between a coordinate and a county in GeoPandas

Its pretty easy with Geopy

``````import geopy.distance
dist = geopy.distance.geodesic((38.8459879,117.454361),(38.846255,117.459880))
dist.meters
``````

• Are the coordinates indise `(lat, lon)` or `(lon, lat)`. Thanks!
– M.O.
Commented Apr 22, 2020 at 3:59
• its (lat, long) Commented Apr 22, 2020 at 14:47

Pythagoras only works on a flat plane and not an sphere. The distance between two points on the surface of a sphere is found using great-circle distance:

where φ's are latitude and λ's are longitudes. To convert the distance to meter you need to know the radius of the sphere (6371km for Earth) and multiply it by Δσ in radians. Here is a code that does that:

``````def haversine(coord1, coord2):
import math
# Coordinates in decimal degrees (e.g. 2.89078, 12.79797)
lon1, lat1 = coord1
lon2, lat2 = coord2
R = 6371000  # radius of Earth in meters

a = math.sin(delta_phi / 2.0) ** 2 + math.cos(phi_1) * math.cos(phi_2) * math.sin(delta_lambda / 2.0) ** 2

c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))

meters = R * c  # output distance in meters
km = meters / 1000.0  # output distance in kilometers

meters = round(meters)
km = round(km, 3)
print(f"Distance: {meters} m")
print(f"Distance: {km} km")
``````

and run it for the coordinates you get 479m.

``````haversine(coord1= (117.454361,38.8459879), coord2= (117.459880,38.846255))
``````

The other way to this is to do what @CJ Xu did and convert the coordinates to UTM.

• There is a library to do that in a vectorized way (with the possibility of using numba optimizations) : - github.com/mapado/haversine Commented May 8, 2023 at 22:26

It's a little unfortunate that this isn't built into GeoPandas. But you can ask Pyproj to do the geodesic calculations for you:

``````import geopandas as gpd
import pyproj

df = gpd.GeoDataFrame({
"lon": [117.454361, 117.459880, 117.462374],
"lat": [38.8459879, 38.846255, 38.847381],
})

df["pt"] = gpd.points_from_xy(df["lon"], df["lat"], crs="epsg:4326")
df.set_geometry("pt", inplace=True)

p1 = df["pt"].iloc[1:]
p2 = df["pt"].shift().dropna()

geod = df.crs.get_geod()
dist = geod.inv(
p1.x.to_numpy(),
p1.y.to_numpy(),
p2.x.to_numpy(),
p2.y.to_numpy(),
)[2]

print(dist)
``````

Result:

``````[480.04157936 250.0043709 ]
``````

Using only pyproj, calculate the length of a geodesic from lat/long coords:

``````from pyproj import Geod

geod = Geod(ellps="WGS84")

lons = [117.454361, 117.459880]
lats = [38.8459879, 38.846255]

total_length = geod.line_length(lons, lats)
print(f"total_length = {total_length:.3f} m")
# total_length = 480.042 m
``````