Python代码将大量遥感数据的值缩放指定倍数的方法

  本文介绍基于Python中的gdal模块，批量读取大量多波段遥感影像文件，分别对各波段数据加以数值处理，并将所得处理后数据保存为新的遥感影像文件的方法。
  起首，看一下本文的具体需求。我们现有一个文件夹，此中含有大量.tif格式的遥感影像文件；此中，这些遥感影像文件均含有4个波段，每1个波段都表示其各自的反射率数值。而对于这些遥感影像文件，有的文件其各波段数值已经处于0至1的区间内（也就是反射率数据的正常数值区间），而有的文件其各波段数值则是还没有乘上缩放系数的（在本文中，缩放系数是0.0001）。
  比方，如下图所示，即为文件夹中某一景遥感影像。可以看到其各波段数值都是大于1的，这是因为其数值都是还没有乘上缩放系数的，即是真实的反射率数值的10000倍。

  我们希望实现的是，对于这些遥感影像中，还没有乘上缩放系数0.0001的遥感影像，将其像元值乘上这个缩放系数；而对于已经缩放过（也就是像元数值已经落在0至1区间内）的遥感影像，则不加以任那边理。最后，将经过上述操纵后的所有图像（无论是否执行缩放）均保存至指定的输出结果文件夹中。
  本文所需代码如下。

1. # -*- coding: utf-8 -*-
2. """
3. Created on Thu Apr 18 12:37:22 2024
5. @author: fkxxgis
6. """
8. import os
9. from osgeo import gdal
11. original_folder = r"E:\04_Reconstruction\99_MODIS\new_data\GF_Original"
12. output_folder = r"E:\04_Reconstruction\99_MODIS\new_data\GF_Rec"
14. for filename in os.listdir(original_folder):
15.     if filename.endswith('.tif'):
16.         dataset = gdal.Open(os.path.join(original_folder, filename), gdal.GA_ReadOnly)
17.         width = dataset.RasterXSize
18.         height = dataset.RasterYSize
19.         
20.         band_count = dataset.RasterCount
21.         driver = gdal.GetDriverByName('GTiff')
22.         output_dataset = driver.Create(os.path.join(output_folder, "New_" + filename), width, height, band_count, gdal.GDT_Float32)
23.         
24.         for band_index in range(1, band_count + 1):
25.             band = dataset.GetRasterBand(band_index)
26.             data = band.ReadAsArray()
27.             
28.             if band_index == 1:
29.                 data = data.astype(float)
30.                 data[data > 1] /= 10000
31.             elif band_index == 2:
32.                 data = data.astype(float)
33.                 data[data > 1] /= 10000
34.             elif band_index == 3:
35.                 data = data.astype(float)
36.                 data[data > 1] /= 10000
37.             elif band_index == 4:
38.                 data = data.astype(float)
39.                 data[data > 1] /= 10000
41.             output_band = output_dataset.GetRasterBand(band_index)
42.             output_band.WriteArray(data)
43.             output_band.FlushCache()
45.         output_dataset.SetGeoTransform(dataset.GetGeoTransform())
46.         output_dataset.SetProjection(dataset.GetProjection())
48.         dataset = None
49.         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的区间内；如下图所示。

  至此，大功告成。

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。