本文介绍基于Python中的gdal模块，批量读取大量多波段遥感影像文件，分别对各波段数据加以数值处理，并将所得处理后数据保存为新的遥感影像文件的方法。
首先，看一下本文的具体需求。我们现有一个文件夹，其中含有大量.tif格式的遥感影像文件；其中，这些遥感影像文件均含有4个波段，每1个波段都表示其各自的反射率数值。而对于这些遥感影像文件，有的文件其各波段数值已经处于0至1的区间内（也就是反射率数据的正常数值区间），而有的文件其各波段数值则是还没有乘上缩放系数的（在本文中，缩放系数是0.0001）。
例如，如下图所示，即为文件夹中某一景遥感影像。可以看到其各波段数值都是大于1的，这是因为其数值都是还没有乘上缩放系数的，即是真实的反射率数值的10000倍。
我们希望实现的是，对于这些遥感影像中，还没有乘上缩放系数0.0001的遥感影像，将其像元值乘上这个缩放系数；而对于已经缩放过（也就是像元数值已经落在0至1区间内）的遥感影像，则不加以任何处理。最后，将经过上述操作后的所有图像（无论是否执行缩放）均保存至指定的输出结果文件夹中。
本文所需代码如下：
# -*- coding: utf-8 -*- """ Created on Thu Apr 18 12:37:22 2024 @author: fkxxgis """ import os from osgeo import gdal original_folder = r"E:\04_Reconstruction\99_MODIS\new_data\GF_Original" output_folder = r"E:\04_Reconstruction\99_MODIS\new_data\GF_Rec" for filename in os.listdir(original_folder): if filename.endswith('.tif'): dataset = gdal.Open(os.path.join(original_folder, filename), gdal.GA_ReadOnly) width = dataset.RasterXSize height = dataset.RasterYSize band_count = dataset.RasterCount driver = gdal.GetDriverByName('GTiff') output_dataset = driver.Create(os.path.join(output_folder, "New_" + filename), width, height, band_count, gdal.GDT_Float32) for band_index in range(1, band_count + 1): band = dataset.GetRasterBand(band_index) data = band.ReadAsArray() if band_index == 1: data = data.astype(float) data[data > 1] /= 10000 elif band_index == 2: data = data.astype(float) data[data > 1] /= 10000 elif band_index == 3: data = data.astype(float) data[data > 1] /= 10000 elif band_index == 4: data = data.astype(float) data[data > 1] /= 10000 output_band = output_dataset.GetRasterBand(band_index) output_band.WriteArray(data) output_band.FlushCache() output_dataset.SetGeoTransform(dataset.GetGeoTransform()) output_dataset.SetProjection(dataset.GetProjection()) dataset = None output_dataset = None
首先，我们使用os.listdir()函数遍历原始数据文件夹中的所有文件，并使用if语句筛选出以.tif结尾的文件；随后，使用gdal.Open()函数打开原始影像数据集，并指定只读模式；接下来，使用dataset.RasterXSize和dataset.RasterYSize获取影像数据集的宽度和高度。
随后，使用dataset.RasterCount获取波段数量，并使用gdal.GetDriverByName()创建输出数据集的驱动程序对象；紧接着，通过Create()方法创建输出数据集，并指定输出文件的路径、宽度、高度、波段数量和数据类型（gdal.GDT_Float32表示浮点型）。
接下来，就可以开始使用循环，对每个文件的每个波段进行处理。首先，使用dataset.GetRasterBand()方法获取当前波段对象，然后使用band.ReadAsArray()将波段数据读取为数组；根据波段索引的不同，对波段数据进行处理。在本文中，对4个波段进行的其实是相同的处理，即将大于1的像素值除以10000。
其次，使用output_dataset.GetRasterBand()方法获取输出数据集中的当前波段对象，并使用output_band.WriteArray()方法将处理后的数据写入输出数据集。
再次，使用dataset.GetGeoTransform()和dataset.GetProjection()分别获取原始数据集的地理转换和投影信息，并使用output_dataset.SetGeoTransform()和output_dataset.SetProjection()设置输出数据集的地理转换和投影信息。
最后一步，关闭数据集对象。至此，代码就完成了对每个.tif文件的处理，并将处理后的数据保存到输出文件夹中。
此时，打开本文开头展示的那1景遥感影像，可以看到其像素数值已经是乘上缩放系数之后的了，也就是落在了0至1的区间内；如下图所示。
至此，大功告成。
作者：疯狂学习GIS
