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

Question

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

os.chdir(r"C:\Users\gis_newbie\Downloads\natural_earth_vector\50m_physical")

# shapefile part of download from here https://www.naturalearthdata.com/downloads/
lines=gpd.read_file('ne_50m_coastline.shp')
# 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

Answer 1

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

import geopandas as gpd

coast = gpd.read_file(r'/home/bera/Desktop/GIStest/sweden_outline_wgs84_4326.shp').to_crs(3006) #Replace 3006 with a crs fitting your needs
points = gpd.read_file(r'/home/bera/Desktop/GIStest/sweden_1k_random_points_4326.shp').to_crs(3006)

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

Answer 2

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

os.chdir(r"C:\Users\gis_newbie\Downloads\natural_earth_vector\50m_physical")

# shapefile part of download from here https://www.naturalearthdata.com/downloads/
lines=gpd.read_file('ne_50m_coastline.shp')
# 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()