python解析GeoTiff檔

import rasterio
from rasterio.plot import show
fp = r'GeoTiff_Image.tif'
img = rasterio.open(fp)
show(img)

import rasterio
fname='giam_28_classes_global.tif'
img = rasterio.open(fname)
nx,ny,nz=img.width,img.height,img.count #格數
dxm,dym=(img.bounds.right-img.bounds.left)/img.width,-(img.bounds.top-img.bounds.bottom)/img.height #間距
x0,y0=img.xy(0,0) #左上方中心點
l,b,r,t=img.bounds[:] #左、下、右、上邊界線座標
lonCent,latCent=img.lnglat() #中心點經緯度
data=img.read() #shape=(1,ny,nx)，南北方向為北向南
data=np.flip(data[0,:,:],[0]) #轉向

import rasterio
from rasterio.transform import Affine

resx,resy=(np.max(lon)-np.min(lon))/(nx+2*M-1),(np.max(lat)-np.min(lat))/(ny+2*M-1)
transform = Affine.translation(np.min(lon), np.max(lat)) * Affine.scale(resx, -resy)
new_dataset = rasterio.open(TIF,'w',driver='GTiff',height=grid_z2.shape[0],width=grid_z2.shape[1],count=1,
  dtype=grid_z2.dtype,crs='+proj=latlong',transform=transform,)  
data=np.flip(grid_z2,[0]) #翻轉南北向
new_dataset.write(data, 1) #寫出數據
new_dataset.close() #關閉檔案

img = rasterio.open(fname)
l,b,r,t=img.bounds[:]
LL=False
if (l+r)/2==img.lnglat()[0]:LL=True
...
if LL:
  lon, lat = np.meshgrid(x, y)
else:
  x_g, y_g = np.meshgrid(x, y)
  Xcent=(x[0]+x[-1])/2
  Ycent=(y[0]+y[-1])/2
  Latitude_Pole, Longitude_Pole=twd97.towgs84(Xcent, Ycent)
  pnyc = Proj(proj='lcc', datum='NAD83', lat_1=10, lat_2=40,
        lat_0=Latitude_Pole, lon_0=Longitude_Pole, x_0=0, y_0=0.0)
  xgl,ygl=x_g-Xcent,  y_g-Ycent
  lon,lat=pnyc(xgl, ygl, inverse=True)
result=cntr_kml(data, lon, lat, fname)