การวาดรูปร่างบนภาพด้วย matplotlib

Question

จากรูปหลายเหลี่ยมฉันสร้างบัฟเฟอร์สี่เหลี่ยมเพื่อสร้างภาพ sattelite ของตำแหน่งของฉัน รูปร่างถูกกำหนดให้เป็นไฟล์. shp ที่ฉันอ่านด้วย geopandas

ฉันต้องการแสดงภาพและรูปร่างบนกราฟเดียวกันโดยใช้ 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

ดูความหมายของมาตราส่วนบนแกน x และ y ของภาพโดยใช้ matplotlib คืออะไร

ไฟล์รูปร่างเวกเตอร์จุด (การฉายภาพคาร์ทีเซียน) :

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

เราจะเห็นว่าเราไม่สามารถพล็อตจุดบนภาพได้ แต่ถ้าเราคำนวณส่วนขยายที่แท้จริงของแรสเตอร์ (โดยใช้ผลลัพธ์ของการเปลี่ยนรูปแบบ gdal) สำหรับขอบเขต 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) ด้วย 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()