Skip to content

Instantly share code, notes, and snippets.

@deeplycloudy
Created November 19, 2021 23:22
Show Gist options
  • Save deeplycloudy/54bb5a43939083ccb8d3dd690e6f9170 to your computer and use it in GitHub Desktop.
Save deeplycloudy/54bb5a43939083ccb8d3dd690e6f9170 to your computer and use it in GitHub Desktop.
Distance from Cartopy feature polygon
import cartopy.feature as cfeature
import pyproj
from shapely.geometry import Point
from shapely.ops import transform as shapely_transform
# Some location of interest. Also used as center of projection.
wgs84_pt = Point(-72.2495, 43.886)
# Set up projections
wgs84 = pyproj.Proj(proj='latlon')
aeqd = pyproj.Proj(proj='aeqd', lon_0=wgs84_pt.x, lat_0=wgs84_pt.y)
project = pyproj.Transformer.from_proj(wgs84, aeqd, always_xy=True).transform
aeqd_point = shapely_transform(project, wgs84_pt)
# Loop through shape components, searching only for polygons,
# and generate Shapely polygons in map projection
def gen_transformed_polys(shape):
for geom in shape.geometries():
typestr = str(type(geom))
if 'Polygon' in typestr:
if 'Multi' in typestr:
for poly in geom:
yield shapely_transform(project, poly)
else:
yield shapely_transform(project, geom)
# Test with the Natural Earth land polygons
land = cfeature.NaturalEarthFeature('physical', 'land', '50m',
edgecolor='face',
facecolor=cfeature.COLORS['land'])
land_distances_km = [int(one_shape.exterior.distance(aeqd_point)/1000.0)
for one_shape in gen_transformed_polys(land)]
print(min(land_distances_km), max(land_distances_km))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment