sf | 空间矢量对象的“聚合”操作

前面介绍空间连接函数时（sf | 空间矢量对象的属性连接方法）提到，对于数值型属性，要想将属性连接后的几何要素在sf对象中整合成一行，需要额外执行“聚合”操作。所谓“聚合”，可以包括“要素融合”和“属性加权”两个部分，前者针对sf对象的几何信息而言，后者针对属性信息而言。空间矢量对象的“聚合”操作类似于普通数据框的分类汇总操作，因此也可以使用dplyr包的summarize()函数。除此之外，sf包还可以通过一些函数实现各式的功能。

“融合”（dissolve）操作在ArcGIS中的示意图如下图。在ArcGIS中，可以通过指定字段，将该字段具有相同值的要素聚合成单个要素。

1 要素整体化

通过前篇sf | 空间矢量对象的属性连接方法的介绍，我们知道sf包中的一些函数是针对sf对象中的各几何要素进行运算的。但有时候，我们需要将sf对象看作一个整体参与运算，这就需要使用以下函数：

st_combine(x)
st_union(x, y, by_feature = FALSE, is_coverage = FALSE)

• x和y是sf、sfc或sfg对象，其中y可以省略;

• 当y不缺省时，st_union()函数的参数by_feature会失效。

当y缺省且st_union()函数的参数by_feature = FALSE时，两个函数的功能都是将x中的多个要素整合成单要素，同时丢弃x的属性信息。

st_combine()不会丢弃内部的几何边界信息，输出对象是“复合要素”，而st_union()函数会融合掉内部的几何边界信息，输出对象一般情况下是“简单要素”

比较下面代码中输出对象a1和a2数据结构的差异：

library(sf)
library(tidyverse)
# 导入示例数据
bj <-  read_sf("G:/temp/中国行政区划（2019）/县.shp") %>%
  filter(省 == "北京市")
hb <- read_sf("G:/temp/中国行政区划（2019）/市.shp") %>%
  filter(省 == "河北省")
dim(bj)
dim(hb)

a1 <- st_combine(bj)
a2 <- st_union(bj)
length(a1)
length(a2)
class(a1)
class(a2)

# 部分输出结果
> dim(bj)
[1] 16  8
> dim(hb)
[1] 11  6

> length(a1)
[1] 1
> length(a2)
[1] 1
> class(a1)
[1] "sfc_MULTIPOLYGON" "sfc"             
> class(a2)
[1] "sfc_POLYGON" "sfc"   

可以发现，a1和a2的长度都是1，但前者的几何构成是“MULTIPOLYGON”（复合要素），后者是“POLYGON”（简单要素）。可以通过绘图函数更直观地展示二者的区别:

par(mfrow = c(1,2))
plot(a1, main = "st_combine(bj)")
plot(a2, main = "st_union(bj)")

当sf对象本身呈破碎状态时，st_union()的输出对象也会是“复杂要素”

以河北省为例：

b1 <- st_combine(hb)
b2 <- st_union(hb)
length(b1)
length(b2)
class(b2)
plot(b2, main = "st_union(hb)")

# 输出结果
> class(b2)
[1] "sfc_MULTIPOLYGON" "sfc" 

sf包中的st_cast()函数用于转换sf对象类型，可以将st_combine()和st_union()输出的单个复合要素重新转换为多个简单要素。

st_cast(x, to, ...)

• x为待转换类型的对象，可以是sf、sfc和sfg对象；

• 参数to为转换后的对象类型，包括POINT、MULTIPOINT、LINESTRING、MULTILINESTRING、POLYGON、MULTIPOLYGON和GEOMETRYCOLLECTION等若干类。

对于包含复杂要素的sf对象，st_combine()的输出对象并不一定能通过st_cast()还原，st_union()的输出对象转换后也不一定是单要素

c1 <- st_cast(a1, "POLYGON")
c2 <- st_cast(b2, "POLYGON")
length(c1)
length(c2)
par(mfrow = c(1,2))
plot(c1)
text(st_coordinates(st_centroid(c1)),
     labels = seq(length(c1)), col = "red")
plot(c2)
text(st_coordinates(st_centroid(c2)),
     labels = seq(length(c2)), col = "red")

由于朝阳区在顺义区有一块飞地（即首都机场，图中标4的位置），使得转换后的c1包含17个简单要素，而实际上北京只有16个区县；由于河北有一些与主体地区分离的飞地和岛屿，使得c2包含多个要素（图中1和2的标记重叠了）。

当y缺省且st_union()函数的参数by_feature = TRUE时，st_union()的功能是融合各个单要素内部的几何边界，同时各要素的属性信息不会丢失

d1 <- st_union(bj, by_feature = T)
class(d1)
dim(d1)

# 部分输出结果
> class(d1)
[1] "sf"         "tbl_df"     "tbl"        "data.frame"
> dim(d1)
[1] 16  8

当y不缺省时，st_union()的功能是将x和y中的各要素两两联合，属性也保留

e1 <- st_union(bj, hb)
dim(e1)

# 部分输出结果
> dim(e1)
[1] 176  13

可以看出，e1的大小是bj和hb的行数相乘，列数相加后减1（几何信息合并为一列）。

2 分类汇总

st_combine()和st_union()并不能实现类似于ArcGIS中的“融合”操作，sf包主要通过继承其他工具包中的分类汇总函数来实现这一功能，如dplyr包的summarise()和stats包中的aggregate()函数。

与空间连接函数连用，用于要素属性汇总

使用空间连接函数进行属性连接后，一个要素在输出对象中常占据多行（详见sf | 空间矢量对象的属性连接方法）

st_read(system.file("shape/nc.shp", package="sf")) %>%
  select(FIPS) %>%
  st_transform(2264) -> nc
st_sf(value = 1:100, geom = st_make_grid(st_as_sfc(st_bbox(nc)))) -> ncgrid
ncjoin <- st_join(nc, ncgrid)
ncjoin

# 部分输出结果
> ncjoin
Simple feature collection with 349 features and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 406265 ymin: 48359.7 xmax: 3052877 ymax: 1044143
projected CRS:  NAD83 / North Carolina (ftUS)
First 10 features:
     FIPS value                       geometry
1   37009    84 MULTIPOLYGON (((1270813 913...
1.1 37009    93 MULTIPOLYGON (((1270813 913...
1.2 37009    94 MULTIPOLYGON (((1270813 913...
2   37005    94 MULTIPOLYGON (((1340553 959...
3   37171    84 MULTIPOLYGON (((1570586 910...
3.1 37171    85 MULTIPOLYGON (((1570586 910...
3.2 37171    94 MULTIPOLYGON (((1570586 910...
3.3 37171    95 MULTIPOLYGON (((1570586 910...
4   37053    90 MULTIPOLYGON (((2881206 948...
4.1 37053    99 MULTIPOLYGON (((2881206 948...

可以看出nc中许多FIPS（县编码）对应多个ncgrid中的value，假设使用平均值作为每个县最终的value属性值：

• 使用summarise()函数，需要与group_by()连用

ncjoin %>%
  group_by(FIPS) %>%
  mutate(value = mean(value))

# 部分输出结果
> ncjoin %>%
+   group_by(FIPS) %>%
+   mutate(value = mean(value))
Simple feature collection with 349 features and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 406265 ymin: 48359.7 xmax: 3052877 ymax: 1044143
projected CRS:  NAD83 / North Carolina (ftUS)
# A tibble: 349 x 3
# Groups:   FIPS [100]
   FIPS  value                                                geometry
 * <fct> <dbl>                         <MULTIPOLYGON [US_survey_foot]>
 1 37009  90.3 (((1270813 913326.4, 1251094 927717.9, 1244873 928273.~
 2 37009  90.3 (((1270813 913326.4, 1251094 927717.9, 1244873 928273.~
 3 37009  90.3 (((1270813 913326.4, 1251094 927717.9, 1244873 928273.~
 4 37005  94   (((1340553 959379.4, 1340433 964500, 1334126 973974.4,~
 5 37171  89.5 (((1570586 910379.4, 1564742 914896.4, 1546917 915898.~
 6 37171  89.5 (((1570586 910379.4, 1564742 914896.4, 1546917 915898.~
 7 37171  89.5 (((1570586 910379.4, 1564742 914896.4, 1546917 915898.~
 8 37171  89.5 (((1570586 910379.4, 1564742 914896.4, 1546917 915898.~
 9 37053  96.3 (((2881206 948550.7, 2878541 955075.8, 2873989 954300.~
10 37053  96.3 (((2881206 948550.7, 2878541 955075.8, 2873989 954300.~
# ... with 339 more rows

• 使用aggregate()函数

aggregate()函数是基础包stats中的函数，语法结构如下：

aggregate(x, by, FUN, ..., simplify = TRUE, drop = TRUE)

ncjoin %>%
  aggregate(by = list(.$FIPS), FUN = mean)

# 部分输出结果
> ncjoin %>%
+   aggregate(by = list(.$FIPS), FUN = mean)
Simple feature collection with 100 features and 3 fields
Attribute-geometry relationship: 0 constant, 2 aggregate, 1 identity
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 406265 ymin: 48359.7 xmax: 3052877 ymax: 1044143
projected CRS:  NAD83 / North Carolina (ftUS)
First 10 features:
   Group.1 FIPS    value                       geometry
1    37001   NA 81.00000 POLYGON ((1927130 771026.9,...
2    37003   NA 69.00000 POLYGON ((1374531 742546.1,...
3    37005   NA 94.00000 POLYGON ((1340553 959379.4,...
4    37007   NA 40.50000 POLYGON ((1723928 386344.2,...
5    37009   NA 90.33333 POLYGON ((1270813 913326.4,...
6    37011   NA 78.00000 POLYGON ((1129642 815436.4,...
7    37013   NA 59.00000 MULTIPOLYGON (((2712011 562...
8    37015   NA 83.50000 POLYGON ((2656952 775025.1,...
9    37017   NA 27.00000 POLYGON ((2222874 235545.1,...
10   37019   NA 12.50000 POLYGON ((2104506 72539.63,...
There were 50 or more warnings (use warnings() to see the first 50)

单独使用，用于某些字段值相同的要素之间的边界融合

比如通过融合同一地级行政区下的县级行政区之间的边界得到地级行政区的矢量数据：

ah <-  read_sf("G:/temp/中国行政区划（2019）/县.shp") %>%
  filter(省 == "安徽省")
# 边界融合
ah %>%
  group_by(市代码) %>%
  summarise() -> ahcity
# 地图可视化
par(mfrow = c(1,2))
plot(st_geometry(ah), main = "ah")
plot(st_geometry(ahcity), main = "ahcity")

3 面积加权

面积加权可以看作是属性连接的一种方法，sf包中提供的st_interpolate_aw()函数可以通过面积加权连接属性值；另外，该函数也可用在summarise()函数后计算融合后的sf对象的属性值。语法结构如下：

st_interpolate_aw(x, to, extensive, ...)

• x为属性来源对象，需为sf格式；

• to为属性接收对象，可以为sf或sfc格式；

• extensive为逻辑参数，TRUE表示x的属性以求和的方式赋值给参数to，FALSE表示属性以平均值的形式赋值给参数to。

使用面积加权函数时只能针对数值型属性，且x中只能保留待连接的属性，无关属性需要提前删去

比如将县级单位的属性聚合到地级单位：

ah %>%
  mutate(value = 1:dim(ah)[1]) %>%
  select(value) -> ahvalue
ahcity2 <- st_interpolate_aw(ahvalue, ahcity, extensive = T)
ahcity2

# 部分输出结果
> ahcity2
Simple feature collection with 16 features and 1 field
Attribute-geometry relationship: 0 constant, 1 aggregate, 0 identity
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 114.8785 ymin: 29.39519 xmax: 119.6452 ymax: 34.65234
geographic CRS: CGCS_2000
First 10 features:
       value                       geometry
1   45.00000 POLYGON ((117.5679 31.26984...
2  108.00000 POLYGON ((118.5703 31.29403...
3  147.00000 POLYGON ((117.2921 32.81864...
4  196.00000 POLYGON ((117.087 32.72772,...
5  207.00000 POLYGON ((118.8781 31.53705...
6  158.00000 MULTIPOLYGON (((116.8852 33...
7  174.00000 MULTIPOLYGON (((117.2939 30...
8  505.00000 POLYGON ((116.9 30.35936, 1...
9   17.69134 POLYGON ((118.3465 30.11127...
10 532.00000 POLYGON ((118.415 32.08747,...

---