使用 R 语言绘制带邻国地区的中国地图

在之前「使用 R 语言绘制历年中国各省市区县地图（小地图版本 + 长版）」课程的基础上，我们今天再来学习下如何把邻国的区域也添加到中国地图上：

小地图版本

1. 在地图上添加散点

首先加载所需 R 包，读取地图数据：

library(tidyverse)
library(sf)

# 中国地图通常使用这样的坐标系
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"

# 1. 带小地图版本

# 绘图区域（数据制作的时候预设好的）
st_bbox(c(xmin = -2725586,
          xmax = 2982768,
          ymax = 6000000,
          ymin = 1800655),
        crs = st_crs(mycrs)) %>%
  st_as_sfc() -> plotbbox

# 读取小地图版本的中国城市地图
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

在地图上添加散点图，这里以瞪羚企业分布为例：

readxl::read_xlsx("瞪羚、独角兽、创新型企业经纬度数据.xlsx") %>%
  filter(!is.na(经度)) -> pointdf

# 转换成 sf 对象
pointdf %>%
  st_as_sf(coords = c("经度", "纬度"), crs = 4326) %>%
  st_transform(mycrs) -> pointdfsf

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

# 提取这个范围的点
pointdfsf %>%
  st_intersection(small_bbox) -> pointdfsf_small

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

# 合并两部分
bind_rows(pointdfsf, pointdfsf_small) -> pointdfsfall

然后就可以绘制地图了：

# 绘图
library(ggspatial)
library(ggnewscale)
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 = pointdfsfall, aes(color = 省),
          size = 0.5, show.legend = F) +
  scale_color_manual(values = paletteer::paletteer_d("ggsci::default_igv", n = 32)) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  labs(title = "瞪羚、独角兽、创新型企业城市分布",
       subtitle = "数据爬取&绘制：微信公众号 RStata",
       caption = "数据来源：瞪羚云网站") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_y = unit(0.1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p1

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

2. 填充地图

下面我们再演示下填充地图的绘制。

首先统计每个城市的公司数量：

# 统计每个城市的公司数量
pointdf %>%
  count(市, 市代码) -> citydf

# 和地图数据合并
citymap %>%
  left_join(citydf) %>%
  mutate(n = if_else(is.na(n), 0, n)) -> citymap2

然后就可以绘图了：

n = n + 1 是为了对数化之后 0 变成 -Inf，在图上会反映为缺失值。

# 绘图
citymap2 %>%
  mutate(n = n + 1) %>%
  ggplot() +
  geom_sf(data = plotbbox, fill = "#BBD1EB") +
  geom_sf(aes(fill = n), 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)
  ) +
  scico::scale_fill_scico(palette = "acton",
                          name = "公司数量",
                          trans = "log10") +
  new_scale_color() +
  geom_sf(data = pointdfsfall, aes(color = 省),
          size = 0.1, show.legend = F) +
  scale_color_manual(values = paletteer::paletteer_d("ggsci::default_igv", n = 32)) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  labs(title = "瞪羚、独角兽、创新型企业地理分布",
       subtitle = "数据爬取&绘制：微信公众号 RStata",
       caption = "数据来源：瞪羚云网站") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_y = unit(0.1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p2

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

也可以使用分段填色：

# 使用分段填色
citymap2$n %>%
  quantile(probs = 1:10/10, digits = 1) %>%
  as.integer() %>%
  unique() -> cutlist

cutlist

#> [1]     0     1     3     9    18    39    92 16712

citymap2 %>%
  mutate(group = cut(n, breaks = unique(cutlist),
                     include.lowest = T,
                     labels = c("<= 1", "1～3", "3～9", "9～18",
                                "19～39", "39～92", ">= 92"))) -> citymap3

citymap3 %>%
  ggplot() +
  geom_sf(data = plotbbox, fill = "#BBD1EB") +
  geom_sf(aes(fill = group), 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)
  ) +
  scico::scale_fill_scico_d(palette = "acton",
                            name = "公司数量",
                            end = 0.95) +
  new_scale_color() +
  geom_sf(data = pointdfsfall, aes(color = 省),
          size = 0.1, show.legend = F) +
  scale_color_manual(values = paletteer::paletteer_d("ggsci::default_igv", n = 32)) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  labs(title = "瞪羚、独角兽、创新型企业地理分布",
       subtitle = "数据爬取&绘制：微信公众号 RStata",
       caption = "数据来源：瞪羚云网站") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_y = unit(0.1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p3

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

长版地图

长版地图的绘制更简单一些，这里仅仅演示散点的绘制：

library(tidyverse)
library(sf)
library(ggspatial)
library(ggnewscale)
# 2. 长版地图

# 中国地图通常使用这样的坐标系
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"

# 绘图区域（数据制作的时候预设好的）
st_bbox(c(xmin = -2925762.0,
          xmax = 2507277.8,
          ymax = 6221888.6,
          ymin = 377031.1),
        crs = st_crs(mycrs)) %>%
  st_as_sfc() -> longbbox

# 读取长版的中国省级地图
read_sf("chinaprov2021long/chinaprov2021long.shp") %>%
  filter(!is.na(省代码)) -> provmap

# 线条
read_sf("chinaprov2021long/chinaprov2021long_line.shp") %>%
  filter(class %in% c("九段线", "海岸线")) %>%
  select(class) -> provlinemap

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

# 读取数据
readxl::read_xlsx("瞪羚、独角兽、创新型企业经纬度数据.xlsx") %>%
  filter(!is.na(经度)) -> pointdf

# 转换成 sf 对象
pointdf %>%
  st_as_sf(coords = c("经度", "纬度"), crs = 4326) %>%
  st_transform(mycrs) -> pointdfsf

pointdfsf

#> Simple feature collection with 37688 features and 8 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -1489994 ymin: 2072412 xmax: 1962404 ymax: 5261103
#> 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,688 × 9
#>    公司编号                      公司名称 省    省代码 市    市代码 县    县代码
#>  * <chr>                         <chr>    <chr>  <dbl> <chr>  <dbl> <chr>  <dbl>
#>  1 00073fd4d5204e1db2c517acb9b2… 北京领…  北京… 110000 北京… 110000 东城… 110101
#>  2 0007fc40037c446c8c73e834f1c1… 北京华…  北京… 110000 北京… 110000 密云… 110118
#>  3 000e2dda3c0f45c8b99dfb515a14… 广州穗…  广东… 440000 广州… 440100 黄埔… 440112
#>  4 0013abf442d84d9caa48eb286aa5… 武汉Ｔ…  湖北… 420000 武汉… 420100 江夏… 420115
#>  5 001434ab5e154531abdae6f3a6cb… 北京启…  北京… 110000 北京… 110000 昌平… 110114
#>  6 00155fb6ee9544c7af936bcaf1c2… 中通铁…  北京… 110000 北京… 110000 西城… 110102
#>  7 0015cb8bb70e4b60b77f1d4de5d4… 广东国…  广东… 440000 广州… 440100 天河… 440106
#>  8 0017ac5d15584cfca9a48e012e24… 思科系…  北京… 110000 北京… 110000 海淀… 110108
#>  9 00190b6fdb964662b7211723bbc2… 合肥英…  安徽… 340000 合肥… 340100 蜀山… 340104
#> 10 001f5df0e4ad43e5b47bed0a622b… 北京好…  北京… 110000 北京… 110000 昌平… 110114
#> # ℹ 37,678 more rows
#> # ℹ 1 more variable: geometry <POINT [m]>

# 汇总统计每个城市的
pointdfsf %>%
  st_drop_geometry() %>%
  count(省, 省代码) -> provdf

# 分组
provmap %>%
  left_join(provdf) %>%
  mutate(n = if_else(is.na(n), 0, n)) -> provdf2

provdf2$n %>%
  quantile(probs = 1:8/8, digits = 1) %>%
  as.integer() %>%
  unique() -> cutlist

cutlist

#> [1]     8    28    58   159   513   870  2446 16712

provdf2 %>%
  mutate(group = cut(n, breaks = c(0, unique(cutlist)),
                     include.lowest = T,
                     labels = c("<= 8", "8～28", "28～58", "58～159",
                                "159～513", "513～870", "870～2459", "> 2450"))) -> provdf3

provdf3 %>%
  ggplot() +
  geom_sf(data = longbbox, fill = "#BBD1EB") +
  geom_sf(aes(fill = group), color = "gray", linewidth = 0.01) +
  geom_sf(data = china_neighboring_long, fill = "#EDEDED") +
  stat_sf_coordinates(data = china_neighboring_long,
                      geom = "text", color = "gray",
                      aes(label = country_cn), family = cnfont,
                      size = 3) +
  geom_sf(data = provlinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  scale_color_manual(
    values = c("九段线" = "#A29AC4",
               "海岸线" = "#0055AA")
  ) +
  scale_linewidth_manual(
    values = c("九段线" = 0.6,
               "海岸线" = 0.3)
  ) +
  scico::scale_fill_scico_d(palette = "acton",
                            name = "公司数量",
                            end = 0.95) +
  new_scale_color() +
  geom_sf(data = pointdfsf, aes(color = 省),
          size = 0.1, show.legend = F) +
  scale_color_manual(values = paletteer::paletteer_d("ggsci::default_igv", n = 32)) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  labs(title = "瞪羚、独角兽、创新型企业地理分布",
       subtitle = "数据爬取&绘制：微信公众号 RStata",
       caption = "数据来源：瞪羚云网站") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_y = unit(0.1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) +
  theme(legend.position = c(0.2, 0.25)) -> p4

ggsave("pic4.png", width = 8, height = 9, device = png)

更多关于地图绘制的内容可以学习下平台上的系列课程「使用 R 语言进行地理计算」。

点击这里跳转到 RStata 短书平台获取附件：使用 R 语言绘制带邻国地区的中国地图