Python基于GDAL镶嵌拼接遥感影像

更新时间：2023年10月23日 14:42:21 作者：RS迷途小书童

这篇文章主要介绍了Python基于GDAL镶嵌拼接遥感影像, 这里有一点需要注意的就是,用这个方法进行镶嵌拼接操作时,影像有一条明显的拼接线,不知道是不是我数据的问题,你们可以自己尝试一下,只要修改主函数中的路径即可,需要的朋友可以参考下

没啥好说的，处理高分辨率影像时，数据高达几十G。用ENVI或者ArcGIS进行影像的拼接时，往往会出现未响应的情况。出现未响应的话，运气好等个一晚上可能会动一动，运气不好就等着强制关闭重做吧。

所以搞了一个Python进行拼接操作的代码，虽然速度不算快，但至少不会未响应。同时如果对代码进行一些改进，还可以进行批量拼接的操作，百利而无一害。

一、导入GDAL库

from osgeo import gdal

二、查看影像信息

为了凑字数的，可以查看影像的投影、宽度、高度、波段数等信息。不过需要注意的是在ENVI中没有投影坐标系，只有地理坐标系是做不了镶嵌拼接的。

这个代码我还不太清楚能不能不要投影坐标系进行拼接，你们可以自己试试。但最好还是用包含投影坐标系的影像进行拼接。所以在拼接之前就可以用这段代码先看一看。

def Get_data(filepath):
    ds = gdal.Open(filepath)  # 打开数据集dataset
    ds_width = ds.RasterXSize  # 获取数据宽度
    ds_height = ds.RasterYSize  # 获取数据高度
    ds_bands = ds.RasterCount  # 获取波段数
    ds_geo = ds.GetGeoTransform()  # 获取仿射地理变换参数
    ds_prj = ds.GetProjection()  # 获取投影信息
    print("影像的宽度为：" + str(ds_width))
    print("影像的高度为：" + str(ds_height))
    print("仿射地理变换参数为：" + str(ds_geo))
    print("投影坐标系为：" + str(ds_prj))
    # data = ds.ReadAsArray(0, 0, ds_width, ds_height)  # 以数组的形式读取整个数据集

三、镶嵌模块

这里用到了GDAL中的Warp函数，是不是有点熟悉。之前裁剪也是用的这个函数，不得不说这个函数是真的强大，之前有文章介绍了其中的函数，大家有兴趣可以去看下，同时给个赞吧！【Python&RS】GDAL批量裁剪遥感影像/栅格数据

代码中的srcSRS,dstSRS分别是输入投影和输出投影，这里用一样的就行了。因为我们做的是镶嵌操作，肯定是不用动原始坐标系的。其他的参数都在代码中表明了，这里就不介绍了，如果大家有什么问题，可以留言交流。

def Mosaic_GDAL(path_image1, path_image2, path_out):
    """
    :param path_image1: 需要镶嵌的影像
    :param path_image2: 需要镶嵌的影像
    :param path_out: 镶嵌后输出的影像路径
    :return: None
    """
    image1 = gdal.Open(path_image1, gdal.GA_ReadOnly)  # 第一幅影像
    input_proj = image1.GetProjection()
    image2 = gdal.Open(path_image2, gdal.GA_ReadOnly)  # 第二幅影像
    options = gdal.WarpOptions(srcSRS=input_proj, dstSRS=input_proj, format='GTiff',
                               resampleAlg=gdal.GRA_NearestNeighbour,callback=Show_Progress)
    # 输入投影，输出投影，输出格式，重采样方法
    gdal.Warp(path_out, [image1, image2], options=options)
    # 输出路径，需要镶嵌的数据，参数配置
    ds = gdal.Open(path_out, gdal.GA_ReadOnly)
    ds.BuildOverviews(overviewlist=[2, 4, 8, 16])
    # 创建金字塔
    del image1, image2, ds

四、回调函数（没啥用）

上一篇博文已经介绍过了，就两点：1.水字数，2.记录进度。省的看着代码发呆。

def Show_Progress(percent, msg, tag):
    """
    :param percent: 进度，0~1
    :param msg:
    :param tag:
    :return:
    """
    if 25 <= percent*100 <= 26:
        print("进度：" + "%.2f" % (percent*100) + "%")
    if 50 <= percent*100 <= 51:
        print("进度：" + "%.2f" % (percent*100) + "%")
    if 75 <= percent*100 <= 76:
        print("进度：" + "%.2f" % (percent*100) + "%")

五、完整代码

# -*- coding: utf-8 -*-
"""
@Time ： 2023/6/25 16:28
@Auth ： RS迷途小书童
@File ：Raster Data Mosaic.py
@IDE ：PyCharm
@Purpose：遥感影像、栅格数据的镶嵌拼接（有明显拼接线）
"""
from osgeo import gdal
def Get_data(filepath):
    ds = gdal.Open(filepath)  # 打开数据集dataset
    ds_width = ds.RasterXSize  # 获取数据宽度
    ds_height = ds.RasterYSize  # 获取数据高度
    ds_bands = ds.RasterCount  # 获取波段数
    ds_geo = ds.GetGeoTransform()  # 获取仿射地理变换参数
    ds_prj = ds.GetProjection()  # 获取投影信息
    print("影像的宽度为：" + str(ds_width))
    print("影像的高度为：" + str(ds_height))
    print("仿射地理变换参数为：" + str(ds_geo))
    print("投影坐标系为：" + str(ds_prj))
    # data = ds.ReadAsArray(0, 0, ds_width, ds_height)  # 以数组的形式读取整个数据集
def Show_Progress(percent, msg, tag):
    """
    :param percent: 进度，0~1
    :param msg:
    :param tag:
    :return:
    """
    if 25 <= percent*100 <= 26:
        print("进度：" + "%.2f" % (percent*100) + "%")
    if 50 <= percent*100 <= 51:
        print("进度：" + "%.2f" % (percent*100) + "%")
    if 75 <= percent*100 <= 76:
        print("进度：" + "%.2f" % (percent*100) + "%")
def Mosaic_GDAL(path_image1, path_image2, path_out):
    """
    :param path_image1: 需要镶嵌的影像
    :param path_image2: 需要镶嵌的影像
    :param path_out: 镶嵌后输出的影像路径
    :return: None
    """
    image1 = gdal.Open(path_image1, gdal.GA_ReadOnly)  # 第一幅影像
    input_proj = image1.GetProjection()
    image2 = gdal.Open(path_image2, gdal.GA_ReadOnly)  # 第二幅影像
    options = gdal.WarpOptions(srcSRS=input_proj, dstSRS=input_proj, format='GTiff',
                               resampleAlg=gdal.GRA_NearestNeighbour, callback=Show_Progress)
    # 输入投影，输出投影，输出格式，重采样方法
    gdal.Warp(path_out, [image1, image2], options=options)
    # 输出路径，需要镶嵌的数据，参数配置
    ds = gdal.Open(path_out, gdal.GA_ReadOnly)
    ds.BuildOverviews(overviewlist=[2, 4, 8, 16])
    # 创建金字塔
    del image1, image2, ds
if __name__ == "__main__":
    path1 = "B:/proj_clip"
    path2 = "B:/2_proj_clip"
    path_output = "warp1.tif"
    print("开始拼接......")
    Mosaic_GDAL(path1, path2, path_output)
    print("拼接完成......")

这里有一点需要注意的就是，用这个方法进行镶嵌拼接操作时，影像有一条明显的拼接线，不知道是不是我数据的问题。你们可以自己尝试一下。只要修改主函数中的路径即可。

如果大家在学习Python或者RS时有什么问题，可以随时留言交流！

到此这篇关于Python基于GDAL镶嵌拼接遥感影像的文章就介绍到这了,更多相关Python遥感影像镶嵌内容请搜索脚本之家以前的文章或继续浏览下面的相关文章希望大家以后多多支持脚本之家！

您可能感兴趣的文章:

OpenCV实现手势虚拟拖拽的使用示例(附demo)
本文主要介绍了OpenCV实现手势虚拟拖拽的使用示例,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧
2023-11-11
pyCharm 实现关闭代码检查
这篇文章主要介绍了pyCharm 实现关闭代码检查，具有很好的参考价值，希望对大家有所帮助。一起跟随小编过来看看吧
2020-06-06
Pycharm如何打断点的方法步骤
这篇文章主要介绍了Pycharm如何打断点的方法步骤，文中通过示例代码介绍的非常详细，对大家的学习或者工作具有一定的参考学习价值，需要的朋友们下面随着小编来一起学习学习吧
2019-06-06
python中numpy 常用操作总结
这篇文章主要介绍了python中numpy常用操作总结，NumPy是Python语言的一个扩充程序库，支持大量高维度数组与矩阵运算，此外也针对数组运算提供大量的数学函数库
2022-09-09
python之pandas用法大全
本文讲解了python的pandas基本用法，大家可以参考下
2018-03-03
Python3.5装饰器典型案例分析
这篇文章主要介绍了Python3.5装饰器,结合实例形式分析了装饰器修饰有参数函数、装饰器修饰函数参数等情况相关使用技巧,需要的朋友可以参考下
2019-04-04
Python上级目录文件导入的几种方法(from.import)
有时候我们可能需要import另一个路径下的python文件,下面这篇文章主要给大家介绍了关于Python上级目录文件导入的几种方法,文中通过代码介绍的非常详细,需要的朋友可以参考下
2023-12-12
Django objects的查询结果转化为json的三种方式的方法
这篇文章主要介绍了Django objects的查询结果转化为json的三种方式的方法，小编觉得挺不错的，现在分享给大家，也给大家做个参考。一起跟随小编过来看看吧
2018-11-11
Python数学建模学习模拟退火算法旅行商问题示例解析
模拟退火算法不仅可以解决连续函数优化问题，KIRKPATRICK在1983年成功将其应用于求解组合优化问题,现已成为求解旅行商问题的常用方法，通常采用反序、移位和交换等操作算子产生新解
2021-10-10
Python使用pymupdf实现PDF加密
这篇文章主要介绍了如何使用 Python 和 wxPython 库创建一个简单的图形用户界面（GUI）应用程序,用于对 PDF 文件进行加密,感兴趣的小伙伴可以了解下
2023-08-08