采用GDAL批量波段运算计算植被指数0基础教程
1. 引言
在传统的遥感数据处理方法中，通常使用ArcGis或ENVI软件进行波段运算。然而，这些软件在处理大量数据时往往效率低下。有没有一种方法可以批量进行波段运算，一下子计算几十个植被指数，并自动保存计算后的指标影像到相应目录下呢？答案是肯定的！本文将介绍如何使用GDAL库实现这一目标。
2. GDAL简介
GDAL（Geospatial Data Abstraction Library）是一个用于处理栅格数据和矢量数据的开源库，是地理信息系统（GIS）和遥感领域的重要工具之一。GDAL的详细介绍及其安装可以查看这篇文章：GDAL介绍及安装。
3. 植被指数计算公式
这里列出了要计算的23种植被指数及其计算公式：
- DVI 为差值植被指数；
- EVI 为增强植被指数；
- GCVI 为绿光归一化差值植被指数；
- GNDVI 为绿光归一化差值植被指数；
- ISR1, ISR2 为短波红外比值植被指数；
- LSWI 为地表水分指数；
- MSAVI 为修正型土壤调节植被指数；
- NDPI 为归一化差异物候指数；
- NDVI 为归一化植被指数；
- NDWI1, NDWI2, NDWI3 为归一化水体指数；
- RENDVI1, RENDVI2 为红边归一化差值植被指数；
- RENDWI1, RENDWI2 为红边归一化水体指数；
- SAVI 为土壤调节植被指数；
- SR 为比值植被指数；
- TR1, TR2 为比值生物学变量植被指数；
- VRE1, VRE2 为红边植被指数。
| 植被指数 | 计算公式 |
|---|---|
| DVI | DVI = nir - red |
| EVI | EVI = 2.5 × (nir - red)/(nir + 6 × red - 7.5 × blue + 1) |
| GCVI | GCVI = nir/green - 1 |
| GNDVI | GNDVI = (nir - green)/(nir + green) |
| ISR1 | ISR1 = nir/swir1 |
| ISR2 | ISR2 = nir/swir2 |
| LSWI | LSWI = (nir - swir1)/(nir + swir1) |
| MSAVI | MSAVI = 1/2 × (2 × nir + 1 - ((2 × nir + 1)^2 - 8 × (nir - red))^(1/2)) |
| NDPI | NDPI = (nir - (0.74 × red + 0.26 × swir1))/(nir + (0.74 × red + 0.26 × swir1)) |
| NDVI | NDVI = (nir - red)/(nir + red) |
| NDWI1 | NDWI1 = (nir - swir1)/(nir + swir1) |
| NDWI2 | NDWI2 = (nir - swir2)/(nir + swir2) |
| NDWI3 | NDWI3 = (green - nir)/(green + nir) |
| RENDVI1 | RENDVI1 = (nir - red1)/(nir + red1) |
| RENDVI2 | RENDVI2 = (nir - red2)/(nir + red2) |
| RENDWI1 | RENDWI1 = (red1 - swir1)/(red1 + swir1) |
| RENDWI2 | RENDWI2 = (red2 - swir1)/(red2 + swir1) |
| SAVI | SAVI = (1 + L) × (nir - red)/(nir + red + 0.5) |
| SR | SR = nir/red |
| TR1 | TR1 = (1/red1 - 1/red1) × nir |
| TR2 | TR2 = (1/red - 1/nir) × red2 |
| VRE1 | VRE1 = nir/red1 |
| VRE2 | VRE2 = nir/red2 |
4. 批量计算植被指数
使用 gdal_calc.py 工具可以批量计算表格中各植被指数。以下是计算这些植被指数的命令示例：
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=4 --outfile=dvi.tif --calc="A-B"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=4 -C ronghe.tif --C_band=2 --outfile=evi.tif --calc="2.5*((A-B)/(A+6*B-7.5*C+1))"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=3 --outfile=gcvi.tif --calc="(A/B)-1"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=3 --outfile=gndvi.tif --calc="(A-B)/(A+B)"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=6 --outfile=isr1.tif --calc="A/B"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=7 --outfile=isr2.tif --calc="A/B"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=6 --outfile=lswi.tif --calc="(A-B)/(A+B)"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=4 --outfile=msavi.tif --calc="0.5*(2*A+1-((2*A+1)**2-8*(A-B))**0.5)"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=4 -C ronghe.tif --C_band=6 --outfile=ndpi.tif --calc="(A-(0.74*B+0.26*C))/(A+(0.74*B+0.26*C))"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=4 --outfile=ndvi.tif --calc="(A-B)/(A+B)"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=6 --outfile=ndwi1.tif --calc="(A-B)/(A+B)"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=7 --outfile=ndwi2.tif --calc="(A-B)/(A+B)"
python gdal_calc.py -A ronghe.tif --A_band=3 -B ronghe.tif --B_band=5 --outfile=ndwi3.tif --calc="(A-B)/(A+B)"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=4 --outfile=rendvi1.tif --calc="(A-B)/(A+B)"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=4 --outfile=rendvi2.tif --calc="(A-B)/(A+B)"
python gdal_calc.py -A ronghe.tif --A_band=4 -B ronghe.tif --B_band=6 --outfile=rendwi1.tif --calc="(A-B)/(A+B)"
python gdal_calc.py -A ronghe.tif --A_band=4 -B ronghe.tif --B_band=6 --outfile=rendwi2.tif --calc="(A-B)/(A+B)"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=4 --outfile=savi.tif --calc="(1+0.5)*(A-B)/(A+B+0.5)"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=4 --outfile=sr.tif --calc="A/B"
python gdal_calc.py -A ronghe.tif --A_band=4 -B ronghe.tif --B_band=4 --C ronghe.tif --C_band=5 --outfile=tr1.tif --calc="(1/A-1/B)*C"
python gdal_calc.py -A ronghe.tif --A_band=4 -B ronghe.tif --B_band=5 -C ronghe.tif --C_band=4 --outfile=tr2.tif --calc="(1/A-1/B)*C"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=4 --outfile=vre1.tif --calc="A/B"
python gdal_calc.py -A ronghe.tif --A_band=5 -B ronghe.tif --B_band=4 --outfile=vre2.tif --calc="A/B"
这里如何快速生成这些指令呢？我将公式截图直接发给给GPT, 然后给了它一个示范指令如何构建，让其构建所有公式的指令，只需要简单对一下即可，一般不会错，省去了手敲指令的过程。
这些命令使用 gdal_calc.py 工具来计算表格中的各个植被指数，其中 A_band、B_band 和 C_band 分别对应 Landsat 8多光谱，这里我们将其改名成ronghe.tif (你可以改为你的影像名字)的不同波段。请根据实际数据文件调整命令中的文件名和波段编号。
5. 批处理脚本自动计算
手动运行每条命令会非常麻烦，因此可以将这些命令写入一个批处理文件（.bat），实现自动计算。
创建批处理文件
-
打开记事本，将上述命令内容复制进去。
-
保存文件，将文件名命名为
run.bat。 -
将gdal_calc.py文件和遥感影像文件放入同一文件夹下。
gdal.py文件已经上传到百度网盘，需要可搜索关注作者微信公众号:Python与遥感 获取
运行批处理文件
- 打开Anaconda Prompt。
- 输入以下命令进入相应的驱动器和目录：
F: cd path\to\your\folder - 输入并运行批处理文件：
run.bat
运行后，批处理文件将自动依次执行每个命令，计算所有植被指数，并将结果保存到相应的文件中。
6. 结果展示
每条命令运行后，对应的植被指数会保存在文件夹中，黑窗里的数字到99就代表计算完了一个植被指数。只需等待一小会儿，即可计算完所有的植被指数计算。
7. 总结
通过本文介绍的方法，使用GDAL库可以高效地批量计算植被指数，大大提高了工作效率。批处理文件的使用使得整个过程更加自动化和简便。希望这篇0基础教程能帮助大家更好地掌握GDAL的使用方法。如果有任何问题或需要进一步的帮助，请在评论区留言。
欢迎关注我的公众号：python与遥感，或许里面有你想要的免费教程。
