Using GeoPandas intersections to fit a grid on a geometry

Question

I'm trying to write a function that creates a grid over a vector shapefile using GeoPandas:

def grid(src_vector, out_vector, size):

    # read the vector data in mercator proj
    gdf = gpd.read_file(src_vector).to_crs('EPSG:3857')
    
    # extract bounds from gdf 
    min_lon, min_lat, max_lon, max_lat = gdf.total_bounds
    
    # compute the longitudes and latitudes top left corner coordinates
    longitudes = np.arange(min_lon, max_lon, size)
    latitudes = np.arange(min_lat, max_lat, size)
    
    # create the grid centroids
    points = []
    for coords in product(longitudes, latitudes):
        points.append(Point(coords[0], coords[1]))

    # create a buffer grid in lat-long
    grid = gpd.GeoDataFrame({'geometry':points}, crs='EPSG:3857') \
        .buffer(size, cap_style=3) \
        .intersection(gdf) \
        .to_crs('EPSG:4326')
    
    # filter empty geometries
    grid = grid[np.invert(grid.is_empty)]
    
    # export as shapefile
    grid.to_file(out_vector)
    
    return

In order to test it I created a fake shape:

nb_cell = 6
grid_size = 30

gdf = gpd.GeoDataFrame({'geometry': [sg.Point(0,0)]}, crs='EPSG:3857').buffer(nb_cell*grid_size, cap_style=3)
dst = tmp_dir.joinpath('shapes.shp')
gdf.to_file(dst)

but when I launch my function grid(dst, 'output.shp', grid_size) I get the following error :

/usr/local/lib/python3.6/dist-packages/geopandas/base.py:39: UserWarning: The indices of the two GeoSeries are different.

Does someone see what I have done wrong?

edit

The objective is to create grid over country geometries. These geometries can include islands so multiple polygons. The idea is to create a grid on the full extends of the initial gdf and then cut all the grid cells that goes outside the country boundaries.

Based on the documentation, intersection seems to be the way to go

Answer 1

You need to use clip not intersection.

def grid(src_vector, out_vector, size):

    # read the vector data in mercator proj
    gdf = gpd.read_file(src_vector).to_crs('EPSG:3857')
    
    # extract bounds from gdf 
    min_lon, min_lat, max_lon, max_lat = gdf.total_bounds
    
    # compute the longitudes and latitudes top left corner coordinates
    longitudes = np.arange(min_lon, max_lon, size)
    latitudes = np.arange(min_lat, max_lat, size)
    
    # create the grid centroids
    points = []
    for coords in product(longitudes, latitudes):
        points.append(Point(coords[0], coords[1]))

    # create a buffer grid in lat-long
    grid = gpd.GeoDataFrame({'geometry':points}, crs='EPSG:3857') \
        .buffer(size, cap_style=3)

    # clip to geometries
    grid = gpd.clip(grid, gdf).to_crs('EPSG:4326')
    
    # export as shapefile
    grid.to_file(out_vector)
    
    return