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()

これが将来誰かに役立つことを願っています!