Find minimum distance to US coastline in Python

distancegeopandaspythonshapefileshapely

I have the following code in Python:

import geopandas as gpd

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
us = world[world['name']=='United States of America'].dissolve(by='name')
coastline = gpd.clip(gpd.read_file('./data/ne_coastline/ne_10m_coastline.shp').to_crs('EPSG:4326'), us)
coastline.distance(Point(30.121735, -81.735374)).min()

# 153.78644123108657

But upon manually examination in Google Maps, this isn't correct. To simplify the problem, I found the coordinate approximating the closest coastline point (30.135054573782618, -81.34892671159952).

import geopandas as gpd
from shapely.geometry import Point

points_df = gpd.GeoDataFrame({'geometry': [Point(30.135054573782618, -81.34892671159952), Point(30.121735, -81.735374)]}, crs='EPSG:4326')
points_df = points_df.to_crs('EPSG:3087') # https://epsg.io/3087
points_df2 = points_df.shift() # shift the dataframe by 1 to align points
points_df.distance(points_df2)

# 0            NaN
# 1    4872.316273
# dtype: float64

But this is also far off from the correct answer in meters (which I double check manually in Google Maps):

from geopy.distance import distance
distance((30.135054573782618, -81.34892671159952), (30.121735, -81.735374)).m

# 37268.01227017112

I imagine my coordinate reference system is incorrect, but I'm unsure why given the reported accuracy of EPSG:3087. How can I correctly output the minimum distance to the coastline in meters given a coordinate?

Best Answer

You're measuring the distance from Antarctica to the US coastline (in decimal degrees in your first example)...

Unless you mean Point(-81.735374, 30.121735)? which is in Florida.

Point(coords) uses x, y ordering so you need Point(longitude, latitude)

import geopandas as gpd
from shapely.geometry import Point

points_df = gpd.GeoDataFrame({
    'geometry': [
        Point(-81.34892671159952, 30.135054573782618), 
        Point(-81.735374, 30.121735)]}, crs='EPSG:4326')

points_df = points_df.to_crs('EPSG:3087') # https://epsg.io/3087
points_df2 = points_df.shift() # shift the dataframe by 1 to align points

points_df.distance(points_df2)

Output:

0             NaN
1    37219.622729
dtype: float64

And based on this answer, you could do something like this:

import geopandas as gpd
from shapely.geometry import Point

def min_distance(point, lines):
    return lines.distance(point).min()

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
us = world[world['name']=='United States of America'].dissolve(by='name')
coastline = gpd.clip(gpd.read_file('ne_10m_coastline.shp'), us).to_crs('EPSG:3087')
points_df = gpd.GeoDataFrame({
    'geometry': [
        Point(-81.34892671159952, 30.135054573782618), 
        Point(-81.735374, 30.121735)]}, crs='EPSG:4326')
points_df = points_df.to_crs('EPSG:3087') # https://epsg.io/3087

points_df['min_dist_to_coast'] = points_df.geometry.apply(min_distance, args=(coastline,))
Related Question