本篇推文翻译自Manuel Gimond的在线教程Intro to GIS and Spatial Analysis的附录F Raster operations in R，网址为https://mgimond.github.io/Spatial/raster-operations-in-r.html。原教程遵循CC BY-NC 4.0共享协议（https://creativecommons.org/licenses/by-nc/4.0/deed.zh）。

原文篇幅较长，共有6小节，这里翻译的是第5、6节。

5 Global operations and functions

全局操作符和函数在计算输出对象的像元值时，会使用输入对象中的所有像元。

全局函数的一个例子是计算欧几里得距离的distance()函数，它的功能是计算每个像元与对象最近位置的距离。为了描述distance()函数的功能，我们首先定义一个具有两个非空像元的栅格。

library(raster)
r1   <- raster(ncol = 100, nrow=100, xmn = 0, xmx = 100,ymn = 0, ymx = 100)
r1[] <- NA              # Assign NoData values to all pixels
r1[c(850, 5650)] <- 1   # Change the pixels #850 and #5650  to 1
crs(r1) <- "+proj=ortho"  # Assign an arbitrary coordinate system (needed for mapping with tmap)library(tmap)
tm_shape(r1) + tm_raster(palette="red") + tm_legend(outside = TRUE, text.size = .8)

接下来，我们来计算到这两个非空像元的距离，输出结果是栅格对象。默认情况下，它的范围与输入栅格对象相同。

r1.d <- distance(r1)

• 译者注：raster工具包中的distance(x, y,...)函数，当输入对应只有一个栅格对象x（y缺省）时，它计算的是x中每个空值像元与最近的非空像元之间的距离。

tm_shape(r1.d) + tm_raster(palette = "Greens", style = "order", title = "Distance") + tm_legend(outside = TRUE, text.size = .8) +tm_shape(r1) + tm_raster(palette = "red", title = "Points")

你也可以使用SpatialPoints类型的对象或简单的x/y表来得到一个距离栅格。接下来的例子中，我们计算的是输出栅格对象每个像元与点(25,30)和(87,80)之间的距离。但是，我们现在要使用的是点对象（而不是前面例子中使用的栅格对象），因此我们需要创建一个新的栅格对象，来定义输出栅格对象的范围。

# Create a blank raster
r2 <- raster(ncol = 100, nrow = 100, xmn = 0, xmx = 100, ymn = 0, ymx = 100)
crs(r2) <- "+proj=ortho"  # Assign an arbitrary coordinate system # Create a point layer
xy <- matrix(c(25,30,87,80), nrow = 2, byrow = T)
p1 <- SpatialPoints(xy)
crs(p1) <- "+proj=ortho"  # Assign an arbitrary coordinate system

现在，让我们使用distanceFromPoints()函数来计算到这两个点的欧几里得距离。

r2.d <- distanceFromPoints(r2, p1)

让我们来绘制输出结果。

tm_shape(r2.d) + tm_raster(palette = "Greens", style = "order") + tm_legend(outside = TRUE, text.size = .8) +tm_shape(p1) + tm_bubbles(col = "red")

6 Computing cumulative distances

这个练习将展示如何使用gdistance工具包中的函数得到累积距离栅格（cumulative distance raster）。目标之一就是说明“邻近像元”（adjacency cells）的定义对计算结果的影响。

加载gdistance工具包。

library(gdistance)

首先，我们定义一个100*100的栅格，并将每个像元的值设定为1。这个值表示穿过该像元的旅行成本（而不是距离本身）。在这个例子中，我们假定成本在整个范围内是相同的。

r   <- raster(nrows = 100, ncols = 100,xmn = 0, ymn = 0,xmx = 100, ymx = 100)
r[] <- rep(1, ncell(r))

如果想包括除距离以外的旅行成本（如海拔），你应该分别定义每个像元的值，而不是使用常量1。

可以使用转移矩阵（translation matrix）来定义邻近像元之间的旅行成本。因为我们假定从一个像元移动到它邻近的像元，除距离外没有其他成本，因此可以使用function(x) {1}来表示将一个像元转移到它邻近的像元上（也就是说，各个方向的转移成本是相同的）。

在transition()函数中，有四种定义“邻近”的方式。接下来的四段代码分别对其进行演示。

在第一个例子中，邻近被定义成四节点（水平和垂直）连接（也就是"rook"型）。

h4 <- transition(r, transitionFunction = function(x){1},directions = 4)

第二个例子中，邻接被定义成八节点连接（也就是“queen”型）。

h8 <- transition(r, transitionFunction = function(x){1},directions = 8)

第三个例子中，邻近是十六节点连接（也就是“queen”和“knight”相结合型）。

h16 <- transition(r, transitionFunction=function(x){1},16, symm=FALSE)

第四个例子中，邻近是四节点对角线连接（也就是“bishop”型）。

hb <- transition(r, transitionFunction = function(x){1},"bishop", symm=FALSE)

transition()函数将所有邻近像元到源像元的距离看作是相同的。geoCorrection()函数则将其纠正为“真实”的局部距离。实际上，它是给一个像元到它邻近的像元添加了额外的成本（原始距离成本是通过transition()定义的）。纠正的必要性将在后文呈现。

注：geoCorrection()函数还可以校正由于地理坐标系统导致的距离失真，但请确保已经使用projection()函数定义了图层的坐标系统。

h4 <- geoCorrection(h4,  scl=FALSE)
h8 <- geoCorrection(h8,  scl=FALSE)
h16 <- geoCorrection(h16, scl=FALSE)
hb <- geoCorrection(hb,  scl=FALSE)

在“queen”型下，对角线方向的邻近像元到源像元的距离是像元宽距的倍。

接下来，我们将使用四种邻近的定义分别计算从中心点A到栅格所有像元的累积距离（accCost()函数）。

A <- c(50,50) # Location of source cell
h4.acc <- accCost(h4, A)
h8.acc <- accCost(h8, A)
h16.acc <- accCost(h16, A)
hb.acc <- accCost(hb, A)

如果在前面的步骤中没有使用geoCorrection()函数，那么累积距离的计算结果会有所不同。下面两幅图展示了在十六节点连接的邻近定义下二者的区别。

未校正（也就是没有对h16使用geoCorrection()函数）：

校正（也就是对h16使用了geoCorrection()函数）：

在“bishop”型情况下还会出现一个特殊的问题。因为只有对角线方向被定义成邻近像元，从而会存在一些与源像元不存在任意阶邻近关系的像元，它们在计算结果中值被标记为无穷（Inf）。我们将把无穷值像元改成空值像元。

hb.acc[hb.acc == Inf] <- NA

接下来，让我们比较以A为中心、7*7范围内的计算结果。

为了强调四种计算结果的差异，我们将与A点距离在20单位以内的像元标记为红色。

很显然，累积距离栅格的精度极大地受邻近定义的影响。红色像元（也就是累积距离在20单位以内）的个数从925到2749不等。