Encuentra el punto más cercano a la costa de shapefile en Python
Tengo una matriz x (674 lats y 488 Lons) y quiero encontrar la distancia más cercana para cada punto a la línea de costa en metros.
Encontré esta solución: Encontrar el punto más cercano al shapefile Coastline Python
que es básicamente lo que quiero hacer. Sin embargo, la distancia se mide en grados y no en metros ( ver aquí ).
Podría convertir grados a metros usando 1deg = 111km, pero esto no sería muy preciso para dominios más grandes y dominios más al sur.
Mi ejemplo de trabajo está a continuación:
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()
Supongo que reemplazaré point['geometry'].distance(coastline)
con algo que use la función haversine, pero no tengo idea de cómo hacerlo, especialmente con algo medio eficiente.
Respuestas
Puede usar el paquete haversine , es bastante fácil de usar. De su documentación:
from haversine import haversine, Unit
lyon = (45.7597, 4.8422) # (lat, lon)
paris = (48.8567, 2.3508)
haversine(lyon, paris) # in kilometers
así que para lo que quieres necesitarías:
haversine(lyon, paris, unit=Unit.METERS) # in meters
Encontré una solución razonablemente rápida combinando las respuestas en https://stackoverflow.com/questions/44681828/efficient-computation-of-minimum-of-haversine-distances
y
Encontrar el punto más cercano al shapefile Coastline Python
El código que funciona ahora se ve así:
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
coords_arr = np.deg2rad(lonlat1)
a = np.deg2rad(lonlat2)
# 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...')
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
#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 = gpd.read_file(c)
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')
print('Type: {} not found'.format(type(type(coastline))))
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()
¡Espero que esto pueda ayudar a alguien en el futuro!