Menggambar bentuk pada gambar dengan matplotlib

Question

Dari bentuk poligon saya membuat buffer persegi untuk membuat gambar satelit dari lokasi saya. Bentuknya didefinisikan sebagai file .shp yang saya baca dengan geopandas.

Saya ingin menampilkan gambar DAN bentuk pada grafik yang sama menggunakan matplotlib, hasil akhirnya akan terlihat seperti ini:

Saya berhasil menampilkan dan meregangkan gambar pada gambar 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

Tapi saya tidak tahu bagaimana menampilkan bentuk di atasnya. Adakah yang tahu bagaimana melakukan trik ini?

Answer 1

Perhatikan apa arti skala pada sumbu x dan y pada citra menggunakan matplotlib

Shapefile vektor titik (proyeksi kartesian) :

 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) dengan GDAL (proyeksi kartesian)

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

Kita dapat melihat bahwa kita tidak dapat memplot titik-titik pada gambar, tetapi jika kita menghitung ekstensi raster yang sebenarnya (menggunakan hasil geotransform gdal) untuk luasan 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) dengan rasterio (langsung menggunakan rasterio transform )

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