(2)R中的地理数据

发表于 更新于 分类于 ， 阅读次数： Waline： 本文字数： 16k 阅读时长 ≈ 57 分钟

本文由SCY翻译自《Geocomputation with R》第二章

R中的地理数据

本章将简要介绍两种基本的地理数据模型：矢量和栅格。我们将介绍每种数据模型背后的理论以及它们在哪些学科中占主导地位,然后演示它们在R中的实现。

前提条件

这是本书的第一章实践内容,因此它带有一些软件需求。

您需要能够访问安装了最新版本R的计算机(R4.2.0 或更高版本)。我们建议不仅阅读文本内容,也运行每章的代码,以建立您的地理计算技能。

为了跟踪您的学习过程,值得在计算机上新建一个文件夹来保存您的R脚本、输出和其他在学习《用R进行地理计算》过程中产生的东西。您也可以下载或克隆本书背后的源代码以支持您的学习。我们强烈建议在编写/运行/测试R代码时安装集成开发环境(IDE),如RStudio(针对大多数人推荐)或VS Code。^[我们建议使用RStudio项目、VS Code工作区或类似系统来管理项目。通过RStudio中的R控制台,可以使用rstudioapi包快速完成此操作。例如,使用以下命令在主目录中打开一个名为“geocompr-learning”的新项目:rstudioapi::openProject("~/geocompr-learning")。]

如果您是R的新手,我们建议您在深入学习《用R进行地理计算》代码之前,首先学习一些R入门资源,比如Garrett Grolemund的《用R学习编程》或Claudia Engel的《R语言入门》 。组织好您的工作(例如使用RStudio项目),并给脚本起合理的名称,如chapter-02-notes.R,以记录您在学习过程中的代码。

安装好环境后,是时候运行一些代码了!除非您已经安装了这些包,否则第一件事是使用以下命令安装本章需要的基础R包:^[spDataLarge 不存在于CRAN上,意味着必须通过r-universe或以下命令安装:remotes::install_github("Nowosad/spDataLarge")。]

1
2
3
4
install.packages("sf")
install.packages("terra")
install.packages("spData")
install.packages("spDataLarge", repos = "https://nowosad.r-universe.dev")

📌我们建议按照CRAN上的R安装说明进行操作。
如果您使用Mac或Linux,上述安装sf的命令第一次可能不会起作用。
这些操作系统 (OSs) 都有一些“系统需求”,在该包的README中有描述。
网上可以找到其他特定操作系统的安装说明,包括rtask.thinkr.fr网站上的文章在Ubuntu 22.04.1 LTS上安装R 4.2及空间包提示。

复现本书第1部分内容所需的包可以通过以下命令安装:remotes::install_github("geocompx/geocompkg")。这个命令使用remotes包中的install_packages()函数来安装托管在GitHub代码托管、版本控制和协作平台上的源代码。以下命令将安装复现整本书所需的所有依赖项(警告:这可能需要几分钟):remotes::install_github("geocompx/geocompkg", dependencies = TRUE)本章中提出的代码所需的包可以通过library()函数“加载”(技术上是附加)如下:

1
2
library(sf)          # classes and functions for vector data
#> Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

library(sf)的输出将报告该包正在使用的关键地理库(如GEOS)的版本,如第2.2.1节所述。

1
library(terra)      # classes and functions for raster data

我们安装的其他R包包含了书中将使用的数据:

1
2
3
4
5
6
7
8
9
10
library(spData)        # load geographic data
#> The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
#> which was just loaded, will retire in October 2023.
#> Please refer to R-spatial evolution reports for details, especially
#> https://r-spatial.org/r/2023/05/15/evolution4.html.
#> It may be desirable to make the sf package available;
#> package maintainers should consider adding sf to Suggests:.
#> The sp package is now running under evolution status 2
#>      (status 2 uses the sf package in place of rgdal)
library(spDataLarge)   # load larger geographic data

引言

本章将简要介绍两种基本的地理数据模型:矢量和栅格。我们将介绍每种数据模型背后的理论以及它们在哪些学科中占主导地位,然后演示它们在R中的实现。

矢量数据模型使用点、线和多边形来表示世界。这些要素具有明确定义的边界,这意味着矢量数据集通常具有很高的精度(但不一定是准确性,我们将在第2.5节中看到)。栅格数据模型将表面划分为大小恒定的单元格。栅格数据集是 web 地图中使用的背景图像的基础,自航空摄影和卫星遥感设备起源以来,它一直是地理数据的重要来源。栅格将空间特定特征聚合到给定分辨率,这意味着它们在空间上是一致的并且可扩展的(许多全球栅格数据集都是可用的)。

使用哪种模型呢?答案可能取决于您的应用领域:

  • 矢量数据在社会科学中占主导地位,因为人类定居点往往有明确的边界
  • 栅格在许多环境科学中占主导地位,因为它们依赖遥感数据

在一些领域中两者有很大的重叠,矢量和栅格数据集可以一起使用:例如,生态学家和人口统计学家通常同时使用矢量和栅格数据。此外,可以在两种形式之间进行转换。无论您的工作更多地使用矢量数据集还是栅格数据集,在使用它们之前理解其基本数据模型都是值得的,这将在后续章节中讨论。本书使用 sf 和 terra 包分别处理矢量数据和栅格数据集。

矢量数据

📌在本书中使用“矢量”一词时要注意它有两层含义:
地理矢量数据和 R 中的vector类(注意等宽字体字体)。
前者是一个数据模型,后者是一个R类,就像 data.frame 和 matrix。
但两者之间还是有联系的:地理矢量数据模型的核心空间坐标可以用R中的vector对象来表示。
所以在阅读时要注意上下文,区分这两层含义。一般来说,没有特别说明时,“矢量”指的都是地理矢量数据模型。

地理矢量数据模型基于坐标参考系统(CRS)内的点。点可以代表独立的要素(例如公交车站的位置),也可以链接在一起形成更复杂的几何形状,如线和多边形。大多数点几何只包含两个维度(三维坐标参考系统包含额外的$z$值,通常代表海拔高度)。

在这种系统中,例如伦敦可以用坐标c(-0.1, 51.5)表示。这意味着它相对于原点向东-0.1度,向北51.5度。在这种情况下,原点位于地理坐标系(“经纬度”)中的0度经线(本初子午线)和0度纬线(赤道)上（图2.1，左面板）。在投影坐标系中,同一点也可以用“东向/北向”值为c(530000,180000)来概括表示,在英国国家网格中,这意味着伦敦位于CRS$原点$以东530公里,以北180公里。这可以从视觉上进行验证:与表示伦敦的点略微相隔超过5个“框”——由宽度为100公里的灰色网格线划定的正方形区域——从原点（图2.1，右面板）。

英国国家网格原点位于西南半岛以外的海域,以确保英国内陆的大多数位置都具有正的东向和北向值。^[我们提到的原点,在图中用蓝色表示,事实上是“伪”原点。 真正的原点,也就是畸变最小的位置,位于2°W 49°N。这是英国测绘局选择的位置,在经度上大致位于不列颠岛的中心。简而言之,坐标由两个数字组成,表示距离一个原点的距离,通常是在$x$然后$y$维度上。]关于坐标参考系统的更多内容,将在第节和第7章详细描述,但就本节而言,知道这些坐标原理就足够了。


图2.1：矢量（点）数据的示例，其中伦敦的位置（红色X）相对于一个原点（蓝色圆圈）进行了表示。左图代表了一个地理坐标参照系统（CRS），其原点位于0°经度和纬度。右图代表了一个投影坐标参照系统（CRS），其原点位于西南半岛西海域。

sf为地理矢量数据提供了类,并为地理计算的重要底层库提供了一致的命令行接口:

  • GDAL,用于读写和操作各种地理数据格式（第8章）。

  • PROJ,一个强大的坐标系统转换库（第7章）。

  • GEOS,一个平面几何引擎,用于计算缓冲区和投影坐标系数据的质心等操作（第5章）。

  • S2,由Google开发的用C++编写的球面几何引擎,通过s2包（见下文2.2.9节和第7章）。

这些接口的信息会在sf包首次加载时打印出来:在本章开头library(sf)命令下方出现的消息告诉我们连接的GEOS、GDAL和PROJ库的版本(这在不同计算机和时间上有所不同),以及是否开启了S2接口。当今我们认为理所应当,但是只有与不同地理库的紧密集成,才使R中的可重现地理计算成为可能。

sf的一个漂亮功能是您可以更改在未投影数据上使用的默认几何引擎:使用sf::sf_use_s2(FALSE)命令可以“关闭”S2,这意味着所有几何操作(包括未投影数据上的几何操作)将默认使用平面几何引擎GEOS，平面几何基于二维空间（见第2.2.9节）。平面几何引擎(如GEOS)假设“平面”(投影)坐标,而球面几何引擎(如S2)假设未投影(经纬度)坐标。

本节将介绍sf类,为后续章节做准备（（第5章和第8章分别介绍了GEOS和GDAL接口））。

简单特征介绍

Simple Features 是一个由开放地理空间联盟（OGC）开发和认可的开放标准，这是一个非营利组织，我们将在第8.5节中再次讨论其活动。Simple Features 是一个分层数据模型，代表了各种几何类型。规范支持的 18 种几何类型中，只有 7 种在大多数地理研究中被使用（如图2.2所示）。这些核心几何类型得到了 R 包 sf 的完全支持。^[完整的 OGC 标准包括一些相当不寻常的几何类型，包括 ‘surface’ 和 ‘curve’ 几何类型，这些在现实世界应用中目前用途有限。所有 18 种类型都可以用sf包来表示，尽管在写作时（2022年），绘图只适用于“核心 7 种”。]


图2.2：Simple feature types fully supported by sf.

sf可以表示所有常见的矢量几何类型（栅格数据类别不受sf支持）：点、线、多边形及其各自的’multi’版本（将同一类型的特征组合成单一特征）。sf 还支持几何集合，可以在单一对象中包含多种几何类型。sf 提供了以前在三个包中提供的相同（甚至更多）功能—— sp 用于数据类别，rgdal用于通过接口到GDAL和PROJ进行数据读写，以及rgeos用于通过接口到GEOS进行空间操作。

为了重申第一章的信息，地理R包有着与低级库接口的悠久历史，而sf通过统一的接口延续了这一传统，该接口适用于GEOS的近期版本进行几何操作，GDAL库用于读取和写入地理数据文件，以及PROJ库用于表示和转换投影坐标参考系统。

通过对 Google’s 球面几何库 s2的R接口，sf还能快速且准确地进行“在非平面几何上的测量和操作”。自从 1.0.0 版的sf在2021年6月发布以来，s2的功能现在默认用于具有地理（经度/纬度）坐标系统的几何体，这是sf的一个独特特性，与仅支持用于几何操作的GEOS的空间库不同，例如Python包GeoPandas。我们将在后续章节中讨论s2。

sf能够将多个强大的地理计算库整合到一个框架中，这是一个值得注意的成就，它降低了使用高性能库进行可重复地理数据分析的“准入门槛”。sf的功能在其网站r-spatial.github.io/sf/ 上有详细的文档，包含7篇专题文章。这些可以按如下方式离线查看：

1
2
vignette(package = "sf") # see which vignettes are available
vignette("sf1")          # an introduction to the package

正如第一篇专题文章解释的，R中的简单特征对象存储在数据框中，地理数据占据一个特殊的列，通常命名为’geom’或’geometry’。我们将使用本章开始时加载的spData提供的world数据集，来展示什么是sf对象以及它们是如何工作的。world是一个’sf数据框’，包含空间和属性列，这些列的名称由函数names()返回（在本例中，最后一列包含地理信息）：

1
2
3
4
5
class(world)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"
names(world)
#>  [1] "iso_a2"    "name_long" "continent" "region_un" "subregion" "type"     
#>  [7] "area_km2"  "pop"       "lifeExp"   "gdpPercap" "geom"

这个geom列中的内容赋予了sf对象它们的空间能力：world$geom是一个’list column’，其中包含了国家多边形的所有坐标。sf对象可以通过函数plot()快速绘制。尽管plot()是R默认安装（基础 R）的一部分，它是一个 generic 函数，可被其他包扩展。sf包含非导出的（大多数时间对用户隐藏的）plot.sf()函数，这是在以下命令中幕后调用的，该命令创建了图2.3。

1
plot(world)


图2.3：A spatial plot of the world using the sf package, with a facet for each attribute.

注意，与大多数GIS程序默认为地理对象创建单一地图不同，对sf对象进行plot()操作会导致数据集中每个变量都有一个地图。这种行为在探索不同变量的空间分布方面可能很有用，并将在节进一步讨论。

更广泛地说，将地理对象视为具有空间能力的常规数据框有很多优点，特别是如果您已经习惯于使用数据框。例如，常用的summary()函数提供了world对象内变量的有用概览。

1
2
3
4
5
6
7
8
9
summary(world["lifeExp"])
#>     lifeExp                geom    
#>  Min.   :50.6   MULTIPOLYGON :177  
#>  1st Qu.:65.0   epsg:4326    :  0  
#>  Median :72.9   +proj=long...:  0  
#>  Mean   :70.9                      
#>  3rd Qu.:76.8                      
#>  Max.   :83.6                      
#>  NA's   :10

尽管我们只为summary()命令选择了一个变量，它还输出了关于几何的报告。这演示了sf对象的几何列的“粘性”行为，意味着除非用户故意移除它们，否则几何体将被保留，正如我们将在第3.2节看到的。该结果提供了world中包含的非空间和空间数据的快速总结：所有国家的平均预期寿命为71岁（范围从不到51岁到超过83岁，中位数为73岁）。

📌上面摘要输出中的MULTIPOLYGON一词指的是world对象中特征（国家）的几何类型。对于像印度尼西亚和希腊这样有岛屿的国家，这种表示是必要的。其他几何类型在2.2.4部分有描述。

这个简单特征对象的基础行为和内容值得更深入地了解，这可以视为一个“spatial data frame（空间数据框）”。

sf对象很容易进行子集操作：下面的代码展示了如何返回一个仅包含world对象的前两行和前三列的对象。输出显示了与常规data.frame相比有两个主要差异：包含额外的地理元数据（Geometry type、Dimension、Bounding box以及以Geodetic CRS开始的线上的坐标参考系统信息），以及存在一个名为geom的’几何列’。

1
2
3
4
5
6
7
8
9
10
11
12
world_mini = world[1:2, 1:3]
world_mini
#> Simple feature collection with 2 features and 3 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -180 ymin: -18.3 xmax: 180 ymax: -0.95
#> Geodetic CRS:  WGS 84
#> # A tibble: 2 × 4
#>   iso_a2 name_long continent                                                geom
#>   <chr>  <chr>     <chr>                                      <MULTIPOLYGON [°]>
#> 1 FJ     Fiji      Oceania   (((-180 -16.6, -180 -16.5, -180 -16, -180 -16.1, -…
#> 2 TZ     Tanzania  Africa    (((33.9 -0.95, 31.9 -1.03, 30.8 -1.01, 30.4 -1.13,…

所有这些可能看起来相当复杂，尤其是对于一个“简单”的类系统！然而，用这种方式组织事物并使用sf处理矢量地理数据集是有充分理由的。

在描述 sf 包支持的每一种几何类型之前，有必要先了解 sf 对象的基础组成元素。节展示了简单特征对象是数据框，具有特殊的几何列。这些空间列通常被称为 geom 或 geometry：world$geom 是指上文中描述的 world 对象中的空间元素。这些几何列是 sfc 类的“列表列”（见第节）。反过来，sfc 对象由一个或多个 sfg 类的对象组成：我们在节中描述的简单特征几何。

要理解简单特征的空间组件是如何工作的，了解简单特征几何体是至关重要的。因此，在继续描述如何用基于sfg和sfc对象的sf对象在R中表示这些几何体之前，我们将在2.2.4节中介绍当前支持的每一种简单特征几何类型。

📌前面的代码块使用 = 在命令 world_mini = world[1:2, 1:3] 中创建了一个名为 world_mini 的新对象。
这被称为赋值。
实现相同结果的等价命令是 world_mini <- world[1:2, 1:3]。
尽管“箭头赋值”更为常用，我们使用“等号赋值”，因为它打字稍微快一些，而且由于与常用语言如 Python 和 JavaScript 的兼容性，更容易教授。
使用哪一个主要是个人偏好的问题，只要你保持一致（像 styler 这样的包可以用来改变风格）。

为什么是简单特征?

简单特征（Simple features）是一种得到广泛支持的数据模型，它是包括 QGIS 和 PostGIS 在内的许多 GIS 应用程序中的数据结构基础。这一模型的一个主要优点是，使用这种数据模型确保了你的工作可以在其他设置中进行交叉转移，例如从空间数据库导入和导出。

从R的角度来看，一个更具体的问题是“为什么要使用sf包，而sp包已经经过了反复测试”？有很多原因（与简单特征模型的优点有关）：

  • 快速地读取和写入数据
  • 提升了绘图性能
  • sf 对象在大多数操作中可以被视为数据框
  • sf 函数名相对一致和直观（所有函数名都以 st_ 开头）
  • sf 函数可以与 |> 操作符结合使用，并且与 tidyverse 的 R 包集合非常兼容。

sf 对 tidyverse 包的支持通过 read_sf() 函数体现，这是一个第8.6.1节中详细介绍的用于导入地理矢量数据的函数。与返回存储在基础 R data.frame 中的属性的函数 st_read() 不同（该函数发出详细的消息，在下面的代码块中未显示），read_sf() 会无声地返回数据作为一个 tidyverse tibble。下面展示了这一点：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
world_dfr = st_read(system.file("shapes/world.shp", package = "spData"))
#> Reading layer `world' from data source 
#>   `/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/spData/shapes/world.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 177 features and 10 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -180 ymin: -89.9 xmax: 180 ymax: 83.6
#> Geodetic CRS:  WGS 84
world_tbl = read_sf(system.file("shapes/world.shp", package = "spData"))
class(world_dfr)
#> [1] "sf"         "data.frame"
class(world_tbl)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"

第3章展示了如何使用tidyverse函数操作sf对象，sf现在是R中进行空间矢量数据分析的首选包。spatstat是一个提供众多空间统计函数的软件包生态系统，而terra也有矢量地理数据类，但它们在处理矢量数据方面都没有达到与sf相同的普及程度。许多受欢迎的包都基于sf构建，正如在上一章节中所示，其在每天的下载次数方面的流行程度正在上升。从传统包rgeos和rgdal迁移到已建立的包和工作流程需要时间，但当它们被加载时，打印出的信息指出它们“将在2023年底之前被退役”赋予了这个过程一种紧迫感。这意味着任何仍在使用这些包的人都应该“尽早转向使用 GDAL 和 PROJ 的 sf/stars/terra 函数”。换句话说，sf是对未来有保证的，而sp则不是。对于依赖传统类系统的工作流程，可以按照如下方式将sf对象从sp包的Spatial类中转换和转换为：

1
2
3
4
library(sp)
world_sp = as(world, "Spatial") # from an sf object to sp
# sp functions ...
world_sf = st_as_sf(world_sp)           # from sp to sf

绘制基础矢量地图

在sf中，使用plot()创建基础地图。默认情况下，这会创建一个多面板图，每个对象的变量各对应一个子图，如图2.4的左侧面板所示。如果要绘制的对象只有一个变量，将生成一个带有连续颜色的图例或“键”（见右侧面板）。颜色也可以通过col = 来设置，尽管这样做不会创建一个连续的调色板或图例。

1
2
plot(world[3:6])
plot(world["pop"])


图2.4：Plotting with sf, with multiple variables (left) and a single variable (right).

通过设置add = TRUE，图层会被添加到现有的图像上。^[在幕后，sf对象的plot()使用的是sf:::plot.sf()。plot() 是一个通用方法，其行为因要绘制的对象的类别而有所不同。]为了演示这一点，并提供关于属性和空间操作（涉及属性和空间数据操作）内容的洞见，以下代码块过滤出亚洲的国家，并将它们合并为一个单一的特征：

1
2
world_asia = world[world$continent == "Asia", ]
asia = st_union(world_asia)

现在我们可以在世界地图上绘制亚洲大陆。请注意，第一个图只能有一个面板，这样add = TRUE才能起作用。如果第一个图有一个图例，必须使用reset = FALSE：

1
2
plot(world["pop"], reset = FALSE)
plot(asia, add = TRUE, col = "red")


图2.5：A plot of Asia added as a layer on top of countries worldwide.

通过这种方式添加图层可以用于验证图层之间的地理对应关系：plot()函数执行速度快，需要的代码行数少，但不会创建具有广泛选项的交互式地图。对于更高级的地图制作，我们推荐使用专用的可视化包，如 tmap（见第9章）。

有多种方式可以通过sf的plot()方法修改地图。因为sf扩展了基础R绘图方法，plot()的参数（如main =，用于指定地图的标题）也适用于sf对象（参见?graphics::plot和?par）。^[注意：当绘制多个 sf 列时，许多绘图参数在面板图中会被忽略。]

图2.6通过在世界地图上叠加圆圈来展示这种灵活性，这些圆圈的直径（由cex =设置）代表各国的人口。这个图的未投影版本可以用以下命令创建（参见本章末尾的练习和脚本02-contplot.R，以复现图2.6：

1
2
3
4
plot(world["continent"], reset = FALSE)
cex = sqrt(world$pop) / 10000
world_cents = st_centroid(world, of_largest = TRUE)
plot(st_geometry(world_cents), add = TRUE, cex = cex)


图2.6：Country continents (represented by fill color) and 2015 populations (represented by circles, with area proportional to population).

上面的代码使用了函数st_centroid()来将一种几何类型（多边形）转换为另一种（点）（参见第5章），其美学效果通过cex参数进行变化。

sf的绘图方法也有针对地理数据的特定参数。例如，expandBB可以用于在上下文中绘制一个sf对象。它接受一个长度为四的数值向量，相对于零在以下顺序扩展绘图的边界框：底部、左侧、顶部、右侧。这用于在下面的代码块中绘制印度及其庞大的亚洲邻国的上下文，特别强调了东边的中国，该代码块生成了图（有关在绘图中添加文本的练习见下文）：^[注意使用st_geometry(india)仅返回与对象关联的几何形状，以防止在简单特征列（sfc）对象中绘制属性。另一种选择是使用india[0]，它返回一个不包含属性数据的sf对象。]

1
2
3
india = world[world$name_long == "India", ]
plot(st_geometry(india), expandBB = c(0, 0.2, 0.1, 1), col = "gray", lwd = 3)
plot(st_geometry(world_asia), add = TRUE)

India in context, demonstrating the expandBB argument.
图2.7：注意在绘图代码中使用lwd来强调印度。了解用于表示各种几何类型的其他可视化技术，这是下一节的主题。

几何类型

几何体是简单特征（Simple Features）的基础构建块。在R中的简单特征可以采用由sf包支持的18种几何类型之一。在本章中，我们将重点介绍七种最常用的类型：POINT、LINESTRING、POLYGON、MULTIPOINT、MULTILINESTRING、MULTIPOLYGON 和 GEOMETRYCOLLECTION。在PostGIS手册中找到可能的特征类型的完整列表。

通常，Well-Known Binary（WKB）或Well-Known Text（WKT）是简单特征几何体的标准编码。WKB表示通常是计算机容易读取的十六进制字符串。这就是为什么GIS和空间数据库使用WKB来传输和存储几何对象。另一方面，WKT是简单特征的人类可读文本标记描述。两种格式都是可交换的，如果我们展示其中一个，我们自然会选择WKT表示。

每种几何类型的基础是点。点只是2D、3D或4D空间中的一个坐标（有关更多信息，请参见vignette("sf1")）（见图2.8，左侧面板）：

  • POINT (5 2)

线串是一系列点与直线相连的序列（见图2.8，中间面板）。

  • LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2)

多边形是形成一个封闭、不相交环的点序列。封闭意味着多边形的第一个点和最后一个点具有相同的坐标（见图2.8，右侧面板）。

  • 没有孔的多边形：POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5))


图2.8：Illustration of point, linestring and polygon geometries.

到目前为止，我们创建的几何体每个特征只包含一个几何实体。sf也允许在单一特征内存在单一类型的多个几何体，即每种几何类型的“多”版本：

  • 多点：MULTIPOINT (5 2, 1 3, 3 4, 3 2)
  • 多线串：MULTILINESTRING ((1 5, 4 4, 4 1, 2 2, 3 2), (1 2, 2 4))
  • 多多边形：MULTIPOLYGON (((1 5, 2 2, 4 1, 4 4, 1 5), (0 2, 1 2, 1 3, 0 3, 0 2)))


图2.9：Illustration of multi* geometries.

最后，几何集合可以包含任何组合的几何体，包括（多）点和线串（见图2.10）：

  • 几何集合：GEOMETRYCOLLECTION (MULTIPOINT (5 2, 1 3, 3 4, 3 2), LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2))


图2.10：Illustration of a geometry collection.

sf类

简单特征由两个主要部分组成：几何体和非地理属性。图2.11展示了如何创建一个sf对象——几何体来自一个sfc对象，而属性则取自data.frame或tibble。要了解更多关于从零开始构建sf几何体的信息，请阅读以下第2.2.6节和2.2.7节。


图2.11：Building blocks of sf objects.

非地理属性代表特征的名称或其他属性，如测量值、组别和其他事项。为了说明属性，我们将表示2017年6月21日在伦敦的25°C温度。这个例子包含一个几何体（坐标），以及三个具有三个不同类别（地点名称、温度和日期）的属性。^[其他属性可能包括城市性质类别（城市或村庄），或者如果测量是通过自动站进行的，则包括一个备注。]sf类的对象通过将属性（data.frame）与简单特征几何列（sfc）组合来表示这样的数据。它们是通过下面所示的st_sf()创建的，该函数创建了上面描述的伦敦示例：

1
2
3
4
5
6
7
8
lnd_point = st_point(c(0.1, 51.5))                 # sfg object
lnd_geom = st_sfc(lnd_point, crs = "EPSG:4326")    # sfc object
lnd_attrib = data.frame(                           # data.frame object
  name = "London",
  temperature = 25,
  date = as.Date("2017-06-21")
  )
lnd_sf = st_sf(lnd_attrib, geometry = lnd_geom)    # sf object

发生了什么？首先，坐标被用来创建简单特征几何体（sfg）。其次，几何体被转换成一个带有坐标参考系统（CRS）的简单特征几何列（sfc）。第三，属性被存储在一个data.frame中，该data.frame与st_sf()一起与sfc对象组合。这导致了一个sf对象，如下面所示（省略了部分输出）：

1
2
3
4
5
lnd_sf
#> Simple feature collection with 1 features and 3 fields
#> ...
#>     name temperature       date         geometry
#> 1 London          25 2017-06-21 POINT (0.1 51.5)
1
2
class(lnd_sf)
#> [1] "sf"         "data.frame"

结果显示sf对象实际上有两个类别，sf和data.frame。简单特征实际上就是数据框（方形表格），但是具有存储在通常称为geometry的列表列中的空间属性（见2.2.1）。这种二元性是简单特征概念的核心：大多数时候，sf可以被视为并且表现得像一个data.frame。简单特征本质上是带有空间扩展的数据框。

简单特征几何 (sfg)

sfg类在R中代表不同的简单特征几何类型：点、线串、多边形（以及它们的“多”等效项，如多点）或几何集合。

通常，您不需要进行繁琐的任务来自己创建几何体，因为您可以简单地导入一个已经存在的空间文件。然而，如果需要，有一组函数可以从零开始创建简单特征几何对象（sfg）。这些函数的名称简单且一致，因为它们都以st_前缀开头，并以几何类型的小写字母名称结尾：

  • 一个点：st_point()
  • 一个线串：st_linestring()
  • 一个多边形：st_polygon()
  • 一个多点：st_multipoint()
  • 一个多线串：st_multilinestring()
  • 一个多多边形：st_multipolygon()
  • 一个几何集合：st_geometrycollection()

sfg对象可以从三种基础的R数据类型创建：

  1. 数值向量：一个单一点
  2. 矩阵：一组点，其中每一行代表一个点、一个多点或线串
  3. 列表：一组对象的集合，如矩阵、多线串或几何集合

函数st_point()从数值向量创建单一点：

1
2
3
4
5
6
7
8
st_point(c(5, 2))                 # XY point
#> POINT (5 2)
st_point(c(5, 2, 3))              # XYZ point
#> POINT Z (5 2 3)
st_point(c(5, 2, 1), dim = "XYM") # XYM point
#> POINT M (5 2 1)
st_point(c(5, 2, 3, 1))           # XYZM point
#> POINT ZM (5 2 3 1)

结果显示，XY（2D坐标）、XYZ（3D坐标）和XYZM（带有额外变量的3D，通常是测量精度）点类型分别从长度为2、3和4的向量创建。XYM类型必须使用dim参数（代表维度）来指定。

相比之下，在多点（st_multipoint()）和线串（st_linestring()）对象的情况下，使用矩阵：

1
2
3
4
5
6
7
8
9
# the rbind function simplifies the creation of matrices
## MULTIPOINT
multipoint_matrix = rbind(c(5, 2), c(1, 3), c(3, 4), c(3, 2))
st_multipoint(multipoint_matrix)
#> MULTIPOINT ((5 2), (1 3), (3 4), (3 2))
## LINESTRING
linestring_matrix = rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2))
st_linestring(linestring_matrix)
#> LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2)

最后，使用列表来创建多线串、（多）多边形和几何集合：

1
2
3
4
## POLYGON
polygon_list = list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5)))
st_polygon(polygon_list)
#> POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5))
1
2
3
4
5
6
## POLYGON with a hole
polygon_border = rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))
polygon_hole = rbind(c(2, 4), c(3, 4), c(3, 3), c(2, 3), c(2, 4))
polygon_with_hole_list = list(polygon_border, polygon_hole)
st_polygon(polygon_with_hole_list)
#> POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5), (2 4, 3 4, 3 3, 2 3, 2 4))
1
2
3
4
5
## MULTILINESTRING
multilinestring_list = list(rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2)), 
                            rbind(c(1, 2), c(2, 4)))
st_multilinestring((multilinestring_list))
#> MULTILINESTRING ((1 5, 4 4, 4 1, 2 2, 3 2), (1 2, 2 4))
1
2
3
4
5
## MULTIPOLYGON
multipolygon_list = list(list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))),
                         list(rbind(c(0, 2), c(1, 2), c(1, 3), c(0, 3), c(0, 2))))
st_multipolygon(multipolygon_list)
#> MULTIPOLYGON (((1 5, 2 2, 4 1, 4 4, 1 5)), ((0 2, 1 2, 1 3, 0 3, 0 2)))
1
2
3
4
5
6
## GEOMETRYCOLLECTION
geometrycollection_list = list(st_multipoint(multipoint_matrix),
                              st_linestring(linestring_matrix))
st_geometrycollection(geometrycollection_list)
#> GEOMETRYCOLLECTION (MULTIPOINT (5 2, 1 3, 3 4, 3 2),
#>   LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2))

简单特征列 (sfc)

一个sfg对象只包含一个单一的简单特征几何。简单特征几何列（sfc）是一个sfg对象的列表，此外还能包含有关正在使用的坐标参考系统的信息。例如，要将两个简单特征组合成一个具有两个特征的对象，我们可以使用st_sfc()函数。这一点很重要，因为sfc代表了sf数据框中的几何列：

1
2
3
4
5
6
7
8
9
10
11
12
# sfc POINT
point1 = st_point(c(5, 2))
point2 = st_point(c(1, 3))
points_sfc = st_sfc(point1, point2)
points_sfc
#> Geometry set for 2 features 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 1 ymin: 2 xmax: 5 ymax: 3
#> CRS:           NA
#> POINT (5 2)
#> POINT (1 3)

在大多数情况下，一个sfc对象包含相同几何类型的对象。因此，当我们将多边形类型的sfg对象转换为一个简单特征几何列时，我们也会得到一个多边形类型的sfc对象，这可以通过st_geometry_type()来验证。同样地，一个多线串的几何列将导致一个多线串类型的sfc对象：

1
2
3
4
5
6
7
8
9
# sfc POLYGON
polygon_list1 = list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5)))
polygon1 = st_polygon(polygon_list1)
polygon_list2 = list(rbind(c(0, 2), c(1, 2), c(1, 3), c(0, 3), c(0, 2)))
polygon2 = st_polygon(polygon_list2)
polygon_sfc = st_sfc(polygon1, polygon2)
st_geometry_type(polygon_sfc)
#> [1] POLYGON POLYGON
#> 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE
1
2
3
4
5
6
7
8
9
10
11
# sfc MULTILINESTRING
multilinestring_list1 = list(rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2)), 
                            rbind(c(1, 2), c(2, 4)))
multilinestring1 = st_multilinestring((multilinestring_list1))
multilinestring_list2 = list(rbind(c(2, 9), c(7, 9), c(5, 6), c(4, 7), c(2, 7)), 
                            rbind(c(1, 7), c(3, 8)))
multilinestring2 = st_multilinestring((multilinestring_list2))
multilinestring_sfc = st_sfc(multilinestring1, multilinestring2)
st_geometry_type(multilinestring_sfc)
#> [1] MULTILINESTRING MULTILINESTRING
#> 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

也可以从具有不同几何类型的sfg对象创建一个sfc对象：

1
2
3
4
5
# sfc GEOMETRY
point_multilinestring_sfc = st_sfc(point1, multilinestring1)
st_geometry_type(point_multilinestring_sfc)
#> [1] POINT           MULTILINESTRING
#> 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

如前所述，sfc对象还可以额外存储有关坐标参考系统（CRS）的信息。默认值是NA（Not Available，不可用），这可以通过st_crs()来验证：

1
2
st_crs(points_sfc)
#> Coordinate Reference System: NA

sfc对象中的所有几何体必须具有相同的CRS。可以使用st_sfc()（或st_sf()）的crs参数来指定CRS，该参数接受一个以文本字符串形式提供的CRS标识符，例如crs="EPSG:4326"。

1
2
3
4
5
6
7
8
# Set the CRS with an identifier referring to an 'EPSG' CRS code:
points_sfc_wgs = st_sfc(point1, point2, crs = "EPSG:4326")
st_crs(points_sfc_wgs) # print CRS (only first 4 lines of output shown)
#> Coordinate Reference System:
#>   User input: EPSG:4326 
#>   wkt:
#> GEOGCRS["WGS 84",
#> ...

sfheaders包

sfheaders是一个R包，用于加速构建、转换和操作sf对象。它专注于从向量、矩阵和数据框快速构建sf对象，而无需依赖sf库；并通过头文件（因此得名sfheaders）公开其底层的C++代码。这种方法使其他人能够使用编译和快速运行的代码来扩展它。每个核心的sfheaders函数都有一个相应的C++实现，如the Cpp vignette中所描述。对于大多数人来说，R函数将足以从该包的计算速度中受益。sfheaders是独立于sf开发的，但旨在完全兼容，创建有效的sf对象，其类型如前面的章节所述。

sfheaders最简单的用例在下面的代码块中有示例，其中展示了构建sfg、sfc和sf对象：

  • 向量转换为sfg_POINT
  • 矩阵转换为sfg_LINESTRING
  • 数据框转换为sfg_POLYGON

我们将从创建最简单可能的sfg对象开始，即一个单一的坐标对，赋值给名为v的向量：

1
2
v = c(1, 1)
v_sfg_sfh = sfheaders::sfg_point(obj = v)
1
2
3
4
5
v_sfg_sfh # printing without sf loaded
#>      [,1] [,2]
#> [1,]    1    1
#> attr(,"class")
#> [1] "XY"    "POINT" "sfg"

上面的示例展示了当sf未加载时，sfg对象v_sfg_sfh是如何打印的，展示了其底层结构。当sf被加载（如这里的情况）时，上述命令的结果与sf对象无法区分：

1
2
3
4
5
v_sfg_sf = st_point(v)
print(v_sfg_sf) == print(v_sfg_sfh)
#> POINT (1 1)
#> POINT (1 1)
#> [1] TRUE

下一个示例展示了sfheaders如何从矩阵和数据框创建sfg对象：

1
2
3
4
5
6
7
8
# matrices
m = matrix(1:8, ncol = 2)
sfheaders::sfg_linestring(obj = m)
#> LINESTRING (1 5, 2 6, 3 7, 4 8)
# data.frames
df = data.frame(x = 1:4, y = 4:1)
sfheaders::sfg_polygon(obj = df)
#> POLYGON ((1 4, 2 3, 3 2, 4 1, 1 4))

重用对象v、m和df，我们也可以如下构建简单特征列（sfc）（输出未显示）：

1
2
3
sfheaders::sfc_point(obj = v)
sfheaders::sfc_linestring(obj = m)
sfheaders::sfc_polygon(obj = df)

同样地，sf对象可以如下创建：

1
2
3
sfheaders::sf_point(obj = v)
sfheaders::sf_linestring(obj = m)
sfheaders::sf_polygon(obj = df)

在这些示例中，CRS（坐标参考系统）都没有定义。如果您计划使用sf函数进行任何计算或几何操作，我们建议您设置CRS（有关详细信息，请参见章节@ref(reproj-geo-data)）：

1
2
df_sf = sfheaders::sf_polygon(obj = df)
st_crs(df_sf) = "EPSG:4326"

sfheaders也擅长“解构”和“重构”sf对象，即将几何列转换为包含每个顶点和几何特征（以及多特征）ID的坐标数据的数据框。它在将几何列“转换”为不同类型方面非常快速和可靠。该包的文档和为本书开发的测试代码中的基准测试显示，它比sf包在这类操作上要快得多。

S2球面几何运算

球面几何引擎基于世界是圆形的这一事实，而用于地理计算的简单数学程序，如计算两点之间的直线或多边形所围成的面积，假设是平面（投影）几何。自sf版本1.0.0以来，R通过s2接口包支持“开箱即用”的球面几何操作，该接口连接到Google的S2球面几何引擎。S2可能最为人所知的是作为一个离散全球网格系统（DGGS）的例子。另一个例子是Uber的H3全球六边形分层空间索引。

尽管用于描述地球上任何地方的位置可能很有用，例如使用诸如e66ef376f790adf8a5af7fca9e6e422c03c9143f这样的字符串，但sf与S2的接口的主要优点是它提供了用于距离、缓冲区和面积计算等的即插即用函数，如在sf的内置文档中所描述，该文档可以通过命令vignette("sf7")打开。

sf可以在两种模式下运行，分别是开启和关闭S2。默认情况下，S2几何引擎是开启的，如下面的命令所验证：

1
2
sf_use_s2()
#> [1] TRUE

通过在本章早些时候创建的india对象周围创建缓冲区，下面展示了关闭几何引擎的后果（注意当S2关闭时发出的警告）：

1
2
3
4
5
6
7
india_buffer_with_s2 = st_buffer(india, 1)
sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off
india_buffer_without_s2 = st_buffer(india, 1)
#> Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle =
#> endCapStyle, : st_buffer does not correctly buffer longitude/latitude data
#> dist is assumed to be in decimal degrees (arc_degrees).


图2.12：Example of the consequences of turning off the S2 geometry engine. Both representations of a buffer around India were created with the same command but the purple polygon object was created with S2 switched on, resulting in a buffer of 1 m. The larger light green polygon was created with S2 switched off, resulting in a buffer with inaccurate units of degrees longitude/latitude.

在整本书中，我们将假设S2是开启的，除非明确说明。用以下命令再次开启它。

1
2
sf_use_s2(TRUE)
#> Spherical geometry (s2) switched on

📌尽管在许多情况下，sf使用S2是有道理的，但在某些情况下，有充分的理由在一个R会话或整个项目的过程中关闭S2。
如sf的GitHub仓库中的问题1771所记录，这种默认行为可能会导致在关闭S2（以及使用旧版本的sf）时可以运行的代码失败。
这些边缘情况包括对不符合S2更严格定义的多边形的操作。
如果你看到像#> Error in s2_geography_from_wkb ...这样的错误消息，关闭S2后再尝试生成错误消息的命令可能是值得的。
要在整个项目中关闭S2，你可以在项目的根目录（主文件夹）中创建一个名为.Rprofile的文件，其中包含命令sf::sf_use_s2(FALSE)。

栅格数据

空间栅格数据模型用连续的单元格网格（通常也称为像素；图2.13-A）来表示世界。这种数据模型通常指的是所谓的规则网格，其中每个单元格具有相同的、恒定的大小——在这本书中，我们只会关注规则网格。然而，还存在几种其他类型的网格，包括旋转的、剪切的、直线的和曲线的网格（参见Pebesma and Bivand（2022）的第1章或Tennekes and Nowosad (2022)的第2章）。

栅格数据模型通常由一个栅格头和一个矩阵（带有行和列）组成，该矩阵表示等间隔的单元格（通常也称为像素，图2.13-A）。^[根据文件格式，头部是实际图像数据文件的一部分，例如GeoTIFF，或存储在一个额外的头部或世界文件中，例如ASCII网格格式。还有无头（平面）二进制栅格格式，应便于导入到各种软件程序中。]栅格头定义了坐标参考系统、范围和原点。原点（或起点）通常是矩阵左下角的坐标（然而，terra包默认使用左上角，图2.13-B）。头通过列数、行数和单元格大小分辨率来定义范围。因此，从原点开始，我们可以通过使用单元格的ID或通过明确指定行和列来轻松访问和修改每个单个单元格（图2.13-B）。这种矩阵表示法避免了明确存储每个单元格角的四个角点的坐标（实际上它只存储一个坐标，即原点），就像矩形矢量多边形的情况那样。这一点和地图代数（见第4.3.2节）使栅格处理比矢量数据处理更高效、更快。然而，与矢量数据相比，一个栅格层的单元格只能容纳一个值。该值可能是数值或分类的（图2.13-C）。


图2.13：Raster data types: (A) cell IDs, (B) cell values, © a colored raster map.

栅格地图通常表示连续的现象，如海拔、温度、人口密度或光谱数据。离散特征，如土壤或土地覆盖类别，也可以在栅格数据模型中表示。栅格数据集的这两种用途都在图2.14中有所说明，该图显示了离散特征的边界在栅格数据集中可能会变得模糊。根据应用的性质，离散特征的矢量表示可能更为合适。


图2.14：Examples of continuous and categorical rasters.

处理栅格数据的R包

在过去的二十年里，已经开发了几个用于读取和处理栅格数据集的包。如第1.5节所示其中最主要的是raster，它在2010年推出后使R的栅格处理能力发生了质的飞跃，并一直是这一领域的主要包，直到terra和stars的开发。这两个较新开发的包都提供了强大和高效的函数，用于处理栅格数据集，并且它们可能的用例之间有很大的重叠。在本书中，我们将重点介绍terra，它取代了较旧且（在大多数情况下）较慢的raster。在了解terra的类系统如何工作之前，本节将描述terra和stars之间的相似性和差异；这些知识将有助于决定在不同情况下哪个更合适。

首先，terra专注于最常见的栅格数据模型（规则网格），而stars还允许存储较不流行的模型（包括规则的、旋转的、剪切的、直线的和曲线的网格）。虽然terra通常处理一个或多层的栅格^[它还有一个额外的类SpatRasterDataset，用于存储许多数据集的集合。]，stars包提供了存储栅格数据立方体的方法——一个具有多层（例如，波段）、多个时间点（例如，月份）和多个属性（例如，传感器类型A和传感器类型B）的栅格对象。重要的是，在这两个包中，所有层或数据立方体的元素都必须具有相同的空间尺寸和范围。其次，这两个包都允许将所有栅格数据读入内存，或者只读取其元数据——这通常是根据输入文件大小自动完成的。然而，它们存储栅格值的方式非常不同。terra基于C++代码，主要使用C++指针。stars将值存储为较小栅格的数组列表或较大栅格的文件路径。第三，stars的函数与sf中的矢量对象和函数密切相关，而terra使用其自己的矢量数据对象类，即SpatVector。第四，这两个包对各种函数如何作用于它们的对象有不同的方法。terra包主要依赖于大量的内置函数，其中每个函数都有一个特定的目的（例如，重新采样或裁剪）。另一方面，stars使用一些内置函数（通常名称以st_开头），有其自己的现有R函数（例如，split()或aggregate()）的方法，以及现有的dplyr函数（例如，filter()或slice()）的方法。

重要的是，将对象从terra转换为stars（使用st_as_stars()）以及反过来（使用rast()）都是非常简单的。我们还鼓励您阅读@pebesma_spatial_2022，以获得对stars包最全面的介绍。

terra包介绍

terra包在R中支持栅格对象。与其前身raster（由同一开发者Robert Hijmans创建）一样，它提供了一整套广泛的函数，用于创建、读取、导出、操作和处理栅格数据集terra的功能在很大程度上与更成熟的raster包相同，但也有一些差异：terra的函数通常比raster的等效函数更具计算效率。另一方面，raster的类系统很受欢迎，并被许多其他包使用。您可以无缝地在这两种类型的对象之间进行转换，以确保与旧脚本和包的向后兼容性，例如，使用raster包中的raster()、stack()和brick()函数（有关使用地理数据的R包演变的更多信息，请参见前一章）。

除了用于栅格数据操作的函数外，terra还提供了许多低级函数，这些函数可以作为开发用于处理栅格数据集的新工具的基础。terra还允许您处理过大以至于无法装入主内存的大型栅格数据集。在这种情况下，terra提供了将栅格划分为更小块的可能性，并迭代地处理这些块，而不是将整个栅格文件加载到RAM中。

为了说明terra的概念，我们将使用spDataLarge中的数据集。它包括几个栅格对象和一个覆盖锡安国家公园（美国犹他州）地区的矢量对象。例如，srtm.tif是该地区的数字高程模型（有关更多详细信息，请参见其文档?srtm）。首先，让我们创建一个名为my_rast的SpatRaster对象：

1
2
3
4
5
6
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
my_rast = rast(raster_filepath)
class(my_rast)
#> [1] "SpatRaster"
#> attr(,"package")
#> [1] "terra"

在控制台中输入栅格的名称，将打印出栅格头部信息（尺寸、分辨率、范围、坐标参考系统）以及一些额外信息（类别、数据源、栅格值的摘要）。

1
2
3
4
5
6
7
8
9
10
my_rast
#> class       : SpatRaster 
#> dimensions  : 457, 465, 1  (nrow, ncol, nlyr)
#> resolution  : 0.000833, 0.000833  (x, y)
#> extent      : -113, -113, 37.1, 37.5  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : srtm.tif 
#> name        : srtm 
#> min value   : 1024 
#> max value   : 2892

专用函数报告每个组件：dim()返回行数、列数和图层数；ncell()返回单元格（像素）数；res()返回空间分辨率；ext()返回其空间范围；crs()返回其坐标参考系统。inMemory()报告栅格数据是存储在内存中还是磁盘上。

help("terra-package")返回所有可用的terra函数的完整列表。

绘制基础栅格地图

与sf包类似，terra也为其自己的类提供了plot()方法。

1
plot(my_rast)


图2.15：Basic raster plot.

还有几种在R中绘制栅格数据的方法，这些方法不在本节的讨论范围内，包括：

  • 来自terra包的plotRGB()函数，用于基于SpatRaster对象中的三个图层创建红-绿-蓝图。
  • 诸如tmap之类的包，用于创建栅格和矢量对象的静态和交互式地图（参见第9章）。
  • 函数，例如来自rasterVis包的levelplot()，用于创建分面图，这是一种常用于可视化随时间变化的技术。

栅格类

SpatRaster类在terra中代表栅格对象。在R中创建栅格对象的最简单方法是从磁盘或服务器中读入一个栅格文件（第8.6.2章）。

1
2
single_raster_file = system.file("raster/srtm.tif", package = "spDataLarge")
single_rast = rast(raster_filepath)

terra 包在 GDAL 库的帮助下支持众多的驱动程序。文件中的栅格通常不会完全读入RAM，除了它们的头部和指向文件本身的指针。

也可以使用相同的 rast() 函数从零开始创建栅格。这在后续的代码块中有所展示，该代码块生成一个新的 SpatRaster 对象。生成的栅格包含36个单元格（由 nrows 和 ncols 指定的6列和6行），围绕本初子午线和赤道居中（参见 xmin、xmax、ymin 和 ymax 参数）。值（vals）被分配给每个单元格：1分配给单元格1，2分配给单元格2，依此类推。记住：rast() 按行填充单元格（与 matrix() 不同），从左上角开始，这意味着顶行包含值1至6，第二行包含值7至12，等等。有关创建栅格对象的其他方法，请参见 ?rast。

1
2
3
new_raster = rast(nrows = 6, ncols = 6, 
                  xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5,
                  vals = 1:36)

给定行数和列数以及范围（xmin、xmax、ymin、ymax），分辨率必须为0.5。分辨率的单位是底层 CRS 的单位。在这里，它是度数，因为栅格对象的默认 CRS 是 WGS84。然而，可以用 crs 参数指定任何其他 CRS。

SpatRaster 类还处理多个图层，这些图层通常对应于一个单一的多光谱卫星文件或一系列的栅格。

1
2
3
4
5
6
7
8
9
10
11
12
multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
multi_rast = rast(multi_raster_file)
multi_rast
#> class       : SpatRaster 
#> dimensions  : 1428, 1128, 4  (nrow, ncol, nlyr)
#> resolution  : 30, 30  (x, y)
#> extent      : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 12N (EPSG:32612) 
#> source      : landsat.tif 
#> names       : landsat_1, landsat_2, landsat_3, landsat_4 
#> min values  :      7550,      6404,      5678,      5252 
#> max values  :     19071,     22051,     25780,     31961

nlyr()用于获取存储在SpatRaster对象中的图层数量：

1
2
nlyr(multi_rast)
#> [1] 4

对于多层栅格对象，可以使用terra::subset()来选择图层。^[[[ 和 $ 操作符也可以用来选择图层，例如通过命令multi_rast$landsat_1和multi_rast[["landsat_1"]]。]它接受图层编号或其名称作为第二个参数：

1
2
multi_rast3 = subset(multi_rast, 3)
multi_rast4 = subset(multi_rast, "landsat_4")

相反的操作，即将多个SpatRaster对象合并为一个，可以使用c函数来完成：

1
multi_rast34 = c(multi_rast3, multi_rast4)

📌大多数SpatRaster对象并不存储栅格值，而是存储指向文件本身的指针。
这有一个显著的副作用——它们不能直接保存为".rds"或".rda"文件，也不能用于集群计算。
在这些情况下，有两种可能的解决方案：（1）使用wrap()函数，该函数创建一种特殊类型的临时对象，该对象可以作为R对象保存或用于集群计算；（2）使用writeRaster()将对象保存为常规栅格。

地理、投影坐标参考系

矢量和栅格空间数据类型共享空间数据本质上的概念。其中最基础的概念便是坐标参照系统（CRS），它定义了数据的空间元素与地球（或其他天体）表面的关系。CRS分为地理或投影两种，正如本章开头所介绍。本节解释每种类型，为后面章节奠定基础，该章将深入讨论设置、转换和查询CRS。

地理坐标系

地理坐标系统使用两个值——经度和纬度来标识地球表面上的任何位置（见图2.17）。经度 是从本初子午线平面到东西方向的角距离。纬度 是从赤道平面向北或向南的角距离。因此，地理CRS中的距离不是以米为单位度量的这有重要的影响，如第7章所示。

在地理坐标系统中，地球表面由球体或椭球体表面表示。球体模型假设地球是一个给定半径的完美球体——它们具有简单性的优点，但同时也是不准确的：地球不是一个球体！椭球体模型由两个参数定义：赤道半径和极半径。这些是合适的，因为地球是压缩的：赤道半径比极半径长约11.5公里。^[压缩的程度通常被称为 扁平度，用赤道半径（$a$）和极半径（$b$）来定义：$f = (a - b) / a$。术语 椭圆度 和 压缩 也可以使用。由于 $f$ 是一个相当小的值，数字椭球体模型使用 ‘逆扁平度’（$rf = 1/f$）来定义地球的压缩。通过执行 sf_proj_info(type = "ellps")，可以看到各种椭球体模型中的 $a$ 和 $rf$ 的值。]

椭球体是CRSs的一个更广泛组成部分：基准。这包含有关使用哪个椭球体以及笛卡尔坐标和地球表面位置之间精确关系的信息。有两种类型的基准——地心（如 WGS84）和局部（如 NAD83）。图2.16中可以看到这两种基准的示例。黑线代表一个 地心基准，其中心位于地球重力中心，没有针对特定位置进行优化。在 局部基准 中，如紫色虚线所示，椭球体表面被移动以与特定位置的表面对齐。这允许局部CRS考虑地球表面的局部变化，例如由于大山脉造成的。这可以在图中看到，其中局部基准适合于菲律宾地区，但与地球表面的大部分其他区域不对齐。图2.16中的两个基准都放置在一个大地水准面（geoid）之上——一个全球平均海平面的模型。^[请注意，图中的大地水准面通过放大因子10,000来夸大了大地水准面的凹凸不平表面，以突出地球的不规则形状。]


图2.16：地心和局部大地基准显示在大地水准面之上（以错误的颜色和垂直夸张的10,000倍比例因子）。大地水准面的图像改编自 @essd-11-647-2019 的作品。

投影坐标系

所有的投影坐标系统（Projected CRSs）都基于前一节中描述的地理坐标系统，并依赖于地图投影将地球的三维表面转换为投影坐标系统中的东坐标和北坐标（x和y值）。投影坐标系统基于笛卡尔坐标，在一个隐含的平面表面上（见图2.17的右侧面板）。它们有一个原点、x和y轴，以及一个如米这样的线性度量单位。

这种转换不能在不添加一些形变的情况下完成。因此，在这个过程中，地球表面的一些属性会被扭曲，比如面积、方向、距离和形状。一个投影坐标系统只能保留这些属性中的一个或两个。投影通常根据它们保留的属性来命名：等面积投影保留面积，方位投影保留方向，等距投影保留距离，和保角投影保留局部形状。

有三大主要的投影类型——圆锥投影、圆柱投影和平面投影（方位投影）。在一个圆锥投影中，地球表面被投影到一个沿着一个或两个切线的圆锥上。在这种投影中，扭曲是沿着切线最小化的，并且随着与这些线的距离增加而增加。因此，它最适合用于中纬度地区的地图。圆柱投影将表面映射到一个圆柱上。这种投影也可以通过沿着一个或两个切线与地球表面接触来创建。圆柱投影最常用于整个世界的制图。平面投影将数据投影到一个与地球在一个点或沿一个切线接触的平面上。它通常用于极地地区的制图。使用sf_proj_info(type = "proj")可以获得PROJ库支持的可用投影的列表。

关于不同投影、它们的类型、属性和适用性的快速概要可以在 “Map Projections”（1993）和 https://www.geo-projections.com/ 上找到。我们将在后面章节中详细讨论坐标参照系统（CRSs）并解释如何从一个CRS投影到另一个CRS（见第7章）。目前，足够了解：

  • 坐标系统是地理对象的一个关键组成部分
  • 知道你的数据处于哪种CRS中，以及它是在地理坐标（经/纬度）还是投影坐标（通常是米）中，对于R如何处理空间和几何操作是重要的
  • sf 对象的CRS可以使用 st_crs() 函数查询，terra 对象的CRS可以使用 crs() 函数查询。


图2.17：Examples of geographic (WGS 84; left) and projected (NAD83 / UTM zone 12N; right) coordinate systems for a vector data type.

单位

坐标参考系统（CRS）的一个重要特点是它们包含有关空间单位的信息。很明显，知道一个房屋的测量单位是英尺还是米是至关重要的，地图也同样如此。在地图上添加比例尺或其他距离指示器是良好的制图实践，以展示页面或屏幕上的距离与地面上的距离之间的关系。同样，正式指定几何数据或单元格以何种单位进行测量也是重要的，以提供上下文，并确保后续的计算是在这一上下文中进行的。

sf对象中的几何数据有一个新颖的特点，即它们对单位有原生支持。这意味着在sf中进行的距离、面积和其他几何计算返回的值附带有一个由units软件包定义的units属性。这是有利的，因为它防止了由于不同单位（大多数CRS使用米，一些使用英尺）引起的混淆，并提供了有关尺寸的信息。下面的代码块演示了如何计算卢森堡的面积：

1
luxembourg = world[world$name_long == "Luxembourg", ]
1
2
st_area(luxembourg) # requires the s2 package in recent versions of sf
#> 2.41e+09 [m^2]

输出的单位是平方米（m^2^），显示该结果表示二维空间。这些信息存储为一个属性（有兴趣的读者可以通过attributes(st_area(luxembourg))来发现这一点），可以用于后续使用单位的计算，例如人口密度（通常以每平方千米的人数来衡量）。报告单位可以防止混淆。以卢森堡为例，如果单位保持未指定，人们可能错误地假设单位是公顷。为了将这个巨大的数字转换为更易消化的大小，人们可能会受到将结果除以一百万（一个平方千米中的平方米数量）的诱惑：

1
2
st_area(luxembourg) / 1000000
#> 2409 [m^2]

然而，结果再次错误地以平方米为单位给出。解决方案是使用 units 包来设置正确的单位：

1
2
units::set_units(st_area(luxembourg), km^2)
#> 2409 [km^2]

在栅格数据的情况下，单位同样非常重要。然而，到目前为止，sf 是唯一支持单位的空间包，这意味着处理栅格数据的人应该谨慎地处理单位变换（例如，将像素宽度从英制单位转换为十进制单位）。my_rast 对象（见上文）使用的是WGS84投影，单位是十进制度。因此，其分辨率也以十进制度给出，但你需要知道它，因为 res() 函数仅返回一个数值向量。

1
2
res(my_rast)
#> [1] 0.000833 0.000833

如果我们使用UTM（通用横轴墨卡托）投影，单位会发生变化。

1
2
3
repr = project(my_rast, "EPSG:26912")
res(repr)
#> [1] 83.5 83.5

再次强调，res() 命令返回的是一个没有任何单位的数值向量，这就迫使我们必须知道UTM投影的单位是米。

练习

E1. 在spData包中包含的world数据对象的几何列上使用summary()。输出告诉我们关于：

  • 它的几何类型是什么？
  • 有多少个国家？
  • 它的坐标参考系统（CRS）是什么？

E2. 运行在2.2.3节（基础地图制作）中“生成”世界地图的代码。
找出你电脑上的图像与书中的图像之间的两个相似点和两个不同点。

  • cex参数是做什么用的（参见?plot）？
  • 为什么cex被设置为sqrt(world$pop) / 10000？
  • 奖励：尝试用不同的方式可视化全球人口。

E3. 使用plot()创建尼日利亚的上下文地图（参见2.2.3节）。

  • 调整plot()的lwd、col和expandBB参数。
  • 挑战：阅读text()的文档并标注地图。

E4. 创建一个名为my_raster的空SpatRaster对象，具有10列和10行。
为新的光栅分配0到10之间的随机值并绘制它。

E5. 从spDataLarge包中读取raster/nlcd.tif文件。
你能从这个文件的属性中获取哪些信息？

E6. 检查来自spDataLarge包的raster/nlcd.tif文件的CRS。
你能从中了解到哪些信息？

-------------已经到底啦-------------