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...')
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
#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 = gpd.read_file(c)
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.