利用 Python 处理遥感影像数据：计算年度平均影像

分手后的思念是犯贱 2024-04-17 06:17 102阅读 0赞

在地球科学、气象学以及环境监测等领域，遥感影像数据是一种重要的信息源，它们可以提供地表的地形、植被覆盖、气候变化等丰富信息。然而，随着观测技术的进步，我们通常会获得大量的遥感影像数据，如何高效地处理和分析这些数据成为了一项挑战。本文将介绍如何利用 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 编程和相关库，我们可以实现对遥感影像数据的高效处理和分析，为科学研究和实际应用提供强大的工具支持。