# Drawing a shape on an image with matplotlib

From a polygon shape I create a square buffer to create a sattelite image of my location. The shape is defined as a .shp file that I read with geopandas.

I'd like to display the image AND the shape on the same graph using matplotlib, the finale result should look like that : I manage to display and stretch the image on a matplotlib figure

with rio.open(file) as f:

bands = []
for i in range(3):
band = data[i]
h_, bin_ = np.histogram(band[np.isfinite(band)].flatten(), 3000, density=True) #remove the NaN from the analysis

cdf = h_.cumsum() # cumulative distribution function
cdf = 3000 * cdf / cdf[-1] # normalize

# use linear interpolation of cdf to find new pixel values
band_equalized = np.interp(band.flatten(), bin_[:-1], cdf)
band_equalized = band_equalized.reshape(band.shape)

bands.append(band_equalized)

data = np.stack( bands, axis=0 )

data = data/3000
data = data.clip(0, 1)
data = np.transpose(data,[1,2,0])

i = year - start_year
ax = axes[getPositionPdf(i), getPositionPdf(i)]
ax.imshow(data, interpolation='nearest')
#[...] unrelevant display customization


But I don't know how to display the shape on top of it. Does anyone know how to perform this trick?

A point vector shapefile (cartesian projection):

 import geopandas as gpd
df['x'] = df.geometry.x
df['y'] = df.geometry.y
id         geometry                    x               y
0  1   POINT (203734.167 89573.589)    203734.166875   89573.588721
1  2   POINT (203981.632 89261.402)    203981.631683   89261.402347

fig, ax = plt.subplots()
ax.scatter(df.x,df.y,c='r') A raster:

a) with GDAL (cartesian projection)

  from osgeo import gdal
ds = gdal.Open(dem)
# plot the raster
fig, ax = plt.subplots()
img = ax.imshow(data)
plt.show() We can see that we can't plot the points on the image, but if we compute the real extension of the raster (using the result of gdal geotransform) for matplotlib extent:

 gt = ds.GetGeoTransform()
extent = (gt, gt + ds.RasterXSize * gt,gt + ds.RasterYSize * gt, gt)
fig, ax = plt.subplots()
img = ax.imshow(data, extent=extent, origin='upper')
ax.scatter(df.x,df.y,c='r')
plt.show() b) with rasterio (using directly rasterio transform)

import rasterio
from rasterio.plot import show
ds = rasterio.open(dem)
fig, ax = plt.subplots() • I extraxted the extents directly from rasterio with extent = [ f.bounds.left, f.bounds.right, f.bounds.bottom, f.bounds.top ]  Sep 24, 2020 at 13:17