[GIS] Mask rasterio raster with geopandas shapefile

geopandaspythonrasterrasterioshapefile

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)
folium.GeoJson(StudyA).add_to(m)
m

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]

Best Answer

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)

Related Solutions

[GIS] Sharing a Folium Interactive Map w/ Embedded Markers

I haven't gotten json files to load with nbviewer yet, but here's two things that work in similar situations:

If you only use map.simple_marker() and similar functions, it's possible to completely embed the map HTML into the notebook. I haven't gotten this to work with json yet, but it looks like the upcoming version of follium allows templates, so it possibly can be inlined in the future. I have an example of this here:

http://nbviewer.ipython.org/gist/bburky/7763555/folium-ipython.ipynb

Also, you're running a SimpleHTTPServer but you don't need to. You can use the magic IPython files/ url to refer to files in the directory that IPython was launched from. So you can use an iframe with a source files/index.html. However, I haven't been able to get that to work in the Notebook Viewer. Theoretically you should be able to include multiple files in the gist and make it work.

[GIS] Plot shapefile on top of raster using plot and imshow from matplotlib

I added that recipe to the rasterio documentation. Since it was such a simple shape, in this case I just unzipped the coords in the single record contained by the shapefile. That is, x, y = zip(*features[0]['coordinates'][0]), and then just plot.

More generally, I use PolygonPatch from descartes, and matplotlib.collections.

import fiona
import rasterio
import rasterio.plot
import matplotlib as mpl
from descartes import PolygonPatch

src = rasterio.open("tests/data/RGB.byte.tif")

with fiona.open("tests/data/box.shp", "r") as shapefile:
    features = [feature["geometry"] for feature in shapefile]

rasterio.plot.show((src, 1))
ax = mpl.pyplot.gca()

patches = [PolygonPatch(feature) for feature in features]
ax.add_collection(mpl.collections.PatchCollection(patches))

The appearance of the shapes can be customized with keywords like edgecolor or facecolor passed to PolygonPatch. To produce a thick red line as in the example, replace the last two lines in the example above with:

patches = [PolygonPatch(feature, edgecolor="red", facecolor="none", linewidth=2) for feature in features]
ax.add_collection(mpl.collections.PatchCollection(patches, match_original=True))

The match_original keyword is necessary in the second example because parameters like facecolor and edgecolor can also be set in PatchCollection, and the default is to ignore settings of the passed in patches in favor of those set by PatchCollection. Doing so using PolygonPatch gives more flexibility to set different colors, widths, etc., for each patch that you add.

Best Answer

Related Solutions

[GIS] Sharing a Folium Interactive Map w/ Embedded Markers

[GIS] Plot shapefile on top of raster using plot and imshow from matplotlib

Related Question