【PythonRS】GDAL计算遥感影像光谱指数（如NDVI、NDWI、EVI等）
GDAL(Geospatial Data Abstraction Library)是一个在X/MIT许可协议下的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式。它还有一系列命令行工具来进行数据转换和处理。
Python的GDAL库作为栅格数据的处理转换库，其支持几百种栅格数据格式，如常见的TIFF、ENVI、HFA、HDF4等。因为遥感影像大部分都是栅格数据，所以GDAL库非常适合处理遥感影像、如光谱指数计算、波段合成、批量下载等。
一、GDAL库的安装
因为博主的GDAL库的安装有些时候了，所以pip install 的方法能不能安装需要自己尝试。这里推荐大家下载对应包，使用本地安装。
1.打开GDAL库的下载链接，选择自己电脑Python对应的版本，别下载错了哦！
2.然后按住Win+R，输入cmd进入命令行。在cmd终端中，使用cd跳转至下载GDAL包的目录。如 cd /d G:\try\，G:\try\是你狭窄GDAL包的文件夹目录。然后在命令行里输入pip GDAL-3.4.3-cp38-cp38-win_amd64.whl，这是你下载的GDAL包的名称，每个人不一样，别傻乎乎地用我的！
3.可以在编译器里查看是否安装成功，博主使用的是PyCharm，在设置中可以看到是否安装成功。当然你也可以先编写程序，如果没安装成功，程序会报错=。=
二、编写程序
1.导入计算所需库
numpy是数组处理的库，GDAL是我们安装的栅格数据处理库。遥感影像本质上是多维数组，如果这都不知道，建议先百度。
import numpy as np
from osgeo import gdal2.查看栅格数据信息，如宽度、高度、波段数、投影等
当然这一步和我们的光谱指数计算没有太大的关系，但后续计算时需要这里面的相关参数，所以大家可以先了解一下。
ds = gdal.Open(filepath) # 打开数据集dataset
ds_width = ds.RasterXSize # 获取数据宽度
ds_height = ds.RasterYSize # 获取数据高度
ds_bands = ds.RasterCount # 获取波段数
ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数
ds_prj = ds.GetProjection() # 获取投影信息
print("影像的宽度为：" + str(ds_width))
print("影像的高度为：" + str(ds_height))
print("仿射地理变换参数为：" + str(ds_geo))
print("投影坐标系为：" + str(ds_prj))3.以数组的形式打开红波段、近红外波段
每种卫星对应的红波段和近红外波段都不一样、博主这里red是3，nir是4。
ds = gdal.Open(filepath) # 打开数据集dataset
ds_width = ds.RasterXSize # 获取数据宽度
ds_height = ds.RasterYSize # 获取数据高度
ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数
ds_prj = ds.GetProjection() # 获取投影信息
array_red = ds.GetRasterBand(red).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float)
# red是红波段对应的波段数
array_nir = ds.GetRasterBand(nir).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float)
# nir是近红波段对应的波段数
# 以数组的形式读取红波段和近红外波段4.计算NDVI（归一化植被指数）
其他光谱指数类似，指数更换red、nir对应的波段即可。
b1 = array_nir-array_red
b2 = array_nir+array_red
NDVI_data = np.divide(b1, b2, out=np.zeros_like(b1), where=b2 != 0) # 计算NDVI5.保存计算好的NDVI
这里需要先创建一个数据驱动，用于存储数组，类型可以替换（如可以保存为TIFF、ENVI、HFA等不同格式），这里我保存的TIFF格式，但需要注意的是GDAL中的类型命名不同，TIFF格式不能直接写TIFF，而是GTIFF，具体格式要求大家可以参考Python与开源GIS。
driver = gdal.GetDriverByName('GTiff') # 载入数据驱动，用于存储内存中的数组随后我们需要创建一个空数组，给这个数组写入各类属性，如高度、宽度、波段数、投影等，一般我们都使用原始输入数据的属性。
s_result = driver.Create(out_path, ds_width, ds_height, bands=1, eType=gdal.GDT_Float64)
# 创建一个数组，宽高为原始尺寸
ds_result.SetGeoTransform(ds_geo) # 导入仿射地理变换参数
ds_result.SetProjection(ds_prj) # 导入投影信息
ds_result.GetRasterBand(1).SetNoDataValue(9999) # 将无效值设为9999最后，我们再将计算好的NDVI写入数组，删除缓存即可。注意缓存一定要删除，不然NDVI不会写入数组。
ds_result.GetRasterBand(1).WriteArray(NDVI_data) # 将NDVI的计算结果写入数组
del ds_result
# 删除内存中的结果，否则结果不会写入图像中三、完整程序
# -*- coding: utf-8 -*-
"""
@Time ： 2023/3/30 9:11
@Auth ： RS迷途小书童
@File ：NDVI计算.py
@IDE ：PyCharm
"""
import numpy as np
from osgeo import gdaldef Get_data(filepath):ds = gdal.Open(filepath) # 打开数据集datasetds_width = ds.RasterXSize # 获取数据宽度ds_height = ds.RasterYSize # 获取数据高度ds_bands = ds.RasterCount # 获取波段数ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数ds_prj = ds.GetProjection() # 获取投影信息print("影像的宽度为：" + str(ds_width))print("影像的高度为：" + str(ds_height))print("仿射地理变换参数为：" + str(ds_geo))print("投影坐标系为：" + str(ds_prj))# data = ds.ReadAsArray(0, 0, ds_width, ds_height) # 以数组的形式读取整个数据集def Write_NDVI(filepath, out_path, red, nir):print("-----正在进行归一化植被指数NDVI计算-----")ds = gdal.Open(filepath) # 打开数据集datasetds_width = ds.RasterXSize # 获取数据宽度ds_height = ds.RasterYSize # 获取数据高度ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数ds_prj = ds.GetProjection() # 获取投影信息array_red = ds.GetRasterBand(red).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float)# red是红波段对应的波段数array_nir = ds.GetRasterBand(nir).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float)# nir是近红波段对应的波段数# 以数组的形式读取红波段和近红外波段b1 = array_nir-array_redb2 = array_nir+array_redNDVI_data = np.divide(b1, b2, out=np.zeros_like(b1), where=b2 != 0) # 计算NDVIdriver = gdal.GetDriverByName('GTiff') # 载入数据驱动，用于存储内存中的数组ds_result = driver.Create(out_path, ds_width, ds_height, bands=1, eType=gdal.GDT_Float64)# 创建一个数组，宽高为原始尺寸ds_result.SetGeoTransform(ds_geo) # 导入仿射地理变换参数ds_result.SetProjection(ds_prj) # 导入投影信息ds_result.GetRasterBand(1).SetNoDataValue(9999) # 将无效值设为9999ds_result.GetRasterBand(1).WriteArray(NDVI_data) # 将NDVI的计算结果写入数组del ds_result# 删除内存中的结果，否则结果不会写入图像中print("计算完成")if __name__ == "__main__":file_path = r"G:/Shanghai_TreeClassicfication/Sen2Cor_Shanghai/Shanghai.dat"# 输入的栅格数据路径out_path1 = r"G:/GDAL_try/Shanghai-1.tif"# 导出的文件路径Get_data(file_path) # 执行函数，获取影像基本信息print("\n")print("--------------NDVI计算--------------")red1 = int(input("请输入红波段："))nir1 = int(input("请输入近红外波段："))Write_NDVI(file_path, out_path1, red1, nir1) # 执行函数，计算NDVI
使用Python去处理遥感影像可以大大减少我们的工作时间，提升工作效率。可能在一幅遥感影像处理时感觉还没有ENVI计算NDVI的速度快，但实际上如果你需要研究长时间序列的NDVI或者大量卫星影像，你就会发现编程给你带来的便利不是一点半点。针对上面分享的程序，只需加入os.listdir和for循环即可遍历文件夹中所有的栅格数据，成百上千幅影像在后台一会就处理好了。同时如果有需要，还可以直接将程序作为函数调用，用在完整的遥感影像预处理、定量分析中。本文章仅作为学习使用，如有侵权请联系作者删除。
如果大家在学习Python或者RS时有什么问题，可以随时留言交流！如果大家对批量处理有兴趣同样可以留言给博主，博主会分享相关代码以供学习！
本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！