Zeichnen einer Form auf einem Bild mit matplotlib

Question

Aus einer Polygonform erstelle ich einen quadratischen Puffer, um ein Satellitenbild meiner Position zu erstellen. Die Form ist als .shp-Datei definiert, die ich mit Geopandas gelesen habe.

Ich möchte das Bild UND die Form mit matplotlib im selben Diagramm anzeigen. Das Endergebnis sollte folgendermaßen aussehen:

Es gelingt mir, das Bild auf einer Matplotlib-Figur anzuzeigen und zu dehnen

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

Aber ich weiß nicht, wie ich die Form darüber anzeigen soll. Weiß jemand, wie man diesen Trick ausführt?

Answer 1

Schauen Sie sich mit matplotlib an, was Skalierung auf der x- und y-Achse des Bildes bedeutet

Ein Punktvektor-Shapefile (kartesische Projektion) :

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

Ein Raster:

a) mit GDAL (kartesische Projektion)

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

Wir können sehen, dass wir die Punkte auf dem Bild nicht zeichnen können, aber wenn wir die tatsächliche Ausdehnung des Rasters (unter Verwendung des Ergebnisses der gdal-Geotransform) für die Matplotlib-Ausdehnung berechnen:

 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) mit Rasterio (mit direkter Rasterio-Transformation )

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