6 地图

地图绘制是一项常见的可视化任务，需要用到专门的工具。主要解决两个方面的问题：使用一个数据源绘制地图，以及将另一个信息的数据添加到地图中。

6.1 矢量边界

使用geom_polygon()或许是绘制地图最简单的方法，它是为不同区域绘制边界。通过使用ggplot2::map_data()数据， 该地图包不是特别准确或最新的包，但它内置在 R 中，它是密歇根州县边界的数据集：

mi_counties <- map_data("county", "michigan") %>% 
  select(lon = long, lat, group, id = subregion)
head(mi_counties)
#>     lon  lat group     id
#> 1 -83.9 44.9     1 alcona
#> 2 -83.4 44.9     1 alcona
#> 3 -83.4 44.9     1 alcona
#> 4 -83.3 44.8     1 alcona
#> 5 -83.3 44.8     1 alcona
#> 6 -83.3 44.8     1 alcona

在这个数据集中有四个变量：lat和long指定边界点的纬度和经度（即边界点的坐标），id指一个区域的名称，group提供了一个区域内的连续区域的独特标识符（例如，如果一个区域由多个岛屿组成）。为了更好地了解数据包含的内容，我们可以使用mi_counties数据绘图geom_point()，如下所示。在此图中，数据框中的每一行都绘制为一个点，从而生成一个散点图，显示每个县的角点。为了将这个散点图变成地图，我们改为使用geom_polygon()，它将每个县绘制为一个不规则的多边形。

ggplot(mi_counties, aes(lon, lat)) + 
  geom_point(size = .25, show.legend = FALSE) +
  coord_quickmap()

ggplot(mi_counties, aes(lon, lat, group = group)) +
  geom_polygon(fill = "white", colour = "grey50") + 
  coord_quickmap()

image

image

在这两个图中，coord_quickmap()`来调整轴以确保经度和纬度以相同的比例呈现。coord中讨论了 ggplot2 中的坐标系，但正如我们将在下面看到的，地图数据通常需要更严格的方法。因此，ggplot2 提供geom_sf()和coord_sf()处理以简单的空间数据。

6.2 简单的特征图

上面概述的方法有一些局限性，尤其是简单的“经纬度”数据格式在实际映射中通常不使用。地图的矢量数据通常使用开放地理空间联盟产生的“简单特征”标准进行编码。Edzer Pebesma https://github.com/r-spatial/sf开发的[sf] (https://github.com/r-spatial/sf)包为处理这些数据提供了一个很好的工具，并且ggplot2中的geom_sf()和coord_sf()函数可以与sf包一起工作。

为了引入这些功能，我们依赖于Michael Sumner https://github.com/mdsumner/ozmaps/的 ozmaps 包，该包提供澳大利亚州边界、地方政府区域、选举边界等的地图，为了说明sf数据集是什么样子的，我们导入了一个数据集，描述了澳大利亚各州和领土的边界:

library(ozmaps)
library(sf)

oz_states <- ozmaps::ozmap_states
oz_states
#> Simple feature collection with 9 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 106 ymin: -43.6 xmax: 168 ymax: -9.23
#> Geodetic CRS:  GDA94
#> # A tibble: 9 × 2
#>   NAME                                                                  geometry
#> * <chr>                                                       <MULTIPOLYGON [°]>
#> 1 New South Wales   (((151 -35.1, 151 -35.1, 151 -35.1, 151 -35.1, 151 -35.2, 1…
#> 2 Victoria          (((147 -38.7, 147 -38.7, 147 -38.7, 147 -38.7, 147 -38.7)),…
#> 3 Queensland        (((149 -20.3, 149 -20.4, 149 -20.4, 149 -20.3)), ((149 -20.…
#> 4 South Australia   (((137 -34.5, 137 -34.5, 137 -34.5, 137 -34.5, 137 -34.5, 1…
#> 5 Western Australia (((126 -14, 126 -14, 126 -14, 126 -14, 126 -14)), ((124 -16…
#> 6 Tasmania          (((148 -40.3, 148 -40.3, 148 -40.3, 148 -40.3)), ((147 -39.…
#> # … with 3 more rows

该输出显示了与数据相关的一些元数据(稍后讨论)，并告诉我们该数据实际上是一个9行2列的tibble。sf数据的一个优势是显而易见的，我们可以很容易地看到数据的整体结构:澳大利亚由6个州和一些地区组成。这里有9个不同的地理单位，所以在这个tibble中有9行(参见mi_counties数据，每个多边形顶点有一行)。

最重要的列是geometry，它指定每个州和地区的空间几何图形。几何列中的每个元素都是多个多边形对象，顾名思义，它包含指定一个或多个多边形顶点的数据，这些顶点划分了区域的边界。对于这种格式的数据，我们可以使用geom_sf()和coord_sf()来绘制一个可用的地图，而不需要指定任何参数，甚至不需要明确声明任何图形属性:

ggplot(oz_states) + 
  geom_sf() + 
  coord_sf()

image

要理解为什么这样做，请注意geom_sf()依赖于ggplot2中其他地方没有使用的几何对象。这种图形属性可以通过以下三种方式指定：

• 在最简单的情况下（如上所示），当用户什么都不做时， geom_sf()将尝试将其映射到名为geometry的列。

• 如果data参数是 sf 对象，则geom_sf()可以自动检测几何列，即使它没有调用geometry。

• 您可以使用 aes(geometry = my_column)以通常的方式手动指定映射。如果您有多个列，这很有用。

coord_sf()函数管理在后续讨论地图投影。

6.2.1 分层地图

在某些情况下，您可能希望将一个地图覆盖在另一个地图之上。ggplot2包支持这一点，它允许您向一个绘图添加多个geom_sf()层。例如，我将使用oz_states数据以不同的颜色绘制澳大利亚各州，并将澳大利亚地区的边界覆盖在此图上。为此，需要执行两个预处理步骤。首先，我将使用dplyr::filter()从州边界删除“Other Territories”。

下面的代码绘制了一个带有两个地图层的图:第一个使用oz_states以不同的颜色填充各州，第二个使用oz_votes绘制选区边界。其次，我将使用rmapshaper包中的ms_simplify()函数以简化的形式提取选区边界如果原始数据集(在这种情况下是ozmaps::abs_ced),存储在比你的绘图要求更高的分辨率，这通常是一个好主意，以减少渲染绘图所花费的时间。

oz_states <- ozmaps::ozmap_states %>% filter(NAME != "Other Territories")
oz_votes <- rmapshaper::ms_simplify(ozmaps::abs_ced)
#> Registered S3 method overwritten by 'geojsonlint':
#>   method         from 
#>   print.location dplyr

现在有了分别代表州边界和选区边界的数据集oz_states和oz_votes，可以通过向图形中添加两个geom_sf()层来构建理想的绘图结果:

ggplot() + 
  geom_sf(data = oz_states, mapping = aes(fill = NAME), show.legend = FALSE) +
  geom_sf(data = oz_votes, fill = NA) + 
  coord_sf()

image

值得注意的是，该图的第一层将fill图形属性映射到数据中的变量上。在这种情况下，NAME变量是一个分类变量，不传达任何附加信息，但可以使用相同的方法来可视化其他类型的区域元数据。例如，如果oz_states有一个额外的列，指定每个州的失业率，我们可以将fill图形属性映射到该变量。

6.2.2 地图标签

向地图添加标签是注释图（第8章）的一个示例，并受[geom_sf_label()](https://ggplot2.tidyverse.org/reference/ggsf.html)和支持[geom_sf_text()](https://ggplot2.tidyverse.org/reference/ggsf.html)。例如，虽然可以合理地期望澳大利亚观众知道澳大利亚各州的名称（并且在上图中没有标记），但很少有澳大利亚人会知道悉尼大都市区不同选区的名称。为了绘制悉尼的选举地图，我们首先需要提取相关选区的地图数据，然后添加标签。下图通过指定xlim和ylimin来放大悉尼地区[coord_sf()](https://ggplot2.tidyverse.org/reference/ggsf.html)，然后使用[geom_sf_label()](https://ggplot2.tidyverse.org/reference/ggsf.html)标签覆盖每个选区：
给地图添加标签是注释图的一个例子，它由geom_sf_label()和geom_sf_text()支持。例如，虽然澳大利亚人可能会知道澳大利亚各州的名字(在上面的图中没有标注)，但很少澳大利亚人会知道悉尼大都市区不同选区的名称。为了绘制悉尼的选区地图，我们首先需要提取相关选区的的地图数据，然后添加标签。下面的图通过在coord_sf()中指定xlim和ylim来放大悉尼地区，然后使用geom_sf_label()用标签覆盖每个选区

# filter electorates in the Sydney metropolitan region
sydney_map <- ozmaps::abs_ced %>% filter(NAME %in% c(
  "Sydney", "Wentworth", "Warringah", "Kingsford Smith", "Grayndler", "Lowe", 
  "North Sydney", "Barton", "Bradfield", "Banks", "Blaxland", "Reid", 
  "Watson", "Fowler", "Werriwa", "Prospect", "Parramatta", "Bennelong", 
  "Mackellar", "Greenway", "Mitchell", "Chifley", "McMahon"
))

# draw the electoral map of Sydney
ggplot(sydney_map) + 
  geom_sf(aes(fill = NAME), show.legend = FALSE) + 
  coord_sf(xlim = c(150.97, 151.3), ylim = c(-33.98, -33.79)) + 
  geom_sf_label(aes(label = NAME), label.padding = unit(1, "mm"))
#> Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
#> give correct results for longitude/latitude data

image

警告信息值得注意。在内部[geom_sf_label()](https://ggplot2.tidyverse.org/reference/ggsf.html)使用[st_point_on_surface()](https://r-spatial.github.io/sf/reference/geos_unary.html)sf 包中的函数 来放置标签，并且出现警告消息是因为 sf 用于计算几何量（例如质心、内部点）的大多数算法都是基于点位于平面二维上的假设表面并用笛卡尔坐标参数化。这种假设不是严格保证的，并且在某些情况下（例如，极地附近的区域）以这种方式处理经度和纬度的计算将给出错误的答案。出于这个原因， sf 包在依赖此近似值时会产生警告消息。
警告信息是值得注意的。在内部geom_sf_label()使用来自sf包的函数st_point_on_surface()来放置标签，并且出现警告消息是因为sf使用的大多数算法来计算几何数量(例如，质心，内点)是基于一个假设，即这些点位于平面二维表面上，并用笛卡尔坐标参数化。这种假设并不严格，在某些情况下(例如，靠近极点的地区)，以这种方式处理经度和纬度的计算将给出错误的答案。因此，sf包在依赖这个近似时产生警告消息。

6.2.3 添加其他几何对象

虽然geom_sf()在某些方面很特殊的，但它的行为与任何其他geom基本相同，允许使用标准geom在地图上绘制额外的数据。例如，我们可能希望使用geom_point()在地图上绘制澳大利亚首都城市的位置。下面的代码演示了如何做到这一点:

oz_capitals <- tibble::tribble( 
  ~city,           ~lat,     ~lon,
  "Sydney",    -33.8688, 151.2093,  
  "Melbourne", -37.8136, 144.9631, 
  "Brisbane",  -27.4698, 153.0251, 
  "Adelaide",  -34.9285, 138.6007, 
  "Perth",     -31.9505, 115.8605, 
  "Hobart",    -42.8821, 147.3272, 
  "Canberra",  -35.2809, 149.1300, 
  "Darwin",    -12.4634, 130.8456, 
)

ggplot() + 
  geom_sf(data = oz_votes) + 
  geom_sf(data = oz_states, colour = "black", fill = NA) + 
  geom_point(data = oz_capitals, mapping = aes(x = lon, y = lat), colour = "red") + 
  coord_sf()

image

在这个示例中，geom_point仅用于指定首都城市的位置，但其基本思想可以扩展到更一般地处理点元数据。例如，如果oz_capitals数据包含一个额外的变量，该变量指定每个大都市区域内的选民数量，那么我们可以使用size图形属性对该数据进行编码。

6.3 光栅图

为制图提供地理空间信息的第二种方法是依赖栅格数据。与简单特征格式不同，在简单特征格式中，地理实体是由一组线、点和多边形指定的，栅格采用图像的形式。在最简单的情况下，光栅数据可能只是一个位图文件，但有许多不同的图像格式。特别是在地理空间方面，有些图像格式包含元数据(例如大地基准面、坐标参考系统)，可用于将图像信息映射到地球表面。例如，一种常见的格式是GeoTIFF，这是一个常规的TIFF文件，并提供了额外的元数据。幸运的是，在GDAL(地理空间数据抽象库，https://gdal.org/)的帮助下，大多数格式都可以轻松地读入R。例如，sf包包含一个函数sf::gdal_read()，它提供从R访问GDAL光栅驱动程序。然而，您很少需要直接调用这个函数，因为有其他高级函数会为您处理这个问题。

例如，假设我们希望绘制澳大利亚气象局 (BOM) 在其 FTP 服务器上公开提供的卫星图像。bomrang 包²⁴为服务器提供了一个方便的接口，包括一个get_available_imagery()返回文件名向量的get_satellite_imagery()函数和一个下载文件并将其直接导入 R的函数。然而，为了说明的目的，我将使用一个更灵活的方法可适配任何FTP服务器，使用[download.file()](https://rdrr.io/r/utils/download.file.html)功能：
举例来说，假设我们希望绘制澳大利亚气象局(BOM)在其FTP服务器上公开提供的卫星图像。bomrang包为服务器提供了一个方便的接口,包括一个get_available_imagery()函数,它返回一个向量的文件名和get_satellite_imagery()函数,下载一个文件,直接进口到R .用于说明的目的,我将使用一个更灵活的方法,可以适应任何FTP服务器,并使用download.file()函数:

# list of all file names with time stamp 2020-01-07 21:00 GMT 
# (BOM images are retained for 24 hours, so this will return an
# empty vector if you run this code without editing the time stamp)
files <- bomrang::get_available_imagery() %>%
  stringr::str_subset("202001072100") 

# use curl_download() to obtain a single file, and purrr to 
# vectorise this operation
purrr::walk2(
  .x = paste0("ftp://ftp.bom.gov.au/anon/gen/gms/", files),
  .y = file.path("raster", files),
  .f = ~ download.file(url = .x, destfile = .y)
)

请注意，如果您想自己运行此代码，您需要将时间戳字符串从"202001072100"当前日期的前一天更改为当前日期的前一天，并且您需要确保在您的工作目录中有一个名为“raster”的文件夹，其中文件将被下载。在本地缓存文件后，我们可以检查已下载的文件列表：

dir("raster")
#> [1] "IDE00421.202001072100.tif" "IDE00422.202001072100.tif"

所有 14 个文件均由日本气象厅运营的 Himawari-8 地球静止卫星拍摄的图像构建，并在 13 个不同波段拍摄图像。澳大利亚 BOM 发布的图像包括可见光谱（通道 3）和红外光谱（通道 13）的数据：

img_vis  <- file.path("raster", "IDE00422.202001072100.tif")
img_inf <- file.path("raster", "IDE00421.202001072100.tif")

要将img_visible文件中的数据导入R，我将使用stars包将数据导入为stars对象：

library(stars)
sat_vis <- read_stars(img_vis, RasterIO = list(nBufXSize = 600, nBufYSize = 600))
sat_inf <- read_stars(img_inf, RasterIO = list(nBufXSize = 600, nBufYSize = 600))

在上面的代码中，第一个参数指定了光栅文件的路径，该RasterIO参数用于将低级参数列表传递给 GDAL。在这种情况下，我使用nBufXSize和nBufYSize确保 R 以低分辨率读取数据（作为 600x600 像素图像）。要查看 R 导入了哪些信息，我们可以检查sat_vis对象：

sat_vis
#> stars object with 3 dimensions and 1 attribute
#> attribute(s), summary of first 1e+05 cells:
#>                            Min. 1st Qu. Median Mean 3rd Qu. Max.
#> IDE00422.202001072100.tif     0       0      0 18.1       0  255
#> dimension(s):
#>      from  to   offset    delta                  refsys point values x/y
#> x       1 600 -5500000  18333.3 Geostationary_Satellite FALSE   NULL [x]
#> y       1 600  5500000 -18333.3 Geostationary_Satellite FALSE   NULL [y]
#> band    1   3       NA       NA                      NA    NA   NULL

这个输出告诉我们一些关于星星对象结构的信息。对于sat_vis对象，底层数据存储为一个三维数组，其中x和y维度指定空间数据。band在这种情况下，维度对应于颜色通道 (RGB)，但由于数据是灰度的，因此对于该图像来说是多余的。在其他数据集中，可能有对应于不同传感器的波段，也可能有时间维度。请注意，空间数据还与坐标参考系统（在输出中称为“refsys”）相关联。

要在 ggplot2 中绘制数据sat_vis，我们可以使用stars 包提供的函数geom_stars()。图形可能是这样的：

ggplot() + 
  geom_stars(data = sat_vis) + 
  coord_equal()

image

该geom_stars()函数要求data参数是一个星星对象，并将栅格数据映射到fill图形属性。因此，上面卫星图像中的蓝色阴影是由 ggplot2 比例决定的，而不是图像本身。也就是说，虽然sat_vis包含三个波段，但上面的图只显示第一个波段，原始数据值（范围从 0 到 255）映射到 ggplot2 用于连续数据的默认蓝色调色板。要查看图像文件“真正”的样子，我们可以使用facet_wrap()以下命令分离波段：

ggplot() + 
  geom_stars(data = sat_vis, show.legend = FALSE) +
  facet_wrap(vars(band)) + 
  coord_equal() + 
  scale_fill_gradient(low = "black", high = "white")

image

仅显示原始图像的一个限制是不容易确定相关陆地的位置，我们可能希望将卫星数据与oz_states矢量图重叠以显示澳大利亚政治实体的轮廓。但是，这样做需要小心，因为两个数据源与不同的坐标参考系统相关联。要正确投影数据oz_states，应使用sf 包中的函数st_transform()转换数据。在下面的代码中，我从sat_vis栅格对象中提取 CRS ，并将oz_states数据转换到相同的系统下使用。

oz_states <- st_transform(oz_states, crs = st_crs(sat_vis))

完成后，我现在可以在光栅图像的顶部绘制矢量图，以使图像更易于读者理解。现在通过检查可以清楚地看到卫星图像是在澳大利亚日出时拍摄的：

ggplot() + 
  geom_stars(data = sat_vis, show.legend = FALSE) +
  geom_sf(data = oz_states, fill = NA, color = "white") + 
  coord_sf() + 
  theme_void() + 
  scale_fill_gradient(low = "black", high = "white")

image

如果我们想在顶部绘制更常规的数据怎么办？一个简单的例子是根据oz_capitals包含纬度和经度数据的数据框绘制澳大利亚首府城市的位置。然而，由于这些数据不与CRS相关联，并且不在同一规模作为栅格数据sat_vis，这将需要进行改造。为此，我们首先需要使用oz_capitals以下数据从数据创建一个 sf 对象[st_as_sf()](https://r-spatial.github.io/stars/reference/st_as_sf.html)：

cities <- oz_capitals %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = FALSE)

此投影是使用 EPSG 代码 4326 设置的，这是一种使用纬度和经度值作为坐标并依赖于 WGS84 数据的椭球投影。完成后，我们现在可以将坐标从经纬度几何转换为匹配我们数据sat_vis的几何对象：

cities <- st_transform(cities, st_crs(sat_vis))

现在可以使用geom_sf()以下方法覆盖转换后的数据：

ggplot() + 
  geom_stars(data = sat_vis, show.legend = FALSE) +
  geom_sf(data = oz_states, fill = NA, color = "white") + 
  geom_sf(data = cities, color = "red") + 
  coord_sf() + 
  theme_void() + 
  scale_fill_gradient(low = "black", high = "white")

image

这个版本的图像更清楚地表明卫星图像是在达尔文大约日出时拍摄的：所有东部城市的太阳都已经升起，但珀斯没有。这可以在数据可视化中使用geom_sf_text()为每个城市添加标签的功能更加清晰。例如，我们可以使用这样的代码向绘图添加另一个图层，

geom_sf_text(data = cities, mapping = aes(label = city)) 

如何更好的设定文本的位置（参见注释）。