使用 R 语言绘制城市间专利合作申请数量网络图(二)

继续上次课的内容，今天我们继续讲解如何使用 ggraph 包的网络图层绘制城市间专利合作网络。

也可以补充学习系列课程「ggplot2 数据可视化」的网络图绘制课时：https://rstata.duanshu.com/#/brief/course/9d8e0cc791644376979d78edc73eb18a

首先还是加载 R 包和读取相关的数据，注意这里加载了 tidygraph 和 ggraph 包。

library(tidyverse)
library(ggspatial)
library(tidygraph)
library(ggraph)
library(sf)
library(ggnewscale)

haven::read_dta("2020年城市间各类型专利合作数量统计.dta") -> df
df

#> # A tibble: 12,316 × 7
#>     年份 城市1    城市2    合作申请专利数量 合作申请外观设计专利数量
#>    <dbl> <chr>    <chr>               <dbl>                    <dbl>
#>  1  2020 七台河市 北京市                  3                        0
#>  2  2020 七台河市 哈尔滨市               12                        0
#>  3  2020 万宁市   大庆市                  1                        0
#>  4  2020 万宁市   江门市                  1                        0
#>  5  2020 万宁市   海口市                  1                        0
#>  6  2020 三亚市   三沙市                  2                        0
#>  7  2020 三亚市   上海市                  1                        0
#>  8  2020 三亚市   九江市                  1                        0
#>  9  2020 三亚市   北京市                 14                        0
#> 10  2020 三亚市   南京市                  2                        0
#> # ℹ 12,306 more rows
#> # ℹ 2 more variables: 合作申请发明专利数量 <dbl>,
#> #   合作申请实用新型专利数量 <dbl>

# 中国地图通常使用这样的坐标系
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

# 线条
read_sf("chinaprov2021mini/chinaprov2021mini_line.shp") %>%
  filter(!str_detect(class, "_") & !class %in% c("胡焕庸线", "秦岭-淮河线")) %>%
  select(class) -> provlinemap
read_sf("chinacity2021mini/chinacity2021mini.shp") %>%
  filter(!is.na(省代码)) -> citymap

# 配色：https://tidyfriday.cn/colors

df %>%
  select(from = 城市1, to = 城市2, value = 合作申请专利数量) %>%
  filter(!from %in% "三沙市" & !to %in% "三沙市") -> df

df

#> # A tibble: 12,312 × 3
#>    from     to       value
#>    <chr>    <chr>    <dbl>
#>  1 七台河市 北京市       3
#>  2 七台河市 哈尔滨市    12
#>  3 万宁市   大庆市       1
#>  4 万宁市   江门市       1
#>  5 万宁市   海口市       1
#>  6 三亚市   上海市       1
#>  7 三亚市   九江市       1
#>  8 三亚市   北京市      14
#>  9 三亚市   南京市       2
#> 10 三亚市   哈尔滨市     3
#> # ℹ 12,302 more rows

计算城市质心位置：

# 城市质心位置
read_sf("2021行政区划/市.shp") %>%
  st_transform(mycrs) -> city

city %>%
  select(省, 市) %>%
  st_point_on_surface() -> city_centroid

bind_cols(
  city_centroid %>%
    st_drop_geometry(),
  city_centroid %>%
    st_coordinates() %>%
    as_tibble()
) -> city_centroiddf

city_centroiddf

#> # A tibble: 371 × 4
#>    省     市            X        Y
#>    <chr>  <chr>     <dbl>    <dbl>
#>  1 安徽省 安庆市 1079679. 3295407.
#>  2 安徽省 蚌埠市 1122942. 3589681.
#>  3 安徽省 亳州市 1028999. 3617694.
#>  4 安徽省 池州市 1170731. 3267887.
#>  5 安徽省 滁州市 1187495. 3536891.
#>  6 安徽省 阜阳市  967020. 3560606.
#>  7 安徽省 合肥市 1147083. 3438594.
#>  8 安徽省 淮北市 1062868. 3651788.
#>  9 安徽省 淮南市 1085199. 3509706.
#> 10 安徽省 黄山市 1248987. 3249687.
#> # ℹ 361 more rows

计算每个 from 出发的总 value，然后按照 sum 排序：

df %>%
  group_by(from) %>%
  mutate(sum = sum(value, na.rm = T)) %>%
  arrange(sum) %>%
  ungroup() -> df

这样按照由低到高的顺序排列的目的是为了让高 sum 的边绘制到上面。在 ggplot2 绘图中，最先被绘制的线条会被绘制在下面，最后绘制的线条会盖住前面的线条。在后面的代码中我会把高 sum 的一些边设置成彩色，通过这样的排序就可以让彩色的线条浮在上面不被遮盖。

把 df 转换网络数据：

df %>%
  as_tbl_graph() %>%
  activate(nodes) %>%
  left_join(city_centroiddf, by = join_by("name" == "市")) -> graphdf

注意，要先给 df 排序再把其转换成网络数据。这样 graphdf 数据依然会保留 df 的顺序，这样 df 的排序才会对后面网络图的绘制生效。

简化地图数据使得图表绘制的更快：

city %>%
  st_simplify(dTolerance = 2000) -> city

在正式绘图前，我们先了解主要用的几种图层。为了绘图代码的简洁，我把部分代码存放在了 basemap.R 中：

basemap <- function(plot) {
  plot +
    geom_sf(data = citymap, fill = NA,
          color = "gray10", linewidth = 0.05, show.legend = F) +
    geom_sf(data = provlinemap, aes(color = class,
                                    linewidth = class),
            fill = NA, show.legend = F) +
    scale_fill_manual(
      name = "",
      values = c("#fed439", "#709ae1", "#8a9197", "#d2af81", "#fd7446", "#d5e4a2", "#197ec0", "#f05c3b")
    ) +
    scale_color_manual(
      name = "",
      values = c("九段线" = "black",
                 "海岸线" = "#0055AA",
                 "小地图框格" = "black",
                 "省份" = "gray30")
    ) +
    scale_linewidth_manual(
      name = "",
      values = c("九段线" = 0.5,
                 "海岸线" = 0.2,
                 "小地图框格" = 0.2,
                 "省份" = 0.3)
    ) +
    new_scale("fill") +
    new_scale("color") +
    new_scale("linewidth") +
    annotation_scale(
      width_hint = 0.2,
      text_family = cnfont, text_face = "plain",
      pad_x = unit(0.3, "cm")
    ) +
    annotation_north_arrow(
      location = "tr", which_north = "false",
      width = unit(1.6, "cm"),
      height = unit(2, "cm"),
      style = north_arrow_fancy_orienteering(
        text_family = cnfont,
        text_face = "plain"
      )
    ) +
    theme(axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          legend.title.position = "top",
          legend.position = "inside",
          legend.position.inside = c(0.12, 0.2))
}

首先是 geom_edge_bundle_force() 图层。该图层是基于力导向算法（force-directed algorithm），模拟物理力（如弹簧力）来将边捆绑在一起。通过模拟节点和边之间的相互作用力，使得边在图中自然地聚集在一起。但是这个计算是非常耗时的，因此它会把计算结果临时存储起来，所以第一次运行的时候会很费时间，但是之后就快多了。

ggraph(graphdf, x = X, y = Y) %>%
  basemap() +
  geom_edge_bundle_force(width = 0.05, color = "black") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  labs(title = "geom_edge_bundle_force, width = 0.05") -> p1
ggsave("pic2020a.png", width = 10, height = 8.5, device = png)

geom_edge_bundle_path() 是基于路径捆绑（path bundling）算法，将边沿着预定义的路径进行捆绑。通常使用层次结构或树状结构来定义路径。这个速度会更快一些：

ggraph(graphdf, x = X, y = Y) %>%
  basemap() +
  geom_edge_bundle_path(width = 0.05, color = "black") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  labs(title = "geom_edge_bundle_path, width = 0.05") -> p2
ggsave("pic2020b.png", width = 10, height = 8.5, device = png)

geom_edge_bundle_minimal()：是最小生成树，最快，但是最粗糙：

ggraph(graphdf, x = X, y = Y) %>%
  basemap() +
  geom_edge_bundle_minimal(width = 0.05, color = "black",
                           max_distortion = 10) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  labs(title = "geom_edge_bundle_minimal, width = 0.05, max_distortion = 10") -> p3
ggsave("pic2020c.png", width = 10, height = 8.5, device = png)

下面我们使用 geom_edge_bundle_path() 绘制，然后给图表添加更多细节。这里细节太多，我还是不在讲义里面一句句解读了，还是推荐大家结合视频讲解观看。

df %>%
  group_by(from) %>%
  summarise(sum = sum(value, na.rm = T)) %>%
  ungroup() %>%
  arrange(desc(sum)) -> sumdf

col_pal <- c("#fed439", "#709ae1", "#8a9197", "#d2af81", "#fd7446", "#d5e4a2", "#197ec0", "#f05c3b", "#46732e", "#71d0f5")

graphdf %>%
  left_join(sumdf, by = join_by("name" == "from")) -> graphdf

as.integer(seq(range(df$value)[1],
               range(df$value)[2],
               length = 6)) -> myrange1
as.integer(seq(range(sumdf$sum)[1],
               range(sumdf$sum)[2],
               length = 6)) -> myrange2

# 添加更多细节

此处代码需下载讲义材料查看～

这样我们就绘制好了这幅地图，感兴趣的小伙伴也可以试试区县和省份的。

点击这里跳转到 RStata 短书平台获取附件：使用 R 语言绘制城市间专利合作申请数量网络图(二)