使用terra::extract处理矢量数据时如何保留XY坐标？

我有一组从栅格提取的XY坐标点，想要把它们和所属行政区划（比如点a在纽约、点b在佛蒙特州）关联起来。用terra包的extract函数可以实现点和栅格/矢量数据的关联，但处理矢量数据时，没法用xy=TRUE参数保留XY坐标（这个参数在处理栅格时有效），会报错unused argument (xy = TRUE)。想问处理矢量数据时怎么保留坐标，或者有没有更优化的关联方法？

library(terra)

# 加载示例数据
r <- rast(system.file("ex/elev.tif", package="terra")) # 栅格数据
v <- vect(system.file("ex/lux.shp", package="terra"))  # 矢量行政区划数据

# 生成区域内的采样点
points <- spatSample(r, 5, as.points=TRUE)    
# 查看点坐标
crds(points)
#>             x        y
#> [1,] 5.929167 49.67083
#> [2,] 6.279167 49.94583
#> [3,] 6.487500 49.84583
#> [4,] 6.004167 49.49583
#> [5,] 6.379167 49.62083

# 关联点与矢量数据，但结果不含XY坐标
terra::extract(v, points)
#>   id.y ID_1       NAME_1 ID_2           NAME_2 AREA    POP
#> 1    1    3   Luxembourg    8         Capellen  185  48187
#> 2    2   NA         <NA>   NA             <NA>   NA     NA
#> 3    3   NA         <NA>   NA             <NA>   NA     NA
#> 4    4    3   Luxembourg    9 Esch-sur-Alzette  251 176820
#> 5    5    2 Grevenmacher   12     Grevenmacher  210  29828

# 处理矢量数据时用xy=TRUE会报错
terra::extract(v, points, xy=TRUE)
#> Error in .local(x, y, ...) : unused argument (xy = TRUE)

# 处理栅格数据时xy=TRUE正常工作
terra::extract(r, points, xy=TRUE)
#>   ID elevation        x        y
#> 1  1       323 5.929167 49.67083
#> 2  2        NA 6.279167 49.94583
#> 3  3        NA 6.487500 49.84583
#> 4  4       358 6.004167 49.49583
#> 5  5       283 6.379167 49.62083

解决方案 #

方法1：手动提取坐标并合并结果 #

extract处理矢量时返回的结果里有id.y字段，对应输入点的序号，我们可以提取点的坐标后，通过这个序号合并：

# 提取矢量关联结果
vec_extract <- terra::extract(v, points)
# 提取点的坐标并转为数据框，添加id列对应id.y
point_crds <- as.data.frame(crds(points))
point_crds$id.y <- seq(nrow(point_crds))
# 合并坐标与矢量属性
result <- merge(vec_extract, point_crds, by = "id.y")
print(result)
#>   id.y ID_1       NAME_1 ID_2           NAME_2 AREA    POP        x        y
#> 1    1    3   Luxembourg    8         Capellen  185  48187 5.929167 49.67083
#> 2    2   NA         <NA>   NA             <NA>   NA     NA 6.279167 49.94583
#> 3    3   NA         <NA>   NA             <NA>   NA     NA 6.487500 49.84583
#> 4    4    3   Luxembourg    9 Esch-sur-Alzette  251 176820 6.004167 49.49583
#> 5    5    2 Grevenmacher   12     Grevenmacher  210  29828 6.379167 49.62083

方法2：使用terra::intersect直接关联 #

intersect函数可以直接将点与矢量图层相交，结果会保留点的几何信息（包含坐标）和矢量的属性：

# 点与矢量相交，结果为SpatVector对象
intersect_result <- terra::intersect(points, v)
# 转为数据框查看
as.data.frame(intersect_result)
#>        x        y ID_1       NAME_1 ID_2           NAME_2 AREA    POP
#> 1 5.929167 49.67083    3   Luxembourg    8         Capellen  185  48187
#> 2 6.004167 49.49583    3   Luxembourg    9 Esch-sur-Alzette  251 176820
#> 3 6.379167 49.62083    2 Grevenmacher   12     Grevenmacher  210  29828

注意：intersect只会返回落在矢量范围内的点，不在范围内的点会被过滤掉。如果需要保留所有点（包括无匹配的），可以用left_join：

# 将点转为带ID的数据框
points_df <- cbind(id = seq(nrow(points)), as.data.frame(crds(points)))
# 将矢量属性转为数据框
v_df <- as.data.frame(v)
# 空间连接（左连接保留所有点）
joined_result <- terra::left_join(points, v, by = "geometry")
as.data.frame(joined_result)
#>        x        y id ID_1       NAME_1 ID_2           NAME_2 AREA    POP
#> 1 5.929167 49.67083  1    3   Luxembourg    8         Capellen  185  48187
#> 2 6.279167 49.94583  2   NA         <NA>   NA             <NA>   NA     NA
#> 3 6.487500 49.84583  3   NA         <NA>   NA             <NA>   NA     NA
#> 4 6.004167 49.49583  4    3   Luxembourg    9 Esch-sur-Alzette  251 176820
#> 5 6.379167 49.62083  5    2 Grevenmacher   12     Grevenmacher  210  29828

方法3：使用terra::join进行属性连接 #

如果已经通过extract得到了关联结果，也可以用join函数将坐标合并进去，逻辑和方法1类似，但更简洁：

# 给points添加id字段
points$id <- seq(nrow(points))
# 提取矢量属性
vec_attr <- terra::extract(v, points)
# 合并点的属性（含坐标）与矢量属性
final_result <- merge(as.data.frame(points), vec_attr, by.x = "id", by.y = "id.y")
print(final_result)
#>   id        x        y ID_1       NAME_1 ID_2           NAME_2 AREA    POP
#> 1  1 5.929167 49.67083    3   Luxembourg    8         Capellen  185  48187
#> 2  2 6.279167 49.94583   NA         <NA>   NA             <NA>   NA     NA
#> 3  3 6.487500 49.84583   NA         <NA>   NA             <NA>   NA     NA
#> 4  4 6.004167 49.49583    3   Luxembourg    9 Esch-sur-Alzette  251 176820
#> 5  5 6.379167 49.62083    2 Grevenmacher   12     Grevenmacher  210  29828

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