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.
A simpler way would be to use geopandas and shapely:
import geopandas as gpd
import shapely.geometry as geom
def polyline_to_polygone(input_file, output_file):
gdf = gpd.read_file(input_file)
gdf['geometry'] = gdf['geometry'].apply(lambda x: geom.Polygon(x.coords))
# in order to have donut geometries, subset polygons from each others
for index in list(gdf.index):
indices = list(gdf.index)
indices.remove(index)
for index_ in indices:
try:
gdf.loc[index, 'geometry'] = gdf.loc[index, 'geometry'].difference(gdf.loc[index_, 'geometry'])
except ValueError:
pass
gdf.to_file(output_file)
Best Answer
I'd recommend using gdal_contour . The results are likely to be much better than any attempt to re-implement it :-)
Having said that, there should be a way to do this in rasterio. I've not tried this, so it may behave differently to how I'd expect it to work in QGIS.
Step 1. Quantize the raster into contour bands using rasterio. Use rasterio as a raster calculator. Here's the formula I use for this. This assumes a float raster. G is the elevation gap between contours. E is the existing elevation value). You might need to tweak this.
That will 'posterize' a float DEM into multiples of the contour elevation (0, G, 2G, 3G...), should look like this...
Step 2. Apply the existing rasterio vectorize example to the quantized raster. See this example in github.