Encontre o ponto mais próximo da linha costeira do shapefile em Python

Question

Eu tenho um xarray (674 lats e 488 Lons) e desejo encontrar a distância mais próxima de cada ponto até a costa em metros.

Eu encontrei esta solução: Encontrar o ponto mais próximo do shapefile litoral Python

que é basicamente o que eu quero fazer. No entanto, a distância é medida em graus e não em metros ( veja aqui ).

Eu poderia converter graus em metros usando 1graus = 111km, mas isso não seria muito preciso para domínios maiores e domínios mais ao sul.

Meu exemplo de trabalho está abaixo:

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

Meu palpite é substituir o point['geometry'].distance(coastline)por algo que use a função haversine, mas não tenho ideia de como fazer isso, especialmente algo meio eficiente.

Answer 1

Você pode usar o pacote haversine , é muito fácil de usar. De sua documentação:

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

então para o que você precisa, você precisa:

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

Answer 2

Eu encontrei uma solução razoavelmente rápida combinando as respostas em https://stackoverflow.com/questions/44681828/efficient-computation-of-minimum-of-haversine-distances

e

Encontrando o ponto mais próximo do shapefile litoral Python

O código que funciona agora se parece com este:

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

Espero que isso possa ajudar alguém no futuro!