本篇推文翻译自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不等。