Tuesday, 20 December 2016

python - Plot shapefile on top of raster using plot and imshow from matplotlib



I need to plot a shapefile on top of a raster. I would like to use the python package rasterio and some combination of fiona or geopandas to load the shapefile.


I found an example in the rasterio documentation but it doesn't provide code for the actual plotting.


enter image description here


I have tried loading the raster with rasterio, loading the shapefile using fiona, calling imshow on the raster, and calling plot on the shapefile. However, the shapefile does not appear on the plot.



Answer



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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...