Unit of GeoPandas distance function after determining closest coastal line of a point

I am using this (I think working) code to find the distance of a lat/long point to some open source coast lines.

``````import os
import geopandas as gpd
import pandas as pd

# force CRS but I think it is already in this format
lines.to_crs("EPSG:4326")
lines.plot()

points_df = pd.DataFrame({'Latitude': [57.58125], 'Longitude': [-3.98848]}) # taken from lines data - so produced distance should be about 0

points = gpd.GeoDataFrame(points_df, geometry=gpd.points_from_xy(points_df.Longitude, points_df.Latitude, crs="EPSG:4326"))

for index, row in lines.iterrows():
lines.at[index, 'distance'] = row['geometry'].distance(points.iloc[0]['geometry'])

lines = lines.sort_values(by=['distance'], ascending=True)
lines = lines.head(1) # minimum distance
print(lines)
lines.plot()
``````

I ensured that the data is in the same coordinate reference system. Despite lots of Google searches I am still not sure about the unit of GeoPandas' distance function. Is it degrees? Could one get it in km or meters without having to resort to using:

``````from math import radians, cos, sin, asin, sqrt

def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance in kilometers between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])

# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a))
r = 6371 # Radius of earth in kilometers. Use 3956 for miles. Determines return value units.
return c * r
``````
• Are you trying to find the distance to nearest coast?
– Bera
Commented Aug 4, 2022 at 7:20
• yes - in other words given a land point I want to estimate distance to sea (as crow flies is good enough). Commented Aug 4, 2022 at 7:45
• EPSG:4326 has length units in degrees, so the distance units are degrees. Your question is essentially how to get geopandas to calculate geodesic distances. Commented Aug 4, 2022 at 7:49
• You can also reproject the coastline to an aeqd projection centred around the point of interest, which can be done on-the-fly with `lines.to_crs(f"+proj=aeqd +lat_0={pt.y} +lon_0={pt.x} +x_0=0 +y_0=0")` (assuming `pt = Point(-3.98848, 57.58125)`) Commented Aug 4, 2022 at 8:24
• @MikeT - you put it very eloquently: "calculate geodesic distances" and also your other comment. would you mind posting an answer pls? Commented Aug 4, 2022 at 8:32

Reproject you data to a coordinate system with units in meters.

``````import geopandas as gpd

def closest_coastline(input_coast_df, input_point):
"""Finds the closest distance from one input point to the lines in an input line dataframe"""
dist = round( min([input_point.distance(coastline) for coastline in input_coast_df.geometry])/1000, 0)
return dist

#closest_coastline(coast, points.at[0, 'geometry'])
#Out[55]: 16.0
``````

• You should not use 3857 for any distance measurement purposes, ever, for any reason. It is always wrong, and increases in wrongness as you get further from the Equator. Commented Aug 4, 2022 at 12:20
• @vince what else to use please? Commented Aug 4, 2022 at 12:23
• it is mainly UK/Europe Commented Aug 4, 2022 at 12:38
• posted an answer - happy to delete - the code uses estimate_utm_crs and appears to produce sensible results in meters ... Commented Aug 4, 2022 at 12:39

I found this in the meantime and adapted my somewhat reproducible code (see below). It appears to produce sensible results. What do you experts think?

``````import os
import geopandas as gpd
import pandas as pd

# force CRS but I think it is already in this format

lines.plot()

points_df = pd.DataFrame({'Latitude': [57.58125], 'Longitude': [-3.98848]}) # taken from lines data - so produced distance should be about 0
#points_df = pd.DataFrame({'Latitude': [51.509865],'Longitude': [-0.118092]}) # London

points = gpd.GeoDataFrame(points_df, geometry=gpd.points_from_xy(points_df.Longitude, points_df.Latitude, crs="epsg:4326"))

lines.to_crs("epsg:4326") #"EPSG:4326"

lines = lines.to_crs(lines.estimate_utm_crs())
points = points.to_crs(lines.crs)

for index, row in lines.iterrows():
lines.at[index, 'distance'] = row['geometry'].distance(points.iloc[0]['geometry'])

lines = lines.sort_values(by=['distance'], ascending=True)
lines = lines.head(1) # minimum distance
print(lines)
lines.plot()
``````