(4)空间操作

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

空间操作，包括矢量数据集之间的空间连接以及栅格数据集上的局部和焦点操作，是地理计算的重要部分。本章展示了如何基于空间对象的位置和形状以多种方式进行修改。许多空间操作有与之对应的非空间（属性）操作，因此上一章中演示的如子集提取和数据集连接等概念在此也适用。

前提条件

• 此章节需要使用和之前章节相同的包。

library(sf)
library(terra)
library(dplyr)
library(spData)

引言

空间操作，包括矢量数据集之间的空间连接以及栅格数据集上的局部和焦点操作，是地理计算的重要部分。本章展示了如何基于空间对象的位置和形状以多种方式进行修改。许多空间操作有与之对应的非空间（属性）操作，因此上一章中演示的如子集提取和数据集连接等概念在此也适用。对于矢量操作来说，这一点尤为适用：矢量属性操作为理解空间对应关系——即空间子集（在空间矢量提取部分中介绍）提供了基础。空间连接和空间聚合也具有非空间相似操作，已在上一章中介绍。

空间操作与非空间操作在许多方面有所不同，例如：空间连接可以通过多种方式进行，包括匹配与目标数据集相交或位于目标数据集一定距离范围内的实体，而上一章节矢量属性连接中讨论的属性连接只能通过一种方式进行（除非使用模糊连接，如 fuzzyjoin 包的文档所述）。对象之间不同类型的空间关系，包括相交和分离，均在章节拓扑关系中描述。空间对象的另一个独特之处在于距离，所有的空间对象都通过空间相关联，距离计算可用于探索这种关联的强度，正如在距离关系节中描述的矢量数据的情况一样。

栅格对象上的空间操作包括子集提取——在空间栅格子集提取节中进行了介绍-以及将多个栅格"瓦片"合并成一个对象，在合并栅格节中演示。地图代数涵盖了一系列修改栅格单元值的操作，可以参考周围单元值，也可以不参考周围单元值。地图代数的概念对许多应用至关重要，本节介绍了地图代数，并分别在地图代数节介绍了局部、焦点和分区地图代数操作。全局地图代数操作会生成代表整个栅格数据集的汇总统计数据，与栅格数据的距离计算在全局操作和距离章节中进行了讨论。在章节中，讨论了合并两个栅格数据集的过程，并结合一个可重现的示例进行演示。

📌需要注意的是，使用两个空间对象进行空间操作需要这两个对象具有相同的坐标参考系统。这个话题在crs介绍中引入，并在地理数据重投影中进行了更详细的讨论。

矢量数据空间操作

本节概述了在sf包中表示为简单要素的矢量地理数据上的空间操作。栅格数据空间操作节使用 terra 包中的类和函数对栅格数据集进行空间操作。

空间矢量提取

空间子集提取是将空间对象进行处理并返回一个包含与另一个对象空间相关特征的新对象的过程。类似于属性子集提取（在矢量属性提取子集）节中介绍），可使用方括号（[）运算符来创建 sf 数据框的子集，其语法为x[y, , op = st_intersects]，其中x是 sf 对象的一个子集，y是"子集对象"，op = st_intersects是一个可选参数，用于指定进行子集提取处理的拓扑关系（也称为二元谓词）。当未提供op参数时，st_intersects()是默认的拓扑关系。命令x[y, ]与上述x[y, , op = st_intersects]完全相同，但不同于x[y, , op = st_disjoint]（这些以及其他拓扑关系的含义在下一节中描述）。tidyverse中的filter()函数也可以使用，但这种方法更冗长，如下例所示。

为了演示空间子集提取，我们将使用 spData 包中的nz和nz_height数据集，它们分别包含新西兰16个主要地区和101个最高点的地理数据，以投影坐标系为基础。以下代码块创建的对象表示Canterbury，然后使用空间子集提取返回该区域内所有的高点：

canterbury = nz |> filter(Name == "Canterbury")
canterbury_height = nz_height[canterbury, ]

就像属性子集提取一样，命令x[y, ]（相当于nz_height[canterbury, ]）使用源对象y的内容对目标x进行特征子集提取。然而，与y是逻辑或整数类的向量不同的是，对于空间子集提取，x和y都必须是地理对象。具体而言，以这种方式用于空间子集提取的对象必须具有类sf或sfc。nz和nz_height都是地理向量数据框，具有类sf，操作的结果将返回另一个sf对象，表示目标nz_height对象中与（在本例中位于）canterbury地区相交的特征（即位于高处的点）。

用于空间子集提取的各种拓扑关系决定了目标对象中的特征，必须与要选择的子集对象具有的空间关系类型。这些关系包括接触、交叉或包含，在拓扑关系部分中我们将很快看到。默认设置st_intersects是一个"全包括"的拓扑关系，它将返回与源"子集"对象接触、交叉或包含的目标中的特征。如上所示，可以用op=参数指定其他空间运算符，如下面的命令所示，该命令返回st_intersects()的相反内容，即与坎特伯雷不相交的点（请参阅部分）：

nz_height[canterbury, , op = st_disjoint]

📌请注意前面代码块中的空参数用 , , 表示是为了突出 op，即 sf 对象 [ 的第三个参数。
可以使用这个参数以多种方式改变子集操作。
例如，nz_height[canterbury, 2, op = st_disjoint] 返回相同的行，但仅包括第二个属性列（详见 sf:::`[.sf` 和 ?sf）。

对于许多应用而言，关于矢量数据的空间子集提取的知识就是你需要了解的全部内容：它只会按预期工作。如果你急于了解更多拓扑关系，超出了 st_intersects()和st_disjoint()的范围，请跳到下一节拓扑关系。如果你对细节感兴趣，包括其他子集方式，请继续阅读。

另一种进行空间子集提取的方法是使用拓扑操作返回的对象。这些对象本身就很有用，例如在探索相邻区域之间的关系图网络时，但它们也可以用于子集操作，如下方的代码块中所演示的：

sel_sgbp = st_intersects(x = nz_height, y = canterbury)
class(sel_sgbp)
#> [1] "sgbp" "list"
sel_sgbp
#> Sparse geometry binary predicate list of length 101, where the
#> predicate was `intersects'
#> first 10 elements:
#>  1: (empty)
#>  2: (empty)
#>  3: (empty)
#>  4: (empty)
#>  5: 1
#>  6: 1
#>  7: 1
#>  8: 1
#>  9: 1
#>  10: 1
sel_logical = lengths(sel_sgbp) > 0
canterbury_height2 = nz_height[sel_logical, ]

以上代码块创建一个 sgbp 类的对象 (一个稀疏几何二元谓词，一个在空间操作中长度为 x 的列表)，然后将其转换为逻辑向量 sel_logical（包含仅为 TRUE 和 FALSE 的值，这也可以由 dplyr 的 filter 函数使用）。lengths()识别nz_height中的哪些要素与 y 中的任何对象相交。在这种情况下，1 是最大可能值，但对于更复杂的操作，可以使用该方法仅对与源对象中的 2 个或更多个要素相交的要素进行子集操作。

📌注意：通过在st_intersects()等运算符中设置sparse = FALSE（意思是“返回密集矩阵而不是稀疏矩阵”）也可以返回逻辑输出。例如，命令st_intersects(x = nz_height, y = canterbury, sparse = FALSE)[, 1]将返回与sel_logical相同的输出。
注意：涉及sgbp对象的解决方案更具普遍性，因为它适用于多对多的操作并且对内存要求较低。

使用 sf 函数st_filter()可以实现相同的结果，该函数是为了增加sf对象与dplyr数据操作代码之间的兼容性而创建的：

canterbury_height3 = nz_height |>
  st_filter(y = canterbury, .predicate = st_intersects)

到目前为止，有三个完全相同（除了行名）的canterbury_height版本，一个使用[操作符创建，一个通过中间选择对象创建，另一个使用sf的便捷函数st_filter()创建。

下一节探讨了不同类型的空间关系，也称为二元谓词，可以用来确定两个特征是否存在空间关系。

拓扑关系

拓扑关系描述了对象之间的空间关系。完整称呼为"二元拓扑关系"，是关于由有序点集（通常形成点、线和多边形）在两个或多个维度中定义的两个对象之间的空间关系的逻辑陈述（答案只能是TRUE或FALSE）。这听起来可能相当抽象，实际上，拓扑关系的定义和分类是基于1966年首次以书籍形式出版的数学基础，代数拓扑领域一直延续到21世纪。

尽管拓扑关系起源于数学，但通过参考常用函数的可视化，可以直观地理解用于测试常见空间关系类型的拓扑关系。显示了各种几何对及其关联关系。图中的第三和第四对（从左到右然后向下）表明，对于某些关系，顺序很重要：虽然等于、相交、交叉、接触和重叠的关系是对称的，意味着如果function(x, y)为真，则function(y, x)也为真，但包含和内部等顺序重要的几何关系则不是。注意，每一对几何图形都有一个"DE-9IM"字符串，例如FF2F11212，将在下一节中描述。

在sf中，测试不同拓扑关系的函数被称为"二元谓词"，如在操作简单要素几何的文献 Manipulating Simple Feature Geometries 中所述，可以使用命令 vignette(“sf3”)查看，也可以在帮助页面中查看?geos_binary_pred。为了更好地理解拓扑关系的实用性，我们将建立一个简单可重现的例子，基于上图中所述的关系，巩固了前一章节几何所介绍的矢量几何图形表示的知识。请注意，为了创建代表多边形顶点坐标（x 和 y）的表格数据，我们使用基本 R 函数 cbind() 创建表示坐标点的矩阵，一个POLYGON，最后是 sfc 对象，如空间类所述：

polygon_matrix = cbind(
  x = c(0, 0, 1, 1,   0),
  y = c(0, 1, 1, 0.5, 0)
)
polygon_sfc = st_sfc(st_polygon(list(polygon_matrix)))

请注意在上面创建的多边形之上，我们将创建额外的几何体来展示它们在空间中的关系。所示的命令在绘制时与该多边形相关联。请注意在转换数据框时使用了函数st_as_sf()和参数coords，以高效地将包含坐标列的数据框转换为包含点的sf对象：

line_sfc = st_sfc(st_linestring(cbind(
  x = c(0.4, 1),
  y = c(0.2, 0.5)
)))
# 创建点
point_df = data.frame(
  x = c(0.2, 0.7, 0.4),
  y = c(0.1, 0.2, 0.8)
)
point_sf = st_as_sf(point_df, coords = c("x", "y"))

一个简单查询是：point_sf 中的哪些点与多边形 polygon_sfc 以某种方式相交？此问题通过检查可得到答案（点1和点3分别与多边形相切和在多边形内）。可以使用空间谓词 st_intersects()，如下所示：

st_intersects(point_sf, polygon_sfc)
#> Sparse geometry binary predicate... `intersects'
#>  1: 1
#>  2: (empty)
#>  3: 1

结果应该符合你的直觉：第一个和第三个点返回正（1）的结果，第二个点在多边形的边界之外返回负结果（用一个空向量表示）。而令人意想不到的是，结果以向量列表的形式呈现。这个稀疏矩阵输出只记录存在关系的部分，减少了对多要素对象进行拓扑操作时的内存需求。正如我们在前面的部分中所看到的，当sparse = FALSE时，返回的是一个由TRUE或FALSE值组成的稠密矩阵：

st_intersects(point_sf, polygon_sfc, sparse = FALSE)
#>       [,1]
#> [1,]  TRUE
#> [2,] FALSE
#> [3,]  TRUE

在上述输出中，每一行代表目标（参数x）对象中的一个特征，每一列代表选择对象（y）中的一个特征。在这个情况下，y对象polygon_sfc中只有一个特征，所以结果只有一列，这个结果可以用于我们在空间子集提取部分所看到的子集提取。

st_intersects()即使在特征只是相切的情况下也返回TRUE，intersects是一种"全能"的拓扑操作，它识别许多类型的空间关系，如下图所示。更有限制的问题包括哪些点位于多边形内，以及哪些特征在y上或包含与y共享的边界？这些问题可以如下回答（结果未显示）：

st_within(point_sf, polygon_sfc)    # 在范围内
st_touches(point_sf, polygon_sfc)   # 相切

请注意，尽管第一个点接触多边形的边界，但它并不在其中；第三个点在多边形内部，但不接触其边界的任何部分。st_intersects()的反义词是st_disjoint()，它只返回与选择对象在空间上完全不相关的对象（注意[, 1]将结果转换为向量）：

st_disjoint(point_sf, polygon_sfc, sparse = FALSE)[, 1]
#> [1] FALSE  TRUE FALSE

函数st_is_within_distance()检测到那些几乎接触选择对象的特征，它有一个额外的dist参数。它可以用来设置目标对象需要多近才能被选择。请注意，尽管点2距离polygon_sfc的最近顶点的距离超过0.2个单位，但当距离设置为0.2时，它仍然被选择。这是因为距离是测量到最近的边缘，在这种情况下是多边形直接在点2上方的部分，如下图所示。（你可以用命令st_distance(point_sf, polygon_sfc)来验证点2和多边形之间的实际距离是0.13。）下面的代码块演示了’is within distance’二元空间谓词，其结果显示每个点都在多边形的0.2个单位内：

st_is_within_distance(point_sf, polygon_sfc, dist = 0.2, sparse = FALSE)[, 1]
#> [1] TRUE TRUE TRUE

📌算拓扑关系的函数使用空间索引来大大提高空间查询性能。
它们使用Sort-Tile-Recursive（STR）算法来实现。
下一节中提到的st_join函数也使用空间索引。
您可以在https://www.r-spatial.org/r/2017/06/22/spatial-index.html中了解更多信息。

DE-9IM strings

在前一节所示的二元谓词之下，隐含的是Dimensionally Extended 9-Intersection Model (DE-9IM)。正如这个神秘的名字所暗示的那样，这并不是一个容易的话题。不过，学习它有可能更好地理解空间关系。此外，DE-9IM的高级用途还包括创建自定义空间谓词。该模型最初被其发明者标记为"DE + 9IM"，指的是"两个要素的边界、内部和外部的交界面的维度" ，但现在被称为DE-9IM 。

为了演示DE-9IM strings的工作原理，让我们看看中第一个几何对之间的各种关系。下图展示了9 intersection model (9IM)，显示了每个对象的内部、边界和外部之间的交点，当第一个对象x的每个组件被排列为列，而y的每个组件被排列为行时，会创建出一个带有每个元素交点突出显示的分面图形。

DE-9IM strings是基于每种关系的不同维度而产生的。在这种情况下，图中的红色交点分别有着0（点）、1（线）和2（多边形）不同的维度，如简单表格所示。

	Interior (x)	Boundary (x)	Exterior (x)
Interior (y)	2	1	2
Boundary (y)	1	1	1
Exterior (y)	2	1	2

按行展开这个矩阵（即按顺序连接第一行、第二行、第三行）得到字符串212111212。另一个例子可以用来说明该系统：下图中展示的关系（第三列和第一行中的第三个多边形对）可以用DE-9IM系统定义如下：

• 较大对象x的内部与y的内部、边界和外部之间的交集分别有维数2、1和2
• 较大对象x的边界与y的内部、边界和外部之间的交叉点分别有F，F 和1的维度，其中F表示false，物体是不相交的
• x的外部与y的内部、边界和外部之间的交集分别具有F、F和2的维度，更大对象的外部不接触y的内部或边界，但更小和更大对象的外部覆盖相同的区域

当这三个组件连接在一起时，创建字符串212、FF1和FF2。这与从函数st_relations()获得的结果相同(参见本章的源代码，查看图) 中的其他几何图形是如何创建的) :

xy2sfc = function(x, y) st_sfc(st_polygon(list(cbind(x, y))))
x = xy2sfc(x = c(0, 0, 1, 1, 0), y = c(0, 1, 1, 0.5, 0))
y = xy2sfc(x = c(0.7, 0.7, 0.9, 0.7), y = c(0.8, 0.5, 0.5, 0.8))
st_relate(x, y)
#>      [,1]       
#> [1,] "212FF1FF2"

理解DE-9IM strings可以推出新的二元空间谓词。?st_relate 帮助页面包含’Queen’关系和’rook’关系的函数定义，其中多边形共享边界或仅共享一个点。'Queen’关系意味着"边界-边界"关系(简单表格中第二列和第二行的单元格或 DE-9IM string的第五个元素)不能是空的，它相当于模式 F***T****，而对于’rook’关系，则相同的元素必须是1（表示线性相交）。这些措施的执行情况如下:

st_queen = function(x, y) st_relate(x, y, pattern = "F***T****")
st_rook = function(x, y) st_relate(x, y, pattern = "F***1****")

在先前创建的对象x的基础上，我们可以使用新创建的函数来找出网格中哪些元素是’queen’和’rook’相对于网格中间的正方形，如下所示:

grid = st_make_grid(x, n = 3)
grid_sf = st_sf(grid)
grid_sf$queens = lengths(st_queen(grid, grid[5])) > 0
plot(grid, col = grid_sf$queens)
grid_sf$rooks = lengths(st_rook(grid, grid[5])) > 0
plot(grid, col = grid_sf$rooks)

空间连接

连接两个非空间数据集依赖于一个共享的’key’变量，如矢量属性连接节中所述。空间数据连接应用了相同的概念，但是依赖于前面部分描述的空间关系。与属性数据一样，连接从源对象(y)向目标对象(参数x在联接函数中)添加新的列。

这个过程可以通过以下例子来说明：假设你在地球表面随机分布了十个点，问这些点中哪些在陆地上，属于哪些国家。在reproducible example中实现这个想法将会提高你的地理数据处理技能，并演示空间连接的运作方式。首先，需要创建随机分布在地球表面的点：

set.seed(2018) # set seed for reproducibility
(bb = st_bbox(world)) # the world's bounds
#>   xmin   ymin   xmax   ymax 
#> -180.0  -89.9  180.0   83.6
random_df = data.frame(
  x = runif(n = 10, min = bb[1], max = bb[3]),
  y = runif(n = 10, min = bb[2], max = bb[4])
)
random_points = random_df |> 
  st_as_sf(coords = c("x", "y"), crs = "EPSG:4326") # set coordinates and CRS

下图中的场景显示，Random_points对象(左上)缺乏属性数据，而world(右上)具有属性，包括图例中显示的国家样本的国家名称。空间连接使用st_join()实现，如下面的代码块所示。输出是Random_join对象，如图(左下)。在创建连接的数据集之前，我们使用空间子集提取来创建world_random，其中只包含含有随机点的国家，以验证在连接的数据集中返回的国家名称数量应该是4(参见下图的右上面板)。

world_random = world[random_points, ]
nrow(world_random)
#> [1] 4
random_joined = st_join(random_points, world["name_long"])

默认情况下，st_join()执行左连接，这意味着结果是一个包含来自 x 的所有行的对象，其中包括没有匹配 y 的行（请参见矢量属性连接），但是也可以通过设置参数 left = FALSE 进行内连接。与空间子集提取相似，st_join() 默认使用的拓扑运算符是 st_intersects()，可以通过设置 join 参数来更改（详见 ?st_join）。上面的示例演示了将多边形图层的一列添加到点图层的方法，但是该方法适用于任何几何类型。在这种情况下，例如当x包含多边形，每个多边形都与y中多个对象匹配时，空间连接将通过为每个y的匹配创建新行而导致重复的要素。

非重叠连接

有时候，两个地理数据集虽然没有触碰，但仍然存在着强烈的地理关系。cycle_hire和cycle_hire_osm这两个数据集已经被附加在spData包中，它们提供了很好的例子。将它们绘制出来可以发现它们经常密切相关，但并不触碰，如下图所示。以下代码创建了基础图：

plot(st_geometry(cycle_hire), col = "blue")
plot(st_geometry(cycle_hire_osm), add = TRUE, pch = 3, col = "red")

我们可以检查是否有任何点与下面所示的st_intersect()相同:

any(st_touches(cycle_hire, cycle_hire_osm, sparse = FALSE))
#> [1] FALSE

想象一下，我们需要将cycle_hire_osm中的capacity变量与官方"目标"数据中的 cycle_hire进行连接，此时需要使用非重叠连接。最简单的方法是使用二元谓词 st_is_within_distance()，如下所示，使用20米的阈值距离。如果启用了球面几何引擎（s2），则可以将度量单位的阈值距离设置为未投影数据（例如 lon/lat CRSs，如 WGS84），因为它在sf中默认启用（请参见s2节）。

sel = st_is_within_distance(cycle_hire, cycle_hire_osm, 
                            dist = units::set_units(20, "m"))
summary(lengths(sel) > 0)
#>    Mode   FALSE    TRUE 
#> logical     304     438

这表明，目标对象cycle_hire内有438个点位于与cycle_hire_osm的阈值距离内。如何检索与各个cycle_hire_osm点相关联的值？解决方案即为使用st_join()，但添加了一个dist参数（设为20 m）：

z = st_join(cycle_hire, cycle_hire_osm, st_is_within_distance, 
            dist = units::set_units(20, "m"))
nrow(cycle_hire)
#> [1] 742
nrow(z)
#> [1] 762

请注意，连接结果中的行数大于目标行数。这是因为在cycle_hire中的某些自行车租赁站在cycle_hire_osm中有多个匹配项。为了聚合重叠点的值并返回平均值，我们可以使用属性章节学习到的聚合方法，得到行数与目标相同的对象：

z = z |> 
  group_by(id) |> 
  summarize(capacity = mean(capacity))
nrow(z) == nrow(cycle_hire)
#> [1] TRUE

附近站点的容量可以通过比较源数据cycle_hire_osm的容量绘图和这个新对象的结果来进行验证(图表未显示)。

plot(cycle_hire_osm["capacity"])
plot(z["capacity"])

这种连接的结果使用了空间操作，以改变与简单要素相关的属性数据；而与每个要素相关的几何体保持不变。

空间聚合

与属性数据聚合相同，空间数据聚合也可以压缩数据：聚合结果的行数比非聚合输入要少。统计学的聚合函数（如平均数或总和）可以总结变量的多个数值，并返回每个分组变量的单个值。在矢量属性聚合节中，我们演示了如何使用aggregate()和group_by() |> summarize()根据属性变量压缩数据，本节介绍了相同的函数如何与空间对象一起工作。

回到新西兰的例子中，假设您想要查找每个地区高峰的平均高度，源对象（在这种情况下是 y 或 nz）定义了如何将目标对象（x 或 nz_height）中的价值分组。使用基础 R 的 aggregate() 方法可以在一行代码中完成此操作：

nz_agg = aggregate(x = nz_height, by = nz, FUN = mean)

前一条命令的结果是一个具有与（空间）聚合对象（nz）相同几何形状的sf对象，您可以使用命令identical(st_geometry(nz),st_geometry(nz_agg))来验证。前一操作的结果如下图所示，该图显示了在新西兰的16个地区中，nz_height每个要素的平均值。同样的结果也可以通过将st_join()的输出导入到"tidy"函数group_by()和summarize()中来生成，具体如下：

nz_agg2 = st_join(x = nz, y = nz_height) |>
  group_by(Name) |>
  summarize(elevation = mean(elevation, na.rm = TRUE))

通过使用函数mean()，nz_agg对象与分组对象nz具有相同的几何形状，但增加了一个新列，该列汇总每个地区中x的值。其他函数也可以用于此处，包括median()、sd()和其他每个组返回单个值的函数。需要注意的是，aggregate()与group_by() |> summarize()方法之间的一个区别是前者在不匹配的区域名称上显示NA值，而后者保留区域名称。因此，"tidy"方法在聚合功能和结果列名称方面更具灵活性。同时，也在合并节涵盖了创建新几何体的聚合操作。

连接不一致图层

空间一致性是与空间聚合相关的重要概念。聚合对象（我们将其称为y）与目标对象（x）是一致的，如果这两个对象有共享的边界。通常情况下，行政边界数据满足此条件，其中大单位——例如英国的中层超输出区（MSOAs）或许多其他欧洲国家的区域——由许多较小单位组成。

相比之下，不一致聚合对象与目标对象没有共同的边界。这对于空间聚合(和其他空间操作)是有问题的，如下图所示。聚合每个子区域的质心将不会返回准确的结果。面积插值通过将值从一组面积单位转移到另一组面积单位来克服这个问题，使用一系列算法，包括简单的面积加权方法和更复杂的方法，如"pycnophylactic"方法。

spData包有一个名为incongruent的数据集（在上图的右面板中带有黑色边框的彩色多边形），以及一个名为aggregating_zones的数据集（在的右面板中带有半透明蓝色边框的两个多边形）。假设incongruent的value列指的是以百万欧元为单位的总区域收入。我们如何将九个基础空间多边形的值转换为aggregating_zones中的两个多边形？

这个最简单、有用的方法是面积加权空间插值方法，它按照重叠面积的比例，将值从不相容对象传递到聚合区域中的一个新列中：输入和输出要素之间的空间交叉越大，相应的值也越大。这在下面的代码片段中通过st_interpolate_aw()实现。

iv = incongruent["value"] # keep only the values to be transferred
agg_aw = st_interpolate_aw(iv, aggregating_zones, extensive = TRUE)
#> Warning in st_interpolate_aw.sf(iv, aggregating_zones, extensive = TRUE):
#> st_interpolate_aw assumes attributes are constant or uniform over areas of x
agg_aw$value
#> [1] 19.6 25.7

案例中，由于总收入是一种所谓的空间广泛变量（随着地区增大而增加），所以将落入聚合区域的交叉值进行总结是有意义的。这里假设收入在较小的区域内均匀分布（所以有上面的警告信息）。但对于空间密集变量intensive（如平均收入或百分比），情况会有所不同，它们不会随着区域的增加而增加。st_interpolate_aw()在处理空间密集变量时同样有效：将extensive参数设置为FALSE，它将在执行聚合时使用平均函数而非求和函数。

距离关系

拓扑关系是二元的──一个要素要么相交，要么不相交──而距离关系是连续的。两个对象之间的距离是用st_distance()函数计算的。下面的代码块说明了这一点，找到了新西兰最高点与坎特伯雷地区地理重心之间的距离，该地理重心是在空间子集提取中创建的：

nz_highest = nz_height |> slice_max(n = 1, order_by = elevation)
canterbury_centroid = st_centroid(canterbury)
st_distance(nz_highest, canterbury_centroid)
#> Units: [m]
#>        [,1]
#> [1,] 115540

这个结果有两个令人惊讶之处:

• 它具有units，告诉我们距离是100,000米，而不是100,000英寸或任何其他距离的度量。
• 即使结果只包含一个单一值，它仍以矩阵形式返回

第二个特性暗示了st_length()的另一个有用的特性，它能够在对象x和y中的所有特征组合之间返回距离矩阵。下面的命令说明了这一点，它找出了nz_height中前三个特征与co所代表的奥塔戈和坎特伯雷地区之间的距离。

co = filter(nz, grepl("Canter|Otag", Name))
st_distance(nz_height[1:3, ], co)
#> Units: [m]
#>        [,1]  [,2]
#> [1,] 123537 15498
#> [2,]  94283     0
#> [3,]  93019     0

请注意，nz_height中第二和第三个要素以及co中第二个要素之间的距离为零。这证明了点与多边形之间的距离是指到多边形的任何部分的距离。nz_height中第二个和第三个点位于奥塔哥地区，这可以通过绘制它们来验证（结果未显示）：

plot(st_geometry(co)[2])
plot(st_geometry(nz_height)[2:3], add = TRUE)

栅格数据空间操作

本节建立在栅格数据操作基础上，该节重点介绍了处理栅格数据的各种基本方法，以便展示更高级和明确的空间栅格操作，并使用在本节中手动创建的elev和grain对象。为了方便读者，这些数据集也可以在spData包中找到。

elev = rast(system.file("raster/elev.tif", package = "spData"))
grain = rast(system.file("raster/grain.tif", package = "spData"))

空间栅格提取

前一章栅格数据操作展示了如何检索与特定单元格ID或行列组合相关联的值。栅格对象也可以通过位置（坐标）和其他空间对象进行提取。要使用坐标进行子集提取，可以使用terra函数cellFromXY()将坐标’translate’为单元格ID。另一种方法是使用terra::extract()（请注意，在tidyverse中还有一个名为extract()的函数）来提取值。下面展示了两种方法来找到覆盖在坐标为0.1，0.1的点上的单元格的值。

id = cellFromXY(elev, xy = matrix(c(0.1, 0.1), ncol = 2))
elev[id]
# the same as
terra::extract(elev, matrix(c(0.1, 0.1), ncol = 2))

栅格对象可以与另一个栅格对象进行子集提取，示例如下：

clip = rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45,
            resolution = 0.3, vals = rep(1, 9))
elev[clip]
# we can also use extract
terra::extract(elev, ext(clip))

这等于获取第一个栅格对象（在本例中为elev）的值，这些值落在第二个栅格（即clip）的范围内，如下图所示。

上述示例返回了特定单元格的值，但在许多情况下，需要对栅格数据集进行子集提取得到空间输出。这可以通过将[运算符的drop参数设置为FALSE来实现。下面的代码返回elev的前两个单元（即顶部行的前两个单元格），并作为栅格对象返回（仅显示输出的前两行）：

elev[1:2, drop = FALSE]    # spatial subsetting with cell IDs
#> class       : SpatRaster 
#> dimensions  : 1, 2, 1  (nrow, ncol, nlyr)
#> ...

另一个常见的空间子集提取例子是当一个具有logical（或NA）值的栅格用于掩膜具有相同范围和分辨率的另一个栅格，如上图所示。在这种情况下，可以使用[和 mask()函数（结果未显示）:

# create raster mask
rmask = elev
values(rmask) = sample(c(NA, TRUE), 36, replace = TRUE)

在上面的代码块中，我们创建了一个名为rmask的掩膜对象，其值被随机分配为NA和TRUE。接下来，我们想要保留那些在rmask中为TRUE的elev值。换句话说，我们想要使用rmask屏蔽elev。

# spatial subsetting
elev[rmask, drop = FALSE]           # with [ operator
mask(elev, rmask)                   # with mask()

以上方法还可用于使用NA替换某些值（例如，预计出现错误）。

elev[elev < 20] = NA

这些操作实际上是布尔局部操作，因为我们逐单元格比较了两个栅格。下一小节将更详细地探讨这些和相关的操作。

地图代数

“地图代数"是在上世纪70年代末提出的，用于描述地理栅格数据和（虽然不那么突出）矢量数据的分析的"一套约定、功能和技术”。在这个背景下，我们更加明确地定义地图代数，作为一种修改或汇总栅格单元值的操作，涉及周围单元、区域或应用于每个单元的统计函数。

地图代数操作往往是快速的，因为栅格数据集只隐式地存储坐标，因此有一句古老格言 “栅格更快，但矢量更正确”。栅格数据集中单元的位置可以通过使用其矩阵位置和数据集的分辨率和原点（存储在标头中）来计算。然而，就处理而言，只要我们确保处理后单元格位置不变，单元格的地理位置就几乎不相关。此外，如果两个或多个栅格数据集具有相同的范围、投影和分辨率，可以将它们视为矩阵进行处理。

这是使用terra包进行地图代数的方法。首先，会查询栅格数据集的标头，并（在需要处理多个数据集的地图代数运算中）检查数据集是否兼容。其次，地图代数保留所谓的一对一定位对应关系，这意味着单元格不能移动。这与矩阵代数不同，矩阵代数中，例如在矩阵乘法或除法时，值会改变位置。

地图代数（或栅格数据的制图建模）将栅格操作分为四个子类，每个子类同时在一个或多个栅格上进行操作：

1. 局部或逐单元格操作
2. 焦点或邻域操作。最常见的输出单元格值是3x3输入单元格块的结果
3. 区域操作与焦点操作相似，但计算新值的周围像素网格可能具有不规则的大小和形状
4. 全局或逐栅格操作。这意味着输出单元格可能从一个或多个整个栅格派生其值

这种分类按照用于每个像素处理步骤的单元格数量和输出类型对地图代数操作进行分类。为了完整起见，我们应该提到栅格操作还可以按学科分类，例如地形、水文分析或图像分类。以下部分解释了每种类型的地图代数操作如何使用，参考了实际示例。

局部操作

局部操作包括在一个或多个层上的所有逐单元格操作。栅格代数是局部操作的典型用例——这包括从栅格中添加或减去值，平方和乘以栅格。栅格代数还允许逻辑操作，例如查找大于特定值的所有栅格单元格（在我们下面的示例中为5）。terra包支持所有这些操作及更多，如下图所示：

elev + elev
elev^2
log(elev)
elev > 5

局部操作的另一个很好的例子是将数字值的区间分为诸如将数字高程模型分为低（类别1）、中（类别2）和高（类别3）海拔的组。使用classify()命令，我们首先需要构建一个重分类矩阵，其中第一列对应于类别的下限，第二列对应于上限。第三列代表第一和第二列中指定范围的新值。

rcl = matrix(c(0, 12, 1, 12, 24, 2, 24, 36, 3), ncol = 3, byrow = TRUE)
rcl
#>      [,1] [,2] [,3]
#> [1,]    0   12    1
#> [2,]   12   24    2
#> [3,]   24   36    3

在这里，我们将范围为0-12、12-24和24-36的栅格值进行了重分类，分别赋值为1、2和3。

recl = classify(elev, rcl = rcl)

classify()函数也可以用于我们想要减少分类栅格中类别数量的情况。我们将在地理营销中执行几个附加的重分类操作。

除了算术运算符外，还可以使用app()、tapp()和lapp()函数。它们更有效率，因此，在大型栅格数据集存在的情况下，它们是首选。此外，它们允许你直接保存输出文件。app()函数将一个函数应用到栅格的每个单元格，并用于将多个图层的值汇总（例如，计算总和）到一个图层。tapp()是app()的扩展，允许我们选择要执行某个操作的图层子集（参见index参数）。最后，lapp()函数允许使用图层作为参数将一个函数应用到每个单元格——下面将介绍lapp()的一个应用。

归一化植被指数（NDVI）的计算是一种众所周知的局部（逐像素）栅格操作。它返回一个值介于-1和1之间的栅格；正值表示存在活植物（通常 > 0.2）。NDVI是从遥感影像的红色和近红外（NIR）波段计算的，通常来自Landsat或Sentinel等卫星系统。植被在可见光光谱中大量吸收光线，特别是在红色通道中，同时反射NIR光线，从而解释了NVDI公式：

$$
\begin{split}
NDVI&= \frac{\text{NIR} - \text{Red}}{\text{NIR} + \text{Red}}\
\end{split}
$$

让我们计算锡安国家公园的多光谱卫星影像的NDVI。

multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
multi_rast = rast(multi_raster_file)

栅格对象具有四个卫星波段：蓝、绿、红和近红外（NIR）。我们下一步应该在 R 函数中使用NDVI公式：

ndvi_fun = function(nir, red){
  (nir - red) / (nir + red)
}

这个函数接受两个数值参数，nir和red，并返回一个带有NDVI值的数值向量。它可以用作lapp()的fun参数。我们只需记住，我们的函数只需要两个波段（不是原始栅格的四个），并且它们需要按NIR、red的顺序排列。这就是为什么我们在进行任何计算之前，使用multi_rast[[c(4, 3)]]来提取输入栅格的子集。

ndvi_rast = lapp(multi_rast[[c(4, 3)]], fun = ndvi_fun)

结果显示在下图的右侧面板上，可以与同一区域的RGB图像（同一图的左侧面板）进行比较。这让我们看到最大的NDVI值与该区域北部的密集森林区域相连，而最低的值则与北部的湖泊和积雪覆盖的山脊有关。

预测映射是局部栅格操作的另一个有趣应用。响应变量对应于空间中测量或观察到的点，例如物种丰富度、滑坡的存在、树木疾病或农作物产量。因此，我们可以轻松从各种栅格（例如海拔、pH值、降水量、温度、土地覆盖、土壤类型等）检索空间或预测变量。随后，我们使用lm()、glm()、gam()或机器学习技术将响应建模为预测因子的函数。因此，将估计系数应用于预测栅格值，并对输出栅格值求和（参见 生态章节），可以对栅格对象进行空间预测。

焦点操作

虽然局部函数可能在多个层上操作一个单元格，但焦点操作会考虑中心（焦点）单元格及其邻居。通常考虑的邻域（也称为内核、滤波器或移动窗口）大小为3x3个单元格（即中心单元格及其周围的八个邻居），但可以根据用户的定义采用任何其他（不一定是矩形的）形状。焦点操作将一个聚合函数应用于指定邻域内的所有单元格，将相应的输出用作中心单元格的新值，然后移至下一个中心单元格，如下图所示。此操作的其他名称是空间滤波和卷积。

在R中，我们可以使用focal()函数来执行空间滤波。我们通过一个matrix定义移动窗口的形状，其值对应于权重（参见下面代码块中的w参数）。其次，fun参数让我们指定我们希望应用于这个邻域的函数。这里，我们选择最小值，但可以使用任何其他汇总函数，包括sum()、mean()或var()。

r_focal = focal(elev, w = matrix(1, nrow = 3, ncol = 3), fun = min)

这个函数还接受其他参数，例如，在过程中是否应删除NA值（na.rm = TRUE）或不删除（na.rm = FALSE）。

我们可以快速检查输出是否符合我们的预期。在我们的示例中，最小值必须始终位于移动窗口的左上角（请记住，我们通过从左上角开始逐行增加单元格值来创建输入栅格）。在此示例中，权重矩阵只包括1，这意味着每个单元格对输出具有相同的权重，但这可以更改。

焦点函数或滤波器在图像处理中起着主导作用。低通或平滑滤波使用平均函数来消除极值。在分类数据的情况下，我们可以用众数（最常见的值）替换平均数。相反，高通滤波强调特征。此处，线检测的拉普拉斯和索贝尔滤波可以作为例子。你可以查看focal()的帮助页面了解如何在 R 中使用它们（这也将在本章节结束时的练习中使用）。

地形处理，计算地形特性如坡度、方向和流向，依赖于焦点函数。虽然terrain()可用于计算这些指标，但一些地形算法（包括用于计算坡度的Zevenbergen和Thorne方法）并未在这个terra函数中实现。许多其他算法——包括曲率、贡献区域和湿度指数 ——是在开源桌面地理信息系统（GIS）软件中实现的。GIS桥梁章介绍了如何从R内部访问这样的GIS功能。

区域操作

与焦点操作一样，区域操作将聚合函数应用于多个栅格单元格。然而，在区域操作的情况下，第二个栅格（通常具有分类值）定义了分区滤波器（或"区域"），与上一节中呈现的焦点操作的预定义邻域窗口相对。因此，定义区域滤波器的栅格单元格不一定要相邻。grain 栅格就是一个很好的例子，不同的粒度大小在整个栅格中不规则地分布。最后，区域操作的结果是按区域分组的汇总表，这就是为什么这个操作在GIS世界中也被称为分区统计。这与返回栅格对象的焦点操作形成了对比。

以下代码块使用zonal()函数来计算与每个粒度类别相关的平均海拔。

z = zonal(elev, grain, fun = "mean")
z
#>   grain elev
#> 1  clay 14.8
#> 2  silt 21.2
#> 3  sand 18.7

这将返回每个类别的统计，在这里是每个粒度大小类别的平均海拔。注意：通过将as.raster参数设置为TRUE，也可以获得每个区域的计算统计数据的栅格。

全局操作和距离

全局操作是分区操作的特殊情况，整个栅格数据集代表单个区域。整个栅格数据集的描述性统计是最常见的全局操作，例如最小值或最大值——我们已经在章节栅格数据汇总中讨论过这些。

除此之外，全局操作还可用于计算距离和权重栅格。在第一种情况下，可以计算每个单元格到特定目标单元格的距离。例如，人们可能想要计算到最近海岸的距离（参见terra::distance()）。我们也可能想要考虑地形，这意味着，我们不仅对纯粹的距离感兴趣，而且还想避免在前往海岸时穿越山脉。为此，我们可以通过海拔为距离赋权，以便每增加一个海拔米就"延长"欧氏距离。可见性和视域计算也属于全局操作的一类。

矢量数据的"地图代数"

许多地图代数操作在矢量处理中有对应操作。在仅考虑最大距离（逻辑焦点操作）的情况下计算距离栅格（全局操作）等同于矢量缓冲操作（裁剪章节）。重分类栅格数据（根据输入是局部还是区域函数）等同于溶解矢量数据（空间连接章节）。将两个栅格叠加（局部操作），其中一个包含表示遮罩的NULL或NA值，类似于矢量裁剪（章节）。与空间裁剪非常相似的是交叉两个图层（空间矢量提取章节）。区别在于这两个图层（矢量或栅格）仅共享重叠区域。不过，要小心措辞。有时相同的词在栅格和矢量数据模型中具有略有不同的含义。虽然聚合多边形几何意味着溶解边界，但对于栅格数据几何，它意味着增加单元格大小，从而降低空间分辨率。区域操作根据另一个栅格数据集的区域（类别）使用聚合函数溶解一个栅格的单元格。

合并栅格

假设我们想计算NDVI（见局部操作章节），并且还想从高程数据中计算研究区域内观测的地形属性。这些计算依赖于遥感信息。相应的图像常常被分割成覆盖特定空间范围的场景，而研究区域通常覆盖多个场景。那么，我们就需要合并研究区域所覆盖的场景。在最简单的情况下，我们可以只是合并这些场景，即将它们并排放置。例如，可以使用数字高程数据（SRTM，ASTER）来实现。在下面的代码块中，我们首先下载了奥地利和瑞士的SRTM高程数据（有关国家代码，请参见geodata函数 country_codes()）。第二步，我们将这两个栅格合并为一个。

aut = geodata::elevation_30s(country = "AUT", path = tempdir())
ch = geodata::elevation_30s(country = "CHE", path = tempdir())
aut_ch = merge(aut, ch)

terra 的 merge()命令结合了两个图像，如果它们重叠，则使用第一个栅格的值。

当重叠的值彼此不对应时，这种合并方法用处不大。当您想合并在不同日期拍摄的场景的光谱图像时，通常会出现这种情况。merge() 命令仍然会工作，但在结果图像中您会看到一个明显的边界。另一方面，mosaic() 命令允许您为重叠区域定义一个函数。例如，我们可以计算平均值 – 这可能会平滑合并结果中的明显边界，但最有可能的是它不会让其消失。

练习

library(sf)
library(dplyr)
library(spData)

E1. 在章节2.2中，我们确认了坎特伯雷是新西兰包含最高的100个点中最多点的地区。坎特伯雷地区包含了多少个这些高点？

加分题: 使用 plot() 函数绘制结果，展示整个新西兰，canterbury 地区以黄色高亮，坎特伯雷的高点用红色的交叉表示（提示: pch = 7），而新西兰其它地区的高点用蓝色的圆表示。参见帮助页面 ?points 以获得不同 pch 值的图示。

E2. 哪个地区拥有第二高数量的 nz_height 点，它有多少个点？

E3. 将问题推广到所有地区：新西兰的16个地区中有多少个包含了国内最高100个点中的点？是哪些地区？

• 加分题: 创建一个表格，按照点的数量和地区名字的顺序列出这些地区。

E4. 通过查找和绘制美国各州之间以及与其他空间对象的关系来测试你的空间谓词知识。

这个练习的起点是创建一个代表美国科罗拉多州的对象。使用命令colorado = us_states[us_states$NAME == "Colorado",] (base R) 或使用 filter() 函数 (tidyverse) 完成，并在美国各州的背景下绘制所得对象。

• 创建一个新对象，代表所有与科罗拉多州地理上相交的州，并绘制结果（提示: 使用子集方法 [ 是最简洁的方式）。
• 创建另一个对象，代表所有与科罗拉多州接触（有共享边界的）的对象，并绘制结果（提示: 记住你可以在 base R 的空间子集操作中使用参数 op = st_intersects 和其他空间关系）。
• 加分题: 从美国东海岸附近的哥伦比亚特区的中心到美国西海岸附近的加利福尼亚州的中心创建一条直线（提示: 第5章描述的函数 st_centroid()、st_union() 和 st_cast() 可能会有帮助），并确定这条长的东西线穿过了哪些州。

E5. 使用 dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))，并将海拔分为三个类别: 低(<300)，中等和高(>500)。其次，读取 NDVI 栅格（ndvi = rast(system.file("raster/ndvi.tif", package = "spDataLarge"))），并计算每个海拔类别的平均 NDVI 和平均海拔。

E6. 对 rast(system.file("ex/logo.tif", package = "terra")) 应用线检测滤波器。绘制结果。提示: 阅读 ?terra::focal()。

E7. 计算Landsat图像的归一化差异水指数（NDWI; (green - nir)/(green + nir)）。使用 spDataLarge 包提供的Landsat图像（system.file("raster/landsat.tif", package = "spDataLarge")）。同时，计算该区域的 NDVI 和 NDWI 之间的相关性（提示: 你可以使用 layerCor() 函数）。

E8. StackOverflow 的一篇帖子展示了如何使用 raster::distance() 计算到最近海岸线的距离。尝试做类似的事情，但使用 terra::distance()：检索西班牙的数字高程模型，并计算表示该国各地到海岸的距离的栅格（提示: 使用 geodata::elevation_30s()）。将得到的距离从米转换为公里。注意: 在此操作期间，增加输入栅格的单元大小以减少计算时间可能是明智的（aggregate()）。

E9. 尝试通过将海拔栅格与距离栅格结合，修改上述练习中使用的方法；每上升100个海拔米，到海岸的距离应增加10公里。接下来，计算和可视化使用欧几里得距离（E7）创建的栅格与按海拔加权的栅格之间的差异。