ToB企服应用市场:ToB评测及商务社交产业平台

标题: 批量计算遥感影像NDVI：Python代码 [打印本页]
作者: 去皮卡多    时间: 2024-11-9 14:24
标题: 批量计算遥感影像NDVI：Python代码
  本文介绍基于Python中的gdal模块，批量基于大量多波段遥感影像文件，计算其每1景图像各自的NDVI数值，并将多景效果依次保存为栅格文件的方法。
  如下图所示，现在有大量.tif格式的遥感影像文件，其中均含有红光波段与近红外波段（别的也可以含有其他光谱波段，有没有都不影响）；我们希望，批量计算其每1景遥感影像的NDVI。

  在之前的文章中，我们多次介绍过在不同软件或平台中计算NDVI的方法；而在本文中，我们就介绍一下基于Python中的gdal模块，实现NDVI批量计算的方法。
  这里所需的代码如下。
  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\2021_48STA_Result\Original"
  12. output_folder = r"E:\04_Reconstruction\99_MODIS\new_data\2021_48STA_Result\NDVI_Original"
  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.         driver = gdal.GetDriverByName('GTiff')
  21.         output_dataset = driver.Create(os.path.join(output_folder, "NDVI_" + filename), width, height, 1, gdal.GDT_Float32)
  22.         
  23.         band_red = dataset.GetRasterBand(3)
  24.         data_red = band_red.ReadAsArray()
  25.         data_red = data_red.astype(float)
  26.         band_nir = dataset.GetRasterBand(4)
  27.         data_nir = band_nir.ReadAsArray()
  28.         data_nir = data_nir.astype(float)
  29.         data_ndvi = (data_nir - data_red) / (data_nir + data_red)
  31.         output_band = output_dataset.GetRasterBand(1)
  32.         output_band.WriteArray(data_ndvi)
  33.         output_band.FlushCache()
  34.         output_dataset.SetGeoTransform(dataset.GetGeoTransform())
  35.         output_dataset.SetProjection(dataset.GetProjection())
  37.         dataset = None
  38.         output_dataset = None
  39.         print(filename, "finished!")
复制代码
  代码团体也非常简朴。起首，我们界说输入文件与输入效果文件的路径，前者就是待计算NDVI的遥感影像文件路径，后者则是NDVI效果的遥感影像文件路径。
  接下来，遍历original_folder文件夹中的文件。其中，os.listdir()用于获取文件夹中的文件列表，其后的endswith('.tif')用于筛选出以.tif扩展名结尾的文件。
  随后，对于每个以.tif结尾的文件，起首利用gdal.Open()打开文件——其中的os.path.join()用于构建完备的文件路径；接下来获取影像数据集的宽度和高度，并利用gdal.GetDriverByName()获取GTiff驱动步伐，用于创建输出影像文件；同时，利用driver.Create()创建一个与原始影像具有相同大小的输出影像文件。
  紧接着，从数据集中获取红光和近红外波段的数据。dataset.GetRasterBand()用以获取指定的栅格波段，而band.ReadAsArray()则将波段数据读取为数组；同时，我这里还用了astype()转换数组的格式，克制本来遥感影像的数据格式带来的题目——比方，如果本来遥感影像是无符号整型的数据格式，那么这里不加astype()计算NDVI就会有题目。
  其次，即可计算NDVI。利用获取的红光和近红外波段数据计算NDVI，并将NDVI数据保存在data_ndvi数组中。
  末了，将NDVI数据写入输出影像文件。output_dataset.GetRasterBand()获取输出影像文件的波段，band.WriteArray()将数据写入波段，band.FlushCache()革新波段缓存。
  别的，记得通过output_dataset.SetGeoTransform()和output_dataset.SetProjection()设置输出影像文件的地理变更和投影信息。
  同时，必要清理和关闭数据集，将数据集和输出数据集设置为None以开释资源。还可以打印文件名和finished!，表示当前文件处置惩罚完成。
  实行上述代码，我们即可在效果文件夹中看到计算得到的NDVI数据；如下图所示。

  至此，大功告成。

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



欢迎光临 ToB企服应用市场:ToB评测及商务社交产业平台 (https://dis.qidao123.com/) Powered by Discuz! X3.4