Найдите ближайшую точку к береговой линии шейп-файла в Python

Aug 21 2020

У меня есть рентгеновский снимок (674 латов и 488 латов), и я хочу найти ближайшее расстояние в метрах для каждой точки до береговой линии.

Я нашел это решение: поиск ближайшей точки к береговой линии шейп-файла Python

что в основном то, чем я хочу заниматься. Однако расстояние измеряется в градусах, а не в метрах ( см. Здесь ).

Я мог бы преобразовать градусы в метры, используя 1 градус = 111 км, но это было бы не очень точно для больших областей и областей южнее.

Мой рабочий пример ниже:

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

Вы можете использовать пакет 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

1 drcrisp Aug 24 2020 at 05:11

Я нашел достаточно быстрое решение, объединив ответы в https://stackoverflow.com/questions/44681828/efficient-computation-of-minimum-of-haversine-distances