首页 > 编程语言 >【Python&RS】基于GDAL镶嵌拼接遥感影像

【Python&RS】基于GDAL镶嵌拼接遥感影像

时间:2023-06-29 14:56:34浏览次数:48  
标签:RS Python 投影 拼接 影像 GDAL path ds 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时有什么问题,可以随时留言交流!如果大家对批量处理有兴趣同样可以留言给博主,博主会分享相关代码以供学习!

标签:RS,Python,投影,拼接,影像,GDAL,path,ds,gdal
From: https://www.cnblogs.com/RSran/p/17514172.html

相关文章

  • 【python爬虫案例】用python爬豆瓣读书TOP250排行榜!
    目录一、爬虫对象-豆瓣读书TOP250二、python爬虫代码讲解三、讲解视频四、完整源码一、爬虫对象-豆瓣读书TOP250今天我们分享一期python爬虫案例讲解。爬取对象是,豆瓣读书TOP250排行榜数据:https://book.douban.com/top250​开发好python爬虫代码后,爬取成功后的csv数据,如下:​......
  • rsync备份任务练习
    06-备份任务实战今天的任务主要以实际备份任务入手,完成综合练习,完成对rsync的综合运用。先看需求再讲解再次动手实践  客户端需求客户端需求:1.客户端每天凌晨1点在服务器本地打包备份(/etc目录和/var/log目录)2.客户端备份的数据必须存放至以"主机名_ip地......
  • Python调试中特殊的属性值: __iter__
    在调试python源码时遇到一件奇怪的事情,某个属性的值为__iter__,并且给它明确的赋别的值,发现并没有作用,如下所示:仔细思考了一下,有点儿明白了,__getattr__方法的作用是当找不到对象的属性时会调用此函数,这时,method的值'deviceInfo'就是此对象不存在的属性,所以进......
  • python将excel中的数据写入固定的word格式中
    fromopenpyxlimportload_workbookfromdocximportDocumentimportdatetimefromdatetimeimportdatetimeasdatetwb=load_workbook("收集表.xlsx")ws=wb['收集表']fromdocx.sharedimportInches,Cm#读取表中某行数据ll=[]foriinws.column......
  • Python控制鼠标操作神器
    我们以前讲过怎样使用Python在浏览器中实现页面自动化操作,不管用哪种方式实现,都是通过定位页面中的元素来进行相应的操作。今天我们来聊一聊如何在桌面实现自动化操作。与浏览器页面自动化操作类似,桌面自动化操作也是需要定位鼠标在桌面的位置,然后根据定位的位置执行对应的操作。......
  • 利用python将数据写入CSV文件中
    利用python将数据写入CSV文件中全部代码如下:importcsv#1.创建文件对象f=open('cav_file.csv','w',encoding='utf-8',newline="")#2.基于文件对象构建csv写入对象csv_write=csv.writer(f)#3.构建列表头csv_write.writerow(['学号',......
  • 关于python刷题的语法要点
    赋值lista,bb=a,itisjustareferenceassignment.a改变,b也改变b=a[:],itiscopy.a改变,b不变b=a+[4]a[::]里::是一种operater,默认是a[::1],1代表从头到尾的取间隔为1的元素组成新的list。所以a[::-1],就是reversed,还可以有a[::-2]排序listxx.sort()ory=sor......
  • Java or Python?测试开发工程师如何选择合适的编程语言?
    很多测试开发工程师尤其是刚入行的同学对编程语言和技术栈选择问题特别关注,毕竟掌握一门编程语言要花不少时间成本,也直接关系到未来的面试和就业(不同企业/项目对技术栈要求也不一样),根据自身情况做一个相对正确的选择确实要比盲目投入更明智也更高效。目前最常见的情况是纠结选择Jav......
  • Course 3 EC2 Fundamentals
    AmazonEC2EC2=ElasticComputeCloud=InfrastructureasServiceItmainlyconsistsinthecapabilityof:Rentingvirtualmachines(EC2)Storingdataonvirtualdrives(EBS)Distributingloadacrossmachines(ELB)Scalingtheservicesusinganauto-sc......
  • python题目:使用python将指定字符串如:"aabbc"中的"ab"串删除,删除过后的数据为:c【杭州多
    定义一个函数,在函数体内完成以下逻辑:使用python将指定字符串如:"aabbc"中的"ab"串删除,删除过后的数据为:cimportsyssys.setrecursionlimit(100000)#设置递归的深度为10万defdel1(str1):#str1='aabbc'while"ab"notinstr1:print(str1)exi......