Find closest point to shapefile coastline in Python

I have a xarray (674 lats & 488 Lons) and want to find the closest distance for each point to the coastline in meters.

I found this solution: Finding closest point to shapefile coastline Python

which is basically what I want to do. However, the distance is measured in degrees and not in meters (see here).

I could convert degrees to meter by using 1deg=111km but this would not be very accurate for larger domains and domains further south.

My working example is below:

``````import geopandas as gpd
from shapely.geometry import Point, box
from random import uniform
from concurrent.futures import ThreadPoolExecutor
from tqdm.notebook import tqdm
import cartopy
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import pandas as pd

lon = np.arange(129.4, 153.75+0.05, 0.05)
lat = np.arange(-43.75, -10.1+0.05, 0.05)

precip = 10 * np.random.rand(len(lat), len(lon))

ds = xr.Dataset({"precip": (["lat", "lon"], precip)},coords={"lon": lon,"lat": lat})

ds['precip'].plot()

def get_distance_to_coast(arr):

def compute_distance(point):
point['dist_to_coastline'] = point['geometry'].distance(coastline)
return point

print('Get shape file...')

#single geom for Norway
aus = world[world["name"]=="Australia"].dissolve(by='name').iloc[0].geometry

#single geom for the coastline
c = cartopy.io.shapereader.natural_earth(resolution='50m', category='physical', name='coastline')

c.crs = 'EPSG:4326'

print('Get coastline...')
coastline = gpd.clip(c.to_crs('EPSG:4326'), aus.buffer(0.25)).iloc[0].geometry

print('Group lat/lon points...')
points = []
i = 0
for ilat in arr['lat']:
for ilon in arr['lon']:
points.append({'id':i, 'geometry':Point(ilon,ilat)})
i+=1

print('Computing distances...')
with ThreadPoolExecutor(max_workers=4) as tpe:
result = list(tqdm(tpe.map(compute_distance, points), desc="computing distances", total=len(points)))

gdf = gpd.GeoDataFrame.from_records(result)

print('Convert to xarray...')
lon = gdf['geometry'].x
lat = gdf['geometry'].y
df1 = pd.DataFrame(gdf)
df1['lat'] = lat
df1['lon'] = lon
df1 = df1.drop(columns=['id','geometry'])
df1 = df1.set_index(['lat', 'lon'])
xarr = df1.to_xarray()

return xarr

dist = get_distance_to_coast(ds['precip'])

plt.figure()
dist['dist_to_coastline'].plot()
plt.show()
``````

My guess is to replace the `point['geometry'].distance(coastline)` with something using the haversine function, but I have no idea how do to this, especially something halfway efficient.

You could use the haversine package, its quite easy to use. From their documentation:

``````from haversine import haversine, Unit
lyon = (45.7597, 4.8422) # (lat, lon)
paris = (48.8567, 2.3508)
haversine(lyon, paris) # in kilometers
``````

so for what you want you would need:

``````haversine(lyon, paris, unit=Unit.METERS) # in meters
``````
• Thank you for pointing out this package. The problem with this is that it only takes a single point. I need the minimum distance of a point to the coastline (which is a list). I modified my function using the haversine_vector function explained on the website. It is currently running and definitively a lot slower than before. – drcrisp Aug 21 '20 at 7:20
• It would make sense for it to be somewhat slower since the haversine formula is more computationally expensive than cartesian distance – Louis Cottereau Aug 21 '20 at 8:33

I found a reasonably fast solution combining the answers in https://stackoverflow.com/questions/44681828/efficient-computation-of-minimum-of-haversine-distances

and

Finding closest point to shapefile coastline Python

The code that works now looks like this:

``````import geopandas as gpd
from shapely.geometry import Point, box
from random import uniform
from concurrent.futures import ThreadPoolExecutor
from tqdm.notebook import tqdm
import cartopy
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import pandas as pd
import shapely

lon = np.arange(129.4, 153.75+0.05, 0.25)
lat = np.arange(-43.75, -10.1+0.05, 0.25)

precip = 10 * np.random.rand(len(lat), len(lon))

ds = xr.Dataset({"precip": (["lat", "lon"], precip)},coords={"lon": lon,"lat": lat})

ds['precip'].plot()

def hv(lonlat1, lonlat2):
AVG_EARTH_RADIUS = 6371000. # Earth radius in meter

# Get array data; convert to radians to simulate 'map(radians,...)' part

# Get the differentiations
lat = coords_arr[:,1] - a[:,1,None]
lng = coords_arr[:,0] - a[:,0,None]

# Compute the "cos(lat1) * cos(lat2) * sin(lng * 0.5) ** 2" part.
# Add into "sin(lat * 0.5) ** 2" part.
add0 = np.cos(a[:,1,None])*np.cos(coords_arr[:,1])* np.sin(lng * 0.5) ** 2
d = np.sin(lat * 0.5) ** 2 +  add0

# Get h and assign into dataframe
h = 2 * AVG_EARTH_RADIUS * np.arcsin(np.sqrt(d))
return {'dist_to_coastline': h.min(1), 'lonlat':lonlat2}

def get_distance_to_coast(arr, country, resolution='50m'):

print('Get shape file...')

#single geom for country
geom = world[world["name"]==country].dissolve(by='name').iloc[0].geometry

#single geom for the coastline
c = cartopy.io.shapereader.natural_earth(resolution=resolution, category='physical', name='coastline')

c.crs = 'EPSG:4326'

print('Group lat/lon points...')
points = []
i = 0
for ilat in arr['lat'].values:
for ilon in arr['lon'].values:
points.append([ilon, ilat])
i+=1

xlist = []
gdpclip = gpd.clip(c.to_crs('EPSG:4326'), geom.buffer(1))
for icoast in range(len(gdpclip)):
print('Get coastline ({}/{})...'.format(icoast+1, len(gdpclip)))
coastline = gdpclip.iloc[icoast].geometry #< This is a linestring

if type(coastline) is shapely.geometry.linestring.LineString:
coastline = [list(i) for i in coastline.coords]
elif type(coastline) is shapely.geometry.multilinestring.MultiLineString:
dummy = []
for line in coastline:
dummy.extend([list(i) for i in line.coords])
coastline = dummy
else:
print('In function: get_distance_to_coast')
exit()

print('Computing distances...')
result = hv(coastline, points)

print('Convert to xarray...')
gdf = gpd.GeoDataFrame.from_records(result)
lon = [i[0] for i in gdf['lonlat']]
lat = [i[1] for i in gdf['lonlat']]
df1 = pd.DataFrame(gdf)
df1['lat'] = lat
df1['lon'] = lon
df1 = df1.set_index(['lat', 'lon'])
xlist.append(df1.to_xarray())

xarr = xr.concat(xlist, dim='icoast').min('icoast')
xarr = xarr.drop('lonlat')

return xr.merge([arr, xarr])

dist = get_distance_to_coast(ds['precip'], 'Australia')

plt.figure()
dist['dist_to_coastline'].plot()
plt.show()
``````

I hope this might help someone in the future!