如何将多年份土地类型栅格加入循环，整合至df.area数据框

我有一份2050年的全球土地类型投影栅格（用颜色区分1-6类），已经裁剪到全球海草分布多边形数据集（cmu）。我通过for循环计算了海草多边形上各土地类型（1.水域、2.森林、3.城市、4.裸地、5.农田、6.草地）的面积。现在想把2060、2070、2100年的额外全球土地类型栅格加入df.area数据框，以便分析/绘制土地覆盖类型随时间的变化。（注：代码仅筛选澳大利亚区域测试，因为全球计算耗时太长）

原测试代码 #

##CALCULATE AREA OF RASTER CLASSES###
library(terra)
library(sf)
library(tidyverse)
cmu <- st_read('units-attributes-wgs84L2.gpkg')
land <- rast('7landtypes/SSP1_RCP19/global_SSP1_RCP19_2050.tif')
plot(land)

cmu2 <- filter(cmu, TERRITORY1 == 'Australia') %>% st_transform("+proj=cea +lat_ts=30 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

land.sub <- crop(land, vect(cmu2))

tmp <- list() # tmp object for storing results of loop as a list

system.time(
  for(i in 1:nrow(cmu2)){#loop through cmus and calculate area of diff. land classes
    s <- cmu2[i,] # subset a polygon from the cmu dataset
    rs <- crop(land.sub, vect(s), mask = TRUE)
    if(is.na(minmax(rs)[2]) == FALSE){
      tmp[[i]] <- data.frame(class = unname(values(rs)),
                             area_ha = unname(values(cellSize(rs, unit="ha")))) %>% 
        filter(!is.na(class)) %>% 
        group_by(class) %>% 
        summarise(area_ha = sum(area_ha, na.rm = T)) %>% 
        mutate(unit_ID = st_drop_geometry(s)[,1],
               raster = names (rs))
    }else{
      next
    } 
  } 
)

df.area <- do.call(rbind, tmp)

head(df.area)
# save
write.csv(df.area, 'land-class-area-cmu-2050.csv', row.names = F)

核心思路是新增年份循环嵌套进现有流程，一次性处理所有年份的栅格，最终生成带时间维度的统一数据框。以下是优化后的代码：

步骤1：整理所有年份的栅格路径 #

先把需要处理的栅格文件路径整理成列表，同时提取年份信息：

# 定义所有年份的栅格文件路径（根据实际存储路径调整）
raster_paths <- c(
  "7landtypes/SSP1_RCP19/global_SSP1_RCP19_2050.tif",
  "7landtypes/SSP1_RCP19/global_SSP1_RCP19_2060.tif",
  "7landtypes/SSP1_RCP19/global_SSP1_RCP19_2070.tif",
  "7landtypes/SSP1_RCP19/global_SSP1_RCP19_2100.tif"
)
# 从文件名中提取年份（确保文件名格式一致，否则需调整提取规则）
years <- stringr::str_extract(raster_paths, "\\d{4}")

步骤2：重构双层循环逻辑 #

外层循环遍历年份，内层循环处理多边形，同时给每条记录添加年份标签：

library(terra)
library(sf)
library(tidyverse)

# 读取并预处理海草多边形数据
cmu <- st_read('units-attributes-wgs84L2.gpkg')
cmu2 <- filter(cmu, TERRITORY1 == 'Australia') %>% 
  st_transform("+proj=cea +lat_ts=30 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

# 创建空列表存储所有年份的计算结果
all_results <- list()
result_idx <- 1

system.time({
  # 外层循环：遍历每个年份的栅格
  for (y in seq_along(raster_paths)) {
    current_year <- years[y]
    land <- rast(raster_paths[y])
    land.sub <- crop(land, vect(cmu2))
    
    # 内层循环：遍历每个海草多边形
    for(i in 1:nrow(cmu2)){
      s <- cmu2[i,]
      rs <- crop(land.sub, vect(s), mask = TRUE)
      
      # 跳过无有效栅格数据的多边形
      if(is.na(minmax(rs)[2])) next
      
      # 计算当前多边形各土地类型面积
      area_df <- data.frame(
        class = unname(values(rs)),
        area_ha = unname(values(cellSize(rs, unit="ha")))
      ) %>% 
        filter(!is.na(class)) %>% 
        group_by(class) %>% 
        summarise(area_ha = sum(area_ha, na.rm = T)) %>% 
        mutate(
          unit_ID = st_drop_geometry(s)[,1],
          year = current_year,  # 新增年份字段，关联时间维度
          raster = names(rs)
        )
      
      all_results[[result_idx]] <- area_df
      result_idx <- result_idx + 1
    }
  }
})

# 合并所有结果为一个完整数据框
df.area <- do.call(rbind, all_results)

# 查看结果示例
head(df.area)

# 保存全年份数据
write.csv(df.area, 'land-class-area-cmu-all-years.csv', row.names = F)

关键优化说明 #

• 新增year字段：确保每条面积记录都绑定对应年份，为后续时间序列分析、可视化提供核心维度
• 双层循环结构：避免重复加载预处理多边形数据，提升计算效率
• 保留无效数据过滤：跳过无栅格数据的多边形，避免报错
• 统一结果存储：用单个列表收集所有年份的计算结果，最后一次性合并，减少内存占用

后续分析示例（绘制土地类型变化） #

以城市用地（类别3）为例，绘制各海草区域的面积变化折线图：

df.area %>% 
  filter(class == 3) %>% 
  ggplot(aes(x = year, y = area_ha, group = unit_ID)) +
  geom_line(alpha = 0.6) +
  labs(title = "澳大利亚海草区域城市用地面积变化", x = "年份", y = "面积(公顷)") +
  theme_minimal()

内容的提问来源于stack exchange，提问作者Alex White