(01) 欢迎（空间数据科学——R语言应用）

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

本文由SCY翻译自《Spatial Data Science With Applications in R》第一章

本章介绍了一系列与处理空间和时空数据相关的概念，并指向后续章节，其中这些概念将会被更详细地讨论。它还介绍了一系列构成所有空间数据科学语言实现基础的开源技术。

第一张地图

绘制空间数据的典型方法是创建一张地图。让我们考虑一张简单的地图，如下图所示。

1
2
3
4
5
6
7
8
library(tidyverse)
library(sf)
system.file("gpkg/nc.gpkg", package="sf") |>
	read_sf() -> nc
nc.32119 <- st_transform(nc, 'EPSG:32119')
nc.32119 |>
	select(BIR74) |>
	plot(graticule = TRUE, axes = TRUE)

图1.1：第一张地图：1974-78年北卡罗来纳州各县出生人数

在这个例子中，有多个图形元素存在：

  • 多边形用黑色轮廓绘制，并根据一个名为BIR74的变量用颜色进行填充，该变量的名字出现在标题中。
  • 图例键解释了颜色的含义，并具有特定的色彩调色板和色彩断点，也就是颜色变化的值。
  • 地图的背景显示了具有恒定纬度或经度的曲线（经纬网）。
  • 坐标轴刻度显示了纬度和经度的值。

多边形是一种特殊形式的几何体；空间几何体（点、线、多边形、像素）将在几何章节中详细讨论。多边形由点的序列组成，通过直线连接。空间数据的点位置如何表示或测量，在坐标系章节中有讨论。从上图可以看出，相等的纬度和经度并没有形成直线，这表明在绘图之前进行了某种投影；投影也在坐标系章节和投影章节中讨论。

在上图中的颜色值来自一个变量BIR74的数值，该变量与每个几何体或特征有一个单一的关联值。属性和支持章节讨论了这样的特征属性，以及它们如何与特征几何体相关联。在这个例子中，BIR74是出生人数，即该区域的人数。这意味着这个计数并不是指多边形内每个点关联的一个值，如连续颜色可能暗示的那样，而是测量了多边形上的积分（总和）。

在绘制上图之前，我们需要读取数据，在这个例子中是从一个文件中读取的（sf包）。打印前三条记录的三个属性变量的数据摘要显示：

1
nc |> select(AREA, BIR74, SID74) |> print(n = 3)
1
2
3
4
5
6
7
8
9
10
11
12
## Simple feature collection with 100 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS:  NAD27
## # A tibble: 100 × 4
##    AREA BIR74 SID74                                                                              geom
##   <dbl> <dbl> <dbl>                                                                <MULTIPOLYGON [°]>
## 1 0.114  1091     1 (((-81.47276 36.23436, -81.54084 36.27251, -81.56198 36.27359, -81.63306 36.3406…
## 2 0.061   487     0 (((-81.23989 36.36536, -81.24069 36.37942, -81.26284 36.40504, -81.26624 36.4372…
## 3 0.143  3188     5 (((-80.45634 36.24256, -80.47639 36.25473, -80.53688 36.25674, -80.54501 36.2766…
## # ℹ 97 more rows

打印输出结果如下所示:

  • （选择）的数据集包含 100 个要素（条记录）和 3 个字段（属性）。
  • 几何类型是 MULTIPOLYGON （详见几何章节）。
  • 其维度为XY，表明每个点由2个坐标值组成。
  • 几何的$x$和$y$值的范围。
  • 坐标参考系统（CRS）是大地坐标系统，其坐标以经度和纬度表示，与NAD27基准相关联（坐标系章节）。
  • 三个选定的属性变量后面是一个名为geom的MULTIPOLYGON类型变量，单位为度，包含多边形信息

更复杂的图表可能涉及到每个方面都有一张地图的分面图，如下图所示。

1
2
3
4
5
6
7
8
year_labels <- c("SID74" = "1974 - 1978", "SID79" = "1979 - 1984")
nc.32119 |> select(SID74, SID79) |>
	pivot_longer(starts_with("SID")) -> nc_longer
ggplot() + geom_sf(data = nc_longer, aes(fill = value), linewidth = 0.4) + 
  facet_wrap(~ name, ncol = 1, labeller = labeller(name = year_labels)) +
  scale_y_continuous(breaks = 34:36) +
  scale_fill_gradientn(colors = sf.colors(20)) +
  theme(panel.grid.major = element_line(color = "white"))

图1.2：1974-78年和1979-84年北卡罗来纳州各县突发性婴儿死亡综合症计数的分面地图

::: {.content-visible when-format=“html”}

leaflet 交互式地图如下所示。

1
2
3
library(mapview) |> suppressPackageStartupMessages()
mapviewOptions(fgb = FALSE)
nc.32119 |> mapview(zcol = "BIR74", legend = TRUE, col.regions = sf.colors)
1
## PhantomJS未找到。你可以通过`webshot::install_phantomjs()`进行安装。如果已经安装，请确保PhantomJS可执行文件可以通过PATH变量找到。

图1.3：使用mapview创建的交互式地图：拖动和缩放可移动地图并更改比例；点击一个县将弹出一个窗口，显示可用的县属性。

坐标参考系

在第一张图中，灰色线表示经纬网，即沿着恒定的纬度或经度的网格。显然，这些线并不是直的，这表明数据使用了一种$x$和$y$轴与经度和纬度不对齐的投影。在**交互式地图中，我们看到北卡罗来纳州的北界再次被绘制成一条直线，这表明使用了另一种投影。

第一张图的椭球坐标与特定的基准（这里是：NAD27）关联，这涉及到一组规则，即地球的形状是什么，以及它如何与地球相连接（与地球的哪个点关联的原点，以及如何指向）。如果有人用GPS设备（如手机）测量坐标，它通常会报告与世界大地测量系统1984年（WGS84）基准关联的坐标，与北美大地测量基准1927年（NAD27）关联的相同坐标值相比，可能有约30米的差距。

投影描述了我们如何在以下两者之间来回转换：

  • 椭球坐标，以经度和纬度的度数表示，指向一个近似于地球形状（椭球或球体）的位置，以及

  • 投影坐标，这是在绘制地图时使用的平面的二维坐标系统。

基准转换与从一个基准转移到另一个基准有关。两个主题都由_空间参考系统_涵盖，并在坐标系章节中更详细地描述。

栅格数据和矢量数据

多边形、点和线几何体是矢量数据的例子：点坐标描述了可以位于任何地方的"精确"位置。另一方面，栅格数据描述了值对齐在一个栅格上，意味着在一个通常是方形像素的规则排列的点阵上。一个例子显示在下图中。

1
2
3
4
5
6
7
8
9
10
11
12
13
library(stars)
par(mfrow = c(2, 2))
par(mar = rep(1, 4))
tif <- system.file("tif/L7_ETMs.tif", package = "stars")
x <- read_stars(tif)[,,,1]
image(x, main = "(a)")
image(x[,1:10,1:10], text_values = TRUE, border = 'grey', main = "(b)")
image(x, main = "(c)")
set.seed(131)
pts <- st_sample(st_as_sfc(st_bbox(x)), 3)
plot(pts, add = TRUE, pch = 3, col = 'blue')
image(x, main = "(d)")
plot(st_buffer(pts, 500), add = TRUE, pch = 3, border = 'blue', col = NA, lwd = 2)

图1.4：巴西大西洋海岸Olinda的栅格地图：Landsat-7蓝色波段，颜色值来源于数据值（a），来自（a）的左上角$10 imes 10$子图像，显示数值（b），并由两种不同类型的矢量数据叠加：三个样本点（c），以及以点为中心的500米半径的多边形（d）

矢量和栅格数据可以以不同的方式组合；例如，我们可以在上图©的三个点处查询栅格，或者计算一个聚合值，比如在上图(d)中显示的圆形区域内的平均值。

1
st_extract(x, pts) # query at points
1
2
3
4
5
6
7
8
9
## Simple feature collection with 3 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 290019 ymin: 9114499 xmax: 291693.3 ymax: 9119219
## Projected CRS: SIRGAS 2000 / UTM zone 25S
##   L7_ETMs.tif                 geometry
## 1          80 POINT (290829.6 9114499)
## 2          58   POINT (290019 9119219)
## 3          63 POINT (291693.3 9116038)
1
aggregate(x, st_buffer(pts, 500), FUN = mean) |> st_as_sf() # aggregate over circles
1
2
3
4
5
6
7
8
9
## Simple feature collection with 3 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 289519 ymin: 9113999 xmax: 292193.3 ymax: 9119719
## Projected CRS: SIRGAS 2000 / UTM zone 25S
##         V1                       geometry
## 1 77.18201 POLYGON ((291329.6 9114499,...
## 2 60.12603 POLYGON ((290519 9119219, 2...
## 3 71.57808 POLYGON ((292193.3 9116038,...

在栅格矢量转化中讨论了其他栅格到矢量的转换，包括：

  • 将栅格像素转换为点值
  • 将栅格像素转换为小多边形，可能合并具有相同值的多边形（“多边形化”）
  • 生成描绘具有某一值范围的连续像素区域的线或多边形（“等高线”）
1
plot(st_rasterize(nc["BIR74"], dx = 0.1), col = sf.colors(), breaks = "equal")

图1.5：通过栅格化1974-1978年期间县出生计数而获得的地图，显示在1.1中

矢量到栅格的转换可以像栅格化多边形一样简单，如上图所示。其他更一般的涉及统计建模的矢量到栅格转换包括：

  • 将点值插值到规则网格上的点（空间差值章节）
  • 估计规则网格上点的密度（点格局分析章节）
  • 将多边形值通过面积加权插值到网格单元（面积加权差值）
  • 点、线或多边形的直接栅格化（栅格矢量转化）

栅格类型

栅格维度描述了行和列如何与空间坐标相关联。下图展示了几种不同的可能性。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
x <- 1:5
y <- 1:4
d <- st_dimensions(x = x, y = y, .raster = c("x", "y"))
m <- matrix(runif(20),5,4)
r1 <- st_as_stars(r = m, dimensions = d)

r <- attr(d, "raster")
r$affine <- c(0.2, -0.2)
attr(d, "raster") = r
r2 <- st_as_stars(r = m, dimensions = d)

r <- attr(d, "raster")
r$affine <- c(0.1, -0.3)
attr(d, "raster") = r
r3 = st_as_stars(r = m, dimensions = d)

x <- c(1, 2, 3.5, 5, 6)
y <- c(1, 1.5, 3, 3.5)
d <- st_dimensions(x = x, y = y, .raster = c("x", "y"))
r4 <- st_as_stars(r = m, dimensions = d)

grd <- st_make_grid(cellsize = c(10,10), offset = c(-130,10), n = c(8,5), crs = st_crs('OGC:CRS84'))
r5 <- st_transform(grd, "+proj=laea +lon_0=-70 +lat_0=35")

par(mfrow = c(2,3), mar = c(0.1, 1, 1.1, 1))
r1 <- st_make_grid(cellsize = c(1,1), n = c(5,4), offset = c(0,0))
plot(r1, main = "regular")
plot(st_geometry(st_as_sf(r2)), main = "rotated")
plot(st_geometry(st_as_sf(r3)), main = "sheared")
plot(st_geometry(st_as_sf(r4, as_points = FALSE)), main = "rectilinear")
plot(st_geometry((r5)), main = "curvilinear")

图1.6：多种栅格几何类型

像上图中展示的规则栅格具有恒定的、不一定是正方形的单元格大小，以及与$x$和$y$（东经和北纬）轴对齐。其他类型的栅格包括那些轴不再与$x$和$y$对齐（旋转）、轴不再垂直（倾斜）或者单元格大小沿一个维度变化（矩形）的情况。最后，曲线栅格具有单元格大小或方向属性，这些属性不再独立于其他栅格维度。

当一个在给定坐标参考系统中是规则的栅格投影到另一个栅格上，同时保持每个栅格单元格完好无损时，它会改变形状，并可能变成矩形（例如，从椭球坐标到墨卡托投影，如图3所示）或曲线（例如，从椭球坐标到兰伯特圆锥保角投影，如在图1中使用）。当撤销这一过程时，人们可以恢复到精确的原始栅格。

在新的投影中创建一个新的、规则的网格被称为栅格（或图像）重新投影或者扭曲（变换和扭曲栅格）。扭曲是有损的、不可逆的，并且需要指明是否应对栅格单元格进行插值、平均或求和，或者是否应使用最近邻居进行重新采样。对于这种选择，单元格值是反映分类变量还是连续变量是有关系的（另见支持章节）。

时间序列、数组、数据立方体

许多空间数据不仅仅是空间性的，还具有时间性。就像任何观察都与一个观察地点相关，同样也与一个观察时间或周期相关。上面显示的关于北卡罗来纳州各县的数据集包含了在两个时间段内计数的疾病案例（如图2所示）。尽管原始数据集中的这些变量分布在两个不同的列中，但为了绘制它们，这些列首先必须被堆叠起来，同时重复关联的几何形状——这种形式被称为整洁（tidy）。当我们有与几何形状关联的更长的时间序列时，两种选项——在多个列上分布时间，或者在重复几何形状的同时堆叠列——都不是很理想，更有效的存储这种数据的方式将是一个矩阵或数组，其中一个维度指代时间，其他维度指代空间。对于图像或栅格数据来说，自然的存储方式已经是矩阵；然后，栅格的时间序列会导致一个三维数组。这种数据的一般术语是（时空的）数据立方体（data cube），其中立方体指代具有任意数量维度的数组。数据立方体可以指代栅格和矢量数据，例子在数据立方体章节中给出。

支持

当我们拥有的空间数据的几何形状不是点，而是点的集合（多点、线、多边形、像素）时，与这些几何形状关联的属性可能有几种不同的关系。一个属性可以具有：

  • 几何形状每个点的常量值
  • 几何形状所有点的**汇总（aggregate）**单一值
  • 仅用于描述此几何形状身份的唯一值

一个常量的例子是多边形的土地利用或基岩类型。汇总的一个例子是在给定时间段内某个县的出生人数。身份的一个例子是县的名称。

属性值所指的空间区域被称为其支持（support）：汇总属性具有“块”（或区域、或线）支持，常数属性具有“点”支持（它们适用于每个点）。在操作数据时，支持很重要。例如，图5源自具有多边形支持的变量：每个县的出生人数。栅格化这些值会得到仍与各县相关联的像素值。栅格化的结果是一个无意义的地图：数值值（“出生总数”）与栅格单元格没有关联，县界也不再存在。全州的出生总数或出生密度无法从像素值中恢复。忽略支持很容易导致无意义的结果。属性和支持章节进一步讨论了这一点。

栅格单元格值可能具有点支持或块支持。点支持的一个例子是高程，当单元格在数字高程模型中记录单元格中心点的高程时。块（或单元格）支持的一个例子是卫星图像像素，该像素给出了（与像素相似区域的）平均颜色值。大多数文件格式并不提供这种信息，但在聚合、重新网格化或扭曲栅格，提取点位置的值时，了解这一点可能很重要。

空间数据科学软件

尽管本书在空间数据科学方面主要使用R语言和R包，但这些包中的许多使用了并非专为R开发的软件库。例如，R包sf对其他R包和系统库的依赖如下图所示。

图1.7：sf 包及其依赖；箭头表示强依赖，虚线箭头表示弱依赖。

使用的C或C++库（GDAL、GEOS、PROJ、liblwgeom、s2geometry、NetCDF、udunits2）都是由（空间）数据科学社群开发、维护和使用的，这些社群大都与R社群不同。通过使用这些库，R用户与这些其他社群共享我们对所做工作的理解。因为R、Python和Julia为这些软件提供交互式接口，许多用户比使用基于这些库的其他软件的用户更接近这些库。本书的第一部分描述了这些库中实现的许多概念，这些概念对空间数据科学有一般的相关性。

GDAL

GDAL（地理空间数据抽象库）可以看作是空间数据的瑞士军刀；除了在R中使用外，还用于Python、QGIS、PostGIS以及100多个其他软件项目。

GDAL是一个“库中的库”——为了读取和写入这些数据，它需要大量的其他库。它通常链接到超过100个其他库，每个库可能提供对特定数据文件格式、数据库、网络服务或特定压缩编解码器的访问。

由CRAN分发的二进制R包仅包含静态链接的代码：CRAN不想对主机系统上第三方库的存在做任何假设。因此，当从CRAN以二进制形式安装sf包时，它包括所有所需的外部库以及它们的依赖项，总量可能达到100 Mb。

PROJ

PROJ（或PR$\phi$J）是一个用于地图投影和基准转换的库：它将空间坐标从一个坐标参考系统转换到另一个坐标参考系统。它配备了大量已知投影的数据库，并可访问基准网格（用于基准转换的高精度、预先计算的值）。它符合坐标参考系统的国际标准。坐标系章节处理坐标系统和PROJ。

GEOS and s2geometry

GEOS（Geometry Engine Open Source）和s2geometry是两个用于几何操作的库。它们用于找到测量值（长度、面积、距离），并计算谓词（两个几何体是否有任何共同点？）或新几何体（这两个几何体有哪些共同点？）。GEOS针对的是平面的二维空间（由$R^2$表示），而s2geometry针对的是球面上的几何体（由$S^2$表示）。@sec-cs介绍了坐标参考系统，而球状几何章节更多地讨论了在这两个空间中工作的区别。

NetCDF, udunits2, liblwgeom

NetCDF 既指一种文件格式，也指用于读取和写入NetCDF文件的C库。它允许定义任何维度的数组，并广泛用于空间和时空信息，特别是在（气候）建模社群中。Udunits2 是一个单位测量的数据库和软件库，允许单位转换，处理派生单位，并支持用户定义的单位。liblwgeom“库”是PostGIS 的一个软件组件，其中包含了GDAL或GEOS缺失的若干例程，包括方便地访问随PROJ一同提供的GeographicLib例程。

练习

  1. 列举五个raster（栅格）和vector（矢量）数据之间的差异；

  2. 除了图1下面列出的以外，常见于地图上的五个其他图形组件；

  3. 为什么在图5中显示的数字信息是误导性的（或无意义的）；

  4. 在哪些情况下，你期望在进行$S^2$上的几何操作时与在$R^2$上进行它们时有很大的不同？

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