@@ -13,15 +13,27 @@ def utm_for_geometry(lat: float, lng: float) -> int:
1313 epsg_code = 32600 + utm_band
1414 return epsg_code
1515
16+ def expand_recurse_geoms (shape : shapely .Geometry , radius : float ) -> shapely .Geometry :
17+ # Calling .buffer directly on a shape is incredibly slow when it contains lots of geometry.
18+ # So we call it on each piece of geometry instead. No reason that should be faster, but it is
19+ if hasattr (shape , 'geoms' ):
20+ expanded = shapely .GeometryCollection ([])
21+ for geom in shape .geoms :
22+ expand_geom = expand_recurse_geoms (geom , radius )
23+ expanded = expanded .union (expand_geom )
24+ return expanded
25+ else :
26+ return shape .buffer (radius )
27+
1628def expand_boundaries (source : gpd .GeoDataFrame , radius : float ) -> gpd .GeoDataFrame :
1729 original_projection = source .crs
1830
1931 utm_codes = np .array ([utm_for_geometry (x .centroid .y , x .centroid .x ) for x in source .geometry ])
2032 utm_code = scipy .stats .mode (utm_codes , keepdims = False ).mode
2133 projected_boundaries = source .to_crs (f"EPSG:{ utm_code } " )
2234
23- expanded = [ x . buffer ( radius ) for x in projected_boundaries .geometry ]
24- simplified = shapely . unary_union ( expanded )
35+ project = projected_boundaries .unary_union
36+ simplified = expand_recurse_geoms ( project , radius )
2537
2638 finished = gpd .GeoSeries (simplified )
2739 utm_result = gpd .GeoDataFrame .from_features (finished , crs = f"EPSG:{ utm_code } " )
0 commit comments