如何在R中提取Landsat影像值并绘制不同地物光谱曲线

嘿，针对你的需求，这里有一套清晰且高效的实现步骤，都是R中处理遥感数据的常用方案：

1. 先确保数据坐标系一致（关键前提） #

首先得确认你的Landsat栅格堆栈和研究区shapefile的投影是匹配的，不然提取会出错。可以用以下代码检查和转换：

# 检查投影信息
crs(landsat)
crs(roi)

# 如果投影不一致，将矢量数据转换为栅格的投影
if (!identical(crs(landsat), crs(roi))) {
  roi <- spTransform(roi, crs(landsat))
}

2. 用extract()函数提取光谱值 #

raster包的extract()是提取栅格值到矢量的标准工具，根据你的需求有两种常见用法：

2.1 提取每个地物多边形的平均光谱值（适合快速获取典型光谱曲线） #

如果你的roi属性表中有区分地物的字段（比如叫land_class，值为"水体"、"道路"、"植被"），可以直接按多边形提取均值：

library(raster)

# 提取每个多边形的各波段均值，忽略NA值（如云、阴影）
spectral_values <- extract(landsat, roi, fun = mean, na.rm = TRUE, df = TRUE)

# 将提取结果和原矢量的属性表合并
spectral_df <- cbind(roi@data, spectral_values)

# 移除extract默认生成的冗余ID列
spectral_df <- spectral_df[, !names(spectral_df) %in% "ID"]

2.2 提取每个像元的光谱值（适合后续精细统计或分类验证） #

如果需要研究区内所有像元的光谱数据，可以去掉fun参数，保留每个像元的原始值：

# 提取所有像元的光谱值，返回结构化数据框
pixel_spectra <- extract(landsat, roi, df = TRUE)

# 合并属性表，让每个像元关联对应的地物类型
pixel_spectra <- merge(pixel_spectra, roi@data, by.x = "ID", by.y = row.names(roi@data))

3. 整理数据并绘制光谱曲线 #

提取后的数据需要转成长格式才能方便绘图，用tidyverse工具包会很顺手：

library(tidyverse)

# 转换为长格式，便于分组绘图
long_spectra <- spectral_df %>%
  pivot_longer(cols = starts_with("B"), # 匹配所有波段列
               names_to = "Band",
               values_to = "Reflectance")

# 绘制不同地物的光谱曲线
ggplot(long_spectra, aes(x = Band, y = Reflectance, color = land_class, group = land_class)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  labs(title = "研究区地物光谱曲线", x = "Landsat波段", y = "反射率") +
  theme_minimal()

4. 额外优化建议 #

• 如果你的Landsat数据未做大气校正，建议先完成辐射定标和大气校正步骤，这样得到的光谱值才更具实际意义。
• 处理超大规模栅格时，可以用extract()的blockSize参数分块提取，避免内存溢出：

spectral_values <- extract(landsat, roi, fun = mean, na.rm = TRUE, df = TRUE, blockSize = blockSize(landsat))

内容的提问来源于stack exchange，提问作者Teri P.