2.1.1 maps

分析和展示地理信息数据是一项常规任务，约 30 年前，Richard A. Becker 等为 S 语言引入地理可视化，特别是本文介绍的地区分布图，以及简单的区域面积和区域中心的计算能力[3, 8]，而后到了 2003年，Ray Brownrigg、Thomas P Minka 和 Alex Deckmyn 等将其引入 R 语言社区并持续维护至今[9]。除了最基础的 maps 包，还有坐标投影 mapproj 包[10]，以及提供更多地图数据的 mapdata 包[11]。那个时候，因缺乏一些基础工具，在地图数据获取、空间数据操作和可视化方面都不太容易，仅用这些能做到下图2.2已属不易。

# 加载数据
data(USCancerRates, package = "latticeExtra")
library(maps)
# 郡地图数据
us_county <- map("county", plot = FALSE, fill = TRUE, projection = "polyconic")
# 奈何 maps 内置的地图数据不全，仅保留部分观察数据
USCancerRates2 <- subset(
  x = USCancerRates,
  subset = rownames(USCancerRates) %in% us_county$names
)
# 调色板
colors <- viridisLite::plasma(13)
# 图例文本
leg.txt <- mapply(paste, 50*0:12, 50*1:13, collapse = " ", sep = "~")
# 癌症死亡率划分区间
USCancerRates2$colorBuckets_male <- as.numeric(cut(USCancerRates2$rate.male, 50*0:13))
# 根据城镇名称给地图上每个区域的癌症死亡率匹配颜色
colorsmatched_male <- USCancerRates2$colorBuckets_male[match(
  us_county$names, rownames(USCancerRates2)
)]
# 对女性癌症死亡率类似操作
USCancerRates2$colorBuckets_female <- as.numeric(cut(USCancerRates2$rate.female, 50*0:13))
colorsmatched_female <- USCancerRates2$colorBuckets_female[match(
  us_county$names, rownames(USCancerRates2)
)]
# 绘图
par(mar = c(0, 0, 3, 0), mfrow = c(2, 1))
# 添加地图背景
map("county",
  col = "grey80", fill = TRUE, resolution = 0,
  lty = 0, projection = "polyconic"
)
# 绘制区县地图
map("county",
  col = colors[colorsmatched_male], fill = TRUE, resolution = 0,
  lty = 0, projection = "polyconic", add = TRUE
)
# 添加州边界线
map("state",
  col = "white", fill = FALSE, add = TRUE,
  lty = 1, lwd = 0.2, projection = "polyconic"
)
map.scale()
# 添加图标题
title("1999-2003 年美国各个郡的年平均癌症死亡率（单位：十万分之一）", line = 2)
mtext(text = "男性", side = 3, adj = 0.5)
par(mar = c(1, 0, 2, 0))
map("county",
  col = "grey80", fill = TRUE, resolution = 0,
  lty = 0, projection = "polyconic"
)
map("county",
  col = colors[colorsmatched_female], fill = TRUE, resolution = 0,
  lty = 0, projection = "polyconic", add = TRUE
)
# 添加州边界线
map("state",
  col = "white", fill = FALSE, add = TRUE,
  lty = 1, lwd = 0.2, projection = "polyconic"
)
map.scale()
mtext(text = "女性", side = 3, adj = 0.5)
mtext(text = "数据源：美国国家癌症研究所", side = 1, adj = 0)
# 添加图例
legend("bottomright",
  legend = c(leg.txt, "NA"), title = "死亡率", box.col = NA, border = NA,
  horiz = FALSE, fill = c(colors, "grey80"), cex = 0.85, xjust = 0.5
)

图 2.2: 1999-2003 年美国各个郡的年平均癌症死亡率分布

maps 包[3]内置的美国郡级地图数据欠缺阿拉斯加、夏威夷和波多黎各。数据集 USCancerRates 没有夏威夷、波多黎各各个郡的数据，阿拉斯加的部分郡有数据，如图2.2所示，不少郡没有收集到癌症死亡率数据，以灰色表示。值得一提的是，关于死亡率分级，不同的分法会带给人不同的印象甚至是错觉，此处是可以有操作空间的，特定的死亡率分割方式可以让男女死亡率的空间分布看起来差异不大或很大[12]，详细介绍见 Marynia Kolak 和 Susan Paykin 在 2021 R/Medicine 大会上的视频和网页材料。

在数据操作方面，麻烦的是建立数据指标和各个地理区域的映射，上面反复用到函数 match()，Base R 有很多这样短小精悍的函数。下面简单介绍一下以助理解，函数 match() 返回一个向量，向量的长度与参数 x 一致，向量的元素是整型的，表示参数 x 中的元素出现在参数 table 中的位置，下面是三个小示例：

match(x = c("A", "B"), table = c("A"))
# [1]  1 NA
match(x = c("A", "B"), table = c("C", "A"))
# [1]  2 NA
match(x = c("A", "B"), table = c("C", "A", "D"))
# [1]  2 NA

另外，maps 包内置的地图数据制作起来比较复杂，也很长时间没有更新了，笔者不推荐读者再从零开始构建 map 数据对象。sp 包[13]发布后，maps 包支持将 Spatial 数据对象转为 map 数据对象，这相当于引入了 sp 包及其生态在空间数据获取和数据操作方面的能力[14]。目前，美国国家人口统计局已提供了历年的州、郡、普查级多个比例尺的行政区划地图数据，配合 sf 包[6]，可以很好地解决地理可视化的背景底图问题。一些以前看起来很难的问题，随着时间变迁，技术革新，已经解决了。

2.1.2 latticeExtra

maps 包是基于 Base R 绘图系统的，以线和多边形为基础，辅以颜色填充，还没有任何分层绘图的概念，Deepayan Sarkar 将其引入新的 Trellis（栅格）图形系统，这就有了 latticeExtra 包 [2]，它是以 lattice 包[15]为基础的，特别适合多元数据可视化，不管数据是「宽格式」还是「长格式」都能轻松应对，既保留了 Base R 那种精细操作图形元素的能力，也引入图层面板的概念。如图2.3，在代码形式上将所有的操作都集于一身，没有严格的分层绘图，保留更大的灵活性的代价是提升了复杂性。绘图主要步骤如下：

1. 准备数据：准备想要展示的数据指标 x ，指标所在的数据集 data，指标所描述的地理区域 map，同时传递给函数 mapplot(x, data, map, ...) 完成图形主体部分。
2. 添加配色：了解死亡率的分布情况，设置合适的分段，选择合适的调色板。
3. 添加边界：函数 mapplot(x, data, map, panel, ...) 中的参数 panel 需要一个函数作为参数值，修改默认的面板函数 panel = panel.mapplot，添加各郡的边界线，再调函数 latticeExtra::layer 叠加各州的边界线，设置颜色深浅和边界线的粗线形成层次。
4. 添加标题：添加横轴、纵轴、图例、分面标题，图形副标题，适当调整位置。
5. 最后调整布局、刻度等细节。除了步骤1必须先做，后续步骤没有严格的顺序，图形中各部分是相对独立的，若遇覆盖情况则需注意顺序。

# 美国州边界数据
us_state <- map("state", plot = FALSE, fill = TRUE, projection = "polyconic")
library(lattice)
# 绘图
latticeExtra::mapplot(rownames(USCancerRates) ~ rate.female + rate.male,
  # 观测数据
  data = USCancerRates,
  # 修改默认的面板函数
  panel = function(x, y, map, ...) {
    # 对未收集到癌症死亡率数据的郡，添加郡边界线
    panel.lines(
      x = us_county$x, y = us_county$y,
      lty = 1, lwd = 0.2, col = "black"
    )
    latticeExtra::panel.mapplot(x, y, map, ...)
  },
  # 郡地图数据为背景
  map = us_county,
  # 填充地图用的调色板
  colramp = viridisLite::plasma,
  # 覆盖郡边界线
  border = NA,
  # 多图布局 2 行 1 列
  layout = c(1, 2),
  # 死亡率分桶数
  # cuts = 14, # 指定 breaks 后就不需要 cuts
  # 死亡率数据分割点
  breaks = 50 * 0:13,
  # 仅展示地图数据中包含的郡死亡率数据
  subset = rownames(USCancerRates) %in% us_county$names,
  # 取消坐标轴刻度
  scales = list(draw = F),
  # 修改分面展示文本
  strip = strip.custom(factor.levels = c("女性", "男性")),
  # 去掉横轴标签
  xlab = "",
  # 添加图例标题
  legend = list(top = list(
    fun = grid::textGrob("死亡率", y = 0.1, x = 1.03)
  )),
  # 常规图形参数列表 trellis.par.get('par.sub.text')
  par.settings = list(
    # 副标题放在左下角
    par.sub.text = list(
      font = 2,
      just = "left",
      x = grid::unit(5, "mm"),
      y = grid::unit(5, "mm")
    ),
    # 分面边界和背景色
    strip.border = list(col = "white"),
    strip.background = list(col = "white"),
    # 轴线设置白色以隐藏
    axis.line = list(col = "white")
  ),
  # 副标题
  sub = "数据源：美国国家癌症研究所",
  # 主标题
  main = "1999-2003 年美国各个郡的年平均癌症死亡率"
) +
  # 添加州边界
  latticeExtra::layer(panel.lines(
    x = us_state$x, y = us_state$y,
    lty = 1, lwd = 0.2, col = "white"
  ))

图 2.3: 1999-2003 年美国各个郡的年平均癌症死亡率分布

用 lattice 包绘图往往可以一个函数搞定，参数非常多，都放在一起，这和 Base R 的绘图函数类似，也同时提供全局图形参数控制，但似乎更加复杂。下面谈几个细节：

1. 只需传给参数 colramp 一个生成颜色值向量的函数即可更改调色板，比如 R 内置的 hcl.colors() 或 terrain.colors() 等，为保持全文配色风格一致，图中配色采用 viridisLite 包提供的 plasma 调色板。

2. 参考 Markus Gesmann 在 2015 年写的一篇文章[16]，设置 Lattice 图形参数 par.settings 对图中的标题做了细致的调节，比如副标题的文本 par.sub.text 所有可调整的细节有：

   trellis.par.get("par.sub.text")
   # $alpha
   # [1] 1
   # 
   # $cex
   # [1] 1
   # 
   # $col
   # [1] "#000000"
   # 
   # $font
   # [1] 2
   # 
   # $lineheight
   # [1] 1

   想必读者已看出其规律，以 R 语言的列表结构来传递各个层级的图形参数值。

3. 参考 SO 论坛帖子设置参数 strip 自定义了分面子图的标题文本，再在 par.settings 里对背景 strip.background 和边界 strip.border 微调，而类似的设置在 ggplot2 包的主题函数 theme() 里也有，在 R 控制台执行 formalArgs(ggplot2::theme) 可获得主题函数的参数列表。

4. 读者下一个疑惑可能是如何知道所有的图形控制参数，以及控制的精细程度，Deepayan Sarkar 在书里以图2.4归纳了，纵轴是图形参数，横轴是参数值名称列表[15]，按图索骥一定有所帮助。

   

   图 2.4: lattice 包常用图形参数一览

   参数控制的效果预览如图 2.5 所示，不难看出，lattice 包[15]可以提供精细化的图形调整，是一个非常成熟的绘图工具。

   

   图 2.5: 图形控制参数

5. 如前所述，lattice 和 ggplot2 同出一脉，既然 ggplot2 有图层的概念，lattice 自然也有，除了 提供丰富的内置图层（约100个，读者可在控制台运行 apropos('panel') 查看当前环境下可用的图层），也支持用户自定义图层，图2.6 展示叠加图层 panel.superpose 的效果，传递自定义的符号给 pch 参数，实际上修改了默认图形参数 "superpose.symbol"。

   p1 <- xyplot(Sepal.Length ~ Sepal.Width, data = iris, groups = Species)
   p2 <- xyplot(Sepal.Length ~ Sepal.Width, data = iris, groups = Species,
     pch = c("L", "M", "N"), panel = panel.superpose
   )
   print(p1, split = c(1, 1, 2, 1), more = TRUE)
   print(p2, split = c(2, 1, 2, 1), more = FALSE)

   

   图 2.6: panel 叠加图层：左图默认散点图层，右图符号图层

6. 此外，latticeExtra 扩展了图层功能，特别是 layer() 函数，可在主体绘图函数之后直接叠加新图层。如上图2.3，在 latticeExtra::mapplot() 完成主体绘图工作后，添加美国各州边界线，有助于识别郡位置。大家应该都有这样一种感觉，将一张只有中国国家边界的地图放在面前，你不一定能清晰地指出每一个省份的位置，但只要画上各个省的边界，你肯定能指出更多的省份位置，类似地，从省到市、乃至区县，边界给了我们很好的参照。

2.1.3 ggplot2

考虑到 maps 包内置的地图数据的缺陷，ggplot2 包的流行度，下面采用 ggplot2 包绘制分面地区分布图，相比于 latticeExtra 包，ggplot2 包更适合「长格式」的数据，因此，需要先重塑 USCancerRates 数据集，将郡死亡率数据和郡地理边界数据配对，根据死亡率的分布设置合适分段，最后恢复地图数据的原始顺序。数据操作过程中有两点强调：

1. 数据重塑：USCancerRates 数据集除了 state 和 county 列，剩余列是由三个原子指标按性别维度衍生出来的，分别是癌症死亡率及其置信区间的上、下限值。在「宽格式」转「长格式」过程中，要注意转化前后各个列名的对应关系。

2. 数据关联：先将地图数据放左边，观测数据放右边，以 LEFT JOIN （左关联） 的方式关联起来，接着将连续性的死亡率数据分段，最后调用绘图函数绘制地区分布图。

# 宽格式转长格式
us_cancer_rates <- reshape(
  data = USCancerRates,
  # 需要转行的列
  varying = c(
    "LCL95.male", "rate.male", "UCL95.male",
    "LCL95.female", "rate.female", "UCL95.female"
  ),
  times = c("男性", "女性"), # 构成新列 sex 的列值
  v.names = c("LCL95", "rate", "UCL95"), # 列转行 列值构成的新列，指定名称
  timevar = "sex", # 列转行 列名构成的新列，指定名称
  idvar = c("state", "county"), # 可识别郡的 ID
  # 原数据有 3041 行，性别字段只有两个取值，转长格式后有 2*3041 行
  new.row.names = 1:(2 * 3041),
  direction = "long"
)

# 数据重塑前
head(USCancerRates) # 宽格式
#                 rate.male LCL95.male UCL95.male rate.female LCL95.female
# alabama,pickens     363.7      311.1      423.2       151.0        123.6
# alabama,bullock     345.7      274.2      431.4       140.5        102.8
# alabama,russell     340.7      304.5      380.9       182.3        161.3
# alabama,barbour     335.9      288.9      389.1       185.3        157.2
# alabama,dallas      330.1      293.4      370.6       172.0        151.4
# alabama,greene      328.1      255.9      416.6       124.1         88.5
#                 UCL95.female   state         county
# alabama,pickens        183.6 Alabama Pickens County
# alabama,bullock        189.7 Alabama Bullock County
# alabama,russell        205.5 Alabama Russell County
# alabama,barbour        217.5 Alabama Barbour County
# alabama,dallas         195.0 Alabama  Dallas County
# alabama,greene         171.7 Alabama  Greene County

# 数据重塑后
head(us_cancer_rates, 12) # 长格式
#      state          county  sex LCL95  rate UCL95
# 1  Alabama  Pickens County 男性 311.1 363.7 423.2
# 2  Alabama  Bullock County 男性 274.2 345.7 431.4
# 3  Alabama  Russell County 男性 304.5 340.7 380.9
# 4  Alabama  Barbour County 男性 288.9 335.9 389.1
# 5  Alabama   Dallas County 男性 293.4 330.1 370.6
# 6  Alabama   Greene County 男性 255.9 328.1 416.6
# 7  Alabama    Perry County 男性 261.4 327.9 408.0
# 8  Alabama   Walker County 男性 299.1 327.4 358.2
# 9  Alabama    Macon County 男性 276.4 323.6 377.5
# 10 Alabama    Lamar County 男性 265.3 321.4 386.3
# 11 Alabama Lawrence County 男性 276.1 318.5 366.9
# 12 Alabama  Marengo County 男性 267.4 315.6 370.7

# 从 usmapdata 包获取地图数据
county_df <- usmapdata::us_map("counties")
# 给每个郡的癌症死亡率数据配上地图数据
dat <- merge(
  x = county_df, y = us_cancer_rates,
  by.x = c("full", "county"),
  by.y = c("state", "county"),
  all.y = T
)
# 准备州边界线数据
state_df <- usmapdata::us_map("states")
# 癌症死亡率分级
dat$rate_d <- cut(dat$rate, breaks = 50 * 0:13)
# 恢复地图数据顺序
dat <- dat[order(dat$order), ]

数据准备工作完成后，绘图过程主要有八个步骤：

1. 底图：绘制美国郡地图，填充灰色背景，以此为底图。
2. 映射：添加每个郡的癌症死亡率数据，填充颜色根据死亡率而定。同时以浅白色绘制郡的边界线，将线调细一些，与后面州的边界线形成层次。
3. 边界：添加美国各个州的边界线，帮助熟悉美国地图的读者从州到郡快速定位。
4. 配色：类似前文设置，采用 plasma 调色板，未采集到死亡率数据的郡填充灰色，和地图背景融为一体。
5. 布局：以函数 facet_wrap() 实现分面，分面标题放在图形上方，布局为两行一列。
6. 标题：添加整个图形的主标题、副标题和图例标题。
7. 主题：设置整个图形主题样式 theme_void() 以符合地图特色，将图形主标题位置居中。
8. 细节：打磨图形配色、长宽尺寸、字体大小等方面，直至满意。

library(ggplot2)
ggplot() +
  geom_polygon(
    data = county_df, aes(x, y, group = group),
    fill = "grey80"
  ) +
  geom_polygon(
    data = dat, aes(x, y, group = group, fill = rate_d),
    colour = alpha("white", 1 / 4), size = 0.1
  ) +
  geom_polygon(
    data = state_df, aes(x, y, group = group),
    colour = "gray80", fill = NA, size = 0.15
  ) +
  scale_fill_viridis_d(option = "plasma", na.value = "grey80") +
  facet_wrap(~sex, ncol = 1, strip.position = "top") +
  labs(
    fill = "死亡率", title = "1999-2003 年美国各个郡的年平均癌症死亡率",
    caption = "数据源：美国国家癌症研究所"
  ) +
  theme_void(base_size = 13) +
  theme(plot.title = element_text(hjust = 0.5))

图 2.7: 1999-2003 年美国各个郡的年平均癌症死亡率分布

尽管已经调整了上图2.7的长宽比例，可一旦和非 ggplot2 绘制的图形对比，就可以看出明显变形了。

2.1.4 tmap

同 sf 包[6]的绘图函数 plot() 一样，tmap [5] 也是基于 Base R 图形系统，但使用语法更加贴合 ggplot2 包[4]，对空间数据可视化有更多专业支持，比如基于比例的符号/气泡图[17] ¹和基于二维核密度估计的热力图[5]。

tmap 支持 sp 包提供的 Spatial 类数据对象，也支持 sf 包提供的 Simple Features 类数据对象，后者是新一代更好的工具，因此接下来的示例都将基于 sf 包。首先从美国人口调查局下载州和郡级别的多边形边界数据，再类似 usmapdata 包对阿拉斯加、夏威夷和波多黎各做适当的调整，借助 tidycensus 包[18]让这个过程变得容易，调整地图数据的代码与源文档放在一起，读者可自取。

# 准备地图数据
library(sf)
us_state_map <- readRDS(file = "data/us_state_map.rds")
us_county_map <- readRDS(file = "data/us_county_map.rds")
# 重命名 fips_codes 的 state 列为 state_abbr
fips_codes <- within(tidycensus::fips_codes, {state_abbr = state})
fips_codes <- subset(x = fips_codes, select = setdiff(colnames(fips_codes), "state"))
# 癌症死亡率数据关联州、郡代码 FIPS
us_cancer_rates <- merge(
  x = us_cancer_rates,
  y = fips_codes,
  by.x = c("state", "county"),
  by.y = c("state_name", "county"),
  all.x = TRUE
)
# 合并地图数据和观测数据
us_county_cancer2 <- merge(
  x = us_county_map, y = us_cancer_rates,
  by.x = c("STATEFP", "COUNTYFP"),
  by.y = c("state_code", "county_code"), all.x = TRUE
)

准备好数据后，绘图过程比较类似之前使用 ggplot2 绘图，所不同的是，添加了指北针、比例尺等地区分布图特有的内容。

# 绘图
tmap::tm_shape(us_county_cancer2) +
  tmap::tm_polygons(border.col = "gray") +
  tmap::tm_shape(us_county_cancer2) +
  tmap::tm_polygons(
    col = "rate",
    palette = "plasma", border.alpha = 0.2,
    legend.reverse = TRUE, title = "死亡率",
    colorNA = "gray",
    textNA = "NA", border.col = "gray",
    breaks = 50 * 0:13
  ) +
  tmap::tm_facets(by = "sex", drop.NA.facets = T, ncol = 1) +
  tmap::tm_shape(us_state_map) +
  tmap::tm_polygons(col = NA, border.col = "gray", alpha = 0) +
  tmap::tm_compass(position = c("right", "bottom")) +
  tmap::tm_scale_bar(position = c("right", "bottom")) +
  tmap::tm_credits(text = "数据源：美国国家癌症研究所", position = "left") +
  tmap::tm_layout(
    legend.outside = TRUE,
    legend.outside.position = "right",
    legend.outside.size = 0.15,
    legend.format = list(text.separator = "~"),
    outer.margins = 0,
    asp = 0
  )

图 2.8: 1999-2003 年美国各个郡的年平均癌症死亡率分布

2.1.5 sf

在上一节已经准备好了美国各个郡的边界数据，接下来合并各个郡的癌症死亡率数据 USCancerRates，调用 sf 包[6]的绘图函数 plot()。除了调用函数 plot()的几行代码，绘制图2.9的代码和前面调用 maps 包绘图大体是类似的，但灵活性高得多，因 sf 包的强大，支持大量的空间数据存储格式，不受限于 maps 包内置的不完整地图数据，也不因给郡区域上色添加额外的数据匹配操作。

# 癌症死亡率数据关联州、郡代码
USCancerRates2 <- merge(
  x = USCancerRates,
  y = fips_codes,
  by.x = c("state", "county"),
  by.y = c("state_name", "county"),
  all.x = TRUE
)
# 空间数据合并观测数据
us_county_cancer <- merge(
  x = us_county_map, y = USCancerRates2,
  by.x = c("STATEFP", "COUNTYFP"),
  by.y = c("state_code", "county_code"), 
  all.x = TRUE
)
par(mar = c(0, 0, 4, 0), mfrow = c(2, 1))
# 添加郡地图
plot(st_geometry(us_county_map),
  reset = FALSE, border = NA, col = "grey80", main = ""
)
# 添加癌症死亡率数据
plot(us_county_cancer["rate.male"],
  pal = viridisLite::plasma, reset = FALSE, key.pos = NULL,
  breaks = 50 * 0:13,
  border = NA, lwd = 0.25, add = TRUE
)
# 美国各个州边界
plot(st_geometry(us_state_map), border = "gray80", lwd = 0.25, add = TRUE)
# 添加主标题
title(main = "1999-2003 年美国各个郡的年平均癌症死亡率", line = 2)
# 添加分面标题
mtext(text = "男性", side = 3, adj = 0.5)
# 调整边空
par(mar = c(1, 0, 2, 0))
plot(st_geometry(us_county_map),
  reset = FALSE, border = NA, col = "grey80", main = ""
)
plot(us_county_cancer["rate.female"],
  pal = viridisLite::plasma, reset = FALSE, key.pos = NULL,
  breaks = 50 * 0:13,
  border = NA, lwd = 0.25, add = TRUE
)
plot(st_geometry(us_state_map), border = "gray80", lwd = 0.25, add = TRUE)
mtext(text = "女性", side = 3, adj = 0.5)
mtext(text = "数据源：美国国家癌症研究所", side = 1, adj = 0)
# 调色板
colors <- viridisLite::plasma(13)
# 图例文本
leg.txt <- mapply(paste, 50 * 0:12, 50 * 1:13, collapse = " ", sep = "~")
# 添加图例
legend("bottomright",
  legend = c(leg.txt, "NA"), title = "死亡率", box.col = NA, border = NA,
  horiz = FALSE, fill = c(colors, "grey80"), cex = 0.85, xjust = 0.5
)

图 2.9: 1999-2003 年美国各个郡的年平均癌症死亡率分布

2.1.6 ggplot2 + sf

ggplot2 截止当前最新版本，ggplot2 在空间数据可视化都是非常不专业的，最糟糕的一个点是它会将地图弄变形了。 尽管ggspatial [19] 在 ggplot2 的基础上补充了一些地图特有的元素，如比例尺、指北针，但并没有解决地图变形的核心问题。正如 Roger Bivand 所言，还有其它可能难以预料的问题：

For visualisation, please avoid using ggplot2 unless you use this package often (daily). For thematic mapping, tmap and mapsf are to be preferred, because they are written for making maps. Do not simplify/generalise coastlines unless you really need to do so. You can use tmap and mapview to view thematic maps interactively - as you zoom in, you see artefacts created by line generalization.

因此，尽管 ggplot2 包已经流行开来，并且在很多方面取得成绩，但严格来说，不推荐使用 ggplot2 包来绘制任何和地图相关的图形，除非很清楚研究区域的情况，即在合适的地理区域采用合适的投影、合适的工具绘制准确的地理图形。ggplot2 在 2018 年发布 3.0.0 版本，开始借助 sf 包支持 Simple Features 数据对象的绘图，这就一定程度上缓解了地图变形的问题。

us_county_cancer2$rate_d <- cut(us_county_cancer2$rate, breaks = 50 * 0:13)
ggplot() +
  geom_sf(
    data = us_county_map,
    fill = I("grey80"), colour = NA
  ) +
  geom_sf(
    data = us_county_cancer2[!is.na(us_county_cancer2$sex), ],
    aes(fill = rate_d), colour = alpha("white", 1 / 4), size = 0.1
  ) +
  geom_sf(
    data = us_state_map,
    colour = alpha("gray80", 1 / 4), fill = NA, size = 0.15
  ) +
  scale_fill_viridis_d(option = "plasma", na.value = "grey80") +
  coord_sf(crs = st_crs("ESRI:102003")) +
  facet_wrap(~sex, ncol = 1) +
  labs(
    fill = "死亡率", title = "1999-2003 年美国各个郡的年平均癌症死亡率",
    caption = "数据源：美国国家癌症研究所"
  ) +
  theme_void(base_size = 13) +
  theme(plot.title = element_text(hjust = 0.5))

图 2.10: 1999-2003 年美国各个郡的年平均癌症死亡率分布

2.1.7 mapsf

如 Roger Bivand 所推荐，下面介绍 mapsf 包 [7]，它在 sf 包的基础上添加更多地理可视化的功能，比如指北针、比例尺、三维特效等，同时也支持比例符号图、地区分布图及其混合地图。mapsf 绘图思路和 sf 是一致的，只是封装了一些更加便捷的绘图函数，如图2.11所示，大体分三步：其一准备主题样式，配色、布局等；其二表达地理数据，选择合适的图形，建立数据到地图的映射；其三补充一些说明，如注释、比例尺、标题等。

library(mapsf)
par(mfrow = c(2, 1))
# 设置主题
mf_theme("default", inner = TRUE, tab = TRUE, pos = "center")
# 绘制地图
mf_map(
  x = us_county_cancer,
  var = "rate.male",
  type = "choro",
  breaks = 50 * 0:13,
  pal = "Plasma",
  border = "gray",
  lwd = 0.5,
  leg_val_rnd = 0,
  leg_pos = "left",
  leg_no_data = "NA",
  leg_title = "死亡率\n(每10万人)"
)
# 指北箭头
mf_arrow(pos = "topleft")
# 比例尺
mf_scale(pos = "bottomleft")
# 标题
mf_title(txt = "男性", bg = "#f7f7f7", fg = "black", cex = 1)
# 出处
mf_credits(txt = "数据源：美国国家癌症研究所", pos = "bottomright", cex = 1)

mf_map(
  x = us_county_cancer,
  var = "rate.female",
  type = "choro",
  breaks = 50 * 0:13,
  pal = "Plasma",
  border = "gray",
  lwd = 0.5,
  leg_val_rnd = 0,
  leg_pos = "left",
  leg_no_data = "NA",
  leg_title = "死亡率\n(每10万人)"
)
mf_arrow(pos = "topleft")
mf_scale(pos = "bottomleft")
mf_title(txt = "女性", bg = "#f7f7f7", fg = "black", cex = 1)
mf_credits(txt = "数据源：美国国家癌症研究所", pos = "bottomright", cex = 1)

图 2.11: 1999-2003 年美国各个郡的年平均癌症死亡率分布

3.1.1 数据获取

受上个癌症死亡率地区分布图的启发，除了郡级别，还有社区普查级别的地区分布图，以描述更加精细的空间区域数据分布，同时，这次直接从美国人口调查局获取数据。数据来自 2015-2019 年度美国社区调查（American Community Survey，简称 ACS），每年约350万家庭的抽样，每年更新一次。本节使用的数据有：

1. 北卡罗来纳州社区级家庭年收入，统计每个家庭过去12个月的收入，这相当于家庭年收入，再取社区内所有家庭年收入的中位数，众所周知，中位数比平均数更加稳定，更有代表性。家庭年收入已根据 2019 年的美国消费价格指数 CPI 调整，这样，各个年份的数据是有可比性的。
2. 北卡罗来纳州社区级白人居民的百分比，美国除了白人，还有其它有色人种，且占比不小，分布广泛。
3. 2019 年北卡罗来纳州郡级、社区级比例尺 1:500,000 多边形边界地图数据，2019 年美国州级、郡级比例尺 1:20,000,000 多边形边界地图数据。

美国人口调查局（US Census Bureau）提供数据 API 作为公共服务，Kyle Walker 为此开发了 R 接口 tidycensus 包[18]，原始数据可借助 tidycensus 包从美国官方网站下载。下载数据之前，需要先注册一个访问令牌，保存到 R 环境变量 CENSUS_API_KEY，可在 .Renviron 文件中存储环境变量 CENSUS_API_KEY，以便后续调用。在网络一切都好的情况下，下载所需数据指标和北卡郡级、社区级地图数据非常简单，如下：

library(tidycensus)
library(sf)
# 郡级数据
nc_income_race_county <- get_acs(
  state = "NC", # 北卡州
  geography = "county", # 郡级
  # 三个数据指标：家庭年收入，白人数量，总人数
  variables = c("B19013_001", "B02001_001", "B02001_002"),
  geometry = TRUE, # 下载北卡郡级地图数据
  output = "wide", # 输出宽格式数据
  year = 2019, # 2015-2019 年度
  survey = "acs5", # 5 年
  moe_level = 90 # 置信水平 90%
)
# 社区级数据
nc_income_race_tract <- get_acs(
  state = "NC",
  geography = "tract", # 社区级
  variables = c("B19013_001", "B02001_001", "B02001_002"),
  geometry = TRUE, output = "wide",
  year = 2019, survey = "acs5",
  moe_level = 90
)
# 北卡州郡级地图数据
nc_county_map <- tigris::counties(state = 37, cb = T, year = 2019, class = "sf")
# 北卡州社区级地图数据
nc_tract_map <- tigris::tracts(state = 37, cb = T, year = 2019, class = "sf")

数据下载后，保存为四个文件：

1. nc_income_race_county.rds 北卡州郡级家庭年收入和白人数量。
2. nc_income_race_tract.rds 北卡州社区级家庭年收入和白人数量。
3. nc_county_map.rds 北卡州郡级地图数据。
4. nc_tract_map.rds 北卡州社区级地图数据。

下面补充几个细节：

1. 不同年份 ACS 调查的郡、社区有变化，美国也会定期更新地图数据，因此，北卡郡级、社区级地图数据的年份和 ACS 调查数据的发布年份要一致。

2. R 软件内置了数据集 state.name 和 state.abb，分别对应美国各个州的名称及其简称，tidycensus 包内置的数据集 fips_codes 还提供各州、郡的 FIPS 代码，用于下载所需州郡的数据。

3. 读者若对其它数据指标感兴趣，可通过如下命令查看：

   acs_vars <- load_variables(year = 2019, dataset = "acs5")
   View(acs_vars)

   一共 27040 个指标，其中很多是衍生指标，比如上面用到的白人数量 B02001_001，人口数按种族分组统计，自然还可以有按性别、年龄分组的衍生指标。

4. Kyle Walker 的著作《Analyzing US Census Data: Methods, Maps, and Models in R》[22] 提供了更多关于美国社区调查数据及其使用的介绍。

5. 若网络不好，查看函数 get_acs() 源码获取数据下载地址，打开浏览器从网页获取，嫌麻烦的话，直接用本文随附的数据吧。

3.1.2 数据探索

The phenomenon external to an area of interest affects what goes on inside.

— Waldo Tobler [23]

在展示北卡州的数据之前，先看下美国各州、郡的收入分布概况，如图3.1 和图3.2 分别是美国州级、郡级的收入分布图，结合这两张图，不难看出除了少量的西岸加州一带城市和东岸波士顿、纽约、华盛顿一带城市，收入的空间分布较为平衡。北卡州的收入水平属于中等偏上，也算典型，境内分布的罗利、达勒姆、教堂山三座城市被称为东岸硅谷。北卡罗来纳大学有16个分校，俨然一个超大规模的大学城，教堂山分校的统计系是北美最早建立的统计系之一。

图 3.1: 2015-2019 年度美国各个州家庭年收入

图 3.2: 2015-2019 年度美国各个郡家庭年收入

接下来，看下2015-2019年度美国北卡罗来纳州社区级家庭年收入的空间分布，如图3.3所示，图中空白处表示该区域的样本量太少或调查区域时间不连续，导致统计模型的估计不可用。结合图3.3，从空间上看，高收入和低收入群体有明显的圈层，在图中高亮的地区（也可以叫城镇），低收入社区穿插在高收入社区中，且有清晰的穿插线。

# 读取数据
nc_income_race_county <- readRDS(file = "data/nc_income_race_county.rds")
nc_income_race_tract <- readRDS(file = "data/nc_income_race_tract.rds")
# 计算家庭年收入和白人占比数据
nc_income_race_tract <- within(nc_income_race_tract, {
  pop <- B02001_001E
  pctWhite <- B02001_002E / B02001_001E
  medInc <- B19013_001E
})

图 3.3: 北卡州社区级家庭年收入的空间分布

如图 3.4 所示，2015-2019 年度 ACS 调查，北卡州一共划分 2195 个社区，除了 36 个没有收集到收入数据的社区外，北卡州超过一半的社区家庭年收入超过 50000 美元，中位数是 51439 美元，图中三条红线分别是 25%，50% 和 75% 分位点。

图 3.4: 北卡州社区级家庭年收入的分布

如图 3.5 所示，全州大部分社区人口数 10000 以内，极少的社区人口在 15000 以上。图中紫黑色表示人口数为 0，因社区人口太少或区域在时间上不连续导致模型无法给出估计。

图 3.5: 北卡州社区级人口的空间分布

如图3.6所示，北卡州超过一半的社区人口数不足 5000，中位数是 4389，图中三条红线分别是 25%，50% 和 75% 分位点。ACS 调查会根据人口密度划分采样区，人口密度越高的城市采样区越多，反之，采样区越少。若某些区域的人口密度不断增大，还要继续拆分，一旦拆分，两次调查的空间划分就不一样了，就会出现局部地区不可比。整体上，每个采样区的人口数应在 1200-8000，最好在 4000 人，以便抽样调查的样本比较均衡，方便后续模型估计，更多内容详见帮助文档?tigris::tracts。

图 3.6: 北卡州社区级人口的分布

如图3.7所示，白人集中分布在东岸城市（北卡州的东部边缘，也是大西洋西岸）和北卡州的西部，整体上，白人占全州比例非常高，再结合图3.5，粗估也在 70%-80% 之间。

图 3.7: 北卡州社区级白人占比的空间分布

如图3.8所示，将家庭年收入和白人占比两个数据的分布联合起来看，融合了上图3.3和图3.7的数据信息。图中以白色块表示缺失数据。

图 3.8: 2015-2019 年度北卡州社区级家庭年收入和白人占比的空间相关性

3.1.3 数据解读

继续社区级普查数据，先通过简单的线性模型粗估一下，家庭年收入和白人占比之间的关系：

\[ \mathrm{medInc} = \alpha + \mathrm{pctWhite} \cdot \beta \]

fit = lm(data = nc_income_race_tract, medInc ~ pctWhite)
summary(fit)
# 
# Call:
# lm(formula = medInc ~ pctWhite, data = nc_income_race_tract)
# 
# Residuals:
#    Min     1Q Median     3Q    Max 
# -43312 -15651  -6136   8213 174442 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)    27808       1652    16.8   <2e-16 ***
# pctWhite       43800       2271    19.3   <2e-16 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Residual standard error: 24900 on 2157 degrees of freedom
#   (36 observations deleted due to missingness)
# Multiple R-squared:  0.147,   Adjusted R-squared:  0.147 
# F-statistic:  372 on 1 and 2157 DF,  p-value: <2e-16

根据调整的 \(R^2 = 0.147\)，各社区家庭年收入和白人占比相关性比较低，结合图3.9也不难看出稳定正向的关系，平均来说，社区的白人占比增加一个百分点，该地区家庭平均年收入增加438美元。

图 3.9: 家庭年收入和白人占比线性相关性

Everything is related to everything else, but near things are more related than distant things.

— Waldo Tobler [24]

接下来，回归到本节主旨：空间相关性。先来看下社区级家庭年收入的空间相关性，这可以用局部莫兰 I 指数刻画，可调 rgeoda包[25]计算，rgeoda 包站在libgeoda 和 GeoDa的肩膀上提供许多空间数据分析算法和模型。如图3.10所示，灰白色区域表示统计上空间聚集现象不显著，红色区域表示高收入显著的空间聚集现象，而大片的蓝色区域表示低收入显著的空间聚集现象。

图 3.10: 北卡社区级家庭年收入的空间相关性

相比于图3.3，此图描述的空间结构清晰得多了，整体上，高收入的聚在一起，低收入的也聚在一起。高收入地区，其内圈是低收入，外圈也是低收入，高收入低收入一圈包一圈。结合图3.1，图3.2和图3.3，从州、郡、社区层层下钻和分析就可以看出收入不平衡的区域，以及精细的空间聚类模式。

图 3.11: 北卡社区级白人占比的空间相关性

类似地，图3.11更清晰地描述白人在北卡的分布区域，结合收入分布图3.3和图3.10来看，收入越高的地区白人占比越高，由此可见，引入地区分布图帮助我们更直观地了解了白人占比和家庭收入的关系。