在地球科学、气象学以及环境监测等领域,遥感影像数据是一种重要的信息源,它们可以提供地表的地形、植被覆盖、气候变化等丰富信息。然而,随着观测技术的进步,我们通常会获得大量的遥感影像数据,如何高效地处理和分析这些数据成为了一项挑战。本文将介绍如何利用 Python 中的 GDAL 库处理遥感影像数据,并通过计算年度平均影像来提取更有意义的信息。
1. 环境准备
在开始之前,确保你已经安装了 Python 和 GDAL 库。如果还没有安装,你可以通过 pip 进行安装:
pip install gdal
2. 处理单个 TIFF 文件
我们首先定义了一个函数 process_tiff_folder,它用于处理一个包含多个 TIFF 文件的文件夹。在这个函数中,我们遍历文件夹中的每个 TIFF 文件,读取其数据并提取地理信息。然后,我们将每个像素的经纬度与高程值一起保存在一个二维数组中,以便后续处理使用。
3. 计算年度平均影像
接下来,我们定义了一个名为 calculate_yearly_mean 的函数,它用于计算给定文件夹中所有影像文件的年度平均影像。在这个函数中,我们首先读取输入文件夹中的所有影像文件,并创建一个字典来存储每年的影像数据。然后,我们遍历每个影像文件,累加每年的像素值和像素计数。最后,我们计算每年的平均影像,并将结果保存为新的 TIFF 文件。
4. 示例代码
下面是一个示例代码,演示了如何使用上述函数处理遥感影像数据:
# 输入文件夹和输出文件夹
input_folder = "path/to/input/folder"
output_folder = "path/to/output/folder"
# 获取栅格数据
cols = process_tiff_folder(input_folder, output_folder)
# 计算年度平均影像
calculate_yearly_mean(input_folder, output_folder)
5. 完整代码
import os
import numpy as np
from osgeo import gdal
def process_tiff_folder(folder_path, output_folder):
for root, dirs, files in os.walk(folder_path):
for file in files:
if file.endswith(".tif"):
tif_path = os.path.join(root, file)
folder_name = os.path.basename(root) # 获取文件夹名称
dataset = gdal.Open(tif_path) # 打开tif
# 获取行数列数和地理信息
geo_information = dataset.GetGeoTransform()
col = dataset.RasterXSize
row = dataset.RasterYSize
dem = dataset.GetRasterBand(1).ReadAsArray()
# 获取行列数,对应其经纬度,j对于x坐标
cols = []
for y in range(row):
rows = []
for x in range(col):
# 有效高程
if dem[y][x] != dataset.GetRasterBand(1).GetNoDataValue():
# 输出经纬度
lon = geo_information[0] + x * geo_information[1] + y * geo_information[2]
lat = geo_information[3] + x * geo_information[4] + y * geo_information[5]
child = [lon, lat, dem[y][x], y, x]
rows.append(child)
cols.append(rows)
return cols
def calculate_yearly_mean(input_folder, output_folder):
# 获取输入文件夹中的所有影像文件路径
input_files = [os.path.join(input_folder, f) for f in os.listdir(input_folder) if f.endswith('.tif')]
# 创建输出文件夹
if not os.path.exists(output_folder):
os.makedirs(output_folder)
# 初始化年度影像字典
yearly_images = {}
# 遍历所有输入影像文件
for file_path in input_files:
# 从文件名中提取年份和月份
year = int(file_path.split('_')[1])
month = int(file_path.split('_')[2].split('.')[0])
# 读取影像数据
dataset = gdal.Open(file_path)
image = dataset.GetRasterBand(1).ReadAsArray()
# 处理无效值
invalid_value = dataset.GetRasterBand(1).GetNoDataValue()
image[image == invalid_value] = np.nan
# 初始化年份数据字典
if year not in yearly_images:
yearly_images[year] = {'sum': np.zeros(image.shape), 'count': np.zeros(image.shape)}
# 累加每年的像素值和计数
yearly_images[year]['sum'] += np.where(np.isnan(image), 0, image)
yearly_images[year]['count'] += np.where(np.isnan(image), 0, 1)
# 遍历年度影像字典,计算每年的平均影像并保存
for year, data in yearly_images.items():
# 计算每年的平均影像
yearly_mean = np.divide(data['sum'], data['count'], out=np.zeros_like(data['sum']), where=data['count'] != 0)
# 获取输入影像的地理转换信息
dataset = gdal.Open(input_files[0])
geotransform = dataset.GetGeoTransform()
projection = dataset.GetProjection()
# 创建输出影像
driver = gdal.GetDriverByName('GTiff')
output_path = os.path.join(output_folder, f'{year}_mean.tif')
output_dataset = driver.Create(output_path, yearly_mean.shape[1], yearly_mean.shape[0], 1, gdal.GDT_Float32)
output_dataset.SetGeoTransform(geotransform)
output_dataset.SetProjection(projection)
output_dataset.GetRasterBand(1).WriteArray(yearly_mean)
# 关闭输出数据集
output_dataset = None
print("年度平均影像计算完成!")
# 输入文件夹和输出文件夹
input_folder = r"D:\lky\person\month"
output_folder = r"D:\lky\person\month_year"
# 获取栅格数据
cols = process_tiff_folder(input_folder, output_folder)
# 计算年度平均影像
calculate_yearly_mean(input_folder, output_folder)
6. 结论
通过本文介绍的方法,我们可以轻松地处理遥感影像数据,并从中提取出更有意义的信息,如年度平均影像。这些信息对于地球科学研究、自然资源管理以及环境监测等领域具有重要意义,帮助我们更好地理解和应对地球上的变化。
通过利用 Python 编程和相关库,我们可以实现对遥感影像数据的高效处理和分析,为科学研究和实际应用提供强大的工具支持。