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 gdaldef extract_pixel_time_series(input_folder, output_csv):tif_files = [file for file in os.listdir(input_folder) if file.endswith('.tif')]target_row = 495target_col = 60time_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_valuedataset = Nonefor 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