from osgeo import gdal

# 讀圖像文件

def read_img(filename):

dataset = gdal.Open(filename) # 打開文件

im_width = dataset.RasterXSize # 柵格矩陣的列數

im_height = dataset.RasterYSize # 柵格矩陣的行數

im_geotrans = dataset.GetGeoTransform() # 仿射矩陣

im_proj = dataset.GetProjection() # 地圖投影信息

im_data = dataset.ReadAsArray(0, 0, im_width, im_height).astype(np.float) # 將數據寫成數組，對應柵格矩陣

del dataset # 關閉對象，文件dataset

return im_proj, im_geotrans, im_data, im_height, im_width

def write_img(filename, im_proj, im_geotrans, im_data):

# gdal數據類型包括

# gdal.GDT_Byte,

# gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,

# gdal.GDT_Float32, gdal.GDT_Float64

# 判斷柵格數據的數據類型

if 'int8' in im_data.dtype.name:

datatype = gdal.GDT_Byte

elif 'int16' in im_data.dtype.name:

datatype = gdal.GDT_UInt16

else:

datatype = gdal.GDT_Float32

# 判讀數組維數

if len(im_data.shape) == 3:

im_bands, im_height, im_width = im_data.shape

else:

im_bands, (im_height, im_width) = 1, im_data.shape

# 創建文件

driver = gdal.GetDriverByName("GTiff") # 數據類型必須有，因為要計算需要多大內存空間

dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

dataset.SetGeoTransform(im_geotrans) # 寫入仿射變換參數

dataset.SetProjection(im_proj) # 寫入投影

if im_bands == 1:

dataset.GetRasterBand(1).WriteArray(im_data) # 寫入數組數據

else:

for i in range(im_bands):

dataset.GetRasterBand(i + 1).WriteArray(im_data[i])

del dataset

if __name__=='__main__':

proj, geotrans, values, row, column = read_img(輸入數據) # 讀數據

write_img(r'輸出地址', proj, geotrans, 輸出影像名稱)#寫數據