Overlay cropped raster with vector layer
I recently faced a problem of having to plot “cropped raster” layer and a vector layer on the same axes. It is known that we first need to identify the spatial extent of each layer, having the same coordinate reference system.
Rasterio does offer a plotting function
show which can plot a raster layer with the correct spatial extent for you when we pass the dataset reader object.
When we pass a reader object, the spatial extent is automatically read by
with rs.open(path_to_file, "r") as src: # import rasterio as rs f, ax = plt.subplots(figsize=(9,9)) _ = show(src, ax=ax) # from rasterio.plot import show _ = vector_layer.plot(ax=ax) # `vector_layer` is a geodataframe (geopandas)
Moreover, if we pass a numpy array to the
show function, the spatial extent of that array has to be explicitly passed using the
transform parameter of the
show function since the numpy array does not know the corner location of the raster and thus the plot would begin with x,y: 0,0 as shown below.
with rs.open(path_to_file, "r") as src: img = src.read(1) # img is a numpy array f, ax = plt.subplots(figsize=(9,9)) _ = show(img, transform = src.transform, ax=ax) _ = vector_layer.plot(ax=ax)
But what if you want to plot a subset of the raster image, in the sense that you would like to slice the image arbitrarily and plot it. When you slice the image, the affine transformation is not the same anymore and thus plotting the sliced image would result in a plot having the spatial extent of the original image while the sliced image being magnified (Fig. 2).
with rs.open(path_to_file, "r") as src: img = src.read(1)[1:-1,1:-1] f, ax = plt.subplots(figsize=(9,9)) _ = show(img, transform = src.transform, ax=ax) _ = vector_layer.plot(ax=ax)
To avert this problem, we need to find the new affine transformation of the cropped image. Luckily rasterio has a
window_transform method on the dataset reader which can compute the new transformation from the old one by passing the bounds of the layer. The
window_transform function can either take a 2D N-D array indexer in the form of a tuple
((row_start, row_stop), (col_start, col_stop)) or provide offset as written in its documentation
Cropped raster and vector overlay
The above method returns the new affine transformation, which can be passed to the
show function for the numpy array through the
transform parameter. We also change the read method instead of slicing the array by window parameter to maintain uniformity
# load raster with rs.open(path_to_file, "r") as src: # window = (((row_start), (row_stop)), ((col_start), (col_stop))) img = src.read(1, window = ((1,-1), (1,-1))) f, ax = plt.subplots(figsize=(9,9)) show(img, transform=src.window_transform(((1,-1), (1,-1))), ax=ax) _ = vector_layer.plot(ax=ax)
show method is helpful for plotting rasters or even RGB images for that matter. One of the differences with matplotlib’s plotting is the order of axes.
show expects it the bands to be the last axis while matplotlib, the first. It can also plot 4-band image, which is almost always the for satellite images.
While there is an
extent paramter in matplotlib’s plotting function,
show function is much tidier and straight-forward to implement cropped raster and overlay vector layer on it.