两幅影像处理为相同地理坐标
from osgeo import osr, gdal
def read_img(img_path):
dataset = gdal.Open(img_path)
width = dataset.RasterXSize
height = dataset.RasterYSize
im_geotrans = dataset.GetGeoTransform()
im_proj = dataset.GetProjection()
src_img = dataset.ReadAsArray()
del dataset
return src_img,im_geotrans,im_proj,width,height
def change_geotrans_proj(img1,img2,img_new):
"""
坐标对齐:两幅相同大小的影像,将其中一幅的geotrans、proj赋值给另一幅
newRasterfn: 输出tif路径
rasterOrigin: 栅格数据左上角的经纬度
xsize: x方向像元大小
ysize: y方向像元大小
array: 计算后的栅格数据
"""
_, im_geotrans, im_proj, width, height = read_img(img1)
array, im_geotrans2, im_proj2, width2, height2 = read_img(img2)
cols = array.shape[1]
rows = array.shape[0]
driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create(img_new, width, height,3, gdal.GDT_Byte)
outRaster.SetGeoTransform(im_geotrans)
for i in range(1, 4):
outband = outRaster.GetRasterBand(i)
outband.WriteArray(array[i-1,:, :])
outRaster.SetProjection(im_proj)
outband.FlushCache()
del outRaster
if __name__ == "__main__":
img1=r'./1_after.tif'
img2=r"./1_before.tif"
img_new=r"./1_1_before.tif"
change_geotrans_proj(img1,img2,img_new)
print('finsh...')