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

前不久给大家分享了「1985～2024年省份间、城市间及区县间各类型专利合作数量面板数据」，在推文中我给大家展示了区县间、城市间和省份间的专利合作网络。今天我们以城市间专利申请合作网络为例讲解如何绘制这种空间网络图。

在本次课及下次课我们将介绍两种方法：

1. 使用 ggplot2 的基础图层绘制空间网络图，这是个很基础和直观的方法，非常有必要掌握；
2. 使用 ggraph 包的网络相关图层绘制，这个可以实现更多功能，也非常值得掌握：

基础 ggplot2 绘制	ggraph 包图层绘制

在学习本课程之前需要预先学习之前的课程：

使用 R 语言绘制历年中国省市区县地图（小地图版本+长版）: https://rstata.duanshu.com/#/brief/course/379d19770956478ebc4d919910fca74c

本课程中使用的地图数据也都是来自该课程。

首先加载所需的 R 包：

library(tidyverse)
library(ggspatial)
library(sf)
library(ggnewscale)

读取 2020 年城市间各类型专利合作数量统计数据：

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
provlinemap

#> Simple feature collection with 37 features and 1 field
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -2625586 ymin: 1825139 xmax: 2983002 ymax: 5921583
#> Projected CRS: +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
#> # A tibble: 37 × 2
#>    class                                                                geometry
#>    <chr>                                                   <MULTILINESTRING [m]>
#>  1 省份  ((1026818 3753771, 1026899 3753311, 1026900 3753103, 1026899 3752894, …
#>  2 省份  ((891889.1 2343649, 891865 2343557, 891648.7 2343372, 890923.8 2343325…
#>  3 省份  ((959992.3 4475646, 960020 4475637, 960065.6 4475631, 960093.2 4475637…
#>  4 省份  ((1298152 2534763, 1298134 2534747, 1298011 2534699, 1297947 2534713, …
#>  5 省份  ((95551.38 3786158, 95502.29 3786112, 95452.29 3786068, 95358.63 37861…
#>  6 省份  ((590293.9 2117685, 589768.4 2117046, 589358.1 2116967, 589057 2117038…
#>  7 省份  ((444678.3 2172854, 444444.2 2172566, 444074.1 2172455, 443583.3 21726…
#>  8 省份  ((9574.094 2605698, 9564.813 2605694, 9555.934 2605691, 9547.356 26056…
#>  9 省份  ((556006.8 1909575, 555989.1 1909520, 555843.4 1909450, 555794.2 19094…
#> 10 省份  ((1126379 4260466, 1126410 4260264, 1126368 4260023, 1126019 4259982, …
#> # ℹ 27 more rows

read_sf("chinacity2021mini/chinacity2021mini.shp") %>%
  filter(!is.na(省代码)) -> citymap
citymap

#> Simple feature collection with 371 features and 6 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -2625586 ymin: 1868655 xmax: 2962768 ymax: 5921583
#> Projected CRS: +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
#> # A tibble: 371 × 7
#>    省     省代码 市     市代码 class objid                              geometry
#>  * <chr>   <dbl> <chr>   <dbl> <chr> <dbl>                    <MULTIPOLYGON [m]>
#>  1 安徽省 340000 安庆市 340800 城市      1 (((1129079 3379225, 1129658 3379349,…
#>  2 安徽省 340000 蚌埠市 340300 城市      2 (((1115455 3634678, 1115645 3634689,…
#>  3 安徽省 340000 亳州市 341600 城市      3 (((973132.7 3682424, 973666.2 368238…
#>  4 安徽省 340000 池州市 341700 城市      4 (((1221172 3343232, 1220927 3342713,…
#>  5 安徽省 340000 滁州市 341100 城市      5 (((1209296 3606762, 1209300 3605694,…
#>  6 安徽省 340000 阜阳市 341200 城市      6 (((969637 3624545, 969637.9 3624505,…
#>  7 安徽省 340000 合肥市 340100 城市      7 (((1122388 3525961, 1122728 3526021,…
#>  8 安徽省 340000 淮北市 340600 城市      8 (((1081361 3696008, 1081801 3696328,…
#>  9 安徽省 340000 淮南市 340400 城市      9 (((1079161 3564117, 1079339 3563864,…
#> 10 安徽省 340000 黄山市 341000 城市     10 (((1238827 3310983, 1239824 3310536,…
#> # ℹ 361 more rows

地图配色可以使用这个网站：https://tidyfriday.cn/colors

仅仅保留需要的变量，三沙市的位置会影响图表的比例，而且专利申请的很少，还是从数据中去除了：

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

统计每个城市与其他城市的总合作量：

df %>%
  group_by(from) %>%
  summarise(value = sum(value, na.rm = T)) -> dfsum

dfsum

#> # A tibble: 360 × 2
#>    from     value
#>    <chr>    <dbl>
#>  1 七台河市    15
#>  2 万宁市       3
#>  3 三亚市      99
#>  4 三明市     239
#>  5 三门峡市   521
#>  6 上海市   16108
#>  7 上饶市     442
#>  8 东方市      30
#>  9 东莞市    4062
#> 10 东营市    1828
#> # ℹ 350 more rows

df %>%
  group_by(from) %>%
  mutate(sum = sum(value)) %>%
  ungroup() -> df

读取未经编辑的城市矢量数据，因为我们需要这个数据来计算各城市的质心：

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

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

把 df 和 city_centroiddf 匹配起来就可以得到每条合作线起点和终点的坐标了：

df %>%
  left_join(city_centroiddf, by = c("from" = "市")) %>%
  rename(fromlon = X, fromlat = Y) %>%
  left_join(city_centroiddf, by = c("to" = "市")) %>%
  rename(tolon = X, tolat = Y, prov = 省.x) %>%
  select(-省.y) -> df

df

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

简化地图数据，这样绘制地图更快一些：

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

然后准备两个数据：first100df 和 first10df，第一个用来绘制城市散点（只绘制 100 个），第二个用来添加文本标签（只添加 10 个）。太多就会密密麻麻的，反而效果很差。

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

first10df

#> # A tibble: 10 × 4
#>    from    fromlon  fromlat   sum
#>    <chr>     <dbl>    <dbl> <dbl>
#>  1 北京市  972929. 4384044. 54652
#>  2 上海市 1545609. 3420743. 16108
#>  3 南京市 1275570. 3476635. 11024
#>  4 深圳市  947535. 2396243. 10807
#>  5 广州市  876561. 2458588.  9926
#>  6 成都市 -106982. 3257911.  8467
#>  7 杭州市 1373275. 3256649.  8392
#>  8 苏州市 1470544. 3443155.  7845
#>  9 天津市 1027425. 4297302.  7410
#> 10 武汉市  885732. 3288717.  6576

然后就可以绘图了。由于细节很多，这里就不在一一解释了，感兴趣的小伙伴可以结合视频讲解学习：

ggplot() +
  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") +
  geom_curve(data = df,
             aes(x = fromlon, y = fromlat,
                 xend = tolon, yend = tolat,
                 color = from,
                 linewidth = value,
                 alpha = from),
             curvature = 0.3) +
  scale_linewidth_continuous(range = c(0.1, 0.5),
                             breaks = as.integer(seq(range(df$value)[1], range(df$value)[2], length = 4)),
                             # limits = c(0, 800),
                             name = "合作申请专利数量",
                             guide = guide_legend(direction = "horizontal",
                                                  nrow = 1,
                                                  order = 1,
                                                  byrow = T)) +
  scale_color_manual(values = rev(c("#fed439", "#709ae1", "#8a9197", "#d2af81", "#fd7446", "#d5e4a2", "#197ec0", "#f05c3b", "#46732e", "#71d0f5", rep("gray", length(unique(df$from)) - 10))),
                     breaks = rev(c(first10df$from, df %>% anti_join(first10df) %>% pull(from) %>% unique()))) +
  scale_alpha_manual(values = rev(c(rep(1, 10),
                                    rep(0.2, length(unique(df$from)) - 10))),
                     breaks = rev(c(first10df$from, df %>% anti_join(first10df) %>% pull(from) %>% unique()))) +
  guides(color = "none", alpha = "none") +
  new_scale(new_aes = "size") +
  geom_point(data = first100df, aes(x = fromlon,
                                   y = fromlat,
                                   size = sum),
             color = "gray10") +
  scale_size_continuous(range = c(0.01, 7),
                        breaks = as.integer(seq(range(first100df$sum)[1], range(first100df$sum)[2], length = 6)),
                        name = "总合作申请数量",
                        # limits = c(0, 8000),
                        guide = guide_legend(direction = "horizontal",
                                             nrow = 2,
                                             order = 2,
                                             byrow = T)) +
  ggrepel::geom_text_repel(data = first100df,
                           aes(x = fromlon,
                               y = fromlat,
                               label = from),
                           color = "gray10",
                           min.segment.length = 0.1,
                           family = cnfont,
                           size = 3) +
  theme_ipsum(base_family = cnfont) +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.title.position = "top") +
  theme(legend.position = "inside",
        legend.position.inside = c(0.12, 0.2),
        plot.background = element_rect(fill = "white", color = "white")) +
  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"
    )
  ) +
  labs(title = "2020 年各城市合作申请专利网络",
       subtitle = "数据处理&绘图：微信公众号 RStata",
       caption = "数据来源：国家知识产权局") -> p

ggsave("pic2020.png", width = 10, height = 8.5, device = png)

附件中也提供了区县和省份的绘图代码，感兴趣的小伙伴也可以拿来参考。

下次课我们再继续讲解如何使用 ggraph 包绘制类似的图表。

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