R语言terra包遥感数据处理进阶：栅格代数与分区统计实战-北京尧图网络

做遥感数据处理前两篇我们把 terra 的基本读写、投影变换、裁剪和重采样都过了一遍。说实话很多人玩到这一步就卡住了数据能读、能裁、能看可一旦需要从一堆栅格里真正“算出点结果”——比如提取植被指数、做变化检测、按行政区统计平均值——就发现光会基础操作远远不够更麻烦的是代码写得不规范几百兆的栅格一跑就内存爆掉或者算出来的结果根本对不上。这一篇我打算换个思路不再按函数功能逐个罗列而是顺着一条完整的实际分析流程往下走假设我们要处理某半干旱流域两期多光谱影像目标是提取植被覆盖变化并把结果按地块汇总成统计表。整个过程会涉及到栅格代数运算、条件掩膜、分区统计、多源数据对齐以及大区域数据处理时的性能问题。这些都是 terra 处理遥感数据时真正绕不开的硬骨头也是我认为系列第三篇最该讲透的内容。1. 栅格代数与波段计算从像元推演到整套指标1.1 app 和 lapp 的分工很多人一开始就搞混了terra 里做栅格代数最常用的两个函数是app()和lapp()。名字长得像用途却完全不同我见过不少人在这一步栽跟头。app()是针对单个 SpatRaster不管它是单波段还是多波段按像元套用同一个函数。什么意思呢比如你有一个多波段影像要对每个像元的所有波段做标准化或者要计算整个栅格的最大值、最小值用app()就对了。它从头到尾只处理“一个栅格对象”函数内部拿到的是该像元位置上的所有波段值构成的向量。lapp()则是针对多个 SpatRaster 对象做逐像元的函数运算。比如你有两期 NDVI 图层想计算它们的差值就把这两个 SpatRaster 放进一个 list用lapp()把函数套用到两套栅格的对应像元上。实际操作中我通常这样把握如果运算逻辑需要的输入都来自同一个栅格的不同波段优先app()如果输入涉及多个独立栅格比如几个不同时间的影像就上lapp()。这个区分不搞清楚后面所有代码都会别扭。多波段 NDVI 的计算就是一个典型例子。terra 读取多光谱影像后是一个多波段的 SpatRaster。NDVI 需要红波段和近红外波段本质上是“同一个影像内部的波段代数”。我把影像读进来按波段名或者波段索引取出子集然后用lapp()把两个波段的 SpatRaster 对象组合运算。library(terra) # 读取多光谱影像波段顺序假设为 B、G、R、NIR img - rast(scene2020.tif) names(img) # 提取红波段和近红外波段 red - img[[3]] nir - img[[4]] # 计算 NDVI ndvi - lapp(list(red, nir), function(r, n) { (n - r) / (n r) }) names(ndvi) - NDVI_2020等一下有细心的朋友会问lapp()为什么传的是 list 而不是两个独立参数这是 terra 的固定要求函数第一个参数必须是 SpatRaster 组成的 list后面的参数才是自定义函数的辅助形参。我在早期经常习惯性把两个栅格放成lapp(red, nir, fun)然后报错一脸懵。1.2 多维数组思维像元位置才是栅格运算的锚点用lapp()写自定义函数时有一个容易忽略的点函数接收的并不是单个数值而是每个像元在所有输入栅格中的值组成的向量。还是拿 NDVI 举例function(r, n)里的r不是红波段整个栅格而是“当前像元”的红波段值对于整个影像terra 会把函数在每个像元上分别执行一次。这个机制保证了所有输入栅格必须拥有相同的行列数、范围和分辨率不然根本无法对齐像元。所以在多栅格运算之前我几乎必做的一步是检查所有输入是否对齐compareGeom(red, nir)如果返回TRUE说明行列数、范围、分辨率一致可以放心运算。如果返回FALSE不可强行塞进lapp()必须先做重采样或对齐后面第4节会细讲。类比对 Excel 的操作lapp()就像你把几张同样大小的表格叠在一起把每个同位置的单元格单独抽出来做计算前提是这几张表的大小、行列顺序完全一致。1.3 把常用指标封装成函数别每次都复制粘贴做遥感数据处理NDVI、EVI、NDWI、RVI 这类指数几乎是家常便饭。我强烈建议把指标计算封装成独立函数放在单独脚本里统一管理。一方面代码可读性大幅提升另一方面同一个项目里不同时期的数据可以反复调用不用每次把公式重写一遍也减少写错波段顺序的概率。calc_ndvi - function(red, nir) { (nir - red) / (nir red) } calc_ndwi - function(green, nir) { (green - nir) / (green nir) } # 使用 ndvi_2020 - lapp(list(img_2020[[3]], img_2020[[4]]), calc_ndvi) ndvi_2023 - lapp(list(img_2023[[3]], img_2023[[4]]), calc_ndvi)这里我刻意把波段名写在[[ ]]里而不是波段号以降低波段顺序记错的概率。我踩过一次真实的坑某个数据源的波段顺序是 BGRN而不是常见的 BGRNIR 布局直接用固定索引取波段算出来的 NDVI 整体偏小成图后色调完全不对。后来我拿到数据的第一件事永远是打印names(img)哪怕多花十秒钟也比后面返工两个小时划算。2. 条件运算与重分类把连续值变成决策结果2.1 ifel 和 classify什么时候用哪个栅格运算算出来的 NDVI 是连续浮点值可实际业务里我们往往需要的是“有没有植被”“水体还是陆地”“高覆盖还是低覆盖”这类离散结论。这个时候就轮到条件运算和重分类出场了。terra 里做条件运算的核心函数是ifel()意思是“if-else”的栅格版本。比如要把 NDVI 中大于 0.3 的像元保留其他置为 NA可以这样写veg_mask - ifel(ndvi_2020 0.3, ndvi_2020, NA)ifel()的三个参数分别是条件、条件为 TRUE 时的输出值、条件为 FALSE 时的输出值。这个函数优势在于逻辑清晰适合做简单的阈值筛选。但一旦阈值多起来比如要把 NDVI 分成 5 个等级逐层嵌套ifel()会写得又臭又长我试过写三层嵌套回头自己都看晕了。多层级的阈值划分正确做法是用classify()。它是栅格版的cut()或reclassify()把连续值映射到离散类别。# 定义分类矩阵: 左闭右开from, to, new_value legends - matrix(c( -1.0, 0.0, 1, 0.0, 0.2, 2, 0.2, 0.4, 3, 0.4, 0.6, 4, 0.6, 1.0, 5 ), ncol 3, byrow TRUE) ndvi_class - classify(ndvi_2020, legends)2.2 阈值选取不要想当然先看数据分布分类阈值怎么定不少教材会直接给一个经验值比如 NDVI 0 算水体0~0.2 算裸地0.2~0.4 算低覆盖植被。但实际项目里不同传感器、不同季节、不同大气校正条件NDVI 分布特征差异非常大。我通常会先看一下全局统计值再决定分级区间global(ndvi_2020, c(min, max, mean, sd), na.rm TRUE)如果最小值只有 0.05最大值也只有 0.6说明这个影像已经经过了某种截断处理固定的经验阈值意义不大。更稳妥的做法是用分位数辅助设计分类区间q - global(ndvi_2020, fun quantile, probs c(0.05, 0.25, 0.5, 0.75, 0.95), na.rm TRUE)基于实际分位数分段能让分类结果更贴合研究区的地表真实状况而不是被一个理论上“通用”的阈值带偏。这里我不建议直接把分位数当成阈值用那样相当于把数据强行均分失去了分类的物理意义但参考它来调整区间端点是合理的。2.3 掩膜把不需要的区域干净利落地清掉条件运算经常和掩膜配合使用。所谓掩膜就是用一个二进制栅格或矢量范围把研究区之外、或者不感兴趣的像元遮掉。terra 的mask()函数默认保留mask参数为 NA 之外的像元# 假设有一个研究区边界 shp aoi - vect(study_area.shp) # 先按 NDVI 阈值筛选植被像元 veg - ifel(ndvi_2020 0.3, 1, NA) # 再用研究区边界掩膜 veg_masked - mask(veg, aoi)注意一个细节mask()的第二个参数可以是 SpatVector也可以是 SpatRaster。如果传 SpatVector它内部会自动转成栅格掩膜层效率稍低但使用方便如果传 SpatRaster则要求掩膜层与目标栅格几何完全对齐。大数据量时我更倾向于直接把研究区矢量栅格化成掩膜层这样后续多次mask()不用重复转换。关于 NA 的语义我一直给项目里的新人强调在栅格运算里NA 不是 0而是“不参与计算”。很多条件判断遇到 NA 会直接返回 NA如果后面还要做统计NA 像元既不计入分子也不计入分母在比例计算时容易造成偏差。为避免这种问题需要明确检查# 查看每个类别像元数排除 NA freq(ndvi_class)freq()会直接展示每个类别的像元数不统计 NA正好用来看分类结果是否合理。3. 分区统计与区域汇总让栅格数据真正“落地”3.1 为什么必须做分区统计栅格结果只是“一张图”业务上很难直接使用真正能被汇报、制表、用于后续决策的往往是“某个区域的平均值、总面积、占比”。把栅格的像元信息按区域地块、行政区、流域子区聚合成统计表这一步叫分区统计也就是zonal()函数干的事。我做一个植被变化分析项目时客户要求输出的不是“NDVI 变化栅格”而是“每个地块里植被覆盖面积增加了多少平方公里退化了多少平方公里”。如果没有分区统计我就只能把变化栅格整张存档毫无交付价值。zonal()的基本逻辑是至少传入两个栅格——被统计的目标栅格和定义分区的分区栅格。分区栅格既可是 SpatRaster也可是 SpatVector。如果用 SpatVectorzonal()会先把矢量栅格化再按区域统计一步到位省不少事。# 统计每个地块 2020 年的平均 NDVI z_ndvi - zonal(ndvi_2020, aoi, fun mean, na.rm TRUE)返回的是 data.frame每一行对应aoi里的一个地块列是fun指定的统计量。最常用的统计函数有mean、sum、min、max、sd、count等需要多个统计量可以同时传多个z_stats - zonal(ndvi_2020, aoi, fun c(min, max, mean, sd), na.rm TRUE)3.2 用面积而不是像元数来计算覆盖比例如果我需要知道“每个地块中高覆盖植被的面积占比”直接数像元个数再乘以像元面积是最朴素的做法但不够稳。尤其是当影像经过重投影或者分辨率不均匀时单个像元的面积可能会随纬度变化这种情况在用经纬度坐标系处理大范围数据时非常常见。terra 提供了cellSize()函数可以逐像元计算真实地表面积结合条件栅格求和能获得更准确的面积统计# 逐像元面积平方米 cell_area - cellSize(ndvi_class, unit m) # 高覆盖植被类别4和5面积单位转换为平方公里 high_veg - ifel(ndvi_class %in% c(4, 5), 1, NA) area_high - zonal(high_veg * cell_area, aoi, fun sum, na.rm TRUE) area_high$area_km2 - area_high$sum / 1e6这里有个小知识点cellSize()返回的是与输入栅格同样规格的栅格每个像元的值为该像元的实际面积。在高纬度地区用普通像元数乘标称分辨率计算的面积会明显偏大而用cellSize()能避免这个问题。3.3 需要注意的统计口径na.rm 和权重分区统计里最容易被忽视的就是na.rm参数。默认情况下不少函数会返回 NA而遥感影像里 NA 几乎是常态云掩膜、范围外、无效值。如果na.rm FALSE只要区域内有一个 NA整个区域统计结果就是 NA这在处理大影像时几乎必然导致结果全是 NA。另外一个容易忽视的问题是面积加权平均。假设你有一个区域的 NDVI 均值区域内像元面积不一致同样是高纬/低纬混排直接用mean得到的其实是“像元平均”而不是“面积平均”。要严谨可以用我上面提到的cellSize()做加权# 面积加权 NDVI ndvi_weighted - ndvi_2020 * cell_area sum_w - zonal(ndvi_weighted, aoi, fun sum, na.rm TRUE) area_by_zone - zonal(cell_area, aoi, fun sum, na.rm TRUE) sum_w$area_weighted_ndvi - sum_w$sum / area_by_zone$sum这种方式在项目汇报时更有说服力尤其是当研究区跨越多个纬度带时。4. 多源数据对齐统整分辨率、投影与边界一次理清4.1 拿到数据的第一件事投影和范围体检做遥感的都知道不同来源的数据很少能直接叠加运算。某系列 30 米分辨率数据、另一类 10 米分辨率数据再加上从第三方平台下载的矢量边界投影五花八门——有的用经纬度有的用投影坐标还有的单位都不是米。我在一个项目里就吃过亏。当时做多源数据集成分析看起来范围差不多就跳过投影检查直接lapp()结果系统报错说“different extents”排查了半天才发现两个影像虽然范围相近但投影坐标系一个用的是以米为单位的多圆锥投影另一个是经纬度行列数差了将近三倍。所以现在我拿到任何栅格数据第一件事永远是做“体检”img_2020 img_2023 crs(aoi) ext(aoi)4.2 project 与 resample别搞混“重投影”和“重采样”terra 里有两个容易混淆的函数project()用于栅格/矢量的投影变换resample()用于把栅格重采样到另一个网格。当目标只是统一到同一投影坐标系时用project()。它的一个优势在于project()默认使用较精确的重采样方法并且你可以在project()时直接指定目标分辨率。# 把影像从当前投影转到目标投影同时设定30米分辨率 img_2023_reproj - project(img_2023, EPSG:32650, res 30)而当两个栅格已经投影一致但网格位置、分辨率不同需要“对齐”时用resample()把其中一个重采到另一个的网格上img_2023_aligned - resample(img_2023, img_2020, method bilinear)method的选择需要想清楚。连续变量NDVI、温度、高程用bilinear效果好平滑合理类别变量土地利用类型、分类结果用near最近邻千万不能用 bilinear否则会插值出不存在的类别。4.3 范围对齐与裁剪先裁剪还是先对齐多源数据处理时我习惯的顺序是先统一投影再统一分辨率最后统一范围。如果数据范围差很多直接resample()会把整个大范围都纳入计算徒增计算量。合理的做法是先用intersect()求交集范围或crop()到研究区再做重采样或投影。# 裁剪到研究对象范围 img_2020_crop - crop(img_2020, aoi) img_2023_crop - crop(img_2023, aoi) # 然后对齐到统一网格 img_2023_final - resample(img_2023_crop, img_2020_crop, method bilinear)很多新人容易犯的错误是顺序颠倒先重采样后裁剪白白计算了研究区之外的数千万个像元数据量大时这一步就可能吃满内存。另外crop()之后建议用trim(img_2020_crop)把外围的 NA 边界裁掉避免后续运算多走无谓的像元。5. 实测中的性能与内存管理大范围栅格运算不卡死的经验5.1 terra 的懒执行和分块处理机制用过raster包的读者一定体会过处理大影像时的绝望一个 10 GB 的影像做局部运算R 直接卡死几小时。换成 terra 之后情况改善明显但前提是你要理解它的工作机制。terra 默认采取“懒执行”加“分块处理”的策略。很多函数并不会立刻把全部数据读入内存而是等到你真正访问结果时才按块chunk逐个计算。这意味着你可以把复杂的多层运算链写在一条代码里而不必担心中间结果瞬间占满内存。比如out - ifel(ndvi_2020 0.3, 1, NA) out - mask(out, aoi) out - project(out, EPSG:32650)如果写成链式terra 会在最后plot()或writeRaster()时才真正执行计算并自动按分块策略逐块处理。这比一步步生成中间栅格文件高效得多。5.2 写文件时要显式控制压缩和分块writeRaster()是执行计算的触发点也是性能瓶颈所在。大数据写入 GeoTIFF 时我强烈建议设置压缩选项和内部分块writeRaster(out, NDVI_change.tif, filetype GTiff, datatype FLT4S, gdal c(COMPRESSDEFLATE, TILEDYES, BIGTIFFYES), overwrite TRUE)COMPRESSDEFLATE能显著减小文件体积TILEDYES让后续读取速度加快不少BIGTIFFYES是为了避免超过 4 GB 时文件格式无法扩展的问题。这几个参数在项目里几乎是标配。5.3 用 SpatRaster 子集代替整幅读入有些操作只需要影像的局部范围比如按图幅处理时直接整幅读取全波段数据很浪费。合理使用rast()配合范围参数或者用crop()先取得一个较小范围再处理能节省大量 I/O 时间。另外如果只是做预处理不需要产生各种中间结果我通常会避免反复writeRaster()中间栅格。每写一个文件就是一次完整遍历项目数据大时写 3 个中间文件就意味着把全部数据读写了 3 遍耗时成倍增长。改为保持“链式”执行只在最终结果才落盘是一种非常直观但被很多人忽略的性能优化。6. 组合成一条完整数据流水线两期影像的植被变化分析6.1 项目背景与数据准备清单讲完上面这些我把它们串成一个可复跑的完整流程。场景是某半干旱流域的两期多光谱影像2020 年和 2023 年外加一版地块矢量边界。最终目标是输出每个地块内植被覆盖变化的统计表以及一张变化专题图。数据清单如下影像1scene2020.tif多光谱含蓝、绿、红、近红外波段影像2scene2023.tif多光谱波段顺序一致矢量边界plots.shp研究区内的地块边界整个流水线无外乎投影体检 → 计算 NDVI → 阈值分类 → 掩膜裁剪 → 定量统计 → 输出结果。6.2 核心代码链路library(terra) # 1. 读取数据与体检 img_2020 - rast(scene2020.tif) img_2023 - rast(scene2023.tif) plots - vect(plots.shp) # 统一投影假设目标投影为 UTM 50N img_2020 - project(img_2020, EPSG:32650, res 30) img_2023 - project(img_2023, EPSG:32650, res 30) plots - project(plots, EPSG:32650) # 裁剪 img_2020 - crop(img_2020, plots) img_2023 - crop(img_2023, plots) # 若两期影像网格有差异统一网格 if (!compareGeom(img_2020, img_2023)) { img_2023 - resample(img_2023, img_2020, method bilinear) } # 2. 计算 NDVI red_2020 - img_2020[[3]]; nir_2020 - img_2020[[4]] red_2023 - img_2023[[3]]; nir_2023 - img_2023[[4]] ndvi_2020 - lapp(list(red_2020, nir_2020), calc_ndvi) ndvi_2023 - lapp(list(red_2023, nir_2023), calc_ndvi) # 3. 变化检测 ndvi_diff - ndvi_2023 - ndvi_2020 names(ndvi_diff) - NDVI_change # 4. 分类改善0.1、稳定-0.1~0.1、退化-0.1 diff_class - classify(ndvi_diff, rbind( c(-1, -0.1, 1), c(-0.1, 0.1, 2), c(0.1, 1, 3) )) # 5. 分区统计各地块变化类别面积 cell_area - cellSize(diff_class, unit m) area_per_class - zonal(diff_class * cell_area, plots, fun sum, na.rm TRUE)这一步出来的area_per_class是一个数据框第一列是plots自带的属性 ID后面是每个类别对应像元面积总和。除以 1e6 就是平方公里再配合plots本身的属性表就能拼出地块尺度的变化统计表了。6.3 结果检查统计表不等于交付结果拿到统计表之后我从来不敢直接交付。先做几个基础检查freq(diff_class)检查分类栅格的类别像元数是否合理有没有出现某个类别占比异常到 99% 的情况。对照global(ndvi_diff, mean, na.rmTRUE)检查全局平均变化的方向是否和分区统计的方向一致。抽查两三个地块把分区统计值和原始 NDVI 栅格目视对比确认没有因为像元对齐问题导致区域错位。另外我通常会把变化栅格写成 GeoTIFF 存档并顺手生成一张简单的专题图cols - c(#d73027, #ffffbf, #1a9850) plot(diff_class, col cols, main NDVI Change 2020-2023)输出 PNG 或 PDF 后用于汇报材料。但要记住制图只是“可视化”背后的统计精度和业务合理性比图好看重要得多。7. 最后再分享一个实用小习惯项目做了多了之后我养成了一个固定习惯把所有中间产物按“处理阶段”命名并且绝不覆盖原始数据。比如scene2020_utm.tif、ndvi2020.tif、ndvi_diff.tif、change_class.tif每一步都留在磁盘上。这样如果某一个环节出了问题不用从头跑起直接从对应阶段的文件继续即可。尤其在处理几百 GB 级别的项目时这个习惯能省下大量的重算时间。另外一个建议是多用sources()和describe()查看栅格对象的来源与元数据。terra 里sources()能告诉你当前 SpatRaster 是哪个文件、哪个波段来的describe()则能打印栅格文件的信息摘要。多看一眼这两条输出能在数据漂移、波段错位等疑难杂症出现时帮你快速定位问题到底出在源头文件还是处理逻辑。三篇写下来从基础读写到进阶运算再到这篇的多源数据统整与统计汇总这套流程基本覆盖了我日常 80% 的遥感数据工作。剩下的异常检测、时序分析、模型预测等方向就要看具体项目需求再继续推进了。如果后面还有机会我会再写一写 terra 搭配其他 R 包做模型预测的内容。

R语言terra包遥感数据处理进阶：栅格代数与分区统计实战

相关推荐

基于Boost.Signals2实现类Qt信号槽：C++模块解耦与线程安全实践

AHP层次分析法Python实现：从原理到技术选型实战

让 Copilot 真正懂你的代码库：4 种上下文投喂方法实战指南

更多资讯

CUA智能体完全指南：从原理到手搓一个极简原型

层次聚类实战：组平均、BIRCH和CURE选型指南

时间序列预测入门：ARIMA模型从原理到实战

数据仓库与数据挖掘考试核心考点：从星型建模到Apriori手算

GeoServer 2.25.3 WAR包部署避坑指南：Tomcat 10+适配与纯净化实战

道路缓和曲线VC计算：回旋线高精度坐标生成与工程落地

今日更新

Cephadm实战指南：Linux系统工程师的生产级Ceph集群部署与排障

不想再被 Token 额度卡脖子？把免费 API 接入 WorkBuddy 的完整教程

ITK中SpatialObject几何体系解析：坐标系、变换与管道集成

本周精选

交通拥堵预测毕设：基于Python+LSTM的车辆流量预测系统全流程解析

MATLAB/Simulink搭建10机39节点电力系统暂态稳定仿真实战

无人机航拍三维重建全流程：从SfM到网格生成的避坑指南

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？