(6)栅格矢量交叉

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

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

本章专注于栅格和矢量地理数据模型之间的相互作用。包括四个主要技术：使用矢量对象进行栅格裁剪和遮罩；使用不同类型的矢量数据提取栅格值；以及栅格与矢量之间的转换。以上概念使用前几章中使用的数据进行演示，以了解其潜在的现实应用。

前提条件

  • 本章要求下列包：
1
2
3
library(sf)
library(terra)
library(dplyr)

引言

本章专注于栅格和矢量地理数据模型之间的相互作用。包括四个主要技术：使用矢量对象进行栅格裁剪和遮罩；使用不同类型的矢量数据提取栅格值；以及栅格与矢量之间的转换。以上概念使用前几章中使用的数据进行演示，以了解其潜在的现实应用。

栅格裁剪

许多地理数据项目涉及从许多不同来源整合数据，例如遥感图像（栅格）和行政边界（矢量）。通常，输入栅格数据集的范围大于感兴趣的区域。在这种情况下，栅格裁剪和掩膜对于统一输入数据的空间范围非常有用。这两个操作减少了对象的内存使用和后续分析步骤的相关计算资源，而且是创建涉及栅格数据地图之前必要的预处理步骤。

我们将使用两个对象来说明栅格裁剪：

  • 一个代表犹他州西南部海拔高度（海平面以上米数）的SpatRaster对象srtm
  • 一个代表锡安国家公园的矢量（sf）对象 zion

目标对象和裁剪对象必须具有相同的投影。因此，下面的代码块不仅从spDataLarge包中读取了章节@ref(spatial-class) 中安装的数据集，还对zion进行了“重投影”。

1
2
3
srtm = rast(system.file("raster/srtm.tif", package = "spDataLarge"))
zion = read_sf(system.file("vector/zion.gpkg", package = "spDataLarge"))
zion = st_transform(zion, crs(srtm))

我们使用terra包的crop()函数来裁剪srtm栅格。该函数根据传递给其第二个参数的对象的范围减小传递给其第一个参数的对象的矩形范围。下面的命令演示了这个功能，生成了图：

1
srtm_cropped = crop(srtm, zion)

与crop()相关的是terra函数mask()，它将传递给第二个参数的对象范围之外的值设置为NA。因此，以下命令会掩膜锡安国家公园边界以外的每一个单元格：

1
srtm_masked = mask(srtm, zion)

重要的是，在大多数情况下，我们想要同时使用crop()和mask() 这两个函数。这组函数的组合将(a)将栅格的范围限制在我们感兴趣的区域内，然后 (b) 将该区域外的所有值替换为NA。

1
2
srtm_cropped = crop(srtm, zion)
srtm_final = mask(srtm_cropped, zion)

更改mask()的设置会产生不同的结果。设置inverse = TRUE将会屏蔽公园范围内的所有内容（详见?mask），而设置 updatevalue = 0 将会将国家公园外的所有像素设置为0。

1
srtm_inv_masked = mask(srtm, zion, inverse = TRUE)
1
2
3
4
5
6
7
8
#> 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)


Illustration of raster cropping and raster masking.

栅格提取

栅格提取是识别并返回基于（通常是矢量）地理“选择器”对象在特定位置的“目标”栅格关联的值的过程。结果取决于所使用的选择器的类型（点、线或多边形）以及传递给terra::extract()函数的参数。栅格提取的反向过程 — 基于矢量对象分配栅格单元格值 — 是栅格化，详见 @ref(rasterization) 部分。

基本示例是在特定点处提取栅格单元格的值。为此，我们将使用 zion_points，其中包含锡安国家公园内的30个位置样本（见图 @ref(fig:pointextr)）。以下命令从srtm中提取高程值，并创建一个包含点的ID（每个矢量的一行一个值）和每个点的相关srtm值的数据框。现在，我们可以使用cbind()函数将结果对象添加到我们的zion_points数据集中：

1
2
3
data("zion_points", package = "spDataLarge")
elevation = terra::extract(srtm, zion_points)
zion_points = cbind(zion_points, elevation)
Locations of points used for raster extraction.

(\#fig:pointextr)Locations of points used for raster extraction.

栅格提取也适用于线选择器。提取线所接触的每个栅格单元格的一个值。然而，线提取方法不推荐用于沿截面获取值，因为很难获得每一对提取的栅格值之间的正确距离。

在这种情况下，更好的方法是将线分割成许多点，然后为这些点提取值。为了演示这一点，下面的代码创建了 zion_transect，一条从锡安国家公园的西北部到东南部的直线，如图 @ref(fig:lineextr)(A)所示（有关矢量数据模型的回顾，请参见 @ref(vector-data) 部分）：

1
2
3
4
zion_transect = cbind(c(-113.2, -112.9), c(37.45, 37.2)) |>
  st_linestring() |> 
  st_sfc(crs = crs(srtm)) |>
  st_sf(geometry = _)

规划徒步旅行时，从线性选择器中提取高度的实用性可见一斑。下面演示的方法提供了路线的“海拔剖面图”（线不需要是直的），这对于估计由于长时间攀爬而需要的时间非常有用。

第一步是为每个截面添加唯一的 id。接下来，我们可以使用 st_segmentize()函数在我们的线（或线）上用提供的密度（dfMaxLength）添加点，并使用 st_cast() 将它们转换为点。

1
2
3
zion_transect$id = 1:nrow(zion_transect)
zion_transect = st_segmentize(zion_transect, dfMaxLength = 250)
zion_transect = st_cast(zion_transect, "POINT")

现在，我们有了一大组点，我们想要推导出我们的截面中第一个点与后续各点之间的距离。在这种情况下，我们只有一个截面，但原则上，代码应该适用于任何数量的截面：

1
2
3
zion_transect = zion_transect |> 
  group_by(id) |> 
  mutate(dist = st_distance(geometry)[, 1])

最后，我们可以提取截面中每个点的高程值，并将此信息与我们的主要对象结合起来。

1
2
zion_elev = terra::extract(srtm, zion_transect)
zion_transect = cbind(zion_transect, zion_elev)

所得到的zion_transect可以用来创建高程剖面图，如图所示。


Location of a line used for raster extraction (left) and the elevation along this line (right).

栅格提取的最后一种地理矢量对象是多边形。与线条一样，多边形每个多边形倾向于返回许多栅格值。以下命令演示了这一点，其结果是一个数据框，列名为ID（多边形的行号）和srtm（相关的海拔值）：

1
zion_srtm_values = terra::extract(x = srtm, y = zion)

这样的结果可以用来生成每个多边形的栅格值的汇总统计信息，例如用来描述单个区域或比较多个区域。下面的代码展示了这一点，它创建了对象zion_srtm_df，包含锡安国家公园内的海拔值的汇总统计信息：

1
2
3
4
5
6
group_by(zion_srtm_values, ID) |> 
  summarize(across(srtm, list(min = min, mean = mean, max = max)))
#> # A tibble: 1 × 4
#>      ID srtm_min srtm_mean srtm_max
#>   <dbl>    <int>     <dbl>    <int>
#> 1     1     1122     1818.     2661

前述代码段使用了dplyr来提供每个多边形ID的单元格值的汇总统计信息。结果提供了有用的总结，例如公园内的最大高度约为海拔2661米（也可以用这种方式计算其他汇总统计信息，如标准差）。因为示例中只有一个多边形，所以返回了一个单行的数据框；然而，当使用多个选择多边形时，该方法同样适用。

相似的方法适用于在多边形内计数分类栅格值的出现次数。这通过一个土地覆盖数据集（nlcd）在上图中进行说明，并在下面的代码中演示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
nlcd = rast(system.file("raster/nlcd.tif", package = "spDataLarge"))
zion2 = st_transform(zion, st_crs(nlcd))
zion_nlcd = terra::extract(nlcd, zion2)
zion_nlcd |> 
  group_by(ID, levels) |>
  count()
#> # A tibble: 7 × 3
#> # Groups:   ID, levels [7]
#>      ID levels         n
#>   <dbl> <fct>      <int>
#> 1     1 Developed   4205
#> 2     1 Barren     98285
#> 3     1 Forest    298299
#> 4     1 Shrubland 203701
#> # ℹ 3 more rows


Area used for continuous (left) and categorical (right) raster extraction.

尽管terra包提供了在多边形内快速提取栅格值的功能，extract()在处理大型多边形数据集时仍可能成为瓶颈。exactextractr包通过exact_extract()函数提供了一个明显更快的替代方案用于提取像素值。exact_extract()函数还默认计算了每个栅格单元被多边形重叠的部分，这更为精确（有关详细信息，请参见下面的注释）。

📌多边形通常具有不规则的形状，因此，多边形可能只与栅格的某些部分重叠。
为了获得更详细的结果，terra::extract() 函数有一个叫做 exact 的参数。
通过设置 exact = TRUE，我们在输出的数据框中多得到一列fraction，代表每个单元格被多边形覆盖的部分。
这可以用于计算例如连续栅格的加权平均值或者分类栅格的更精确覆盖率。
默认情况下，它的值为 FALSE，因为这个操作需要更多的计算。
exactextractr::exact_extract() 函数始终计算多边形在每个单元格中的覆盖部分。

栅格化

栅格化是将矢量对象转换为栅格对象表示的过程。通常，输出栅格随后用于定量分析（例如，地形分析）或建模。栅格数据模型具有一些特性，使其有助于某些方法。此外，栅格化过程可以帮助简化数据集，因为所得的值都具有相同的空间分辨率：栅格化可以看作是地理数据聚合的一种特殊类型。

terra 包中包含了用于执行此工作的 rasterize() 函数。其前两个参数是 x（要栅格化的矢量对象）和y（定义输出范围、分辨率和坐标参考系统的“模板栅格”对象）。输入栅格的地理分辨率对结果有重大影响：如果分辨率过低（单元格大小过大），结果可能会错过矢量数据的全部地理变化；如果分辨率过高，则计算时间可能过长。在决定适当的地理分辨率时没有简单的规则可遵循，这在很大程度上取决于结果的预期用途。例如，当栅格化的输出需要与某些其他现有栅格对齐时，目标分辨率通常由用户决定。

为了演示栅格化的实际效果，我们将使用一个模板栅格，其范围和坐标参考系统与输入矢量数据cycle_hire_osm_projected相同，并且空间分辨率为1000米：

1
2
3
4
cycle_hire_osm = spData::cycle_hire_osm
cycle_hire_osm_projected = st_transform(cycle_hire_osm, "EPSG:27700")
raster_template = rast(ext(cycle_hire_osm_projected), resolution = 1000,
                       crs = st_crs(cycle_hire_osm_projected)$wkt)

栅格化是一项非常灵活的操作：结果不仅取决于模板栅格的性质，还取决于输入矢量的类型（例如，点、多边形）以及 rasterize() 函数所采用的各种参数。

为了说明这种灵活性，我们将尝试三种不同的栅格化方法。首先，我们创建一个表示自行车租赁点存在或不存在的栅格（称为存在/不存在栅格）。在这种情况下，rasterize() 除了上述的矢量和栅格对象 x 和 y 之外，不需要任何参数。

1
ch_raster1 = rasterize(cycle_hire_osm_projected, raster_template)

fun参数用于指定将紧密相邻的多个观测值转换为栅格对象中的关联单元格所使用的汇总统计信息。默认情况下使用 fun = "last"，但也可以使用其他选项，例如 fun = "length"，在这种情况下，用于计算每个网格单元格中自行车租赁点的数量。

1
2
ch_raster2 = rasterize(cycle_hire_osm_projected, raster_template, 
                       fun = "length")

新的输出结果，ch_raster2，展示了每个网格单元中的自行车租赁点数量。自行车租赁地点拥有不同数量的自行车，由capacity变量描述，这引发了一个问题：每个网格单元的容量是多少？为了计算这个值，我们必须对字段（"capacity"）进行求和，从而得到图所示的输出结果（也可以使用其他汇总函数，例如 mean）。

1
2
ch_raster3 = rasterize(cycle_hire_osm_projected, raster_template, 
                       field = "capacity", fun = sum, na.rm = TRUE)


Examples of point rasterization.

另一个基于加利福尼亚州的多边形和边界的数据集（如下所示）用于说明线的栅格化。在将多边形对象转换为多段线后，创建了一个具有0.5度分辨率的模板栅格。

1
2
3
4
california = dplyr::filter(us_states, NAME == "California")
california_borders = st_cast(california, "MULTILINESTRING")
raster_template2 = rast(ext(california), resolution = 0.5,
                        crs = st_crs(california)$wkt)

在考虑线或多边形栅格化时，一个有用的附加参数是touches。默认情况下，它为FALSE，但当更改为TRUE时，所有被线或多边形边界触及的单元格都会得到一个值。下面的代码展示了使用touches = TRUE的线栅格化。

1
2
california_raster1 = rasterize(california_borders, raster_template2,
                               touches = TRUE)

与多边形栅格化相比，默认情况下touches = FALSE，只选择其质心在选择多边形内的栅格单元格，如图所示。

1
california_raster2 = rasterize(california, raster_template2)


Examples of line and polygon rasterizations.

空间矢量化

空间矢量化是栅格化的对应部分，但方向相反。它涉及将空间连续的栅格数据转换为空间离散的矢量数据，例如点、线或多边形。

📌确实，我们应当小心措辞以避免混淆！
在R语言中，矢量化通常指的是通过类似1:10/2这样的操作来替代for循环等结构的可能性（详见@wickham_advanced_2019）。所以，当谈到将栅格数据转换为矢量数据时，我们应避免使用“矢量化”一词，以免与R中的这一特定概念混淆。

矢量化的最简单形式是将栅格单元格的中心点转换为点。as.points()正是针对所有非NA的栅格网格单元格进行此操作。注意，这里我们还使用了st_as_sf()来将结果对象转换为sf类。

1
2
3
elev = rast(system.file("raster/elev.tif", package = "spData"))
elev_point = as.points(elev) |> 
  st_as_sf()


Raster and point representation of the elev object.

另一种常见的空间矢量化类型是创建代表连续高度或温度（例如等温线）的等高线。我们将使用现实世界的数字高程模型（DEM），因为人造栅格elev会产生平行线（读者的任务：验证并解释为什么会发生这种情况）。可以使用 terra 函数 as.contour()创建等高线，该函数本身是R内置函数filled.contour()的封装，如下所示（未显示）：

1
2
3
4
5
dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))
cl = as.contour(dem) |> 
  st_as_sf()
plot(dem, axes = FALSE)
plot(cl, add = TRUE)

等高线也可以通过诸如contour()、rasterVis::contourplot()或tmap::tm_iso()等函数添加到现有的绘图中。如图所示，等值线可以加标签。


DEM with hillshading, showing the southern flank of Mt. Mongón overlaid with contour lines.

最后一种矢量化涉及将栅格转换为多边形。这可以通过 terra::as.polygons()来完成，该函数将每个栅格单元格转换为由五个坐标组成的多边形，所有这些都存储在内存中（这解释了为什么与矢量相比，栅格通常更快！）。

下面通过将grain对象转换为多边形，并随后消除具有相同属性值的多边形之间的边界（还可以参见as.polygons()中的dissolve参数）来进行说明。

1
2
3
grain = rast(system.file("raster/grain.tif", package = "spData"))
grain_poly = as.polygons(grain) |> 
  st_as_sf()


Illustration of vectorization of raster (left) into polygons (dissolve = FALSE; center) and aggregated polygons (dissolve = TRUE; right).

grain数据集的聚合多边形具有直线边界，这是由于通过连接矩形像素来定义所产生的。smoothr包可用于平滑多边形的边缘。由于平滑会消除多边形边界中的锐角，平滑后的多边形将不具有与原始像素完全相同的空间覆盖（有关示例，请参见 smoothrwebsite ）。因此，在使用平滑后的多边形进行进一步分析时应谨慎。

练习

以下部分练习使用来自 spDataLarge 包的矢量（zion_points）和栅格数据集（srtm）。它们还使用从矢量数据集（ch）派生的多边形 ‘凸包’ 来表示感兴趣的区域:

1
2
3
4
5
6
7
8
9
library(sf)
library(terra)
library(spData)
zion_points_path = system.file("vector/zion_points.gpkg", package = "spDataLarge")
zion_points = read_sf(zion_points_path)
srtm = rast(system.file("raster/srtm.tif", package = "spDataLarge"))
ch = st_combine(zion_points) |>
  st_convex_hull() |> 
  st_as_sf()

E1. 根据（1）zion_points数据集和（2）ch数据集裁剪srtm栅格。输出地图之间有任何区别吗？接下来，使用这两个数据集对srtm进行遮罩处理。现在能看到任何区别吗？你如何解释这种情况？

  • 对srtm栅格进行裁剪和遮罩处理可能会产生不同的结果，因为zion_points和ch表示的区域可能有所不同。裁剪操作基于数据集的边界，而遮罩操作则基于数据集的几何形状。如果两个数据集的几何形状不同，结果也会不同。

E2. 首先，从srtm中提取zion_points中表示的点的值。接下来，使用zion_points中每个点周围的90缓冲区提取srtm的平均值，并比较这两组值。什么时候通过缓冲区提取值比仅通过点提取值更合适？

  • 当需要考虑点周围的区域或当点可能不准确表示感兴趣的区域时，使用缓冲区提取值可能更合适。例如，在考虑地形或其他空间变量时，使用缓冲区可能会提供更多的上下文信息。

E3. 提取新西兰（nz_height对象）海拔高度超过3100米的点，并为新点数据集的范围创建一个分辨率为3 km的模板栅格。使用这两个新对象：

  • 计算每个网格单元中最高点的数量。
  • 找到每个网格单元中的最大海拔。

E4. 聚合统计新西兰高点的栅格（在前一个练习中创建），将其地理分辨率减半（因此单元格为6 x 6 km），并绘制结果。

  • 将低分辨率的栅格重新采样回原始的3 km分辨率。结果有何变化？
  • 降低栅格分辨率的两个优点和两个缺点是什么？

E5. 将grain数据集多边形化，并过滤所有代表粘土的方块。

  • 列举矢量数据相对于栅格数据的两个优点和两个缺点。

  • 在你的工作中，什么时候将栅格转换为矢量会很有用？

  • 矢量数据通常提供更高的精度和更清晰的边界表示，同时允许复杂的几何和拓扑关系。但它们可能需要更多的存储空间和计算时间，特别是对于大型或复杂的数据集。

  • 将栅格转换为矢量可能在需要精确边界、复杂的空间分析或与其他矢量数据集交互时有用。

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