Temukan titik terdekat ke garis pantai shapefile dengan Python
Saya memiliki xarray (674 lats & 488 Lons) dan ingin menemukan jarak terdekat untuk setiap titik ke garis pantai dalam meter.
Saya menemukan solusi ini: Menemukan titik terdekat ke garis pantai shapefile Python
yang pada dasarnya ingin saya lakukan. Namun, jarak tersebut diukur dalam derajat dan bukan dalam meter ( lihat di sini ).
Saya dapat mengonversi derajat ke meter dengan menggunakan 1deg = 111km tetapi ini tidak akan terlalu akurat untuk domain dan domain yang lebih besar di selatan.
Contoh kerja saya di bawah ini:
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()
Dugaan saya adalah mengganti point['geometry'].distance(coastline)
dengan sesuatu yang menggunakan fungsi haversine, tetapi saya tidak tahu bagaimana melakukannya, terutama sesuatu yang setengah efisien.
Jawaban
Anda bisa menggunakan paket haversine , ini cukup mudah digunakan. Dari dokumentasi mereka:
from haversine import haversine, Unit
lyon = (45.7597, 4.8422) # (lat, lon)
paris = (48.8567, 2.3508)
haversine(lyon, paris) # in kilometers
jadi untuk apa yang Anda inginkan, Anda akan membutuhkan:
haversine(lyon, paris, unit=Unit.METERS) # in meters
Saya menemukan solusi yang cukup cepat dengan menggabungkan jawaban dalam https://stackoverflow.com/questions/44681828/efficient-computation-of-minimum-of-haversine-distances
dan
Menemukan titik terdekat ke garis pantai shapefile Python
Kode yang berfungsi sekarang terlihat seperti ini:
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()
Saya harap ini dapat membantu seseorang di masa depan!