Python使用gdal实现影像镶嵌

最新推荐文章于 2024-04-15 14:39:01 发布

原创

最新推荐文章于 2024-04-15 14:39:01 发布 · 3.8k 阅读

38 ·

CC 4.0 BY-SA版权

本文详细介绍了如何使用Python和GDAL库对多个遥感影像进行镶嵌处理，包括获取影像坐标、计算行列数、创建镶嵌影像及写入像素值等关键步骤。

Python使用gdal实现影像镶嵌

如果要对某个文件夹下的多景影像进行镶嵌，我们需要知道镶嵌后影像的行列数目，以及该影像的六个地理变换参数，（值得特别注意的是，无论是影像的重采样还是镶嵌，都需要特别关注影像的6个地理变换参数），关于这六个地理变换参数，请参考我的另一篇博文：https://blog.youkuaiyun.com/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