Рисование фигуры на изображении с помощью matplotlib

Question

Из многоугольника я создаю квадратный буфер, чтобы создать спутниковое изображение моего местоположения. Фигура определяется как файл .shp, который я читаю с помощью геопанд.

Я хотел бы отобразить изображение И фигуру на одном графике с помощью matplotlib, финальный результат должен выглядеть так:

Мне удается отобразить и растянуть изображение на фигуре matplotlib

with rio.open(file) as f:
    data = f.read([1, 2, 3], masked=True)
                
    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)[0], getPositionPdf(i)[1]]
    ax.imshow(data, interpolation='nearest')
    #[...] unrelevant display customization

Но я не знаю, как отобразить фигуру поверх нее. Кто-нибудь знает, как проделать этот фокус?

Answer 1

Посмотрите, что означает масштаб по оси x и y изображения с помощью matplotlib

Шейп-файл точечного вектора (декартова проекция) :

 import geopandas as gpd
 df = gpd.read_file("points.shp")
 df['x'] = df.geometry.x
 df['y'] = df.geometry.y
 df.head(2)
    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')

Растр:

а) с GDAL (декартовой проекцией)

  from osgeo import gdal 
  ds = gdal.Open(dem)
  data = ds.ReadAsArray()
  # plot the raster
  fig, ax = plt.subplots()
  img = ax.imshow(data)
  plt.show()

Мы можем видеть, что мы не можем нанести точки на изображение, но если мы вычислим реальное расширение растра (используя результат геотрансформации gdal) для экстента matplotlib:

 gt = ds.GetGeoTransform()     
 extent = (gt[0], gt[0] + ds.RasterXSize * gt[1],gt[3] + ds.RasterYSize * gt[5], gt[3])
 fig, ax = plt.subplots()
 img = ax.imshow(data, extent=extent, origin='upper')
 ax.scatter(df.x,df.y,c='r')
 plt.show()

б) с растерио (используя прямое растерио преобразование )

import rasterio
from rasterio.plot import show
ds = rasterio.open(dem)
fig, ax = plt.subplots()
show(ds2.read(), transform=ds2.transform, ax=ax)
ax.scatter(df.x,df.y,c='r')
plt.show()