Python读取多个栅格文件并提取像元的各波段时间序列数据与变化值

  本文介绍基于Python语言，读取文件夹下大量栅格遥感影像文件，并基于给定的一个像元，提取该像元对应的全部遥感影像文件中，指定多个波段的数值；修改其中不在给定范围内的异常值，并计算像元数值在每一景遥感影像中变化的差值；最终将这些数据保存为一个新的Excel表格文件的方法。

  首先，我们来看一下本文需要具体实现的需求。现在有一个文件夹，如下图所示；其中，存放了大量的遥感影像文件，且每一景遥感影像都是同一个空间位置、不同成像时间对应的遥感影像，因此其空间参考信息、栅格的行数与列数等都是一致的。此外，每一景遥感影像都具有5个不同的波段。

  我们现在希望，给定一个像元（也就是给定了这个像元在遥感影像中的行号与列号），提取出在指定的波段中（我们这里就提取全部的5个波段），该像元对应的每一景遥感影像的数值（也就是提取了该像元在每一景遥感影像、每一个波段的数值）；随后，将提取到的大于1的数值修改为1，并计算像素值在每一景遥感影像中数值的差值；最后，将提取到的数据保存为一个Excel表格文件。

  明确了需求，我们就可以撰写代码；具体如下。

# -*- coding: utf-8 -*-
"""
Created on Thu Jul 27 11:25:55 2023
@author: fkxxgis
"""
import os
import pandas as pd
from osgeo import gdal
def extract_pixel_time_series(input_folder, output_csv):
    tif_files = [file for file in os.listdir(input_folder) if file.endswith('.tif')]
    target_row = 495
    target_col = 60
    time_series_df = pd.DataFrame()
    for tif_file in tif_files:
        file_path = os.path.join(input_folder, tif_file)
        dataset = gdal.Open(file_path)
        for band in range(dataset.RasterCount):
            band_data = dataset.GetRasterBand(band + 1).ReadAsArray()
            pixel_value = band_data[target_row, target_col]
            date = tif_file[10: 24]
            time_series_df.at[date, f'Band_{band + 1}'] = pixel_value
        dataset = None
    for index in range(len(time_series_df.columns)):
        time_series_df = time_series_df.apply(lambda x: x.clip(upper = 1))
        new_col_name = time_series_df.columns[index] + "_diff"
        time_series_df[new_col_name] = time_series_df.iloc[:, index].diff()
    time_series_df.to_csv(output_csv)
# 示例用法
input_folder = r"E:\01_Reflectivity\FiveBands"
output_csv = r"E:\01_Reflectivity\Data.csv"
extract_pixel_time_series(input_folder, output_csv)

  首先，我们需要导入必要的模块和库。其中os用于操作文件和文件夹，pandas用于处理数据和创建DataFrame格式数据，而gdal则用于读取栅格数据；关于gdal库的配置方法，大家可以参考文章Anaconda环境GDAL库基于whl文件的配置方法（https://blog.csdn.net/zhebushibiaoshifu/article/details/128320388）。

  随后，我们对extract_pixel_time_series这个函数加以定义。这个函数接收两个参数input_folder和output_csv，分别表示存储栅格数据的文件夹路径和输出的Excel文件的路径。随后，列出input_folder文件夹下所有以.tif结尾的文件，并存储在列表中。其次，循环遍历每个栅格文件，构建完整的文件路径，用于后面的数据读取，并使用gdal.Open()打开栅格文件，获取数据集对象。

  接下来，通过循环遍历每个波段。读取当前波段的数据，并存储在band_data变量中。随后基于我们给定的像元位置，提取目标像元的数值（位置就是这个[target_row, target_col]）。此外，为了使得我们保存结果时可以记录每一个数值对应的成像日期，因此需要从文件名中提取日期，并存储在date变量中。

  接下来，通过time_series_df.at[date, f'Band_{band + 1}']，将像元值存储在DataFrame中，行索引为日期，列名为Band_1、Band_2等；随后，将数据集对象dataset设为None，释放内存资源。

  接下来，我们将大于1的数值加以处理，并计算每个波段随时间变化的数值之差。遍历time_series_df的每一列，并对于每一列使用clip(upper=1)将超过1的值截断为1；随后，为每一列创建新列，列名为原列名加上_diff，存储该列差值。

  最后，我们将处理后的时间序列数据保存为Excel表格文件即可。

  运行上述代码，我们即可获得多个遥感影像文件中，给定像元位置处，像元数值的时间变化序列，并可以获得其变化值。

  至此，大功告成。

欢迎关注：疯狂学习GIS