如何从土地利用比例数据框创建年度土地利用变化桑基图

原代码核心问题 #

你的代码逻辑错误：将年份作为源节点、土地类型作为目标节点，展示的是各年份的土地类型占比，而非年度间土地利用的流转关系。流转桑基图需要的是前一年的土地类型节点连接后一年的土地类型节点，以此体现土地类型的跨年度转移。

修正后的完整代码 #

以下代码基于你提供的年度比例数据，模拟土地类型的流转关系（减少类型的流失量按比例分配给增加类型），生成符合预期的年度流转桑基图：

library(tidyverse)
library(networkD3)

# 原始数据
lu_dt <- structure(list(LU_class = c("Cropland", "Forest", "Grassland", 
"Other Land", "Settlement", "Water", "Wetland"), lu2000p = c(27.79, 
22.92, 0.78, 0.05, 47.66, 0.34, 0.46), lu2005p = c(27.86, 22.51, 
0.78, 0.05, 48, 0.34, 0.46), lu2010p = c(23.29, 17.37, 0.69, 
0.03, 57.86, 0.34, 0.42), lu2015p = c(21.36, 16.95, 0.66, 0.03, 
60.24, 0.34, 0.42), lu2020p = c(21.07, 16.81, 0.65, 0.03, 60.68, 
0.34, 0.41)), row.names = c(NA, -7L), class = c("tbl_df", "tbl", 
"data.frame"))

# 转换为长格式并提取年份
lu_long <- lu_dt %>%
  pivot_longer(cols = starts_with("lu"), names_to = "year_str", values_to = "proportion") %>%
  mutate(year = as.integer(str_extract(year_str, "\\d{4}"))) %>%
  select(LU_class, year, proportion)

# 获取排序后的年份列表
years <- unique(lu_long$year) %>% sort()

# 初始化链路数据框
links <- tibble()

# 遍历每对相邻年份，构造流转链路
for (i in 1:(length(years)-1)) {
  year_current <- years[i]
  year_next <- years[i+1]
  
  # 获取当前年份和下一年份的土地比例数据
  current_data <- lu_long %>% filter(year == year_current)
  next_data <- lu_long %>% filter(year == year_next)
  
  # 合并数据并计算类型变化量
  merged <- current_data %>%
    left_join(next_data, by = "LU_class", suffix = c("_current", "_next")) %>%
    mutate(change = proportion_next - proportion_current)
  
  # 构造同类型"保持"链路（取两年中较小的比例值作为存量流量）
  keep_links <- merged %>%
    mutate(
      keep_flow = ifelse(change >= 0, proportion_current, proportion_next),
      source = paste(year_current, LU_class, sep = "_"),
      target = paste(year_next, LU_class, sep = "_")
    ) %>%
    select(source, target, value = keep_flow)
  
  # 构造跨类型"转移"链路（将减少类型的流失量按比例分配给增加类型）
  losers <- merged %>% filter(change < 0) %>% mutate(loss = -change)
  gainers <- merged %>% filter(change > 0) %>% mutate(gain = change)
  
  if(nrow(losers) > 0 && nrow(gainers) > 0) {
    total_gain <- sum(gainers$gain)
    transfer_links <- crossing(losers, gainers) %>%
      mutate(
        value = loss * (gain / total_gain),
        source = paste(year_current, LU_class.x, sep = "_"),
        target = paste(year_next, LU_class.y, sep = "_")
      ) %>%
      select(source, target, value)
    
    links <- bind_rows(links, keep_links, transfer_links)
  } else {
    links <- bind_rows(links, keep_links)
  }
}

# 创建节点数据（每个节点为"年份_土地类型"格式）
nodes <- links %>%
  select(source, target) %>%
  pivot_longer(cols = everything(), names_to = "type", values_to = "name") %>%
  distinct(name) %>%
  mutate(id = row_number() - 1) # networkD3使用0起始索引

# 将链路中的源/目标转换为节点ID
links <- links %>%
  left_join(nodes, by = c("source" = "name")) %>%
  rename(IDsource = id) %>%
  left_join(nodes, by = c("target" = "name")) %>%
  rename(IDtarget = id) %>%
  select(IDsource, IDtarget, value)

# 绘制桑基图
sankeyNetwork(
  Links = links,
  Nodes = nodes,
  Source = "IDsource",
  Target = "IDtarget",
  Value = "value",
  NodeID = "name",
  fontSize = 12,
  nodeWidth = 30,
  margin = list(left = 150, right = 150) # 调整边距避免名称截断
)

代码说明 #

1. 数据转换：将宽格式数据转为长格式，提取年份数字便于后续处理。
2. 链路构造：
   • 保持链路：同一土地类型在相邻年份的存量流转，取两年中较小的比例值。
   • 转移链路：模拟土地类型的流转，将减少类型的流失量按增加类型的增益比例分配。
3. 节点处理：为每个年份的土地类型创建唯一节点（格式如2000_Cropland），确保桑基图按年份层级展示。
4. 绘图调整：设置边距避免节点名称被截断，优化可视化效果。

内容的提问来源于stack exchange，提问作者Lily Nature