1

I am trying to calculate the distance to coastline for some points. Using Cartopy and Shapely, I can calculate the distance in degrees, but that is not very helpful. Perhaps if I could transform the geometries, this would not be difficult, but I'm not sure how to do that in Cartopy. I'm trying to wean myself off the depreciated Basemap package. Does anyone know how to do this where the distance is output in meters? I have a dataframe with lat and lon for the points I want to calculate.

import numpy as np
import pandas as pd
import shapely
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature

land = shapereader.gshhs(scale='h', level=1)
                                  
geoms = list(itertools.chain.from_iterable(geom.geoms for geom in shapereader.Reader(land).geometries()))
geometries = shapely.geometry.MultiPolygon(geoms)

df['dist'] = np.nan

for i in df.index:
    point = shapely.geometry.Point(df.iloc[i]['lon'], df.iloc[i]['lat'])
    df['dist'][i] = geometries.distance(point)

print(df)

1 Answer 1

2

I figured it out. I was able to use pyproj to transform the shapely geometries.

land = shapereader.gshhs(scale='h', level=1)
                                  
geoms = list(itertools.chain.from_iterable(geom.geoms for geom in shapereader.Reader(land).geometries()))
geometries = shapely.geometry.MultiPolygon(geoms)

src_crs = pyproj.CRS('EPSG:4326')
tgt_crs = pyproj.CRS('EPSG:32616')

project = pyproj.Transformer.from_crs(src_crs, tgt_crs, always_xy=True).transform
xy_geometries = transform(project, geometries)

df['dist'] = np.nan

for i in df.index:
    xy = transform(project, shapely.geometry.Point(df.iloc[i]['lon'],
                                                   df.iloc[i]['lat']))
    x, y = xy.xy[0][0], xy.xy[1][0]
    point = shapely.geometry.Point(x, y)
    df['dist'][i] = xy_geometries.distance(point)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.