Rysowanie kształtu na obrazie za pomocą matplotlib

Question

Z kształtu wielokąta tworzę kwadratowy bufor, aby utworzyć satelitarny obraz mojej lokalizacji. Kształt jest definiowany jako plik .shp, który czytałem za pomocą geopand.

Chciałbym wyświetlić obraz ORAZ kształt na tym samym wykresie za pomocą matplotlib, końcowy wynik powinien wyglądać następująco:

Udaje mi się wyświetlić i rozciągnąć obraz na figurze 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

Ale nie wiem, jak wyświetlić kształt na wierzchu. Czy ktoś wie, jak wykonać tę sztuczkę?

Answer 1

Spójrz na znaczenie skali na osi X i Y obrazu za pomocą matplotlib

Plik kształtu wektora punktowego (rzut kartezjański) :

 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')

Raster:

a) z GDAL (rzut kartezjański)

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

Widzimy, że nie możemy wykreślić punktów na obrazie, ale jeśli obliczymy rzeczywiste rozszerzenie rastra (używając wyniku gdal geotransform) dla zasięgu 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()

b) z rasterio (używając bezpośrednio transformacji rasterio )

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()