跳转到内容
AI Agent
QGIS
教程

QGIS 栅格计算器教程:实测 NDVI 公式与 4 个不报错的陷阱

发布于

从 Raster > Raster Calculator(栅格 > 栅格计算器)打开 QGIS 栅格计算器,也可以从 Processing Toolbox 打开。计算 Sentinel-2 NDVI 时使用 ("B08@1" - "B04@1") / ("B08@1" + "B04@1"):菜单里的计算器会自动以 Float32 计算。需要留意三件事:GDAL 计算器会让无符号整数溢出回绕,* 掩膜写入的是 0 而不是 nodata,坡度工具在经纬度网格上会失效。

2026 年 11 月 7 日,我在 QGIS 3.44.15 中逐一测试了这些情况。测试区是巴西米纳斯吉拉斯州 Machado 以南一块 5 km × 5 km 的丘陵农田(2,500 公顷,250,000 个 10 m 像素,UTM 23S,EPSG:32723)。数据是 2026 年 11 月 3 日的 Sentinel-2 L2A 影像和 30 m 的 Copernicus GLO-30 DEM,每个输出都在 QGIS 之外用 64 位浮点 numpy 参考结果做了核对。

QGIS 栅格计算器是什么,你打开的是哪一个?

栅格计算器会在一个或多个栅格波段上逐像素计算一个数学或布尔表达式,然后写出一个新网格。在 QGIS 中,波段写作 "layername@band"(例如 "B08@1" 表示图层 B08 的第 1 波段),可以用 +、-、*、/、<、>、AND、OR 或 if(condition, true, false) 组合波段,结果可以写成 GeoTIFF,也可以写成内存中的虚拟图层。

容易让人困惑的是,QGIS 自带五个都叫 "Raster calculator" 的工具,而且它们的计算引擎并不相同。其中四个使用 QGIS 原生表达式解析器,以浮点数计算。GDAL 那个(底层是 gdal_calc.py)使用 A/B/C 字母变量,并按输入波段自身的数据类型计算 numpy 表达式。

QGIS 3.44 中的工具波段写法输出类型uint16 NDVI 中偏差 > 0.01 的像素(共 250,000 个)
Raster > Raster Calculator (QgsRasterCalculator)"B08@1"Float320
Processing: Raster calculator (native:rastercalc)"B08@1"Float320
Processing: Raster calculator (virtual) (native:virtualrastercalc)"B08@1"Float64(内存中)0
Processing: Raster calculator (legacy) (qgis:rastercalculator)"B08@1"Float320
Processing: GDAL Raster calculator (gdal:rastercalculator)A、BFloat32 文件,按输入类型计算703(0.28%,数值最高达 126.7)

如果只想要一个不用转换类型就能用的公式,请使用顶部菜单的 Raster > Raster Calculator,或 Processing 中 QGIS 原生分组下的第一个 Raster calculator。

分步教程:在 QGIS 栅格计算器中计算 NDVI

如果还需要先获取 Sentinel-2 波段,可以看我们的 在 QGIS 中下载 Sentinel-2 影像 指南,里面讲了 STAC 目录的方法。等 B04(红光,10 m)和 B08(近红外,10 m)出现在 Layers 面板后,计算 NDVI 衡量的是什么 只需四步。

选中决定输出网格的图层

打开对话框之前,先在 Layers 面板中点击 B08。Raster > Raster Calculator 打开时,会从 Layers 面板中当前选中的栅格图层读取输出范围、坐标系(CRS)、宽度和高度。确认标题显示 Extent (current: B08)。

输入 NDVI 表达式

在 Raster > Raster Calculator 中,双击 Raster Bands 列表里的 B08@1 和 B04@1(或者用双引号手动输入),组成:

("B08@1" - "B04@1") / ("B08@1" + "B04@1")

确认表达式框正下方的提示是 Expression valid。在这个对话框里,不需要把 Sentinel-2 的无符号 16 位波段转换成浮点数。

选择输出文件或虚拟图层

把 Output layer 设为一个 .tif 路径,Output format 保持 GeoTIFF。如果只想快速预览、不想在磁盘上生成文件,就勾选 Create on-the-fly raster instead of writing layer to disk(它会在项目文件中保存一个公式 URI,其中包含输入数据的绝对路径,结果与写出的文件相差在 3e-8 以内)。

运行并设置结果样式

点击 OK。在我们 500 × 500 像素的窗口上,计算用时 0.12 秒,数值范围从 -0.5507 到 0.9547,与 64 位 numpy 的最大差值为 2.98e-8。打开 Layer Properties > Symbology(图层属性 > 符号化),切换到 Singleband pseudocolor,选择 RdYlGn 色带,最小值设为 0,最大值设为 0.9。

下面是点击 OK 之前的对话框,范围取自 B08:

QGIS Raster Calculator 对话框:Raster Bands 列表中有 B04@1 和 B08@1,Output layer 路径已填写,Extent (current: B08) 为 EPSG:32723 下的 500 列 500 行,NDVI 公式下方显示 Expression valid。
Raster > Raster Calculator 对话框的初始范围、CRS 和像素尺寸取自 Layers 面板中选中的图层。

用 RdYlGn 设置样式后,Machado 周围茂密的植被显示为深绿色,裸土和道路显示为红色,红光高于近红外的 703 个像素(水体和暗斑)变成负值。整个窗口的 NDVI 均值为 0.727。如果你的冠层茂密到让 NDVI 饱和,或者行间露出裸土,可以看看我们对 NDVI 与 EVI、SAVI 和 NDRE 的比较,了解公式会怎样变化。

QGIS 地图画布显示 5 km × 5 km 的 Machado 农田 NDVI 栅格,使用红黄绿色带从 0 到 0.9 设置样式,左侧是 Layers 面板。
Machado 的 Float32 NDVI 输出,从 0(红)到 0.9(绿)设置样式。每个像素与外部 64 位参考结果的差值都在 3e-8 以内。

QGIS 栅格计算器公式大全

下面每个表达式都在 QGIS 3.44.15(QgsRasterCalculator)中用 Machado 数据集运行过,并与 numpy 逐像素核对。复制时请保留 "layer@band" 两侧的双引号。

任务实测表达式(Raster > Raster Calculator)在 Machado 测试窗口上的结果
NDVI(Sentinel-2 10 m)("B08@1" - "B04@1") / ("B08@1" + "B04@1")Float32,250,000 个像素的范围为 -0.5507 到 0.9547;最大误差 2.98e-8。
NDWI(McFeeters 水体指数)("B03@1" - "B08@1") / ("B03@1" + "B08@1")Float32;"ndwi@1" > 0 提取出 637 个水体像素(6.4 公顷池塘和溪流)。
把低于阈值的像素设为 nodata"ndvi@1" / ("ndvi@1" > 0.3)条件为假时除以 0,写入真正的 nodata(排除 11,343 个像素;均值 0.753)。
条件替换(if)if("ndvi@1" > 0.3, "ndvi@1", 0)能通过解析,保留 0.3 以上的值,但其他地方写入 0,与 * 掩膜相同。
DEM 重分类为 5 类(布尔求和)("dem@1" < 800) * 1 + ("dem@1" >= 800 AND "dem@1" < 850) * 2 + ("dem@1" >= 850 AND "dem@1" < 900) * 3 + ("dem@1" >= 900 AND "dem@1" < 1000) * 4 + ("dem@1" >= 1000) * 554,901 个 DEM 像素得到 1 到 5 的编码;与 native:reclassifybytable 逐像素完全一致。
DEM 重分类为 5 类(嵌套 if)if("dem@1" < 800, 1, if("dem@1" < 850, 2, if("dem@1" < 900, 3, if("dem@1" < 1000, 4, 5))))同样的 5 类(96.9 公顷、1,089.4 公顷、1,323.1 公顷、1,759.7 公顷、582.7 公顷)。
坡度阈值"slope@1" > 15大于 15 度写 1,小于等于写 0(52,979 个有效 UTM DEM 像素中有 17,327 个,32.61%)。
高程米转英尺"dem@1" * 3.28084如果用 QGIS 手册里的 3.28 而不是 3.28084,最大误差为 0.98 ft(1,168 m 处为 0.026%)。

四个会给出错误数值的无声陷阱(以及如何发现)

下面四种错误在 QGIS 中都不会弹出任何警告。每一种都会正常运行完毕,往项目里添加一个图层,给你一组看起来合理的数字,直到你去核对。

1. GDAL 栅格计算器在负 NDVI 上让无符号整数回绕

QGIS 原生栅格计算器在 Sentinel-2 波段上不会遇到整数除法问题,但 Processing 中的 GDAL Raster calculator(gdal:rastercalculator)会按输入文件自身的 numpy 数据类型计算 (A-B)/(A+B),于是红光高于近红外的每个像素都会被算错。

Sentinel-2 Level-2A 波段以无符号 16 位整数(uint16)发布。在我们 250,000 个像素的窗口中,有 703 个像素(0.28%,水体和暗斑)的 B04 大于 B08。当 gdal:rastercalculator 在 uint16 下用 B08 减去 B04 时,负数结果会绕回 65,536,工具把这 703 个像素以 0 到 126.7 之间的值写进 Float32 文件,而不是负的 NDVI。

QGIS 地图显示 GDAL 栅格计算器在 uint16 波段上算出的 NDVI,品红色像素标出红光高于近红外、数值大于 1 的位置。
GDAL 栅格计算器直接作用于原始 uint16 波段。品红色表示大于 1 的值:本应为负的像素绕回了 65,536,最高达到 126.7。

如果要用 GDAL 计算器,要么在表达式里用 (A.astype(numpy.float32)-B)/(A+B) 做类型转换(在我们的测试中错误像素降到了 0),要么改用 QGIS 原生栅格计算器。另外,计算指数时千万不要把 GDAL 输出类型设为 Int16 或 Byte:强制 Int16 会把我们 250,000 个像素中的 249,973 个(99.99%)四舍五入成 -1、0 或 1。

2. 乘以条件(*)写入的是 0 而不是 nodata,会拉低均值

把栅格乘以一个布尔条件,例如 ("ndvi@1" > 0.3) * "ndvi@1",会把每个被排除的像素设为 0 而不是 nodata,这让保留区域的 NDVI 均值从 0.753 降到了 0.719。

QGIS 官方手册把 ("elevation@1" >= 0) * "elevation@1" 作为掩膜栅格的标准写法。在我们的 Machado NDVI 图层上,有 11,343 个像素小于等于 0.3。乘以 ("ndvi@1" > 0.3) 后,这 11,343 个像素全部被写成有效数值 0.0,于是图层统计仍然计入全部 250,000 个像素,你接下来做的每个均值、直方图和分区统计里都会混进 11,343 个零。

掩膜表达式被排除的像素变成图层统计计入的像素保留区域的 NDVI 均值
("ndvi@1" > 0.3) * "ndvi@1"0.0(被当作数据)250,0000.719(误差 -0.034)
"ndvi@1" / ("ndvi@1" > 0.3)nodata(-3.40282e+38)238,6570.753(准确)

要一步写出真正的 nodata,就用除法代替乘法:"ndvi@1" / ("ndvi@1" > 0.3)。条件为真(1)时值不变;条件为假(0)时,除以零会得到 nodata。

3. 在经纬度 DEM 上先算坡度再用计算器,会把整张图都标成陡坡

在地理坐标 DEM(例如 EPSG:4326 下的 Copernicus GLO-30)上用默认设置运行 QGIS Slope 工具,再用 "slope@1" > 15 取阈值,结果把我们的丘陵农田 100% 标成了大于 15 度,平均坡度为 89.99 度,而重投影到 UTM 后的真实比例是 32.61%。

Copernicus GLO-30 和其他免费全球 DEM(参见我们的免费 DEM 下载指南)都以地理坐标度数(EPSG:4326)提供,高程单位却是米。在默认 Z 因子下,Slope 工具会把以米为单位的高差除以 0.000278 的水平步长(以度为单位的像素宽度),于是每面山坡都像一堵垂直的悬崖。

在 "slope@1" > 15 之前坡度的计算方式大于 15 度的像素平均坡度
在 EPSG:4326 DEM 上运行原生 Slope,默认 Z 因子 1100.00%89.99°
在 EPSG:4326 DEM 上运行 GDAL Slope,scale 1100.00%89.99°
在 EPSG:4326 DEM 上运行 GDAL Slope,scale 11112031.61%12.54°
把 DEM 重投影到 UTM 23S 30 m 后运行原生 Slope(参考值)32.61%12.68°
把 DEM 重投影到 UTM 23S 30 m 后运行 GDAL Slope32.71%12.71°

即使用 111120 米每度这个换算捷径,结果也低了整整一个百分点(31.61% 对 32.61%),因为在南纬 22° 处,一度经度比一度纬度短约 7%(cos 22° = 0.93)。在运行 Slope 和栅格计算器之前,务必先把 DEM 重投影到当地 UTM 分带(Raster > Projections > Warp (Reproject),双线性重采样,并设置 -9999 为 nodata 值,避免边缘出现虚假的 87° 坡度)。

4. 不检查影像就减去 Sentinel-2 的 +1000 偏移,会毁掉 90.6% 的像素

在计算 NDVI 之前从我们的 Sentinel-2 L2A 波段中减去 1000,导致 250,000 个像素中有 226,413 个(90.6%)落在有效范围 [-1, 1] 之外,数值从 -1273 到 1217,均值为 1.39 而不是 0.727。

ESA 从 Sentinel-2 处理基线 04.00 起加入了 +1000 的数字量偏移,以免暗像素在零处被截断,STAC 元数据中也列出了 scale: 0.0001 和 offset: -0.1。但如果你从 Earth Search STAC API 获取 Cloud-Optimized GeoTIFF(我们获取的是影像 S2A_23KLR_20261103_0_L2A,基线 05.13),数据提供方已经把数值平移过,并在条目元数据中设置了 earthsearch:boa_offset_applied: true。

在我们的 B04 文件中,250,000 个像素里有 226,773 个(90.7%)原始数字量低于 1000(最小值 59,中位数 491)。动公式之前,先在 Layer Properties > Information 里查看 B04 的最小值:如果暗像素在 50 到 300 左右,说明偏移已经去除,直接用 (B08 - B04) / (B08 + B04) 就是准确的(除以 10,000 的缩放因子在比值中会抵消)。只有当最暗的晴空像素高于 1000 时,才需要减去 1000。

附加:「Expression valid」和那条通用红色错误栏的真正含义

Raster > Raster Calculator 还有两个行为会浪费大量排错时间:

  • 在粗网格上混合不同分辨率时,不会取平均值。 组合 10 m 波段(B08)和 20 m 波段(B11)时,原生计算器会用最近邻法即时重采样(GDAL 计算器则会直接停止,报错 Error! Dimensions of file ... (250, 250) are different from other files (500, 500). Cannot proceed)。在 10 m 输出网格上,每个 20 m 像素按 2 × 2 重复,结果没有问题。在 20 m 输出网格上,计算器不会对四个 10 m 像素取平均:在 62,500 个单元中有 62,250 个,它只是取了每个 2 × 2 块右上角的那个 10 m 像素,在 (B08 - B11) / (B08 + B11) 中与四像素均值最多相差 0.71。
  • 「Expression valid」只检查语法,不检查波段名称。 如果末尾留了一个运算符,比如 "B08@1" /,对话框能发现(显示 Expression invalid,OK 按钮变灰)。
QGIS Raster Calculator 对话框中有一个以斜杠结尾的不完整公式,左下角显示 Expression invalid,OK 按钮不可用。
语法错误会立即让 OK 按钮变为不可用。

但如果你输入了项目中不存在的图层或波段(例如 "B04@2",或在没有加载 B02 时输入 "B02@1"),提示仍然显示 Expression valid。点击 OK 后,QGIS 会写出一个 500 × 500、100% 填满 nodata 的 GeoTIFF,地图上什么也不添加,只显示一条红色横幅:Raster calculator: An error occurred while performing the calculation. 每次看到这条横幅,就把每个 "layer@band" 名称和 Raster Bands 列表逐一对照。

QGIS 地图画布顶部的红色错误栏,内容为 Raster calculator: An error occurred while performing the calculation.
波段名称拼错能通过语法检查,却以这条通用错误栏失败,并在磁盘上留下一个全是 nodata 的 GeoTIFF。

跳过对话框:在 QGIS 中让 AI 智能体完成栅格计算

AI Agent 是我们开发的,所以这一节讲的是我们自己的插件。它在 QGIS 里停靠一个聊天面板,可以搜索 STAC 目录、重投影 DEM、用 Python 或 Processing 进行栅格计算,并在地图上为输出设置样式。

为了看它在完全相同的 Machado 数据上能否避开上面的陷阱,我在一个只有 OpenStreetMap 的空项目上画了一块 120.07 公顷的农场多边形(UTM 23S 下 E 405300-406300, N 7594400-7595600),然后发送了三条单条消息提示,使用 Medium effort 和 Balanced permissions,每个任务开一个新对话。

测试 1:NDVI 分 5 级并给出每级公顷数

我输入了:

Compute NDVI for this field from the latest cloud-free Sentinel-2 and class it in 5 classes with hectares per class

智能体查询了 STAC 目录,注意到有一景 2026 年 11 月 6 日的影像但云量很大,于是选择了与我手动测试相同的 2026 年 11 月 3 日影像(瓦片云量 1.65%)。它获取了 10 m 的 B08 和 B04,计算 NDVI,分成五个区间,设置图层样式,并在 10 m 公制网格上报告了每一级的面积。

NDVI 分级AI Agent 报告我的独立 numpy 核对(同样 12,000 个像素)差异
-1 到 < 00.76 ha0.76 ha0.0%
0 到 < 0.223.85 ha23.84 ha< 0.1%
0.2 到 < 0.49.27 ha9.27 ha0.0%
0.4 到 < 0.68.70 ha8.69 ha0.1%
0.6 到 177.42 ha77.44 ha< 0.1%
有效总面积120.00 ha120.00 ha(多边形共 120.07 ha)0.0%

对它导出的 NDVI 栅格逐像素比较,与我的参考结果的平均绝对差为 5e-5,它的分级栅格在 12,000 个像素中只有 4 个与我的不同。

QGIS 界面显示 120 公顷农场多边形按五个 NDVI 等级从红到深绿设置样式,Layers 面板显示每级的公顷数,右侧 AI Agent 面板列出 2026 年 11 月 3 日的 Sentinel-2 数据源和各级面积。
一条消息:智能体跳过了多云的 11 月 6 日过境影像,使用 11 月 3 日的影像,各级面积与我们的参考值相差在 0.1% 以内。

测试 2:坡度大于 15 度(项目中没有 DEM)

在一个只加载了农场多边形的新对话中,我输入了:

Show where slope is above 15 degrees on the farm layer

智能体获取了 Copernicus 30 m DEM,先重投影到公制网格,而不是在度数下计算坡度,然后设置了一个两级样式的图层(15 degrees or less 为浅黄色,above 15 degrees 为红色)。它报告农场 1,320 个像素中有 331 个大于 15 度(25.08%,最大 29.15°)。我在同一矩形上用 UTM 重投影后的原生 Slope 得到的参考值是 1,360 个像素中有 341 个(25.07%,最大 29.29°),相差 0.01 个百分点。

QGIS 地图在 OpenStreetMap 上显示农场多边形,坡度大于 15 度的单元为红色,右侧 AI Agent 聊天面板报告根据 30 m Copernicus DEM,农场 1,320 个像素中有 331 个超过阈值,最大坡度 29.1 度。
智能体在取阈值之前把 Copernicus DEM 重投影到了公制网格,与我们的参考值相差在 0.01 个百分点以内。

测试 3:仍然需要你检查智能体的地方

在第三个对话中,我问:

Mask out everything with NDVI below 0.3 on the farm layer and tell me the mean NDVI of what is left

注意,这条提示里我没有写 "cloud-free" 或 "latest"。智能体避开了 0 掩膜陷阱:它给每个被掩膜的像素写入 -9999 nodata,并返回 2,277 个保留像素的 NDVI 均值 0.552。它的输出文件均值为 0.5515;我用同一影像和同一规则算出的是 0.5513。

不过,它用的是 2026 年 10 月 27 日的 Sentinel-2 影像,而那天农场有 80.5%(Scene Classification 波段中 12,000 个像素里的 9,656 个)被云、云影或卷云覆盖。它在回答中写明了日期,却从没提到 0.552 只代表田块中晴空的那五分之一。在晴朗的 11 月 3 日影像上,同样的规则得到的均值是 0.766,覆盖 8,785 个像素。

我的经验法则:无论是手动使用栅格计算器,还是在 QGIS AI 生态 中让智能体来做,在相信一个汇总数字之前,都要先查看影像日期,并加载真彩色或 SCL 云波段。

我们是如何测量的。 所有测试于 2026 年 11 月 7 日在 macOS 上的 QGIS 3.44.15-Solothurn 中进行(GDAL 3.13.3、numpy 2.5.2),区域是巴西米纳斯吉拉斯州 Machado 以南一个 5 km × 5 km 的窗口(EPSG:32723,E 402000-407000,N 7591000-7596000),使用来自 Earth Search 的 Sentinel-2 L2A 影像 S2A_23KLR_20261103_0_L2A(测试 3 另用 2026 年 10 月 27 日的影像)和 Copernicus GLO-30 瓦片 S22_00_W046_00。手动测试以及智能体测试 1 和 2 的成功标准都是在任何栅格计算之前写好的;测试 3 的标准是在看过测试 1 和 2 之后写的。参考网格在 QGIS 之外用 64 位浮点 numpy 计算。AI Agent 1.9.4 使用 Medium effort 和 Balanced permissions,每个对话一条消息。一景影像、一个地点只是单一案例的基准测试,不是跨处理基线的统计调查。

常见问题

为什么我在 QGIS 中算出的 NDVI 全是 0、-1 或 1?

你在 GDAL 栅格计算器中强制使用了整数输出类型(Int16 或 Byte),它会把 -1 到 1 之间的每个小数四舍五入成 -1、0 或 1(在我们的测试中占 99.99% 的像素)。把输出类型保持为 Float32,或者使用顶部菜单的 Raster > Raster Calculator,它会自动输出 Float32。

为什么 QGIS 栅格计算器报错 'An error occurred while performing the calculation' 并输出 nodata?

公式中的某个 "layer@band" 名称与已加载的图层或波段编号不匹配(例如在单波段文件上输入 "B04@2",或引用了已经移除的图层)。公式框下方的 Expression valid 只检查语法,不检查波段是否存在。在 Raster Bands 列表中双击每个波段,重新插入即可。

Raster > Raster Calculator 和 Processing 中的 GDAL Raster Calculator 有什么区别?

Raster > Raster Calculator(以及 Processing 中原生的 Raster calculator)使用 QGIS 自己的解析器和 "layer@band" 名称,始终以浮点数计算。GDAL Raster calculator(gdal_calc.py)使用字母变量(A、B),并按输入波段的数据类型计算 numpy 数学表达式,所以无符号 16 位整数波段相减时,负数会绕回 65,536,除非用 A.astype(numpy.float32) 做类型转换。

如何在 QGIS 栅格计算器中写 IF 条件或把数值设为 nodata?

可以写 if("ndvi@1" > 0.3, "ndvi@1", 0),把低于 0.3 的值替换为 0。如果希望被排除的像素变成真正的 nodata,不拉低图层均值,就改为除以条件:"ndvi@1" / ("ndvi@1" > 0.3)。条件为假(0)时,除以零会写入 nodata。

免费试用 AI Agent

在 QGIS 里,说一句就能完成

AI Agent 会在你自己的 QGIS 项目中,完成你描述的任务。免费使用,无需绑卡。

免费试用 AI Agent
由以下作者撰写Yvann Barbot
·9 min read