Pythonでシェープファイルの海岸線に最も近いポイントを見つける
Aug 21 2020
xarray(674 lats&488 Lons)があり、海岸線までの各ポイントの最も近い距離をメートル単位で見つけたいと思っています。
私はこの解決策を見つけました:シェープファイルの海岸線Pythonに最も近いポイントを見つける
それが基本的に私がやりたいことです。ただし、距離はメートルではなく度で測定されます(ここを参照)。
1deg = 111kmを使用して度をメートルに変換することはできますが、これは、より大きなドメインやさらに南のドメインではあまり正確ではありません。
私の実際の例は以下のとおりです。
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()
私の推測ではpoint['geometry'].distance(coastline)
、をハバーシン関数を使用して何かに置き換えることですが、これをどのように行うのか、特に中途半端な効率の何かがわかりません。
回答
1 LouisCottereau Aug 21 2020 at 06:57
あなたは使用することができます半正矢、その非常に使いやすいパッケージを。彼らのドキュメントから:
from haversine import haversine, Unit
lyon = (45.7597, 4.8422) # (lat, lon)
paris = (48.8567, 2.3508)
haversine(lyon, paris) # in kilometers
だからあなたが望むもののためにあなたは必要になるでしょう:
haversine(lyon, paris, unit=Unit.METERS) # in meters
1 drcrisp Aug 24 2020 at 05:11
私はの答えを組み合わせた適度に速い解決策を見つけました https://stackoverflow.com/questions/44681828/efficient-computation-of-minimum-of-haversine-distances
そして
シェープファイルの海岸線Pythonに最も近いポイントを見つける
現在機能するコードは次のようになります。
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()
これが将来誰かに役立つことを願っています!