Split polygon

split_me.py

import geopandas as gpd
from shapely.geometry import box
from typing import Union

def split_polygon_into_squares(
    gdf: gpd.GeoDataFrame,
    square_size_km: float,
    epsg_code: int
) -> gpd.GeoDataFrame:
    """
    Split a GeoDataFrame polygon into squares of a specified size in kilometers.
    Useful because you can do `.bounds` on the output to get bounding boxes that
    overlap.

    Args:
        gdf (gpd.GeoDataFrame): Input GeoDataFrame containing a single polygon feature.
        square_size_km (float): Size of the squares in kilometers.
        epsg_code (int): EPSG code for the desired projected coordinate system.
        
    Returns:
        gpd.GeoDataFrame: GeoDataFrame containing the split squares.
    """
    # Reproject the GeoDataFrame to the specified EPSG code
    gdf_reprojected = gdf.to_crs(epsg=epsg_code)
    
    # Extract the polygon from the reprojected GeoDataFrame
    polygon = gdf_reprojected.geometry.iloc[0]
    
    # Calculate the size of the squares in the reprojected coordinate system
    square_size_meters = square_size_km * 1000  # Convert square size from km to meters
    
    # Split the polygon into squares
    squares = []
    minx, miny, maxx, maxy = polygon.bounds
    for x in range(int(minx), int(maxx), int(square_size_meters)):
        for y in range(int(miny), int(maxy), int(square_size_meters)):
            square = box(x, y, x + square_size_meters, y + square_size_meters)
            # Check if the square intersects with the original polygon
            if square.intersects(polygon):
                # If there is an intersection, take the intersection area as the square
                square_intersection = square.intersection(polygon)
                squares.append(square_intersection)
    
    # Create a GeoDataFrame for the resulting squares
    squares_gdf = gpd.GeoDataFrame(geometry=squares, crs=epsg)
    
    # Reproject the squares back to lat/lon
    return squares_gdf.to_crs(4326)