matplotlib로 이미지에 모양 그리기

Question

다각형 모양에서 정사각형 버퍼를 만들어 내 위치의 위성 이미지를 만듭니다. 모양은 geopandas로 읽은 .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

matplotlib를 사용하여 이미지의 x 및 y 축에서 스케일의 의미를 살펴보십시오.

점 벡터 shapefile (직교 투영) :

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

래스터 :

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

이미지에 점을 그릴 수는 없지만 matplotlib 범위에 대해 래스터의 실제 확장 (gdal geotransform의 결과 사용)을 계산하면 다음과 같습니다.

 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 사용 (직접 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()