如何基于bfastlite识别的时间序列断点绘制空间分布图

嗨，我完全懂你现在的需求——用bfastlite找出时间序列的断点后，想把这些断点信息和研究区的像素/采样点绑定，画出像你描述的那种带标记的空间分布图对吧？我给你梳理一套可行的方案，结合代码示例来讲解~

核心思路 #

bfastlite本身是针对单条时间序列做断点检测的，要实现空间可视化，核心是两步：

1. 批量处理研究区每个像素/采样点的时间序列，提取断点信息（比如是否有断点、断点数量、断点年份等）
2. 把提取到的断点信息和对应的空间对象（栅格/矢量点）结合，最后做空间绘图

具体实现步骤（以栅格时间序列为例） #

假设你手里的是遥感时间序列栅格（比如NDVI月度/年度数据），我们用模拟数据来演示整个流程：

1. 加载所需包 #

library(bfast)
library(terra)  # 处理栅格数据
library(tidyverse)  # 数据整理
library(sf)  # 空间矢量处理

2. 模拟栅格时间序列 #

我们模拟一个3×3的小栅格，包含10年的时间序列，并给其中几个像素手动加入断点，模拟真实场景：

set.seed(123)  # 设置随机种子保证结果可复现

# 生成10个时间步长的栅格
r_list <- lapply(1:10, function(i) {
  rast(nrow=3, ncol=3, vals=rnorm(9, mean=0, sd=0.5))
})

# 给第6个时间步的部分像素加入突变（模拟断点）
r_list[[6]] <- r_list[[6]] + c(0, 2, 0, 2, 0, 0, 0, 0, 2)

# 合并成栅格栈，并设置时间属性
ts_rast <- rast(r_list)
time(ts_rast) <- seq.Date(as.Date("2010-01-01"), by="year", length.out=10)

3. 批量检测每个像素的断点 #

我们把栅格转换成每个像素的时间序列，再循环用bfastlite检测断点：

# 提取每个像素的时间序列并整理成嵌套数据框
pixel_ts <- as.data.frame(ts_rast) %>% 
  pivot_longer(cols=everything(), names_to="date", values_to="value") %>% 
  mutate(date = as.Date(date)) %>% 
  group_by(pixel = row_number() %% ncell(ts_rast)) %>% 
  arrange(date) %>% 
  nest()

# 定义断点检测函数：输入单条时间序列，返回断点数量
detect_breaks <- function(df) {
  # 把数据转换成ts对象
  ts_obj <- ts(df$value, start=c(year(min(df$date)), 1), frequency=1)
  # 用bfastlite检测断点
  bp <- bfastlite(ts_obj)
  # 返回断点数量（如果是NA表示时间序列无效）
  return(list(break_count = ifelse(all(is.na(df$value)), NA, length(bp$breakpoints))))
}

# 给每个像素应用断点检测函数
pixel_breaks <- pixel_ts %>% 
  mutate(break_info = map(data, detect_breaks)) %>% 
  unnest_wider(break_info) %>% 
  mutate(has_break = ifelse(break_count > 0, 1, 0))  # 标记是否有断点

4. 绑定空间信息并可视化 #

把断点结果和栅格的空间位置绑定，然后用ggplot2绘制空间分布图，用红色叉号标记有断点的像素：

# 把栅格转换成sf矢量点对象，获取每个像素的空间坐标
pixel_sf <- as.points(ts_rast[[1]]) %>% 
  st_as_sf() %>% 
  mutate(pixel = row_number())

# 合并断点信息到空间对象
pixel_sf_breaks <- pixel_sf %>% 
  left_join(pixel_breaks %>% select(pixel, has_break, break_count), by="pixel")

# 绘制空间分布图
ggplot() +
  # 先画所有像素的基础点，用颜色区分是否有断点
  geom_sf(data=pixel_sf_breaks, aes(color=factor(has_break)), size=3) +
  # 给有断点的像素叠加红色叉号标记
  geom_sf(data=filter(pixel_sf_breaks, has_break == 1), shape=4, size=4, color="red") +
  # 设置颜色标签和标题
  scale_color_manual(values=c("0"="gray", "1"="blue"), labels=c("无断点", "有断点")) +
  labs(title="时间序列断点空间分布", color="断点状态") +
  theme_minimal()

优化与拓展 #

• 大栅格提速：如果你的研究区很大，循环处理会很慢，可以用terra::app函数做批量并行处理，比如：

# 用terra::app直接对栅格栈做逐像素处理
break_count_rast <- app(ts_rast, function(x) {
  if(all(is.na(x))) return(NA)
  ts_obj <- ts(x, start=c(2010,1), frequency=1)
  bp <- bfastlite(ts_obj)
  return(length(bp$breakpoints))
})

# 直接绘制栅格图
plot(break_count_rast, main="每个像素的断点数量", col=viridis::viridis(5))

• 更细致的可视化：除了标记是否有断点，你还可以提取断点的年份信息，用不同颜色表示不同年份发生的断点，让图的信息更丰富。
• 矢量采样点场景：如果你的数据是采样点（不是栅格），逻辑完全一致——每个点对应一条时间序列，检测断点后把结果和点的空间信息合并，再绘图即可。

备注：内容来源于stack exchange，提问作者Nikos