Plot polygon on rasterio window

Question

I'm trying to plot several polygons on a raster. Several posts suggest to use matplotlib and descartes like this one: Plot shapefile on top of raster using plot and imshow from matplotlib

My polygons do not appear on the map although they have the same coordinate system and that the polygon coordinates are within the bbox of my raster. The plots below show the raster and the polygons as separate plots (a and b) and one can see that a) has "window coordinates" in pixels and b) in EPSG: 4326.

How can I plot my polygons on the rasterio window?

from descartes import PolygonPatch
import matplotlib as mpl
import rasterio


# Import .tif image (RGB)
src = rasterio.open(base_tif)  # src.crs = CRS.from_epsg(4326)
# Select window based on 4 coordinates
rst = src.read(window=windows.from_bounds(minx, miny, maxx, maxy, src.transform))  # 7.493877446690874, 7.545662625573129, 45.95562555558752, 45.99882089053693, Affine(0.00011273472018020788, 0.0,6.7805624095004235,0.0, -0.00011273472018020788, 46.66875202450472)

# My polygons
lakes_for_gl = [{'type': 'Polygon',
  'coordinates': [[[7.559165, 45.965006],
    [7.558827, 45.965006],
    [7.558827, 45.965118],
    [7.559277, 45.965231],
    [7.559165, 45.965006],
    [7.559165, 45.965006]]]},
 {'type': 'Polygon',
  'coordinates': [[[7.51362, 45.958918],
    [7.513056, 45.958918],
    [7.513056, 45.959256],
    [7.513733, 45.960271],
    [7.514296, 45.960496],
    [7.514409, 45.959594],
    [7.51362, 45.958918],
    [7.51362, 45.958918]]]}]


fig, ax = plt.subplots()
rasterio.plot.show(rst)
ax = plt.gca()

patches = [PolygonPatch(feature, edgecolor="red", facecolor="red", linewidth=2) for feature in lakes_for_gl]
ax.add_collection(mpl.collections.PatchCollection(patches, match_original=True))
ax.axis('scaled')
plt.show()

Answer 1

There's a couple of issues with your code:

1. Notice the values in the x and y axes of your plot "a"? They're pixel coordinates, not map coordinates as you didn't pass a transform to rasterio.plot.show

2. Your polygon isn't plotting on the same ax as rst.

Here's a working example:

from descartes import PolygonPatch
import matplotlib as mpl
import rasterio
from rasterio.plot import show
import rasterio.windows as windows
from matplotlib import pyplot as plt
Some dummy data and values
base_tif = "https://github.com/rasterio/rasterio/raw/main/tests/data/RGB2.byte.tif"
minx, miny, maxx, maxy = 161317.5, 2665342.5, 279982.5, 2773057.5
simple_box_poly = [[190983.75, 2746128.75], [250316.25, 2746128.75], [250316.25, 2692271.25], [190983.75, 2692271.25]]
with rasterio.open(base_tif) as src:
    # Create a window and get it's transform
    window = windows.from_bounds(minx, miny, maxx, maxy, src.transform)
    transform = windows.transform(window, src.transform)
# read data in window
rst = src.read(window=window)


Features
features = [{'type': 'Polygon', 'coordinates': [simple_box_poly]}]
get fig and ax
fig, ax = plt.subplots()
plot the data using transform on above ax
show(rst, transform=transform, ax=ax)
plot features on data ax
patches = [PolygonPatch(feature, edgecolor="red", facecolor="red", linewidth=2) for feature in features]
ax.add_collection(mpl.collections.PatchCollection(patches, match_original=True))
ax.axis('scaled')
plt.show()