企业官网建设流程全解析

做遥感数据处理，前两篇我们把 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，而不是常见的 BGR＋NIR 布局，直接用固定索引取波段，算出来的 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。如果用 SpatVector，zonal()会先把矢量栅格化，再按区域统计，一步到位，省不少事。

# 统计每个地块 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("COMPRESS=DEFLATE", "TILED=YES", "BIGTIFF=YES"), overwrite = TRUE)

COMPRESS=DEFLATE能显著减小文件体积，TILED=YES让后续读取速度加快不少，BIGTIFF=YES是为了避免超过 4 GB 时文件格式无法扩展的问题。这几个参数在项目里几乎是标配。

5.3 用 SpatRaster 子集代替整幅读入

有些操作只需要影像的局部范围，比如按图幅处理时，直接整幅读取全波段数据很浪费。合理使用rast()配合范围参数，或者用crop()先取得一个较小范围再处理，能节省大量 I/O 时间。

另外，如果只是做预处理，不需要产生各种中间结果，我通常会避免反复writeRaster()中间栅格。每写一个文件就是一次完整遍历，项目数据大时，写 3 个中间文件就意味着把全部数据读写了 3 遍，耗时成倍增长。改为保持“链式”执行，只在最终结果才落盘，是一种非常直观但被很多人忽略的性能优化。

6. 组合成一条完整数据流水线：两期影像的植被变化分析

6.1 项目背景与数据准备清单

讲完上面这些，我把它们串成一个可复跑的完整流程。场景是某半干旱流域的两期多光谱影像（2020 年和 2023 年），外加一版地块矢量边界。最终目标是输出每个地块内植被覆盖变化的统计表，以及一张变化专题图。

数据清单如下：

• 影像1：scene2020.tif，多光谱，含蓝、绿、红、近红外波段
• 影像2：scene2023.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 结果检查：统计表不等于交付结果

拿到统计表之后，我从来不敢直接交付。先做几个基础检查：

1. freq(diff_class)检查分类栅格的类别像元数是否合理，有没有出现某个类别占比异常到 99% 的情况。
2. 对照global(ndvi_diff, "mean", na.rm=TRUE)检查全局平均变化的方向是否和分区统计的方向一致。
3. 抽查两三个地块，把分区统计值和原始 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 包做模型预测的内容。