Trova il punto più vicino alla costa shapefile in Python

Question

Ho una radiografia (674 lats e 488 lons) e desidero trovare la distanza più vicina per ogni punto dalla costa in metri.

Ho trovato questa soluzione: trovare il punto più vicino allo shapefile litorale Python

che è fondamentalmente quello che voglio fare. Tuttavia, la distanza è misurata in gradi e non in metri ( vedi qui ).

Potrei convertire i gradi in metri usando 1deg = 111km ma questo non sarebbe molto preciso per domini più grandi e domini più a sud.

Il mio esempio di lavoro è di seguito:

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

La mia ipotesi è di sostituire point['geometry'].distance(coastline)con qualcosa che utilizza la funzione haversine, ma non ho idea di come farlo, specialmente qualcosa di semi efficiente.

Answer 1

Potresti usare il pacchetto haversine , è abbastanza facile da usare. Dalla loro documentazione:

from haversine import haversine, Unit
lyon = (45.7597, 4.8422) # (lat, lon)
paris = (48.8567, 2.3508)
haversine(lyon, paris) # in kilometers

quindi per quello che vuoi avresti bisogno di:

haversine(lyon, paris, unit=Unit.METERS) # in meters

Answer 2

Ho trovato una soluzione ragionevolmente veloce combinando le risposte in https://stackoverflow.com/questions/44681828/efficient-computation-of-minimum-of-haversine-distances

e

Trovare il punto più vicino allo shapefile litorale Python

Il codice che ora funziona è simile a questo:

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

Spero che questo possa aiutare qualcuno in futuro!