Python自动缩放或扩增栅格文件的多波段数据

  本文介绍基于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的区间内；如下图所示。

  至此，大功告成。