如何使用terra为提取的数据框添加shapefile的state列？

问题描述：
我通过terra::extract函数，利用shapefile中的多边形从tif文件提取数据，代码如下：

library(terra)

tif = rast('data/snow.tif')
shape = vect('data/canada_provinces.shp')

x = extract(tif, shape, na.rm = T)
head(x)
# ID       values
# 1        100
# 1        100
# 1        100
# 1        101
# 1        101
# 2        150
# 2        150
# 2        150
# ...
# 90       101
# 90       101
# 90       101

希望为结果数据框添加shapefile中的state列，得到如下格式：

# ID       values     states
# 1        100        AK
# 1        100        AK
# 1        100        AK
# 1        101        AK
# 1        101        AK
# 2        150        NT
# 2        150        NT
# 2        150        NT
# ...
# 90       101        BC
# 90       101        BC
# 90       101        BC

但使用terra::extract的layer参数无法实现该需求，请问如何为提取得到的数据框添加实际的state列？

解决方案 #

terra::extract返回的ID列对应输入矢量shape中的行序号，我们可以基于这个序号关联shapefile的属性列，具体步骤如下：

1. 从shapefile中提取ID（行号）与state的映射关系：

# 创建映射表，ID为shape的行索引，states为对应的省份/州名称
state_map <- data.frame(
  ID = seq_len(nrow(shape)),
  states = shape$state
)

2. 将映射表与提取结果合并：

# 使用merge函数完成关联，保留所有提取结果
x_with_state <- merge(x, state_map, by = "ID", all.x = TRUE)

# 查看结果
head(x_with_state)

如果习惯使用dplyr，也可以用left_join实现：

library(dplyr)

state_map <- tibble(ID = 1:nrow(shape), states = shape$state)
x_with_state <- x %>% left_join(state_map, by = "ID")

合并后的结果就会包含states列，每个栅格值都会对应所属的省份/州名称。

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