I am loading a shapefile using geopandas as follows:

import geopandas as gpd
StudyA = gpd.read_file('.../Study_Area_Polygon.shp')

Next, I am using folium to plot it interactively:

m = folium.Map([36.43, 43.06], zoom_start=12)

My objective is now to use rasterio to mask a raster using the shapefile I have loaded before with the following code:

img, out_transform = rasterio.mask.mask(img, StudyA, crop=True, all_touched=True)

However, I am facing an issue as the information of the mask should be provided as a list of GeoJSON-like dicts according to rasterio documentation. How can I produce that GeoJSON-like dicts when reading the shapefile with geopandas.

Previously, I was reading the shp using the following code and it worked well. The issue is that I want/need to use geopandas for interactive plotting with folium

with fiona.open(".../Study_Area_Polygon.shp", "r") as StudyArea:
    shape = [feature["geometry"] for feature in StudyArea]
  • You can open it twice: once with Fiona to supply a mask, and one with GeoPandas for plotting with Folium. – alphabetasoup Oct 7 '19 at 20:24
  • Indeed, thank you for the suggestion. However, I would like to avoid repetition and importing twice the same file. I would prefer to extract the GeoJSON information required by rasterio directly from the geopandas import – GCGM Oct 8 '19 at 7:04

In case someone is interested in the future, I managed to make it work with the following code.

First, open an empty list. Then import the shp with geopandas. Next, we access the geometric information and extract it. The output of that is a dict containing the information required by rasterio.mask.mask

geom = []

StudyA = gpd.read_file('.../Study_Area_Polygon.shp')
coord = shapely.geometry.mapping(StudyA)["features"][0]["geometry"]

Assuming you have a geodataframe (gdf), you can just pass the geometry values flattened to rasterio:

masked_raster, masked_raster_transform = rasterio.mask.mask(raster, gdf[['geometry']].values.flatten())

If your geodataframe has a column named "geometry" and geometry consists of Shapely Polygons, the following works for me

from rasterio.mask import mask

out, _ = mask(data, gdf.geometry, invert=False)

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.