执行地图代数表达

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() 函数）。

image20image21image22

海拔> 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 。

14.8.2. 植被指数计算#

目标是根据Landsat 8多光谱图像计算归一化差异植被指数（NDVI） [ROUSE74]， [TUCKER79] 。 该指数的计算方法如下：

NIR – R)/(NIR + R)

其中NIR是近红外线，R是红色。

[ROUSE74]

ROUSE JR J., HAAS R.H., SCHELL J.A. 等人，“用ERTS监测大平原植被系统”， 第三届地球资源技术卫星研讨会论文集，第1卷，第309–313页，1974年。

[TUCKER79]

“用于监测植被的红色和摄影红外线性组合”，环境遥感，第8卷，第2期，第127-150页，1979年。

QGIS处理步骤#

1.打开Landsat 8图像，在QGIS中，仅打开以下文件"LC81990262017019LGN00_MS_L93.TIF"。 提醒您，此文件的组织方式如下：

– 波段1：蓝色，

– 波段2：绿色，

– 波段3：红色，

– 波段4：近红外。

  1. NDVI计算

在 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之间的小数值的栅格编码。 图 14.28 。高值通常对应于活跃植被，图 14.29 。

image74_xay

图 14.28 红外假色#

image75_xt3

图 14.29 NDVI#

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
有关更多信息，请参阅以下网页：

http://www.gdal.org/gdal_calc.html 。

目录