目录
一、多波段合成
本案例使用的是大疆多光谱无人机拍摄的冬小麦多光谱数据,共有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图层,选择数据——导出数据,设置保存位置、名称、格式
继续使用代码查看导出数据的地理信息,发现地理信息已更新
至此,遥感图像的波段合成及地理配准操作已完成