Python使用gdal实现影像镶嵌

Python使用gdal实现影像镶嵌

如果要对某个文件夹下的多景影像进行镶嵌，我们需要知道镶嵌后影像的行列数目，以及该影像的六个地理变换参数，（值得特别注意的是，无论是影像的重采样还是镶嵌，都需要特别关注影像的6个地理变换参数），关于这六个地理变换参数，请参考我的另一篇博文：https://blog.csdn.net/SunStrongInChina/article/details/104262949
为获取影像的行列数和6个地理变换参数，我们需要知道镶嵌后影像的左上角坐标，右下角坐标，像元宽度和像元高度，则待镶嵌影像的行列数可以通过坐标差值除以像元高度或宽度求取。由于待镶嵌影像的左上角的x坐标值为所有影像的最小x坐标值，而左上角的y坐标值为最大的y坐标值，右下角的x坐标值为最大的x坐标值，右下角的y坐标值为最小的y坐标值。
首先，定义获取影像的左上角和右下角坐标值的函数GetExtent：

from osgeo import gdal
import os
import glob
import math
#获取影像的左上角和右下角坐标
def GetExtent(in_fn):
    ds=gdal.Open(in_fn)
    geotrans=list(ds.GetGeoTransform())
    xsize=ds.RasterXSize 
    ysize=ds.RasterYSize
    min_x=geotrans[0]
    max_y=geotrans[3]
    max_x=geotrans[0]+xsize*geotrans[1]
    min_y=geotrans[3]+ysize*geotrans[5]
    ds=None
    return min_x,max_y,max_x,min_y

第二步，获取目录下所有影像的左上角坐标，计算得到镶嵌后影像的左上角和右下角坐标，

path=r"C:\Users\Sun Strong\Desktop\mosaic"
os.chdir(path)
#如果存在同名影像则先删除
if os.path.exists('mosaiced_image.tif'):
    os.remove('mosaiced_image.tif')
in_files=glob.glob("*.tif")#得到该目录下所有的影像名
in_fn=in_files[0]
#获取待镶嵌栅格的最大最小的坐标值
min_x,max_y,max_x,min_y=GetExtent(in_fn)
for in_fn in in_files[1:]:
    minx,maxy,maxx,miny=GetExtent(in_fn)
    min_x=min(min_x,minx)
    min_y=min(min_y,miny)
    max_x=max(max_x,maxx)
    max_y=max(max_y,maxy)

第三步，计算镶嵌后影像的行列数目，并创建镶嵌后的影像句柄

#计算镶嵌后影像的行列号
in_ds=gdal.Open(in_files[0])
geotrans=list(in_ds.GetGeoTransform())
width=geotrans[1]
height=geotrans[5]

columns=math.ceil((max_x-min_x)/width)#列数
rows=math.ceil((max_y-min_y)/(-height))#行数
in_band=in_ds.GetRasterBand(1)

driver=gdal.GetDriverByName('GTiff')
out_ds=driver.Create('mosaiced_image.tif',columns,rows,1,in_band.DataType)
out_ds.SetProjection(in_ds.GetProjection())
geotrans[0]=min_x#更正左上角坐标
geotrans[3]=max_y
out_ds.SetGeoTransform(geotrans)
out_band=out_ds.GetRasterBand(1)
#定义仿射逆变换
inv_geotrans=gdal.InvGeoTransform(geotrans)#用于给定像元实际坐标，求取像元行列号位置

最后，将各幅影像的像元值读取出来，写入到创建好的影像中

#开始逐渐写入
for in_fn in in_files:
    in_ds=gdal.Open(in_fn)
    in_gt=in_ds.GetGeoTransform()
    #将左上角坐标转换为待写入影像的行列号
    offset=gdal.ApplyGeoTransform(inv_geotrans,in_gt[0],in_gt[3])
    x,y=map(int,offset)
    print(x,y)
    #第二种求取像元行列号的方法
    trans=gdal.Transformer(in_ds,out_ds,[])#in_ds是源栅格，out_ds是目标栅格
    success,xyz=trans.TransformPoint(False,0,0)#计算in_ds中左上角像元对应out_ds中的行列号
    x,y,z=map(int,xyz)
    print(x,y,z)
    data=in_ds.GetRasterBand(1).ReadAsArray()
    out_band.WriteArray(data,x,y)#x，y是开始写入时左上角像元行列号
del in_ds,out_band,out_ds

总结

根据像元实际坐标反求像元的行列号时，有两种方法可以使用，第一种方法：

inv_geotrans=gdal.InvGeoTransform(geotrans)#f仿射逆变换参数
in_ds=gdal.Open(in_fn)
in_gt=in_ds.GetGeoTransform()
#将左上角坐标转换为待写入影像的行列号
offset=gdal.ApplyGeoTransform(inv_geotrans,in_gt[0],in_gt[3])
x,y=map(int,offset)

第二种方法：

 trans=gdal.Transformer(in_ds,out_ds,[])#in_ds是源栅格，out_ds是目标栅格
 success,xyz=trans.TransformPoint(False,0,0)#Fasle代表计算in_ds中左上角像元对应out_ds中的行列号，Ture则相反
 x,y,z=map(int,xyz)#x,y为像元所在的列号和行号

代码汇总

'''
Created on 2020年2月11日

@author: Sun Strong
'''
from osgeo import gdal
import os
import glob
import math
#获取影像的左上角和右下角坐标
def GetExtent(in_fn):
    ds=gdal.Open(in_fn)
    geotrans=list(ds.GetGeoTransform())
    xsize=ds.RasterXSize 
    ysize=ds.RasterYSize
    min_x=geotrans[0]
    max_y=geotrans[3]
    max_x=geotrans[0]+xsize*geotrans[1]
    min_y=geotrans[3]+ysize*geotrans[5]
    ds=None
    return min_x,max_y,max_x,min_y
path=r"C:\Users\Sun Strong\Desktop\mosaic"
os.chdir(path)
#如果存在同名影像则先删除
if os.path.exists('mosaiced_image.tif'):
    os.remove('mosaiced_image.tif')
in_files=glob.glob("*.tif")
in_fn=in_files[0]
#获取待镶嵌栅格的最大最小的坐标值
min_x,max_y,max_x,min_y=GetExtent(in_fn)
for in_fn in in_files[1:]:
    minx,maxy,maxx,miny=GetExtent(in_fn)
    min_x=min(min_x,minx)
    min_y=min(min_y,miny)
    max_x=max(max_x,maxx)
    max_y=max(max_y,maxy)
#计算镶嵌后影像的行列号
in_ds=gdal.Open(in_files[0])
geotrans=list(in_ds.GetGeoTransform())
width=geotrans[1]
height=geotrans[5]

columns=math.ceil((max_x-min_x)/width)
rows=math.ceil((max_y-min_y)/(-height))
in_band=in_ds.GetRasterBand(1)

driver=gdal.GetDriverByName('GTiff')
out_ds=driver.Create('mosaiced_image.tif',columns,rows,1,in_band.DataType)
out_ds.SetProjection(in_ds.GetProjection())
geotrans[0]=min_x
geotrans[3]=max_y
out_ds.SetGeoTransform(geotrans)
out_band=out_ds.GetRasterBand(1)
#定义仿射逆变换
inv_geotrans=gdal.InvGeoTransform(geotrans)
#开始逐渐写入
for in_fn in in_files:
    in_ds=gdal.Open(in_fn)
    in_gt=in_ds.GetGeoTransform()
    #仿射逆变换
    offset=gdal.ApplyGeoTransform(inv_geotrans,in_gt[0],in_gt[3])
    x,y=map(int,offset)
    print(x,y)
    trans=gdal.Transformer(in_ds,out_ds,[])#in_ds是源栅格，out_ds是目标栅格
    success,xyz=trans.TransformPoint(False,0,0)#计算in_ds中左上角像元对应out_ds中的行列号
    x,y,z=map(int,xyz)
    print(x,y,z)
    data=in_ds.GetRasterBand(1).ReadAsArray()
    out_band.WriteArray(data,x,y)#x，y是开始写入时左上角像元行列号
del in_ds,out_band,out_ds