2

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

1
  • .intersection is a row-wise operation which aligns both GeoSeries. What are you trying to do with it? Maybe gpd.clip or gpd.overlay would suit the needs? Jan 12 at 14:20
2

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.