使用 R 语言绘制中国地图 + 空间网络图

前不久给大家分享过「上市公司前5大供应商和客户工商注册数据匹配结果（含经纬度和所处的省市区县）」数据，在数据介绍中我展示了「2022 年上市公司前 5 大客户与供应商地理分布」：

这幅图是使用 Stata 绘制，今天的课程中我们将一起学习下该如何在 Stata 中绘制。

在学习该课程前，需要预先学习下面两个课程：

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

今天的课程也将在之前的这两个课程的基础上进行讲解。

绘制地图和在地图上添加散点图的方法之前我们都已经讲解过了，因此我们今天就直接来看如何在地图上添加连接线。

显然，连接线需要的数据至少应该包含起始点的经纬度坐标，下面我们先构造绘制连接线所需的数据。

加载所需的 R 包：

library(tidyverse)
library(sf)

读取上市公司经纬度数据：

haven::read_dta("上市公司注册地址经纬度数据.dta") -> df0

中国地图通常使用这样的坐标系：

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"

由于绘图使用的底图是 mycrs 坐标系下的，所以下面我们需要把所有的经纬度都转换成 mycrs 坐标系下的坐标：

df0 %>%
  st_as_sf(coords = c("注册地址经度", "注册地址纬度"), crs = 4326) %>%
  st_transform(mycrs) -> df0

bind_cols(
  st_drop_geometry(df0),
  st_coordinates(df0) %>%
    as_tibble() %>%
    set_names(c("注册地址经度", "注册地址纬度"))
) -> df0

构造上市公司-供应商连接线数据：

# 构造上市公司-供应商连接线数据
haven::read_dta("2022年上市公司供应商信息.dta") %>%
  select(经度, 纬度, 省, 股票代码) %>%
  rename(供应商经度 = 经度,
         供应商纬度 = 纬度,
         供应商所在省 = 省) -> dfa

dfa %>%
  st_as_sf(coords = c("供应商经度", "供应商纬度"), crs = 4326) %>%
  st_transform(mycrs) -> dfa

bind_cols(
  st_drop_geometry(dfa),
  st_coordinates(dfa) %>%
    as_tibble() %>%
    set_names(c("供应商经度", "供应商纬度"))
) -> dfa

# 连接两个数据
dfa %>%
  inner_join(df0) %>%
  select(股票代码, 供应商经度, 供应商纬度, 供应商所在省, 注册地址经度, 注册地址纬度, 注册地址所在省份) %>%
  rename(上市公司经度 = 注册地址经度,
         上市公司纬度 = 注册地址纬度,
         上市公司所在省 = 注册地址所在省份) -> df1

df1

#> # A tibble: 2,262 × 7
#>    股票代码 供应商经度 供应商纬度 供应商所在省 上市公司经度 上市公司纬度
#>    <chr>         <dbl>      <dbl> <chr>               <dbl>        <dbl>
#>  1 000005     1608294.   3299398. 浙江省            944512.     2386745.
#>  2 000005      867836.   2441568. 广东省            944512.     2386745.
#>  3 000006      356692.   3660158. 陕西省            942786.     2386569.
#>  4 000006      935506.   2385749. 广东省            942786.     2386569.
#>  5 000006     1487424.   3526677. 江苏省            942786.     2386569.
#>  6 000006      768901.   3004543. 湖南省            942786.     2386569.
#>  7 000007     1552862.   3445405. 上海市            940278.     2386306.
#>  8 000007      823579.   2374488. 广东省            940278.     2386306.
#>  9 000007     1536324.   3437444. 上海市            940278.     2386306.
#> 10 000007     1315762.   3482064. 江苏省            940278.     2386306.
#> # ℹ 2,252 more rows
#> # ℹ 1 more variable: 上市公司所在省 <chr>

构造上市公司-客户连接线数据：

haven::read_dta("2022年上市公司客户信息.dta") %>%
  select(经度, 纬度, 省, 股票代码) %>%
  rename(客户经度 = 经度,
         客户纬度 = 纬度,
         客户所在省 = 省) -> dfb

dfb %>%
  st_as_sf(coords = c("客户经度", "客户纬度"), crs = 4326) %>%
  st_transform(mycrs) -> dfb

bind_cols(
  st_drop_geometry(dfb),
  st_coordinates(dfb) %>%
    as_tibble() %>%
    set_names(c("客户经度", "客户纬度"))
) -> dfb

dfb %>%
  inner_join(df0) %>%
  select(股票代码, 客户经度, 客户纬度, 客户所在省, 注册地址经度, 注册地址纬度, 注册地址所在省份) %>%
  rename(上市公司经度 = 注册地址经度,
         上市公司纬度 = 注册地址纬度,
         上市公司所在省 = 注册地址所在省份) -> df2

df2

#> # A tibble: 1,993 × 7
#>    股票代码 客户经度 客户纬度 客户所在省 上市公司经度 上市公司纬度
#>    <chr>       <dbl>    <dbl> <chr>             <dbl>        <dbl>
#>  1 000005   1431886. 3317730. 浙江省          944512.     2386745.
#>  2 000005   1513723. 3162952. 浙江省          944512.     2386745.
#>  3 000006    938641. 2386290. 广东省          942786.     2386569.
#>  4 000007    940338. 2386288. 广东省          940278.     2386306.
#>  5 000007    937929. 2385388. 广东省          940278.     2386306.
#>  6 000007    940290. 2386210. 广东省          940278.     2386306.
#>  7 000010   -256449. 2607587. 云南省          918984.     2384949.
#>  8 000010    855612. 2435194. 广东省          918984.     2384949.
#>  9 000010    145654. 2774150. 贵州省          918984.     2384949.
#> 10 000010   -229982. 2620613. 云南省          918984.     2384949.
#> # ℹ 1,983 more rows
#> # ℹ 1 more variable: 上市公司所在省 <chr>

构造散点数据：

bind_rows(
  df1 %>%
    select(上市公司经度, 上市公司纬度) %>%
    set_names("经度", "纬度") %>%
    mutate(class = "上市公司"),
  df1 %>%
    select(供应商经度, 供应商纬度) %>%
    set_names("经度", "纬度") %>%
    mutate(class = "供应商/客户")
) -> pointdf1

bind_rows(
  df2 %>%
    select(上市公司经度, 上市公司纬度) %>%
    set_names("经度", "纬度") %>%
    mutate(class = "上市公司"),
  df2 %>%
    select(客户经度, 客户纬度) %>%
    set_names("经度", "纬度") %>%
    mutate(class = "供应商/客户")
) -> pointdf2

# 转换成 sf 对象
library(sf)
pointdf1 %>%
  distinct() %>%
  st_as_sf(coords = c("经度", "纬度"), crs = mycrs) -> pointdfsf1

pointdf2 %>%
  distinct() %>%
  st_as_sf(coords = c("经度", "纬度"), crs = mycrs) -> pointdfsf2

在读取相关的地图数据：

# 读取小地图版本的中国城市地图
read_sf("chinacity2019mini/chinacity2019mini.shp") %>%
  filter(!is.na(省代码)) -> citymap

# 线条
read_sf("chinacity2019mini/chinacity2019mini_line.shp") %>%
  filter(class %in% c("九段线", "海岸线", "小地图框格")) %>%
  select(class) -> citylinemap

# 邻国地图
read_sf("china_neighboring/china_neighboring.shp") -> china_neighboring

为了让图表更美观，我们可以提取小地图部分的散点也绘制到图上：

# 小地图的范围
small_bbox <- st_bbox(c(xmin = 120000,
                        xmax = 1766004.1,
                        ymax = 2557786.0,
                        ymin = 320000),
                      crs = st_crs(mycrs)) %>%
  st_as_sfc()

# 提取小地图上的点
pointdfsf1 %>%
  st_intersection(small_bbox) -> pointdfsf1_small
pointdfsf2 %>%
  st_intersection(small_bbox) -> pointdfsf2_small

# 把这些点移动到小地图的位置
pointdfsf1_small %>%
  mutate(geometry = geometry * 0.5 + c(2100000, 1665139)) %>%
  sf::st_set_crs(mycrs) -> pointdfsf1_small
pointdfsf2_small %>%
  mutate(geometry = geometry * 0.5 + c(2100000, 1665139)) %>%
  sf::st_set_crs(mycrs) -> pointdfsf2_small

# 合并两部分
bind_rows(pointdfsf1, pointdfsf1_small) -> pointdfsfall1
bind_rows(pointdfsf2, pointdfsf2_small) -> pointdfsfall2

绘图区域（数据制作的时候预设好的）：

st_bbox(c(xmin = -2725586,
          xmax = 2982768,
          ymax = 6000000,
          ymin = 1800655),
        crs = st_crs(mycrs)) %>%
  st_as_sfc() -> plotbbox

然后就可以绘图了：

# 绘图
library(ggspatial)
library(ggnewscale)

# 设置绘图主题
theme_set(
  hrbrthemes::theme_ipsum(base_family = cnfont) +
    theme(axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(hjust = 0.5),
          plot.subtitle = element_text(hjust = 0.5),
          legend.position = "bottom")
)

# 供应商——上市公司
ggplot(citymap) +
  geom_sf(data = plotbbox, fill = "#BBD1EB") +
  geom_sf(fill = "white", color = "gray", linewidth = 0.01) +
  geom_sf(data = china_neighboring, fill = "#EDEDED") +
  stat_sf_coordinates(data = china_neighboring,
                      geom = "text", color = "gray",
                      aes(label = country_cn), family = cnfont,
                      size = 3) +
  geom_sf(data = citylinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  scale_color_manual(
    values = c("九段线" = "#A29AC4",
               "海岸线" = "#0055AA",
               "小地图框格" = "black")
  ) +
  scale_linewidth_manual(
    values = c("九段线" = 0.6,
               "海岸线" = 0.3,
               "小地图框格" = 0.3)
  ) +
  new_scale_color() +
  geom_sf(data = pointdfsfall1, aes(color = class),
          size = 0.5) +
  scale_color_manual(values = c("#fed439", "#709ae1"),
                     name = "") +
  new_scale_color() +
  geom_segment(data = df1,
             aes(x = 供应商经度, y = 供应商纬度,
                 xend = 上市公司经度, yend = 上市公司纬度,
                 color = 上市公司所在省), linewidth = 0.05,
             arrow = arrow(length = unit(0.01, "npc"), type = "closed"),
             show.legend = F) +
  scale_color_manual(values = paletteer::paletteer_d("ggsci::default_igv", n = 31)) +
  labs(title = "上市公司前 5 大供应商",
       subtitle = "绘图：微信公众号 RStata") +
  guides(color = guide_legend(nrow = 1)) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont,
                   pad_x = unit(1.2, "cm"),
                   pad_y = unit(1, "cm"),) +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_x = unit(1, "cm"),
    pad_y = unit(1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p1

ggsave("pic1.png", width = 9, height = 7.5, device = png)

上市公司和客户的连接线绘制方法类似：

# 上市公司——客户
ggplot(citymap) +
  geom_sf(data = plotbbox, fill = "#BBD1EB") +
  geom_sf(fill = "white", color = "gray", linewidth = 0.01) +
  geom_sf(data = china_neighboring, fill = "#EDEDED") +
  stat_sf_coordinates(data = china_neighboring,
                      geom = "text", color = "gray",
                      aes(label = country_cn), family = cnfont,
                      size = 3) +
  geom_sf(data = citylinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  scale_color_manual(
    values = c("九段线" = "#A29AC4",
               "海岸线" = "#0055AA",
               "小地图框格" = "black")
  ) +
  scale_linewidth_manual(
    values = c("九段线" = 0.6,
               "海岸线" = 0.3,
               "小地图框格" = 0.3)
  ) +
  new_scale_color() +
  geom_sf(data = pointdfsfall2, aes(color = class),
          size = 0.5) +
  scale_color_manual(values = c("#fed439", "#709ae1"),
                     name = "") +
  new_scale_color() +
  geom_segment(data = df2,
               aes(xend = 客户经度, yend = 客户纬度,
                   x = 上市公司经度, y = 上市公司纬度,
                   color = 上市公司所在省), linewidth = 0.05,
               arrow = arrow(length = unit(0.01, "npc"), type = "closed"),
               show.legend = F) +
  scale_color_manual(values = paletteer::paletteer_d("ggsci::default_igv", n = 31)) +
  labs(title = "上市公司前 5 大客户",
       subtitle = "绘图：微信公众号 RStata") +
  guides(color = guide_legend(nrow = 1)) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont,
                   pad_x = unit(1.2, "cm"),
                   pad_y = unit(1, "cm"),) +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_x = unit(1, "cm"),
    pad_y = unit(1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p2

ggsave("pic2.png", width = 9, height = 7.5, device = png)

最后还可以合并两幅图并使用相同的图例：

# 合并两幅图
library(patchwork)

p1 +
  p2 +
  plot_layout(guides = 'collect') &
  theme(legend.position = 'bottom')-> p

ggsave("pic3.png", width = 15, height = 7.5, device = png)

点击这里跳转到 RStata 短书平台获取附件：使用 R 语言绘制中国地图 + 空间网络图