14.8. 执行地图代数表达¶
QGIS软件及其API提供了一个栅格计算器。GDAL 通过 gdal_calc.py 命令具有类似的功能，
该命令仅在 Processing Toolbox 中接口。图像必须具有相同的维度（行和列中的像素数）才能执行数学运算。
对于不知情的用户来说，用 Pythonnumpy 语法编写的表达输入并不总是直观的。
这里提供了两个例子：对海拔和斜坡数据的查询，以及根据Landsat 8多光谱图像计算植被指数。
14.8.1. 执行多准则查询¶
目标是划定海拔和斜坡分别大于400 m和20°的区域。
QGIS处理步骤
1.打开 DTM
和斜坡栅格
在QGIS中：
- 仅打开以下文件：
srtm_37_03_l93.tif
srtm_37_03_l93_slope.tif
2.处理工具箱激活在菜单栏中：
单击处理，
单击“Toolbox”。
“Processing Toolbox”面板出现在屏幕右侧。
3.在“Processing Toolbox”中查询海拔和斜坡数据：
双击 GDAL /OGR > Miscellaneous > Raster Calculator
在“栅格计算器”窗口中：
设置如下：
输入图层 A: srtm_37_03_l93- Input layer B: srtm_37_03_l93_slope
gdalnumeric语法中的计算: logical_and(A > 400,B > 20)
设置输出nodata值： 0
输出光栅类型：字节
计算： alti-sup400_slope-sup20.tif
单击“OK”
4.结果分析
结果是一个二进制栅格数据，其中满足查询的像素被编码为1（白色），其他像素被编码为0（黑色）。 数学表达必须尊重gdalnumeric语法， 例如通过使用+ -/* 运算符或Python numpy库的任何函数（例如 thelogical_and() 函数）。
海拔> 400 m斜坡> 20°结果
6.相应的GDAL 命令
上面使用的工具相当于在终端中启动以下命令：
gdal_calc --calc "logical_and(A>400,B>20)" --format GTiff --type Byte
--NoDataValue 0
-A srtm_37_03_l93.tif --A_band 1 -B srtm_37_03_l93_pente.tif
--B_band 1
--outfile alti_sup400_pente_sup20.tif
其中：
[--calc]: 数学表达式，
[--format]: 设置输出格式，
[--type]: 设置图像编码，
[--NoDataValue]: 设置输出的nodata值，
[--A]: 第一个输入栅格（此处为srtm_37_03_l93.tif），
[-A_band]: 文件A的栅格带数，
[-B]: 第二个输入光栅（此处为srtm_37_03_l93_slope.tif），
[--B_band]: 文件B的栅格带数，
[--outfile]: 设置要生成的输出文件。
有关更多信息，请参阅以下网页： http://www.gdal.org/gdal_calc.html 。
表2.18.使用GDAL 栅格计算器进行多标准查询
14.8.2. 植被指数计算¶
目标是根据Landsat 8多光谱图像计算归一化差异植被指数（NDVI） [ROUSE74] ， [TUCKER79] 。 该指数的计算方法如下：
NIR – R)/(NIR + R)
其中NIR是近红外线，R是红色。
ROUSE JR J., HAAS R.H., SCHELL J.A. 等人，“用ERTS监测大平原植被系统”， 第三届地球资源技术卫星研讨会论文集，第1卷，第309–313页，1974年。
“用于监测植被的红色和摄影红外线性组合”，环境遥感，第8卷，第2期，第127-150页，1979年。
QGIS处理步骤
1.打开Landsat 8图像
在QGIS中：
仅打开以下文件：
LC81990262017019LGN00_MS_L93.TIF
提醒您，此文件的组织方式如下：
– 波段1：蓝色，
– 波段2：绿色，
– 波段3：红色，
– 波段4：近红外。
NDVI计算
In the :
在“Processing Toolbox”中：
双击 GDAL /OGR > Miscellaneous > Raster Calculator
在“Raster Calculator”窗口中：
设置如下：
输入图层 A: LC81990262017019LGN00_MS_L93
栅格A的栅格带数量: 4
输入图层 B: LC81990262017019LGN00_MS_L93
栅格B的栅格带数量： 3
gdalnumeric语法中的计算：
(A.astype(float) - B)/(A.astype(float) + B)
（A.astype（float）- B）/（A.astype（float）+ B）
注意：在上面的公式中，必须使用 .astype (float) 函数将两个带中的至少一个输入为 float ，
才能获得小数值的结果。
Set output nodata value: -9999
Output raster type: Float32
Calculated: LC81990262017019LGN00_NDVI_L93.TIF
单击“OK”
3.结果分析
NDVI是一个以-1到1之间的小数值的栅格编码。高值通常对应于活跃植被。
红外假色
NDVI (RGB = NIR, R, G)
RDX（RB =近红外线、R、G）
4.对应
GDAL命令上面使用的工具相当于在终端中启动以下命令：
gdal_calc --calc "(A.astype(float)-B)/(A.astype(float)+B) " --format
GTiff --type Float32
--NoDataValue -9999 -A LC81990262017019LGN00_MS_L93.tif
--A_band 4
-B LC81990262017019LGN00_MS_L93.tif --B_band 3
--outfile LC81990262017019LGN00_NDVI_L93.tif
有关更多信息，请参阅以下网页： mhttp://www.gdal.org/gdal_calc.html 。
表2.19.使用GDAL 栅格计算器计算NDVI