4

Say I have a GeoPandas GeoDataFrame of polygons (i.e. the "geometry" column contains polygons of GPS coordinates only).

import geopandas as gpd
from shapely.geometry import Polygon

lat_point_list = [50.854457, 52.518172, 50.072651, 48.853033, 50.854457]
lon_point_list = [4.377184, 13.407759, 14.435935, 2.349553, 4.377184]

polygon_geom = Polygon(zip(lon_point_list, lat_point_list))
polygon = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[polygon_geom]) 


#Visualizing this polygon
import folium
m = folium.Map([50.854457, 4.377184], zoom_start=5, tiles='cartodbpositron')
folium.GeoJson(polygon).add_to(m)
folium.LatLngPopup().add_to(m)
m

enter image description here

I want to compute a new GeoDataFrame with a new "geometry" column of polygons where the new polygons are parallel to the old ones with a given perpendicular distance d.

If d is positive, the new polygons are larger than the old ones and encompass them. If d is negative, the new polygons are smaller than the old ones and are contained within them.

d = 100.  #let's assume the unit of kilometres

#construct new bigger polygons that have all sides parallel to and 100 kilometres outside of the above smaller polygon


d = -50
#construct new smaller polygons that have all sides parallel to and 50 kilometres inside of the above bigger polygon

How would I do this efficiently, with a solution generalized for any kind of polygon?

0

1 Answer 1

5

You can use buffer method with the join_style=2 option (2 refers mitre) for visualization purposes.

import geopandas as gpd
from math import pi
from shapely.geometry import Polygon

lat_point_list = [50.854457, 52.518172, 50.072651, 48.853033, 50.854457]
lon_point_list = [4.377184, 13.407759, 14.435935, 2.349553, 4.377184]

polygon_geom = Polygon(zip(lon_point_list, lat_point_list))
polygon_df = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[polygon_geom]) 

# positive distance
d = 100 # km
d = d * 180 / (pi * 6371) # approx. degree equivalent on Earth sphere
polygon1_df = polygon_df.buffer(d, join_style=2)

# negative distance
d = -30 # km
d = d * 180 / (pi * 6371) # approx. degree equivalent on Earth sphere
polygon2_df = polygon_df.buffer(d, join_style=2)

#Visualizing this polygon
import folium
m = folium.Map([50.854457, 4.377184], zoom_start=5, tiles='cartodbpositron')

folium.GeoJson(polygon_df).add_to(m)

style_function = lambda x: {"fillOpacity": 0, "color": "#ff0000"}
folium.GeoJson(polygon1_df, style_function).add_to(m)
folium.GeoJson(polygon2_df, style_function).add_to(m)

folium.LatLngPopup().add_to(m)
m

The result looks a bit incorrect because of the transformation from EPSG:4326 (polygon's CRS) to EPSG 3857 (basemap's CRS).

enter image description here

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.