使用arcgis对农业多光谱数据多波段合成后进行地里配准并保存

一、多波段合成

本案例使用的是大疆多光谱无人机拍摄的冬小麦多光谱数据,共有red、green、blue、rededge、nir五个波段。

打开arcgis,点击工具栏的arctoolbox,选择数据管理工具

然后选择栅格——栅格处理——波段合成

选择文件夹中的五个波段的tif文件,并设置输出路径及名称,点击保存和确定

等待一段时间后,左侧内容列表出现合成结果

二、地理信息配准

由于拍摄多光谱的无人机只有误差较大的GPS地理信息,没有精准的RTK信息,需要手动对其进行地理配准,使用实地打固定点的方法进行配准,固定点的RTK信息存储于照片中,将照片的地理信息转化为shp文件。具体代码如下:

提取照片中的RTK地理坐标信息,转化为十进制并将定点信息保存至表格中:

# 本文件用来读取大疆精灵4拍摄照片中的RTK地理坐标数据,并转化为度、分、秒的格式,最后转化为十进制的地理信息格式
import exifread
import os
import re
import xlwt
# 定义地理信息转换类
class Degree(object):
    def __init__(self):
        None
    @staticmethod
    def dd_to_dms(dd):
        """
        十进制度转为度分秒
        Paramaters:
            dd : 十进制度
        Return:
            dms : 度分秒
        """
        degree=(int)(float(dd))
        minute=(int)((float(dd)-degree)*60)
        second=round((float(dd)-degree-float(minute)/60)*3600,2)
        dms=str(degree)+'°'+str(minute)+ '′'+str(second)+"″"
        return dms
    @staticmethod
    def dms_to_dd(degree,minute,second):
        """
        度分秒转为十进制度
        Paramater:
            degree : 度
            minute : 分
            second : 秒
        Return:
            dd : 十进制度
        """
        dd=degree+minute/60+second/60/60
        return dd
    @staticmethod
    def parse_dms(dms):
        """
        解析度分秒字符串
        Paramater:
            dms : 度分秒字符串
        Returns:
            degree : 度
            minute : 分
            second : 秒
        """
        parts = re.split('[°′″]', dms)
        degree=float(parts[0])
        minute=float(parts[1])
        second=float(parts[2])
        return {"degree":degree,"minute":minute,"second":second}

#十进制度转为度分秒
# Degree.dd_to_dms(104.25456666666666)
#度分秒转为十进制度
# Degree.dms_to_dd(104,15,16.44)

#将每张照片的所有信息分别保存至txt文件中
def read_img(loadfile_root,loadsavepath,filetype):
    """
    loadfile_root:照片文件夹路径
    loadsavepath:输出文件的路径
    filetype:文件类型(区分大小写)如JPG
    """
    file_root =loadfile_root.replace('/','//')#路径转换
    file_list = os.listdir(file_root)#得到该文件夹下的所有文件的名称
    #遍历该文件夹下的所有文件
    for img_name in file_list:
        str1 = loadsavepath.split("\\")
        if img_name != str1[-1]:
            savepath=loadsavepath.replace('/','//')#路径转换
            txtrealpath=savepath+'/'+img_name.replace('.'+filetype,'')+'.txt'#得到.txt的文件路径
            txt = open(txtrealpath, 'w')#新建txt文件
            realpath=file_root+'/'+img_name#得到图片的文件路径
            f = open(realpath, 'rb')#打开文件''中输入文件路径
            tags = exifread.process_file(f)#读取tags,其为dict类型的文件0
        # 遍历字典的元素,将每项元素的key和value分拆组成字符串,注意添加分隔符和换行符
            for k, v in tags.items():
                 txt.write(str(k) + ':' + str(v) + '\n')
            txt.close()#    注意关闭txt
        else:
            continue

# 此函数将读取指定文件夹下的所有txt文件,寻找到地理信息字段,保存为度分秒格式
def get_lon_lat(folder_path, file_path):
    data1 = []  # 存放经度
    data2 = []  # 存放纬度
    # 遍历文件夹下所有txt文件
    for filename in os.listdir(folder_path):
        if filename.endswith(".txt"):
            file_path1 = os.path.join(folder_path, filename)
            with open(file_path1, "r", encoding="utf-8") as file:
                content = file.read()
                # 使用正则表达式找到指定字段后面的字段
                # matches = re.findall(r'Image_Make\s*:\s*(.*?)\s*', content)
                matches1 = re.findall(r'GPS GPSLongitude:\[(.*?)]', content)
                matches2 = re.findall(r'GPS GPSLatitude:\[(.*?)]', content)
                # 使用split()方法分割字符串
                matches1 = [match.strip() for match in matches1[0].split(",")]
                matches1[2] = eval(matches1[2])
                matches2 = [match.strip() for match in matches2[0].split(",")]
                matches2[2] = eval(matches2[2])
                # print(matches)
                data1.append(matches1)
                data2.append(matches2)

    result1 = []
    result2 = []
    combined_strings1 = []
    combined_strings2 = []
    with open(file_path, "w") as file:
        # 经度数据
        for i in data1:
            # 字符串合并,并用空格分开
            numbers_as_strings = [str(num) for num in i]
            combined_string = ' '.join(numbers_as_strings)
            # print(combined_string)
            combined_strings1.append(combined_string)
        # print(combined_strings1)

        # 将列表中的每个字符串进行处理,使用不同的替换字符替换空格
        for j in combined_strings1:
            # 定义要用于替换空格的字符列表
            str_data = ""
            replacement_chars = ['°', '′', '″']
            for k in j:
                if k == ' ':
                    # 如果字符是空格,则从替换字符列表中取出一个字符进行替换,并添加到结果中
                    str_data += replacement_chars.pop(0)
                else:
                    # 如果字符不是空格,则直接添加到结果中
                    str_data += k
            result1.append(str_data + '″')

        # 纬度数据
        for i in data2:
            # 字符串合并,并用空格分开
            numbers_as_strings = [str(num) for num in i]
            combined_string = ' '.join(numbers_as_strings)
            # print(combined_string)
            combined_strings2.append(combined_string)

        # 将列表中的每个字符串进行处理,使用不同的替换字符替换空格
        for j in combined_strings2:
            # 定义要用于替换空格的字符列表
            str_data = ""
            replacement_chars = ['°', '′', '″']
            for k in j:
                if k == ' ':
                    # 如果字符是空格,则从替换字符列表中取出一个字符进行替换,并添加到结果中
                    str_data += replacement_chars.pop(0)
                else:
                    # 如果字符不是空格,则直接添加到结果中
                    str_data += k
            result2.append(str_data + '″')
            # 将列表内容写入到txt文件中
        for i in range(0, len(data1)):
            file.write(result1[i] + ' ' + result2[i] + "\n")
        print(result1)
        print(result2)

#读取txt文件
def openreadtxt(file_name):
    data = []
    with open(file_name, 'r') as file:  # 打开文件
        file_data = file.readlines()  # 读取所有行
        print(file_data)
    for row in file_data:
        tmp_list = row.split(' ')  # 按‘ ’切分每行的数据
        # tmp_list[-1] = tmp_list[-1].replace('\n',',') #去掉换行符
        tmp_list[-1] = tmp_list[-1].replace('\n', '')  # 去掉换行符
        tmp_list[-1] = tmp_list[-1].replace(' ', '')  # 去掉空格
        data.append(tmp_list)  # 将每行数据插入data中
    return data

if __name__ == "__main__":
    # 带有RTK坐标信息图像文件夹
    loadfile_root = r"C:\Users\10755\Desktop\多波段地理配准\0.地面采样点\20240427配准定点"
    # 保存图片全部信息的txt文档路径
    loadsavepath = r"C:\Users\10755\Desktop\多波段地理配准\0.地面采样点\20240427配准定点\1"
    # 最终输出的十进制格式excel表格存放路径
    final_path= r"C:\Users\10755\Desktop\多波段地理配准\0.地面采样点\20240427配准定点\1\final.xls"
    # 图像格式
    filetype = 'JPG'
    # 提取所有图像信息并保存至txt文件
    read_img(loadfile_root, loadsavepath, filetype)

    # txt文档路径
    folder_path = loadsavepath
    # 定义要写入的文件路径
    # file_path = r"C:\Users\10755\Desktop\多波段地理配准\定点\1.txt"
    file_path = loadfile_root + '\\' + 'info.txt'
    # 读取指定文件夹下的所有txt文件,寻找到地理信息字段,保存为度分秒格式
    get_lon_lat(folder_path, file_path)

    # 将度分秒格式的地理信息转化为十进制格式并保存至exxcel表格中
    data = []
    primary_data = openreadtxt(file_path)
    print(primary_data)
    #将txt中的数据转化为十进制坐标并存入列表
    for i in primary_data:
        # print(i[0],i[1])
        dms0 = Degree.parse_dms(i[0])
        dms1 = Degree.parse_dms(i[1])
        # 经度
        lon = Degree.dms_to_dd(dms0["degree"], dms0["minute"], dms0["second"])
        # 纬度
        lat = Degree.dms_to_dd(dms1["degree"], dms1["minute"], dms1["second"])
        data.append([lon,lat])
    # print(data)
    #将列表中的数据存入excel
    workbook = xlwt.Workbook(encoding='utf-8')  # 设置一个workbook,其编码是utf-8
    worksheet = workbook.add_sheet("test_sheet")  # 新增一个sheet
    # a = [1, 2, 3, 4, 5]  # 列1
    # b = ['a', 'b', 'c', 'd', 'e']  # 列2
    worksheet.write(0, 0, label='经度')  # 将‘经度’作为标题
    worksheet.write(0, 1, label='纬度')  # 将‘纬度’作为标题
    for i in range(len(data)):  # 循环将a和b列表的数据插入至excel
        worksheet.write(i + 1, 0, label=data[i][0])
        worksheet.write(i + 1, 1, label=data[i][1])
    workbook.save(final_path)  # 这里save需要特别注意,文件格式只能是xls,不能是xlsx,不然会报错

将表格中的经纬度坐标转换为shp文件:

import pandas as pd
import geopandas as gpd

# 读取Excel文件中的数据
excel_file = r'C:\Users\10755\Desktop\多波段地理配准\0.地面采样点\20240427配准定点\\final.xls'
output_file = r'C:\Users\10755\Desktop\多波段地理配准\0.地面采样点\20240427配准定点\shp\coordinates.shp'
df = pd.read_excel(excel_file, sheet_name='test_sheet')

# 打印数据检查
print(df.head())

# 创建包含坐标点的GeoDataFrame
gdf = gpd.GeoDataFrame(
    df,
    geometry=gpd.points_from_xy(df['经度'], df['纬度']),
    crs="EPSG:4326"  # WGS84坐标系
)

随后将shp文件导入arcgis中,具体步骤为右击图层——添加数据——选择shp文件,导入结果如下:

取消地理配准中的自动校正,随后点击添加配准点进行地理配准,注意绿色为十字坐标应为对应固定点的真实位置。

随后点击链接表,选择RMS总误差最小​​​​​且不会扭曲结果的变换

在地理配准菜单栏处点击更新地理配准,发现配准点已移至真实位置使用如下代码查看合成.tif的地理信息,发现地理配准并未成功,猜测可能是此时合成.tif为工程文件,地理配准并未生效,可以导出数据再查看。

from osgeo import gdal, osr

#查看正确的地理信息
filename1 = r'C:\Users\10755\Desktop\多波段地理配准\blog\【20240528】多光谱\1.tif'
filename2 = r'C:\Users\10755\Desktop\多波段地理配准\blog\【20240528】多光谱\新建文件夹\合成.tif'

file = [filename1,filename2]
# file1 = [filename1,filename11,filename12,filename13,filename14,filename15]
for i in file:
    dataset = gdal.Open(i)
    geo_transform = dataset.GetGeoTransform()
    print("GeoTransform:", geo_transform)

右击合成.tif图层,选择数据——导出数据,设置保存位置、名称、格式

继续使用代码查看导出数据的地理信息,发现地理信息已更新

至此,遥感图像的波段合成及地理配准操作已完成

  • 3
    点赞
  • 4
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值