(5)几何操作

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

本章重点介绍如何操作地理对象的地理元素，例如通过简化和转换矢量几何图形、裁剪栅格数据集以及将矢量对象转换为栅格并从栅格转换为矢量。阅读本章并尝试章末的练习后，你应该能理解并控制sf对象中的几何列，以及与其他地理对象相关的栅格中所表示像素的范围和地理位置。

前提条件

• 本章使用了与R中的地理数据章相同的包，但额外增加了spDataLarge，该包在属性操作章中已安装。

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

引言

迄今为止，本书已经解释了地理数据集的结构、如何基于它们的非地理属性和空间关系进行操作。本章重点介绍如何操作地理对象的地理元素，例如通过简化和转换矢量几何图形、裁剪栅格数据集以及将矢量对象转换为栅格并从栅格转换为矢量。阅读本章并尝试章末的练习后，你应该能理解并控制sf对象中的几何列，以及与其他地理对象相关的栅格中所表示像素的范围和地理位置。

矢量数据几何操作节介绍了使用’一元’和’二元’操作转换矢量几何图形。一元操作在单个几何图形上独立工作，包括简化（线和多边形）、缓冲区和质心的创建，以及使用’仿射变换’移动/缩放/旋转单个几何图形。二元变换基于另一个形状修改一个几何图形，包括剪裁和几何联合节中介绍。类型转换（例如从多边形到线）。

节涵盖了栅格对象上的几何变换。这涉及改变底层像素的大小和数量，并赋予它们新的值。它教你如何更改分辨率（也称为栅格聚合和分解）、范围和栅格的原点。如果想要对齐来自不同来源的栅格数据集，这些操作特别有用。对齐的栅格对象在像素之间共享一对一对应关系，允许它们使用地图代数操作进行处理。栅格和矢量对象之间的交互在后面章节中介绍。它介绍了如何通过矢量几何图形“掩膜”和“提取”栅格值。重要的是，它还展示了如何将栅格“矢量化”和将矢量数据集“栅格化”，使两种数据模型更易于互换。

矢量数据几何操作

这一部分涉及以某种方式改变矢量（sf）对象几何形状的操作。与上一章节中介绍的空间数据操作相比，这里更为高级，因为我们深入到几何结构中：本节讨论的函数不仅作用于 sf 类对象，还作用于 sfc 类对象。

简化

简化是一种矢量对象（线和多边形）的概括过程，通常用于较小比例的地图中。简化对象的另一个原因是减少它们消耗的内存、磁盘空间和网络带宽。在将复杂几何体发布为交互式地图之前，简化它们可能是明智的选择。sf 包提供了 st_simplify() 函数，使用 GEOS 实现的 Douglas-Peucker 算法来减少顶点数量。st_simplify() 使用 dTolerance 来以地图单位控制概括的级别 [具体细节请参见 @douglas_algorithms_1973]。下图展示了代表塞纳河及其支流的 LINESTRING 几何体的简化过程。以下命令创建了简化的几何体：

seine_simp = st_simplify(seine, dTolerance = 2000)  # 2000 m

Comparison of the original and simplified geometry of the seine object.

生成的 seine_simp 对象是原始 seine 的副本，但顶点更少。这一点很明显，结果在视觉上更简单（上图，右侧），并且比原始对象占用的内存更少，如下所验证：

object.size(seine)
#> 18096 bytes
object.size(seine_simp)
#> 9112 bytes

简化也适用于多边形。这通过表示美国本土的 us_states 来说明。正如我们在后面章节中所展示的，GEOS 假设数据在投影的坐标参考系统（CRS）中，而当使用地理 CRS 时，这可能会导致意外的结果。因此，第一步是将数据投影到某个适当的投影 CRS 中，例如美国国家地图等面积（EPSG = 2163）（如上图左上角所示）：

us_states2163 = st_transform(us_states, "EPSG:2163")
us_states2163 = us_states2163

st_simplify() 也同样很好的处理投影面要素。

us_states_simp1 = st_simplify(us_states2163, dTolerance = 100000)  # 100 km

st_simplify()的一个限制是它基于每个几何体来简化对象。这意味着失去了“拓扑”，导致图@ref(fig:us-simp)（右上面板）中所示的重叠和“多洞”区域单位。 rmapshaper 的 ms_simplify() 提供了克服这个问题的替代方案。默认情况下，它使用 Visvalingam 算法，克服了 Douglas-Peucker算法的一些。 以下代码块使用此函数来简化us_states2163。结果只有输入顶点的1%（通过参数 keep设置），但其对象数量保持不变，因为我们设置了 keep_shapes = TRUE：[^1]

[^1]：即使将 keep_shapes 参数设置为TRUE，也可以通过简化多边形对象来删除小的内部多边形。要防止这种情况，您需要设置 explode = TRUE。此选项在简化之前将所有多边形转换为单独的多边形。

# proportion of points to retain (0-1; default 0.05)
us_states_simp2 = rmapshaper::ms_simplify(us_states2163, keep = 0.01,
                                          keep_shapes = TRUE)

对于简化的另一种选择是平滑多边形和线条几何图形的边界，这在 smoothr 包中实现。平滑插值几何图形的边缘，并不一定会导致较少的顶点，但在使用从空间向量化栅格产生的几何图形时尤为有用。

smoothr 包实现了三种平滑技术：高斯核回归、Chaikin的角剪切算法和样条插值。这些都在website.中进行了描述。请注意，与 st_simplify() 类似，平滑算法不保留“拓扑结构”。smoothr 的主要功能是 smooth()，其中 method 参数指定要使用的平滑技术。这些平滑技术可能在进行地理空间分析和图形表示时很有用，特别是在需要减少几何形状的锯齿状外观或复杂度时。通过合理选择平滑参数，可以在保留地理特征的同时，提高几何形状的视觉美感。以下是一个使用高斯核回归（Gaussian kernel regression）来通过method=ksmooth平滑美国各州边界的例子。smoothness 参数控制用于平滑几何形状的高斯带宽，其默认值为1。

us_states_simp3 = smoothr::smooth(us_states2163, 
                                  method = "ksmooth", smoothness = 6)

最后，原始数据集与简化版和平滑版的视觉比较显示在（下图）中。可以观察到Douglas-Peucker（st_simplify）、Visvalingam（ms_simplify）和高斯核回归（smooth(method=ksmooth)）算法的输出之间的差异。

Polygon simplification in action, comparing the original geometry of the contiguous United States with simplified versions, generated with functions from sf (top-right), rmapshaper (bottom-left), and smoothr (bottom-right) packages.

质心

质心运算确定地理对象的中心。就像统计的集中趋势度量（包括平均数和中位数定义的“平均”）一样，有许多方法可以定义对象的地理中心。它们都为更复杂的矢量对象创建单一点表示。

最常用的质心运算是地理质心。这种类型的质心运算（通常称为“质心”）代表空间对象中的质量中心（想象一下用你的手指平衡一个盘子）。地理质心有许多用途，例如创建复杂几何图形的简单点表示，或估计多边形之间的距离。可以使用sf函数st_centroid()来计算，如下面的代码所示，该代码生成了新西兰地区和塞纳河支流的地理质心，以黑点在图中示出。

nz_centroid = st_centroid(nz)
seine_centroid = st_centroid(seine)

有时，地理质心会落在其父对象的边界之外（想象一个甜甜圈）。在这种情况下，可以使用表面上的点运算来确保点位于父对象内（例如，用于标记不规则的多边形对象，如岛屿状态），如图中的红点所示。请注意，这些红点总是位于其父对象上。它们是使用st_point_on_surface()创建的，如下所示：[^2]

[^2]: 网址 https://gis.stackexchange.com/a/76563/20955提供了st_point_on_surface()如何工作？

nz_pos = st_point_on_surface(nz)
seine_pos = st_point_on_surface(seine)

Centroids (black points) and ‘points on surface’ (red points) of New Zealand’s regions (left) and the Seine (right) datasets.

其他类型的质心也存在，包括切比雪夫中心和视觉中心。我们不在此处探讨这些内容，但可以使用R来计算它们，我们将在算法章中看到。

缓冲区

缓冲区是代表距离几何特征一定距离内的区域的多边形：无论输入是点、线还是多边形，输出都是多边形。与简化（通常用于可视化和减小文件大小）不同，缓冲区倾向于用于地理数据分析。
有多少个点在这条线的一定距离内？哪些人口群体在这个新商店的旅行距离内？通过围绕感兴趣的地理实体创建缓冲区，可以回答并可视化这些问题。

下图描述了围绕塞纳河及其支流的不同大小（5和50公里）的缓冲区。这些缓冲区是通过下面的命令创建的，这些命令显示 st_buffer() 命令至少需要两个参数：一个输入几何和一个距离，以CRS的单位提供（在本例中为米）：

seine_buff_5km = st_buffer(seine, dist = 5000)
seine_buff_50km = st_buffer(seine, dist = 50000)

Buffers around the Seine dataset of 5 km (left) and 50 km (right). Note the colors, which reflect the fact that one buffer is created per geometry feature.

📌st_buffer() 的第三个也是最后一个参数是 nQuadSegs，表示“每个象限的线段数”，默认设置为30（意味着缓冲区创建的圆由 $4 \times 30 = 120$ 条线组成）。
这个参数很少需要设置。当缓冲操作的输出占用的内存成为主要问题（在这种情况下应该减少）或者需要非常高的精度（在这种情况下应该增加）时，它可能会有用。

仿射变换

仿射变换是保留线条和平行性的任何变换。然而，角度或长度并不一定会保留。仿射变换包括平移（移位）、缩放和旋转等。此外，可以使用这些操作的任何组合。仿射变换是地理计算的重要组成部分。例如，标签放置需要移位，非连续区域制图用到缩放，而在重新投影或改进基于失真或错误投影的地图创建的几何图形时，会应用许多仿射变换。sf 包为 sfg 和 sfc 类的对象实现了仿射变换。

nz_sfc = st_geometry(nz)

平移是将每个点以地图单位的相同距离移动。它可以通过将数值向量添加到矢量对象来完成。例如，下面的代码将所有y坐标向北移动100,000米，但不改变x坐标。

nz_shift = nz_sfc + c(0, 100000)

缩放通过一个因子放大或缩小对象。它可以全局或局部应用。全局缩放增加或减小与原点坐标有关的所有坐标值，同时保持所有几何图形的拓扑关系不变。它可以通过减法或乘法对sfg或sfc对象进行操作。

局部缩放独立处理几何体，并需要确定几何体将围绕哪些点进行缩放，例如，质心。在下面的示例中，每个几何体都围绕质心缩小了一半。为实现这一目的，首先要以一种方式移动每个对象，使其中心的坐标为 0,0（(nz_sfc - nz_centroid_sfc)）。接下来，减小几何体的尺寸一半（* 0.5）。最后，将每个对象的质心移回到输入数据的坐标（+ nz_centroid_sfc）。

nz_centroid_sfc = st_centroid(nz_sfc)
nz_scale = (nz_sfc - nz_centroid_sfc) * 0.5 + nz_centroid_sfc

对二维坐标的旋转需要一个旋转矩阵：

$$
R =
\begin{bmatrix}
\cos \theta & -\sin \theta \
\sin \theta & \cos \theta \
\end{bmatrix}
$$
它以顺时针方向旋转点。旋转矩阵可以在R中实现如下：

rotation = function(a){
  r = a * pi / 180 #degrees to radians
  matrix(c(cos(r), sin(r), -sin(r), cos(r)), nrow = 2, ncol = 2)
} 

rotation函数接受一个参数a - 以度为单位的旋转角度。旋转可以围绕选定的点（例如质心）进行。有关更多示例，请参见vignette("sf3")。

nz_rotate = (nz_sfc - nz_centroid_sfc) * rotation(30) + nz_centroid_sfc

Illustrations of affine transformations: shift, scale and rotate.

最后，可以使用st_set_geometry()函数将新创建的几何图形替换旧的几何图形。

nz_scale_sf = st_set_geometry(nz, nz_scale)

裁剪

空间裁剪是一种空间子集处理方式，至少涉及一些受影响特征的 geometry 列的更改。

裁剪只能应用于比点更复杂的特征：线、多边形及其“多”的要素。为了说明这个概念，我们将从一个简单的例子开始：两个重叠的圆，中心点彼此相距一单位，半径为一。

b = st_sfc(st_point(c(0, 1)), st_point(c(1, 1))) # create 2 points
b = st_buffer(b, dist = 1) # convert points to circles
plot(b, border = "grey")
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y"), cex = 3) # add text

Overlapping circles.

假设你想选择的不是一个圆或另一个圆，而是由x和y都覆盖的空间。这可以使用函数st_intersection()来实现，使用名为x和y的对象来表示左侧和右侧的圆。

x = b[1]
y = b[2]
x_and_y = st_intersection(x, y)
plot(b, border = "grey")
plot(x_and_y, col = "lightgrey", border = "grey", add = TRUE) # intersecting area

Overlapping circles with a gray color indicating intersection between them.

以下代码块演示了这如何适用于代表x和y的“Venn”图的所有组合，灵感来自*《R for Data Science》*一书的Figure 5.1。

Spatial equivalents of logical operators.

子集提取和裁剪

裁剪对象可以改变其几何形状，但也可以对对象进行子集提取，只返回与裁剪/子集对象相交（或部分相交）的特征。为了说明这一点，我们将对覆盖图中的圆圈x和y的边界框的点进行子集提取。

有些点将仅在一个圆圈内，有些将在两者内，有些则在两者之外。下面使用st_sample()在圆圈x和y的范围内生成简单随机分布的点，得到图中所示的输出，提出了这样的问题：如何对点进行子集提取，以便仅返回与两者x和y相交的点？

bb = st_bbox(st_union(x, y))
box = st_as_sfc(bb)
set.seed(2017)
p = st_sample(x = box, size = 10)
p_xy1 = p[x_and_y]
plot(box, border = "grey", lty = 2)
plot(x, add = TRUE, border = "grey")
plot(y, add = TRUE, border = "grey")
plot(p, add = TRUE)
plot(p_xy1, cex = 3, col = "red", add = TRUE)
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y"), cex = 2)

Randomly distributed points within the bounding box enclosing circles x and y. The point that intersects with both objects x and y is highlighted.

bb = st_bbox(st_union(x, y))
box = st_as_sfc(bb)
set.seed(2017)
p = st_sample(x = box, size = 10)
x_and_y = st_intersection(x, y)

下面的代码块演示了实现相同结果的三种方法。我们可以直接使用x和y的交集（由上一个代码块中的x_and_y表示）作为子集对象，如下面代码块的第一行所示。我们还可以找到输入点（由p表示）与子集/裁剪对象x_and_y之间的交集，如下面代码块的第二行所示。这第二种方法将返回部分与x_and_y相交的特征，但是对于跨越子集对象边界的空间广泛特征，几何形状将被修改。第三种方法是使用二元空间谓词st_intersects()创建子集对象，该谓词在前一章中介绍。

结果是相同的（除了属性名称的表面差异），但实现差异很大：

p_xy1 = p[x_and_y] # way #1
p_xy2 = st_intersection(p, x_and_y) # way #2
sel_p_xy = st_intersects(p, x, sparse = FALSE)[, 1] & 
  st_intersects(p, y, sparse = FALSE)[, 1] # way #3
p_xy3 = p[sel_p_xy]

尽管上面的示例相当做作，主要用于教育目的而不是应用目的，我们鼓励读者复现结果以加深您对在R中处理地理矢量对象的理解，但它提出了一个重要问题：应该使用哪种实现？通常，应优先选择更简洁的实现，即上述的第一种方法。我们将在算法章中回到选择相同技术或算法的不同实现之间的问题。

合并

正如我们在节中所看到的，空间聚合可以无声无息地溶解同一组中相接触的多边形的几何形状。下面的代码块演示了使用基础和dplyr函数将49个us_states聚合成四个区域的过程（见图中的结果）：

regions = aggregate(x = us_states[, "total_pop_15"], by = list(us_states$REGION),
                    FUN = sum, na.rm = TRUE)
regions2 = us_states |> 
  group_by(REGION) |>
  summarize(pop = sum(total_pop_15, na.rm = TRUE))

Spatial aggregation on contiguous polygons, illustrated by aggregating the population of US states into regions, with population represented by color. Note the operation automatically dissolves boundaries between states.

在几何方面发生了什么情况呢？在幕后，aggregate()和summarize()都使用st_union()组合几何体并溶解它们之间的边界。下面的代码块展示了如何创建一个统一的美国西部：

us_west = us_states[us_states$REGION == "West", ]
us_west_union = st_union(us_west)

该函数可以接受两个几何体并将它们联合起来，如下面的代码块所示，其中创建了一个包括德克萨斯州的统一的西部区块（挑战：重现并绘制结果）：

texas = us_states[us_states$NAME == "Texas", ]
texas_union = st_union(us_west_union, texas)

类型转化

几何类型转换是一项强大的操作，可实现几何类型的转换。它在 sf 包的 st_cast() 函数中实现。重要的是，st_cast() 在单一简单要素几何（sfg）对象、简单要素几何列（sfc）和简单要素对象上的行为不同。

让我们创建一个多点对象来说明几何类型转换如何在简单要素几何（sfg）对象上工作：

multipoint = st_multipoint(matrix(c(1, 3, 5, 1, 3, 1), ncol = 2))

在这种情况下，st_cast()可以用于将新对象转换为线或多边形。

linestring = st_cast(multipoint, "LINESTRING")
polyg = st_cast(multipoint, "POLYGON")

Examples of a linestring and a polygon casted from a multipoint geometry.

从多点到线串的转换是一种常见的操作，它通过有序的点观测（例如 GPS 测量或地理标签媒体）创建线对象。反过来，这允许执行空间操作，例如计算所走路径的长度。从多点或线串到多边形的转换通常用于计算面积，例如从围绕湖泊的 GPS 测量集或从建筑用地的角落。

使用 st_cast() 还可以反转转换过程。

multipoint_2 = st_cast(linestring, "MULTIPOINT")
multipoint_3 = st_cast(polyg, "MULTIPOINT")
all.equal(multipoint, multipoint_2)
#> [1] TRUE
all.equal(multipoint, multipoint_3)
#> [1] TRUE

📌对于单个简单要素几何体（sfg），st_cast()还提供了从非多重类型到多重类型（例如，从POINT到MULTIPOINT）以及从多重类型到非多重类型的几何转换。然而，当从多重类型转换到非多重类型时，旧对象中的第一个元素将保留在输出对象中。

几何转换对于简单特征几何列（sfc）和简单特征对象在大多数情况下与sfg的工作方式相同。其中一个重要的区别是多重类型与非多重类型之间的转换。因此，sfc或sf的多重对象被拆分成许多非多重对象。

下表显示了对简单特征对象的可能几何类型转换。单个简单特征几何（由表格的第一列表示）可以转换为多个几何类型，由下表的列表示。有些转换是不可能的：例如，您不能将单个点转换为多段线或多边形，这就解释了表格中的单元格[1，4:5]包含NA的原因。有些转换将单个要素输入拆分为多个子要素，对sf对象进行“扩展”（添加具有重复属性值的新行）。例如，当由五对坐标组成的多点几何体转换为’POINT’几何体时，输出将包含五个特征。

Geometry casting on simple feature geometries (see Section 2.1) with input type by row and output type by column

	POI	MPOI	LIN	MLIN	POL	MPOL	GC
POI(1)	1	1	1	NA	NA	NA	NA
MPOI(1)	4	1	1	1	1	NA	NA
LIN(1)	5	1	1	1	1	NA	NA
MLIN(1)	7	2	2	1	NA	NA	NA
POL(1)	5	1	1	1	1	1	NA
MPOL(1)	10	1	NA	1	2	1	1
GC(1)	9	1	NA	NA	NA	NA	1
注意：像（1）这样的数值代表特征的数量；NA意味着该操作不可能。缩写：POI、LIN、POL 和 GC 分别指 POINT、LINESTRING、POLYGON 和 GEOMETRYCOLLECTION。这些几何类型的 MULTI 版本由一个前置的 M 表示，例如，MPOI 是 MULTIPOINT 的缩写。

让我们尝试在新对象multilinestring_sf上应用几何类型转换作为示例：

multilinestring_list = list(matrix(c(1, 4, 5, 3), ncol = 2), 
                            matrix(c(4, 4, 4, 1), ncol = 2),
                            matrix(c(2, 4, 2, 2), ncol = 2))
multilinestring = st_multilinestring(multilinestring_list)
multilinestring_sf = st_sf(geom = st_sfc(multilinestring))
multilinestring_sf
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 1 ymin: 1 xmax: 4 ymax: 5
#> CRS:           NA
#>                             geom
#> 1 MULTILINESTRING ((1 5, 4 3)...

你可以将其想象成一条道路或河流网络。新对象只有一行，定义了所有线条。这限制了可以执行的操作数量，例如阻止为每个线段添加名称或计算单条线的长度。在这种情况下，可以使用st_cast()函数，因为它将一条多段线分成三条线：

linestring_sf2 = st_cast(multilinestring_sf, "LINESTRING")
linestring_sf2
#> Simple feature collection with 3 features and 0 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 1 ymin: 1 xmax: 4 ymax: 5
#> CRS:           NA
#>                    geom
#> 1 LINESTRING (1 5, 4 3)
#> 2 LINESTRING (4 4, 4 1)
#> 3 LINESTRING (2 2, 4 2)

Examples of type casting between MULTILINESTRING (left) and LINESTRING (right).

新创建的对象允许创建属性和测量长度：

linestring_sf2$name = c("Riddle Rd", "Marshall Ave", "Foulke St")
linestring_sf2$length = st_length(linestring_sf2)
linestring_sf2
#> Simple feature collection with 3 features and 2 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 1 ymin: 1 xmax: 4 ymax: 5
#> CRS:           NA
#>                    geom         name length
#> 1 LINESTRING (1 5, 4 3)    Riddle Rd   3.61
#> 2 LINESTRING (4 4, 4 1) Marshall Ave   3.00
#> 3 LINESTRING (2 2, 4 2)    Foulke St   2.00

栅格数据几何操作

几何栅格操作包括图像的平移、翻转、镜像、缩放、旋转或扭曲。这些操作在许多应用中是必要的，包括地理参考，用于让图像能够叠加在具有已知坐标参考系统（CRS）的准确地图上。存在多种地理参考技术，包括：

• 基于已知 ground control points的地理校正
• 正射校正，还考虑了局部地形
• 图像registration 用于通过使一幅图像与另一幅图像对齐（在坐标系统和分辨率方面）来组合拍摄同一物体但使用不同传感器拍摄的图像

对于前两点，R 语言相当不适合，因为这些操作通常需要手动干预，这就是它们通常是在专用 GIS 软件的帮助下完成的原因。另一方面，在R中对几幅图像进行对齐是可能的，本节将展示如何执行此操作。这通常包括更改图像的范围、分辨率和原点。当然，也需要匹配的投影。

无论如何，有其他原因需要对单个栅格图像执行几何操作。例如，我们将德国的都市区定义为面积超过 20 平方千米、人口超过500,000的像素。然而，原始居民栅格的分辨率为1平方千米，这就是我们将分辨率降低（聚合）20倍的原因。聚合栅格的另一个原因就是简单地减少运行时间或节省磁盘空间。当然，只有当手头的任务允许栅格数据的粗分辨率时，才推荐使用此方法。

几何相交

我们展示了如何从由其他空间对象覆盖的栅格中提取值。为了检索空间输出，我们可以使用几乎相同的子集提取语法。唯一的区别是必须明确表示我们希望把 drop 参数设置为 FALSE 来保持矩阵结构。这将返回一个包含与clip重叠中点的单元格的栅格对象。

elev = rast(system.file("raster/elev.tif", package = "spData"))
clip = rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45,
            resolution = 0.3, vals = rep(1, 9))
elev[clip, drop = FALSE]
#> class       : SpatRaster 
#> dimensions  : 2, 1, 1  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : 1, 1.5, -0.5, 0.5  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> name        : elev 
#> min value   :   18 
#> max value   :   24

对于相同的操作，我们还可以使用intersect()和crop()命令。

范围和原点

当合并栅格或在栅格上执行地图代数时，它们的分辨率、投影、原点和/或范围必须匹配。否则，我们如何将一个分辨率为0.2（十进制）度的栅格的值与分辨率为1（十进制）度的第二个栅格的值相加呢？当我们想要合并来自具有不同投影和分辨率的不同传感器的卫星图像时，会出现相同的问题。我们可以通过对齐栅格来处理不匹配的数据。

在最简单的情况下，两个图像只在其范围上有所不同。下面的代码在栅格的每一侧添加一行和两列，同时将所有新值设置为NA。

elev = rast(system.file("raster/elev.tif", package = "spData"))
elev_2 = extend(elev, c(1, 2))

Original raster (left) and the same raster (right) extended by one row on the top and bottom and two columns on the left and right.

在R中对具有不同范围的两个对象执行代数运算时，terra包会返回一个错误。

elev_3 = elev + elev_2
#> Error: [+] extents do not match

当然，我们可以利用extend()来使两个栅格的范围相匹配。我们不是告诉函数要添加多少行或列，而是通过使用另一个栅格对象让它自行判断。在这里，我们将 elev 对象扩展到与 elev_2 的范围相同。新添加的行和列的值会被设置为 NA。

elev_4 = extend(elev, elev_2)

栅格的起始点是离坐标（0,0）最近的单元格角。origin() 函数返回起始点的坐标。在下面的例子中，存在一个坐标为（0,0）的单元格角，但这并不一定是起始点。

origin(elev_4)
#> [1] 0 0

如果两个栅格具有不同的起始点，它们的单元格将不会完全重叠，从而使地图代数运算变得不可能。要更改起始点，请使用 origin()。[^3]图揭示了以这种方式更改起始点的效果。

[^3]: 如果两个栅格数据集的起始点仅略有不同，有时简单地增加 terra::terraOptions() 的 tolerance 参数就足够了。

# change the origin
origin(elev_4) = c(0.25, 0.25)

Rasters with identical values but different origins.

注意，频繁改变分辨率，（下一节）也会改变原点。

聚合和离散

栅格数据集在分辨率方面也可能有所不同。为了匹配分辨率，可以减小（使用aggregate()）或增加（使用disagg()）其中一个栅格的分辨率。^4 作为示例，我们在此通过一个因子5改变dem（在spDataLarge 包中找到）的空间分辨率。此外，输出单元格值应与输入单元格的平均值相对应（注意也可以使用其他函数，例如 median()、sum() 等）。

dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))
dem_agg = aggregate(dem, fact = 5, fun = mean)

Original raster (left). Aggregated raster (right).

disagg() 函数用于增加栅格对象的分辨率。它提供了两种计算新创建单元格值的方法：默认方法（method = "near"）仅将输入单元格的值赋给所有输出单元格，从而复制值，导致输出呈“块状”。bilinear 方法使用输入图像的四个最近像素中心（图中的鲑鱼色点）来按距离计算加权平均值（图中的箭头）。输出单元格的值由图左上角的方块表示。

dem_disagg = disagg(dem_agg, fact = 5, method = "bilinear")
identical(dem, dem_disagg)
#> [1] FALSE

The distance-weighted average of the four closest input cells determine the output when using the bilinear method for disaggregation.

比较dem和dem_disagg的值，我们可以发现它们并不相同（您也可以使用compareGeom() 或all.equal()）。然而，这几乎是意料之中的，因为细分是一种简单的插值技术。重要的是要记住，细分会导致更精细的分辨率；然而，相应的值只能与其较低分辨率的源数据一样准确。

重采样

上述聚合和细分的方法仅适用于我们想要通过聚合/细分因子来改变栅格的分辨率的情况。然而，当我们有两个或更多具有不同分辨率和原点的栅格时该怎么办呢？这就是重采样的作用 – 一个计算新像素位置值的过程。简而言之，这个过程采取我们原始栅格的值，并为具有自定义分辨率和原点的目标栅格重新计算新值。

Resampling of an original (input) raster into a target raster with custom resolution and origin.

有几种方法可以用来估计具有不同分辨率/原点的栅格的值。主要的重采样方法包括：

• 最近邻：将原始栅格的最近单元格的值分配给目标栅格的单元格。这是一种快速简单的技术，通常适用于对分类栅格进行重采样。
• 双线性插值：将原始栅格的四个最近单元格的加权平均值分配给目标栅格的单元格（见图@ref(fig:bilinear)）。这是适用于连续栅格的最快方法。
• 三次插值：使用原始栅格的16个最近单元格的值来确定输出单元格的值，应用三次多项式函数。用于连续栅格，并且与双线性插值相比，结果表面更平滑，但计算上更耗费资源。
• 三次样条插值：也使用原始栅格的16个最近单元格的值来确定输出单元格的值，但应用三次样条（分段三次多项式函数）。用于连续栅格。
• Lanczos窗口化sinc重采样：使用原始栅格的36个最近单元格的值来确定输出单元格的值。用于连续栅格。[^5]

[^5]: 关于这种方法的更详细解释可以在 https://gis.stackexchange.com/a/14361/20955 找到。

上述解释强调，只有最近邻重采样适用于分类栅格，而所有方法都可以用于连续栅格（结果不同）。还请注意，这些方法从上到下在复杂性和处理时间上都有所增加。

要应用重采样，terra 包提供了一个 resample() 函数。它接受一个输入栅格（x），一个具有目标空间属性的栅格（y），以及一个重采样方法（method）。

我们需要一个具有目标空间属性的栅格来了解 resample() 函数是如何工作的。在这个例子中，我们创建了target_rast，但通常你会使用一个已经存在的栅格对象。

target_rast = rast(xmin = 794650, xmax = 798250, 
                   ymin = 8931750, ymax = 8935350,
                   resolution = 300, crs = "EPSG:32717")

接下来，我们需要提供两个栅格对象作为前两个参数，并选择上述重采样方法中的一个。

dem_resampl = resample(dem, y = target_rast, method = "bilinear")

下图展示了对dem对象应用不同重采样方法的比较。

Visual comparison of the original raster and five different resampling methods.

resample() 函数还有一些额外的重采样方法，包括sum、min、q1、med、q3、max、average、mode和rms。所有这些方法都基于所有非NA贡献网格单元的值计算给定的统计量。例如，当每个栅格单元代表一个空间广泛的变量（例如，人数）时，sum是有用的。使用sum的效果是，重采样后的栅格应与原始栅格具有相同的人口总数。

当我们的目标栅格与原始栅格具有不同的坐标参照系统（CRS）时，栅格重投影是重采样的特殊情况。

📌terra中的大多数几何操作都是用户友好的，相当快速，并可用于大型栅格对象。然而，可能存在一些情况，当涉及到广泛的栅格或许多栅格文件时，terra可能不是最高效的，应该考虑一些替代方案。

最成熟的替代方案来自GDAL库。它包括几个实用功能，包括：

• gdalinfo - 列出有关栅格文件的各种信息，包括其分辨率、CRS、边界框等。
• gdal_translate - 在不同的文件格式之间转换栅格数据。
• gdal_rasterize - 将矢量数据转换为栅格文件。
• gdalwarp - 允许栅格镶嵌、重采样、裁剪和重投影。

以上所有函数都是用C++编写的，但可以使用sf::gdal_utils()、gdalUtilities包或通过系统命令在R中调用。重要的是，所有这些函数都期望栅格文件路径作为输入，并经常将其输出作为栅格文件返回（例如，gdalUtilities::gdal_translate("my_file.tif", "new_file.tif", t_srs = "EPSG:4326")）。这与通常的terra方法非常不同，后者期望SpatRaster对象作为输入。

练习

E1. 生成并绘制 nz 数据集的简化版本。使用 ms_simplify() 的不同 keep 值（范围从0.5到0.00005）和 st_simplify() 的不同 dTolerance 值（从100到100,000）进行实验。

• 对于每种方法，从什么值开始，结果的形态开始崩溃，使新西兰变得不可识别？
• 高级: st_simplify() 的结果几何类型与 ms_simplify() 的几何类型有什么不同？这会造成什么问题，如何解决？

E2. 在空间数据操作章节的第一个练习中，已经确定坎特伯雷地区拥有新西兰101个最高点中的70个。使用 st_buffer()，nz_height 中有多少点在坎特伯雷的100 km内？

E3. 找到新西兰的地理中心。它距离坎特伯雷的地理中心有多远？

E4. 大多数世界地图都是朝北的。通过反射 world 对象的几何体（本章未提及的仿射变换之一），可以创建一个朝南的世界地图。编写代码来实现这一点。提示: 你需要使用一个两元素向量来进行此转换。加分题: 创建一个倒立的国家地图。

E5. 运行5.2.6章节中的代码。参考该章节中创建的对象，找出包含在 x 和 y 中的 p 中的点。

• 使用基本子集运算符。
• 使用 st_intersection() 创建的中间对象。

E6. 以米为单位计算美国各州的边界线长度。哪个州的边界最长，哪个州的边界最短？提示: st_length 函数计算 LINESTRING 或 MULTILINESTRING 几何的长度。

E7. 将 srtm.tif 文件读入 R (srtm = rast(system.file("raster/srtm.tif", package = "spDataLarge")))。此栅格的分辨率为0.00083乘以0.00083度。使用 terra 包中提供的所有方法，将其分辨率更改为0.01乘以0.01度。可视化结果。你能注意到这些重采样方法的结果之间有任何区别吗？