Matplotlib ile bir görüntü üzerine şekil çizme

Question

Çokgen şeklinden, konumumun uydu görüntüsünü oluşturmak için kare bir tampon oluşturuyorum. Şekil, jeopandalarla okuduğum bir .shp dosyası olarak tanımlandı.

Matplotlib kullanarak resmi VE şekli aynı grafikte görüntülemek istiyorum, sonuç şu şekilde görünmelidir:

Görüntüyü bir matplotlib figüründe göstermeyi ve uzatmayı başardım

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

Ama üstündeki şekli nasıl sergileyeceğimi bilmiyorum. Bu numarayı nasıl yapacağını bilen var mı?

Answer 1

Bak x ve görüntünün y ekseni matplotlib kullanımıyla ilgili ölçek anlam nedir

Bir nokta vektörü şekil dosyası (kartezyen projeksiyon) :

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

Bir raster:

a) GDAL ile (kartezyen projeksiyon)

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

Görüntüdeki noktaları çizemediğimizi görebiliriz, ancak matplotlib kapsamı için rasterin gerçek uzantısını (gdal geotransform sonucunu kullanarak) hesaplarsak:

 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) rasterio ile (doğrudan rasterio dönüşümü kullanarak )

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