Python에서 shapefile 해안선에 가장 가까운 지점 찾기

Question

나는 xarray (674 lats & 488 Lons)를 가지고 있으며 각 지점에서 해안선까지 가장 가까운 거리를 미터 단위로 찾고 싶습니다.

이 해결책을 찾았습니다 : shapefile 해안선 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)엔을 haversine 함수를 사용하는 것으로 대체하는 것이지만, 어떻게 해야할지 모르겠습니다. 특히 중간 정도 효율적인 것입니다.

Answer 1

사용하기 매우 쉬운 haversine 패키지를 사용할 수 있습니다. 그들의 문서에서 :

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

Answer 2

나는 대답을 결합하는 합리적으로 빠른 솔루션을 찾았습니다. https://stackoverflow.com/questions/44681828/efficient-computation-of-minimum-of-haversine-distances

과

shapefile 해안선에 가장 가까운 지점 찾기 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()

나는 이것이 미래에 누군가를 도울 수 있기를 바랍니다!