R语言实战（2）-制作地区分布图及其动画

跟着xiangyun练习R语言

R
ggplot2
Author

Lee

Published

November 16, 2022

因为在工作工程中，经常要接触到一些地图的数据，并需要对其进行可视化操作，所以 我希望学习ggplot2包制作地图的详细过程。原文为xiangyun的博文。xiangyun主要使用的是baseR的技术，我则尝试使用tidy风格对文章代码进行复现。

1 地图制作-地区分布图

这是最常用也最常见的一个地图的种类。

1.1 数据准备

使用maps包中内置的vote.repub数据集。

library(tidyverse)
library(maps)
library(sf)
data(votes.repub)
str(votes.repub)
 num [1:50, 1:31] NA NA NA NA 18.8 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:50] "Alabama" "Alaska" "Arizona" "Arkansas" ...
  ..$ : chr [1:31] "1856" "1860" "1864" "1868" ...

数据集为51行31列的矩阵，行代表各州，列代表大选年份。进一步查看votes.repub的维度信息：

attributes(votes.repub)
$dim
[1] 50 31

$dimnames
$dimnames[[1]]
 [1] "Alabama"        "Alaska"         "Arizona"        "Arkansas"      
 [5] "California"     "Colorado"       "Connecticut"    "Delaware"      
 [9] "Florida"        "Georgia"        "Hawaii"         "Idaho"         
[13] "Illinois"       "Indiana"        "Iowa"           "Kansas"        
[17] "Kentucky"       "Louisiana"      "Maine"          "Maryland"      
[21] "Massachusetts"  "Michigan"       "Minnesota"      "Mississippi"   
[25] "Missouri"       "Montana"        "Nebraska"       "Nevada"        
[29] "New Hampshire"  "New Jersey"     "New Mexico"     "New York"      
[33] "North Carolina" "North Dakota"   "Ohio"           "Oklahoma"      
[37] "Oregon"         "Pennsylvania"   "Rhode Island"   "South Carolina"
[41] "South Dakota"   "Tennessee"      "Texas"          "Utah"          
[45] "Vermont"        "Virginia"       "Washington"     "West Virginia" 
[49] "Wisconsin"      "Wyoming"       

$dimnames[[2]]
 [1] "1856" "1860" "1864" "1868" "1872" "1876" "1880" "1884" "1888" "1892"
[11] "1896" "1900" "1904" "1908" "1912" "1916" "1920" "1924" "1928" "1932"
[21] "1936" "1940" "1944" "1948" "1952" "1956" "1960" "1964" "1968" "1972"
[31] "1976"

为方便后续操作，将原数据的矩阵格式更改为数据框格式。

rownames <- rownames(votes.repub)
voted_repub <- as_tibble(votes.repub,
  # 保留原数据的行名
  rownames = pkgconfig::get_config("tibble::rownames", NA)
)

再增加两列state和state_name列，为后续与地图数据集关联，从而将数据映射到提图上。

state列可与 ggplot2::map_data() 生成的地图数据合并，state_name 列可与 tigris 包获取的数据合并。
voted_repub <-
  voted_repub %>%
  mutate(
    state_names = rownames(voted_repub),
    region = str_to_lower(rownames(voted_repub))
  ) %>%
  arrange(state_names)
Note

需要什么格式的州名取决于所能获得的地图数据的形式，一般需要对地图数据有所了解才能知晓。

接下来为方便后续的操作(如ggplot2绘图)，就是数据重塑。

us_voted_repub <-
  voted_repub %>%
  pivot_longer(-c(region, state_names),
    names_to = "year",
    values_to = "repub_prop"
  )
us_voted_repub
# A tibble: 1,550 × 4
   state_names region  year  repub_prop
   <chr>       <chr>   <chr>      <dbl>
 1 Alabama     alabama 1856       NA   
 2 Alabama     alabama 1860       NA   
 3 Alabama     alabama 1864       NA   
 4 Alabama     alabama 1868       51.4 
 5 Alabama     alabama 1872       53.2 
 6 Alabama     alabama 1876       40.0 
 7 Alabama     alabama 1880       37.0 
 8 Alabama     alabama 1884       38.4 
 9 Alabama     alabama 1888       32.3 
10 Alabama     alabama 1892        3.95
# ℹ 1,540 more rows

可以看到，us_vote_repub数据集存在很多缺失值，这是由于美国大选有些历史数据没有记录，我们需要关注从何时起，大选数据基本全覆盖了。

us_voted_repub %>%
  group_by(year) %>%
  summarise(n = sum(is.na(repub_prop))) %>%
  ggplot(aes(x = year, y = n, group = 1)) +
  geom_line() +
  scale_x_discrete(breaks = seq(1856, 1976, by = 4)) +
  scale_y_continuous(
    expand = c(0, 1),
    limits = c(0, 30)
  ) +
  theme(
    axis.text.x = element_text(angle = 45)
  )

图 1: 历届大选未收集到得票率数据的州有多少

如 图 1 所示，大选在20世纪初就基本覆盖了全美国。

接下来，以某一年的数据为基础，绘制一张静态地区分布图：

us_voted_repub_1976 <- us_voted_repub %>%
  filter(year == 1976)

# 获取美国洲际地图数据
usa_map <- map_data(map = "state")
usa_map %>%
  slice_head(n = 6)
       long      lat group order  region subregion
1 -87.46201 30.38968     1     1 alabama      <NA>
2 -87.48493 30.37249     1     2 alabama      <NA>
3 -87.52503 30.37249     1     3 alabama      <NA>
4 -87.53076 30.33239     1     4 alabama      <NA>
5 -87.57087 30.32665     1     5 alabama      <NA>
6 -87.58806 30.32665     1     6 alabama      <NA>

usa_map中，每一对经度long和纬度lat对应州边界多边形的一个点，按一定顺序连接起来，就可以形成每个州的行政区划边界。注意，由于usa_map是州级地图数据，不包含各郡的边界，故而 subregion列都为缺失值。

下面我们需要将观测数据合并到usa_map数据中，注意，地图数据在左边，观测数据在右侧，因为右侧数据是不完整的，而绘制地图时最好保证地图数据的完整。

usa_voting_df <- usa_map %>%
  left_join(us_voted_repub_1976, by = "region")
usa_voting_df %>%
  slice_head(n = 6)
       long      lat group order  region subregion state_names year repub_prop
1 -87.46201 30.38968     1     1 alabama      <NA>     Alabama 1976      43.48
2 -87.48493 30.37249     1     2 alabama      <NA>     Alabama 1976      43.48
3 -87.52503 30.37249     1     3 alabama      <NA>     Alabama 1976      43.48
4 -87.53076 30.33239     1     4 alabama      <NA>     Alabama 1976      43.48
5 -87.57087 30.32665     1     5 alabama      <NA>     Alabama 1976      43.48
6 -87.58806 30.32665     1     6 alabama      <NA>     Alabama 1976      43.48

由于地图上的每个区域都有编号顺序，不能随意打乱，原文中使用baseR需要对数据进行排序。而使用tidy风格的代码，order列的顺序没有改变，故不需要进行重新排序的操作。

1.2 数据展示

与普通的图形不同，绘制地图需要注意底图和投影坐标，ggplot2提供了geom_map()和coord_map()两个函数实现对应配置。如 图 2，是根据maps包的北美州级地图数据绘制而成的州级边界地图。

ggplot(
  usa_voting_df,
  aes(x = long, y = lat, group = subregion)
) +
  geom_map(
    map = usa_map, aes(map_id = region),
    color = "grey80", fill = "grey30", linewidth = 0.3
  )

图 2: 北美州级地图

接着需要设置坐标系。maps包内置的 state 数据集，其经纬度坐标来自常见 WGS 84 参考系，即所谓的世界大地测量系统 World Geodetic System，于 1984 制定，全球定位系统 GPS 即采用此坐标参考系。

为了在二维的地图上更加形象的展示地球的球形特点，使用mapproj包，以北极为中心的方位角进行投影：纬度为同心圆弧线，经度则是等距径向线，详细描述可参看?mapproj::mapproject具体如下：

ggplot(usa_voting_df, aes(long, lat, group = subregion)) +
  geom_map(
    map = usa_map, aes(map_id = region),
    fill = "grey30", color = "grey80", linewidth = 0.3
  ) +
  coord_map(
    projection = "orthographic",
    orientation = c(
      latitude = 39,
      longitude = -98, ortation = 0
    )
  )

接着，利用多边形图层geom_polygon()添加得票数据，采用不同颜色对得票多少进行区分，并对图形进行一定美化。

library(maps)
ggplot(usa_voting_df, aes(long, lat, group = subregion)) +
  geom_map(
    map = usa_map, aes(map_id = region),
    fill = "grey30", color = "grey80", linewidth = 0.3
  ) +
  coord_map(
    projection = "orthographic",
    orientation = c(
      latitude = 39,
      longitude = -98, ortation = 0
    )
  ) +
  geom_polygon(
    aes(group = group, fill = repub_prop / 100),
    color = "black"
  ) +
  scale_fill_gradient2(
    low = "blue", mid = "white", high = "red",
    midpoint = 0.5, labels = scales::label_percent()
  ) +
  theme_minimal() +
  labs(
    title = "1976 年共和党在各州的得票率",
    x = "", y = "", fill = "得票率"
  ) +
  theme(
    plot.title = element_text(
      size = 16, face = "bold", color = "red3"
    ),
    legend.position = "bottom"
  )

图 3: 1976 年共和党在各州的得票率

图 3 中的地图已经初具雏形，但还有几个缺陷：

  • maps包的state数据缺少阿拉斯加和夏威夷两个州。
  • 图中缺少南北纬，东西经的方向指示。

我们还需进一步对图形进行优化：

  • 使用tigris包下载美国的完整地图。(由于网络原因无法顺利下载，这里直接使用现成的rds数据文件)
  • 使用sf包对地图进行绘制。
us_state_map <- read_rds("d:/Myblog/datas/us_state_map.rds")
us_state_map
Simple feature collection with 52 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -3112200 ymin: -1697728 xmax: 2258154 ymax: 1558935
Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic
First 10 features:
   STATEFP  STATENS    AFFGEOID GEOID STUSPS         NAME LSAD        ALAND
1       53 01779804 0400000US53    53     WA   Washington   00 172117954267
2       46 01785534 0400000US46    46     SD South Dakota   00 196346195316
3       39 01085497 0400000US39    39     OH         Ohio   00 105823701267
4       01 01779775 0400000US01    01     AL      Alabama   00 131174192284
5       05 00068085 0400000US05    05     AR     Arkansas   00 134776580080
6       35 00897535 0400000US35    35     NM   New Mexico   00 314197253999
7       48 01779801 0400000US48    48     TX        Texas   00 676668210823
8       06 01779778 0400000US06    06     CA   California   00 403660088482
9       21 01779786 0400000US21    21     KY     Kentucky   00 102282218059
10      13 01705317 0400000US13    13     GA      Georgia   00 149484604230
        AWATER                       geometry
1  12549727444 MULTIPOLYGON (((-2000238 15...
2   3383460688 MULTIPOLYGON (((-633753 865...
3  10274036690 MULTIPOLYGON (((1081987 544...
4   4593183334 MULTIPOLYGON (((708460.5 -5...
5   2956395922 MULTIPOLYGON (((122656.3 -1...
6    727781442 MULTIPOLYGON (((-1226428 -5...
7  18991880422 MULTIPOLYGON (((-998043.8 -...
8  20305454540 MULTIPOLYGON (((-2066285 -2...
9   2372611005 MULTIPOLYGON (((571924.5 -8...
10  4420380296 MULTIPOLYGON (((939223.1 -2...

同样的，我们将地图数据与观测数据进行连接。注意，此处使用的是NAME。

library(maps)
library(sf)
us_voted_repub_1976 <-
  us_voted_repub_1976 %>%
  rename(NAME = "state_names")

usa_voting_map_1976 <- us_state_map %>%
  left_join(us_voted_repub_1976, by = "NAME")

使用sf包，可以大大简化梯度的绘制过程。

library(sf)
ggplot() +
  geom_sf(
    data = usa_voting_map_1976,
    aes(fill = repub_prop / 100)
  )

图 4: 1976年共和党在各州的得票率

如 图 4 所示，所有的州都显示了，且把阿拉斯加、夏威夷和波多黎各适当移动到了全美地图的下方，方面查看。

那么下一步，我们就要将具体的得票率映射在地图上。

ggplot() +
  geom_sf(
    data = usa_voting_map_1976,
    aes(fill = repub_prop / 100)
  ) +
  scale_fill_gradient2(
    low = "blue", mid = "white", high = "red",
    midpoint = 0.5, labels = scales::label_percent()
  ) +
  theme_minimal() +
  labs(
    title = "1976 年共和党在各州的得票率",
    x = "", y = "", fill = "得票率"
  ) +
  theme(
    plot.title = element_text(size = 16, face = "bold", color = "red3"),
    legend.position = "bottom"
  )

图 5: 1976年共和党在各州的得票率-美化

2 动画制作

之前的章节仅展示了其中一个年份的数据分布情况，不便于对比。本部分介绍GIF和WEB两种动画的形式。

GIF图，参看这里。本部分我们主要关注WEB类型的动图，即可交互式的图形。我们要用到的是echarts4r包。

一篇写给初学者的入门文章参见《echarts4r: 从入门到应用》。但是这部分地图的动图上文并未提到。

library(echarts4r)
library(echarts4r.maps)
us_voted_repub %>%
  group_by(year) %>%
  e_charts(x = state_names, timeline = TRUE) %>%
  em_map(map = "USA") %>%
  e_map(serie = repub_prop, map = "USA", name = "得票率（%") %>%
  e_visual_map(serie = repub_prop) %>%
  e_tooltip()
图 6: 1856-1976 年美国历届大选中共和党在各州的得票率

在上 图 6 的基础上，还可以进一步对图像进行定制。

  • 视觉映射：修改视觉映射上下两端文本，使用红蓝对撞的调色板，置于左下方。有点体现共和党和民主党对垒的意思。

  • 主副标题：在图形上方添加动画标题、左下方说明数据来源。在动画场景下，用 e_timeline_serie() 实现，非动画场景下，可用 e_title() 实现。

  • 地图底图：echarts4r.maps 包内置的 GeoJSON 格式美国州级地图。

  • 坐标投影：相比于前面展示的动画，最终效果图 11 采用的坐标投影不一样，它是反映实际区域面积大小的。目前 echarts4r 还不支持其它投影，读者定制需要引入第三方 JavaScript 库 d3-geo 和 d3-array。

  • 数据映射：注意地图数据中州名 name 字段格式，阿拉斯加州是 Alaska 而不是 alaska，这个用于匹配观测数据中的 state_name 列。

  • 动画播放：设置自动播放动画，播放间隔为 1s。

us_voted_repub %>%
  group_by(year) %>%
  e_charts(x = state_names, timeline = TRUE) %>%
  em_map(map = "USA") %>%
  e_map(serie = repub_prop, map = "USA", name = "得票率（%)") %>%
  # 设置图例条
  e_visual_map(
    serie = repub_prop,
    left = 80,
    bottom = 40,
    text = c("高", "低"),
    inRange = list(
      color = RColorBrewer::brewer.pal(n = 11, name = "RdBu")
    )
  ) %>%
  e_tooltip() %>%
  # 设置时间轴条
  e_timeline_opts(
    autoPlay = TRUE,
    show = TRUE,
    bottom = 20,
    playInterval = 1000
  ) %>%
  e_timeline_serie(
    title = purrr::map(
      rep("1856-1976 年美国历届大选中共和党在各州的得票率", 31),
      function(title) {
        list(
          text = title,
          left = "center"
        )
      }
    ), index = 1
  ) %>%
  e_timeline_serie(
    title = purrr::map(
      rep("数据源：maps 包内置数据集 votes.repub", 31),
      function(subtitle) {
        list(
          subtext = subtitle,
          left = 80,
          bottom = 0,
          sublink = "https://cran.r-project.org/package=maps/"
        )
      }
    ), index = 2
  ) %>%
  e_timeline_serie(
    title = purrr::map(
      colnames(votes.repub),
      function(year) {
        list(
          text = year,
          right = 150,
          top = 50,
          textStyle = list(
            fontSize = 100,
            color = "rgb(170, 170, 170, 0.5)",
            fontWeight = "bolder"
          )
        )
      }
    ), index = 3
  )

3 小结

最后，总结一下制作动态地区分布图的步骤，主要分数据操作、地图绘制和动画制作三部分。

  1. 数据操作部分：这部分比较熟悉，不赘述了。

  2. 地图绘制部分：

  • 地图图层：从权威的渠道获取专业的地图数据，推荐使用 geom_sf() 而不是 geom_map()，前者更加流行和通用。
  • 投影图层：推荐使用与 geom_sf() 配套的 coord_sf()，根据地理区域的情况选择合适的投影。
  • 颜色图层：scale_fill_gradient2() 对撞型的调色板衬托美国大选中的两党之争。
  • 主题层：theme_minimal() 简洁美观即可。
  • 标题层：labs() 适当添加一些图例标题等。
  1. 动画制作部分：
  • GIF 动画：在制作地区分布图方面，笔者推荐 gganimate 和 tmap 包，后者在地理可视化方面有很多积累，文档丰富。二者的使用上和大家熟悉的 ggplot2 一脉相承，学习成本也比较低。
  • Web 动画：echarts4r 包的动画示例比较简单，定制一些内容需要多参考时间轴配置项手册，echarts4r 暂不支持自定义坐标投影。