import os
import gdal
import numpy as np
import pandas as pd
import datetime
import glob
def readTif(fileName):
dataset = gdal.Open(fileName)
if dataset == None:
print(fileName + "文件无法打开")
width = dataset.RasterXSize
height = dataset.RasterYSize
bands = dataset.RasterCount
data = dataset.ReadAsArray()
geotrans = dataset.GetGeoTransform()
proj = dataset.GetProjection()
del dataset
return width, height, bands, data, geotrans, proj
def writeTiff(im_data,input_data,im_geotrans,xOffset,yOffset,im_proj, path):
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
im_bands, im_height, im_width = im_data.shape
if len(input_data.shape) == 3:
input_bands, input_height, input_width = input_data.shape
elif len(input_data.shape) == 2:
input_data = np.array([input_data])
input_bands, input_height, input_width = input_data.shape
else:
print("image is error")
#创建文件
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(path, int(im_width), int(im_height), int(input_bands), datatype)
if(dataset!= None):
dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数
dataset.SetProjection(im_proj) #写入投影
for i in range(input_bands):
dataset.GetRasterBand(i + 1).WriteArray(input_data[i],xOffset,yOffset)
del dataset
if __name__ == '__main__':
start = datetime.datetime.now()
ref_tif = r"G:\xdn\class_result\sy_index_result_pro.tif"
samplefile = r'G:\xdn\sample\sy_sample_pro211.tif'
path = r'G:\xdn\sample\extecd_sample5.tif'
width_ref, height_ref, bands_ref, data_ref, geotrans_ref, proj_ref = readTif(ref_tif)
width_s, height_s, bands_s, data_s, geotrans_s, proj_s = readTif(samplefile)
print("开始扩充影像。。。。。。。。")
#计算待扩充影像在扩充影像左上角的偏移量
xOffset=int((geotrans_s[0]-geotrans_ref[0])/geotrans_ref[1])
yOffset=int((geotrans_s[3]-geotrans_ref[3])/geotrans_ref[5])
im_data = data_ref
input_data = data_s
im_geotrans = geotrans_ref
im_proj = proj_ref
print(xOffset,yOffset)
print(width_ref,height_ref)
print(width_s,height_s)
writeTiff(im_data,input_data,im_geotrans,xOffset,yOffset,im_proj, path)
end = datetime.datetime.now()
print(end-start)
python如何使用gdal扩充影像
最新推荐文章于 2024-06-14 17:09:53 发布