ggplot2地图

ggplot2地图

install.packages(c(
  "colorBlindness", "directlabels", "dplyr", "gameofthrones", "ggforce", 
  "gghighlight", "ggnewscale", "ggplot2", "ggraph", "ggrepel", 
  "ggtext", "ggthemes", "hexbin", "Hmisc", "mapproj", "maps", 
  "munsell", "ozmaps", "paletteer", "patchwork", "rmapshaper", "scico", 
  "seriation", "sf", "stars", "tidygraph", "tidyr", "wesanderson" 
))

地图
绘制地理空间数据是一项常见的可视化任务，需要专门的工具。通常，问题可以分解为两个问题：使用一个数据源绘制地图，以及将来自另一个信息源的元数据添加到地图中。本章将帮助您解决这两个问题。我按照以下方式组织了这一章：第6.1节概述了一种使用. 接下来，第6.3节和第6.4节讨论如何使用地图投影和底层 sf 数据结构。最后，第6.5节讨论了如何根据栅格数据绘制地图。geom_polygon()geom_sf()
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()

在这两个图中，我使用coord_quickmap()调整轴以确保经度和纬度以相同的比例呈现。第16章以更一般的术语讨论 ggplot2 中的坐标系，但正如我们将在下面看到的，地理空间数据通常需要更严格的方法。出于这个原因，ggplot2 提供geom_sf()并coord_sf()处理以简单特征格式指定的空间数据。

6.2简单的特征图

上述方法存在一些限制，其中最重要的是简单的“经纬度”数据格式通常不用于现实世界的制图。地图的矢量数据通常使用开放地理空间联盟制定的“简单特征”标准进行编码。Edzer Pebesma https://github.com/r-spatial/sf开发的 sf 包19为处理此类数据提供了出色的工具集，并且 ggplot2 中的和函数旨在与 sf 包一起使用。geom_sf()coord_sf()

为了介绍这些功能，我们依靠 Michael Sumner 的 ozmaps 包https://github.com/mdsumner/ozmaps/提供澳大利亚州边界、地方政府区域、选举边界等的地图。20为了说明 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 列的小标题。sf 数据的一个优势是显而易见的，我们可以很容易地看到数据的整体结构：澳大利亚由六个州和一些地区组成。有 9 个不同的地理单位，所以这个 tibble 中有 9 行（参见 mi_counties data每个多边形顶点有一行）。

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

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

要理解为什么会这样，请注意它geom_sf()依赖geometry于 ggplot2 中其他地方未使用的美学。这种美学可以通过以下三种方式之一来指定：

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

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

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

该coord_sf()函数控制地图投影，在第6.3节中讨论。

6.2.1分层地图
在某些情况下，您可能希望将一张地图叠加在另一张地图上。ggplot2 包通过允许您向geom_sf()绘图添加多个图层来支持这一点。例如，我将使用这些oz_states数据以不同颜色绘制澳大利亚各州，并将此图与澳大利亚选举区的边界重叠。为此，需要执行两个预处理步骤。首先，我将使用dplyr::filter()从州边界中删除“其他领土”。

下面的代码绘制了一个包含两个地图图层的图：第一个用于oz_states填充不同颜色的州，第二个用于oz_votes绘制选举边界。ms_simplify()其次，我将使用rmapshaper 包中的函数以简化形式提取选举边界。21如果原始数据集（在这种情况下）以比您的绘图所需的分辨率更高的分辨率存储，这通常是一个好主意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()

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

6.2.2标注地图

向地图添加标签是注释图（第8章）的一个示例，并且由geom_sf_label()和支持geom_sf_text()。例如，虽然可以合理地期望澳大利亚观众知道澳大利亚各州的名称（并且在上图中未标记），但很少有澳大利亚人会知道悉尼大都市区不同选民的名称。为了绘制悉尼的选举地图，我们首先需要提取相关选区的地图数据，然后添加标签。下图通过指定xlim和ylimin放大悉尼地区coord_sf()，然后使用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

警告信息值得注意。内部geom_sf_label()使用 sf 包中的函数 st_point_on_surface()来放置标签，并且出现警告消息是因为 sf 用于计算几何量（例如质心、内部点）的大多数算法都基于点位于平面二维上的假设表面并用笛卡尔坐标参数化。这个假设没有严格的保证，并且在某些情况下（例如，两极附近的区域）以这种方式处理经度和纬度的计算会给出错误的答案。出于这个原因，当 sf 包依赖这个近似值时，它会产生警告消息。

6.2.3添加其他geom

尽管geom_sf()在某些方面很特别，但它的行为方式与任何其他几何图形大致相同，允许在具有标准几何图形的地图上绘制附加数据。例如，我们可能希望使用 . 在地图上绘制澳大利亚首府城市的位置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()

在此示例geom_point中仅用于指定首府城市的位置，但可以扩展基本思想以更普遍地处理点元数据。例如，如果 oz_capitals 数据要包含一个附加变量，指定每个大都市区内的选民数量，我们可以使用size美学对该数据进行编码。

6.3地图投影

在本章的开头，我通过在笛卡尔平面上绘制经度和​​纬度来绘制地图，就好像地理空间数据与人们可能想要绘制的其他类型的数据没有什么不同。对于第一个近似值，这是可以的，但如果您关心准确性，这还不够好。这种方法有两个基本问题。

第一个问题是行星的形状。地球既不是一个平面，也不是一个完美的球体。因此，要将坐标值（经度和纬度）映射到某个位置，我们需要对各种事物做出假设。地球有多椭圆？地球的中心在哪里？经度和纬度的原点在哪里？海平面在哪里？构造板块如何运动？所有这些事情都是相关的，并且根据一个人所做的假设，相同的坐标可以映射到相距数米的位置。关于地球形状的一组假设被称为大地基准虽然对于某些数据可视化可能无关紧要，但对于其他数据可视化来说却很重要。人们可能会考虑几种不同的选择：如果您的重点是北美，那么“北美基准面”（NAD83）是一个不错的选择，而如果您的观点是全球性的，那么“世界大地测量系统”（WGS84）可能会更好。

第二个问题是地图的形状。地球大致呈椭圆形，但在大多数情况下，您的空间数据需要绘制在二维平面上。将椭球的表面映射到平面上是不可能的，而没有一些变形或切割，您必须选择在绘制地图时准备接受哪些变形。这是地图投影的工作。

地图投影通常根据它们保留的几何属性进行分类，例如

保留面积投影确保地球上面积相等的区域在地图上以相等的面积绘制。

保形（或保形）投影可确保保留区域的局部形状。

不幸的是，任何投影都不可能保形保面积。这使得详细讨论地图投影有点超出了本书的范围，除了要注意简单的功能规范允许您指示要使用的地图投影。有关地图投影的更多信息，请参阅使用 R 的地理计算https://geocompr.robinlovelace.net/。22

综合起来，大地基准（例如 WGS84）、地图投影的类型（例如墨卡托）和投影参数（例如原点位置）指定了一个坐标参考系统或 CRS，这是一套完整的假设用于将经纬度信息转换成二维地图。一个 sf 对象通常包含一个默认的 CRS，如下图所示：

st_crs(oz_votes)
#> Coordinate Reference System:
#>   User input: EPSG:4283 
#>   wkt:
#> GEOGCRS["GDA94",
#>     DATUM["Geocentric Datum of Australia 1994",
#>         ELLIPSOID["GRS 1980",6378137,298.257222101,
#>             LENGTHUNIT["metre",1]]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433]],
#>     CS[ellipsoidal,2],
#>         AXIS["geodetic latitude (Lat)",north,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         AXIS["geodetic longitude (Lon)",east,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>     USAGE[
#>         SCOPE["Horizontal component of 3D system."],
#>         AREA["Australia including Lord Howe Island, Macquarie Islands, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore."],
#>         BBOX[-60.56,93.41,-8.47,173.35]],
#>     ID["EPSG",4283]]

此输出的大部分对应于明确描述 CRS的众所周知的文本(WKT) 字符串。sf 在内部使用这种详细的 WKT 表示，但是有几种方法可以提供 sf 理解的用户输入。一种这样的方法是以EPSG 代码的形式提供数字输入（参见http://www.epsg.org/）。数据中的默认CRSoz_votes对应EPSG代码4283：

st_crs(oz_votes) == st_crs(4283)
#> [1] TRUE
在 ggplot2 中，CRS 由 控制coord_sf()，这确保图中的每一层都使用相同的投影。默认情况下，使用与数据23coord_sf()的几何列关联的 CRS 。因为 sf 数据通常会提供一个明智的 CRS 选择，所以这个过程通常是不可见的，不需要用户干预。但是，如果您需要自己设置 CRS，您可以通过将有效的用户输入传递给. 下面的示例说明了如何从默认 CRS 切换到 EPSG 代码 3112：crsst_crs()

ggplot(oz_votes) + geom_sf()
ggplot(oz_votes) + geom_sf() + coord_sf(crs = st_crs(3112))

6.4使用 sf 数据

如前所述，使用geom_sf()并coord_sf()严重依赖 sf 包提供的工具创建的地图，实际上 sf 包包含许多用于操作简单要素数据的有用工具。在本节中，我将介绍一些此类工具；更详细的介绍可以在 sf 包网站 https://r-spatial.github.io/sf/上找到。

首先，回想一下简单特征相对于空间数据的其他表示的一个优势是地理单位可以具有复杂的结构。澳大利亚地图数据中的一个很好的例子是 Eden-Monaro 的选区，如下图所示：

edenmonaro <- ozmaps::abs_ced %>% filter(NAME == "Eden-Monaro")

p <- ggplot(edenmonaro) + geom_sf()
p + coord_sf(xlim = c(147.75, 150.25), ylim = c(-37.5, -34.5)) 
p + coord_sf(xlim = c(150, 150.25), ylim = c(-36.3, -36)) 

如图所示，Eden-Monaro 被定义为两个不同的多边形，一个位于澳大利亚大陆的大多边形和一个小岛。然而，大区域的中间有一个洞（这个洞的存在是因为澳大利亚首都地区是一个完全包含在伊甸园-莫纳罗内的独特政治单位，并且如上所示，澳大利亚的选举边界不跨越州界）。MULTIPOLYGON在科幻术语中，这是几何的一个例子。在本节中，我将讨论这些对象的结构以及如何使用它们。

首先，让我们使用 dplyr 只抓取几何对象：

edenmonaro <- edenmonaro %>% pull(geometry)
edenmonaro 对象的元数据可以使用辅助函数访问。例如，st_geometry_type()提取几何类型（例如MULTIPOLYGON），st_dimension()提取维数（XY 数据为 2，XYZ 为 3），st_bbox()将边界框提取为数字向量，并将st_crs()CRS 提取为具有两个分量的列表，一个用于EPSG 代码和另一个用于 proj4string。例如：

st_bbox(edenmonaro)
#>  xmin  ymin  xmax  ymax 
#> 147.7 -37.5 150.2 -34.5

通常，当我们打印edenmonaro对象时，输出会显示所有附加信息（尺寸、边界框、大地基准等），但对于本节的其余部分，我将仅显示输出的相关行。在这种情况下，edenmonaro 由包含一个特征的 MULTIPOLYGON 几何定义：

edenmonaro 
#> Geometry set for 1 feature 
#> Geometry type: MULTIPOLYGON
#> MULTIPOLYGON (((150 -36.2, 150 -36.2, 150 -36.3...

但是，我们可以将 MULTIPOLYGON “投射”到两个不同的 POLYGON 几何形状中，使用st_cast()：

st_cast(edenmonaro, "POLYGON")
#> Geometry set for 2 features 
#> Geometry type: POLYGON
#> POLYGON ((150 -36.2, 150 -36.2, 150 -36.3, 150 ...
#> POLYGON ((148 -36.7, 148 -36.7, 148 -36.7, 148 ...

为了说明这在什么时候可能有用，请考虑道森选区，该选区由 69 个岛屿和澳大利亚大陆的一个沿海地区组成。

dawson <- ozmaps::abs_ced %>% 
  filter(NAME == "Dawson") %>% 
  pull(geometry)
dawson
#> Geometry set for 1 feature 
#> Geometry type: MULTIPOLYGON
#> MULTIPOLYGON (((148 -19.9, 148 -19.8, 148 -19.8...

ggplot(dawson) + 
  geom_sf() +
  coord_sf()

然而，假设我们的兴趣只在于绘制岛屿地图。如果是这样，我们可以首先使用该st_cast()函数将 Dawson 选民划分为组成多边形。这样做之后，我们可以使用st_area()计算每个多边形的面积并which.max()找到面积最大的多边形：

dawson <- st_cast(dawson, "POLYGON")
which.max(st_area(dawson))
#> [1] 69

大大陆区域对应于道森内的第 69 个多边形。有了这些知识，我们可以绘制一张只显示岛屿的地图：

ggplot(dawson[-69]) + 
  geom_sf() + 
  coord_sf()

6.5栅格地图

为制图提供地理空间信息的第二种方法是依靠栅格数据。与简单的要素格式不同，其中地理实体是根据一组线、点和多边形来指定的，栅格采用图像的形式。在最简单的情况下，光栅数据可能只不过是一个位图文件，但那里有许多不同的图像格式。具体而言，在地理空间环境中，存在包含元数据（例如，大地基准、坐标参考系统）的图像格式，可用于将图像信息映射到地球表面。例如，一种常见的格式是 GeoTIFF，它是一个带有附加元数据的常规 TIFF 文件。令人高兴的是，在 GDAL（地理空间数据抽象库，https://gdal.org/）。例如， sf 包包含一个函数sf::gdal_read()，它提供从 R 访问 GDAL 光栅驱动程序的功能。但是，您很少需要直接调用此函数，因为还有其他高级函数可以为您处理这个问题。

例如，假设我们希望绘制由澳大利亚气象局 (BOM) 在其 FTP 服务器上公开提供的卫星图像。bomrang 包24为服务器提供了一个方便的接口，包括一个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 包25将数据导入为 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。在这种情况下，我使用nBufXSizeandnBufYSize来确保 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”）。

要绘制sat_visggplot2 中的数据，我们可以使用geom_stars()stars 包提供的函数。一个最小的情节可能如下所示：

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

该geom_stars()函数要求data参数是星星对象，并将栅格数据映射到填充美学。因此，上面卫星图像中的蓝色阴影是由 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")

仅显示原始图像的一个限制是不容易确定相关陆地的位置，我们可能希望将卫星数据与oz_states矢量图叠加以显示澳大利亚政治实体的轮廓。但是，这样做需要一些小心，因为这两个数据源与不同的坐标参考系统相关联。要正确投影oz_states数据，应使用st_transform()sf 包中的函数转换数据。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")

如果我们想在顶部绘制更多常规数据怎么办？一个简单的例子是根据oz_capitals包含纬度和经度数据的数据框绘制澳大利亚首府城市的位置。但是，由于这些数据与 CRS 无关，并且与 中的栅格数据的比例不同，因此也需要对这些数据进行转换sat_vis。为此，我们首先需要使用以下数据从数据中创建一个 sf 对象：oz_capitalsst_as_sf()

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")

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

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

尽管需要注意确保文本定位良好（参见第8章）。

6.6数据来源

USAboundaries 包https://github.com/ropensci/USAboundaries包含美国的州、县和邮政编码数据。26除了目前的边界，它还有可以追溯到 1600 年代的州和县边界。

tigris 软件包https://github.com/walkerke/tigris可以轻松访问美国人口普查 TIGRIS shapefile。27它包含州、县、邮政编码和人口普查区边界，以及许多其他有用的数据集。

rnaturalearth 包28捆绑了来自http://naturalearthdata.com/的免费、高质量数据。它包含国家边界和每个国家内顶级区域的边界（例如美国的州、法国的地区、英国的县）。

osmar 包https://cran.r-project.org/package=osmar封装了 OpenStreetMap API，因此您可以访问广泛的矢量数据，包括单个街道和建筑物29

如果您有自己的形状文件 ( .shp)，您可以将它们加载到 R 中sf::read_sf()