sf | 空间矢量对象的属性连接方法

在处理空间矢量对象时，有时需要为它们增加新的属性数据。属性数据可以有两个来源：一是根据矢量对象各要素已有的数据进行信息匹配；二是利用空间位置关系把其他矢量对象的属性数据传递过来。前者是非空间方法，后者是空间方法。

在ArcGIS中，Join Data的下拉框提供了这两种属性连接方式，如下图：

在R语言的sf工具包中：

• 非空间方法的思路与普通数据框的连接操作一样，主要使用dplyr包的join系列函数。相关函数方法在之前的推文中已经有所介绍：dplyr | 数据导入和预处理的常用函数。

• sf包中的st_join()函数则使用的是空间方法。这个函数依赖于一系列用于判断几何要素空间位置关系的函数，详见之前的推文：sf | 判断点线面等几何对象的空间位置关系。

1 非空间方法

前面说过，sf对象的数据结构是一个特殊的数据框，许多针对数据框的函数同样可以应用于它。dplyr的join系列函数可以将普通数据框内的数据通过共同变量传递给sf对象。

*_join.sf(x, y, by = NULL, copy = FALSE, suffix = c(".x", ".y"), ...)

• join系列函数包括left_join()、right_join()、full_join()、inner_join()、semi_join()和anti_join()等，函数功能见前推文dplyr | 数据导入和预处理的常用函数；

• 该系列函数不能用于两个矢量对象的属性连接，也就是说x可以是sf对象，但y只能为普通数据框；

• sf包的说明文档里为这些函数添加了后缀.sf，但在使用时需要去掉后缀，并且仍然需要提前加载tidyverse包或dplyr包。

基础包base中的merge()函数也可以实现类似的功能。语法结构如下：

merge(x, y, by = intersect(names(x), names(y)), by.x = by, by.y = by,
      all = FALSE, all.x = all, all.y = all,
      sort = TRUE, suffixes = c(".x",".y"), no.dups = TRUE, incomparables = NULL, ...)

• 与join系列函数一样，merge()函数也要求y是普通数据框；

• 当连接变量在x和y同名时可以使用by参数，该参数默认值是全部同名变量；如果连接变量不同名，则使用by.x和by.y分别指定变量名；

• all默认为FALSE，表示只保留x中能够与y匹配的几何要素（sf对象每行代表一个几何要素），这一点与join系列函数有所区别；当设置为TRUE时，则保留x中的所有几何要素，缺失值记为NA；all.x和all.y分别针对x和y做类似的处理，默认与参数all一致；

• sort控制按照匹配变量排序。

y必须为普通数据框，不能是sf对象，否则会报错

library(tidyverse)
library(sf)

nc <- st_read(system.file("shape/nc.shp", package = "sf"))
ync <- select(nc, FIPS) %>% mutate(var = rpois(100, 10))
class(ync)
ync2 <- st_drop_geometry(ync)  # 删去几何信息，变成普通数据框
class(ync2)

a1 <- left_join(nc, ync, by = "FIPS")
a2 <- merge(nc, ync)

# 部分输出结果
> class(ync)
[1] "sf"         "data.frame"

> class(ync2)
[1] "data.frame"

当连接变量在x和y中同名时，merge()的by参数有默认值

# 连接变量在nc和ync2中都叫FIPS
b1 <- left_join(nc, ync2, by = "FIPS")
b2 <- merge(nc, ync2)

# 修改ync2的FIPS的列名后再连接
ync3 <- rename(ync2, NewFIPS = FIPS)
b3 <- left_join(nc, ync3, by = c("FIPS" = "NewFIPS"))
b4 <- merge(nc, ync3, by.x = "FIPS", by.y = "NewFIPS")

除了连接函数，列拼接函数也可以用于属性连接，但这要求x和y对应行的顺序必须一致

• 当x和y均为sf对象时，推荐使用基础包base中的cbind()函数；

• 当y为普通数据框时，推荐使用dplyr中的bind_cols()函数；

• sf包自建的st_bind_cols()函数是通用函数，但该函数已被停止维护了。

# y为sf对象
c1 <- bind_cols(nc, ync)
c2 <- cbind(nc, ync)  # 推荐
c3 <- st_bind_cols(nc, ync) 

# y为普通数据框
d1 <- bind_cols(nc, ync2)  # 推荐
d2 <- rbind(nc, ync2)
d3 <- st_bind_cols(nc, ync2) 

使用行拼接函数，可以实现矢量数据的合并操作

合并(append)操作在ArcGIS中的示意图如下：

在R中，可以使用行拼接函数实现这一操作，并且使用dplyr中的bind_rows()函数和base中的rbind()函数效果一样。

library(tidyverse)
library(sf)

nc <- st_read(system.file("shape/nc.shp", package = "sf"))
min_nc <- filter(nc, AREA < 0.1)
max_nc <- filter(nc, AREA > 0.2)

e1 <- bind_rows(min_nc, max_nc)
e2 <- rbind(min_nc, max_nc)
identical(e1,e2)  # e1和e2完全相同

par(omi = c(0,0,0,0))
par(mfrow = c(1,2))
plot(st_geometry(min_nc), main = "输入1")
plot(st_geometry(max_nc), main = "输入2")
par(mfrow = c(1,1))
plot(st_geometry(e1), main = "输出")

2 空间方法

2.1 st_join()

sf包中的st_join()使用空间位置关系在两个矢量对象之间传递属性信息。其语法结构如下：

st_join(x, y, join = st_intersects, ...,
  suffix = c(".x", ".y"), left = TRUE, largest = FALSE
)

参数说明：

• 输入对象x和y均需是sf格式；

• join参数为判断点线面位置关系的函数，详见推文sf | 判断点线面等几何对象的空间位置关系，...为继承自这些函数的参数；

• 工作原理：以x和y中的几何要素为最小单位，按照参数join引用的空间位置判断函数判断x的各几何要素对y的各几何要素的位置关系，当符合该位置关系时就将y中的属性数据传递给对应x的几何要素；

• suffix参数用于x和y中有同名变量时区分来源；

• left为TRUE时，输出结果保留x的所有几何要素；当设置为FASLE时，输出结果只保留能够连接的几何要素；

• largest设置为TRUE时，若x中的要素能与多个y内的几何要素连接，则只将面积最大的几何要素信息传递给相应x的几何要素。

生成示例文件（参考官方说明文档）：

library(tidyverse)
library(sf)

st_read(system.file("shape/nc.shp", package="sf")) %>%
  select(FIPS) %>%
  st_transform(2264) -> nc

st_sf(label = apply(expand.grid(1:10, LETTERS[10:1])[,2:1], 1, paste0, collapse = " "),
    geom = st_make_grid(st_as_sfc(st_bbox(nc)))) -> ncgrid
ncgrid$col = sf.colors(10, categorical = TRUE, alpha = 0.3)

# 绘图
plot(st_geometry(nc), lwd = 1.5, key.pos = NULL, main = "nc")
plot(st_geometry(ncgrid), col = ncgrid$col, main = "ncgrid")
text(st_coordinates(st_centroid(ncgrid)), labels = ncgrid$label)

nc的属性变量记录的是美国北卡罗来纳州(North Carolina)各县的FIPS编码；ncgrid是根据nc的范围创建的方格型矢量对象，其包含属性变量labels和颜色col；nc和ncgrid分别有100行数据。

nc

ncgrid

x中某几何要素能与几个y中的几何要素相连接，输出结果中该几何要素就占据几行，每行对应一个y中几何要素的属性数据，这样输出结果的行数可能会大于x的行数

比如，nc和ncgrid均包含100行数据，但下面代码的输出结果远大于100行数据：

nc1 <- st_join(nc, ncgrid)
dim(nc1)

# 部分输出结果
> dim(nc1)
[1] 349   4

为了使输出结果中每个几何要素只占据一行，需要进行“聚合”操作：

• 对于数值型变量，可以使用summarise()或aggregate()函数，这会在后续文章中单独介绍；

• 对于非数值型变量，在使用st_join()时可以将largest参数设置为TRUE。

nc2 <- st_join(nc, ncgrid, largest = T)

# 颜色col
plot(st_geometry(nc2), lwd = 1.5, col = nc2$col)
plot(st_geometry(ncgrid), add = T)

left参数设置为FASLE可以将没有连接关系的几何要素删去

比如ncgrid有的质心不能落到任何一个nc的要素中：

# 获取ncgrid的质心
ncpoints <- st_centroid(ncgrid)

nc3 <- st_join(ncpoints, nc, left = F)
dim(nc3)

# 部分输出结果
> dim(nc3)
[1] 55  4

可以看出，输出的nc4只剩下55行数据了。

join参数可以根据需要选择相离、相切、相交和包含函数

# 使用包含函数
nc4 <- st_join(ncpoints, nc, join = st_within, left = F)
nc5 <- st_join(ncgrid, nc4, join = st_contains)
dim(nc5)

plot(st_geometry(nc), lwd = 1.5)
plot(st_geometry(ncgrid), col = nc5$col.y, add = T)

# 部分输出结果
> dim(nc5)
[1] 100   6

试比较和上幅图的区别

2.2 st_filter()

st_filter(x, y, ..., .predicate = st_intersects)

• 该函数是空间筛选函数，它的.predicate参数相当于st_join()的join参数；

• 该函数主要是将x中符合相应位置关系的几何要素筛选出来，而不会进行属性传递。

nc6 <- st_filter(ncpoints, nc, .predicate = st_within)
dim(ncpoints)
dim(nc6)

# 部分输出结果
> dim(ncpoints)
[1] 100   3
> dim(nc6)
[1] 55  3

可以看出nc5中行数相比于ncpoints变少了，但列数（变量数）未变。

st_filter()一般不会直接引用相离函数

由于一般x中的几何要素总会与y中的一些几何要素具有相离关系，因此直接使用相离函数达不到筛选的目的。

为了筛选出x中与y整体相离的几何要素，可以使用st_combine()函数将y转换成具有整体意义的对象。需要注意的是，转换后的sf对象会变成一个list，不再具备数据框的特征。

比如筛选出没有落入nc范围内的质心：

# 不能达到筛选目的
nc7 <- st_filter(ncpoints, nc, .predicate = st_disjoint)
dim(nc7)

# 可以达到筛选目的
nc8 <- st_filter(ncpoints, st_combine(nc), .predicate = st_disjoint)
dim(nc8)

 

# 部分输出结果
> dim(nc7)
[1] 100   3

> dim(nc8)
[1] 45  3