使用 R 测量每个补丁之间的距离
Measuring distances between every patch using R
我一直在尝试找到一个快速函数来使用 R 同时测量栅格中多个补丁之间的距离。特别是，我想测量到每个方向上最近的补丁的距离（不仅是最近的一）。由于在每个方向上识别最近的一个可能很耗时，因此到每个补丁的距离也可以解决这个问题。
首先我使用了 gDistance，但结果并不直观（请参见下面的示例）。特别是，很难将电导与光栅中正方形之间的实际距离联系起来。
然后我尝试使用栅格包对每个补丁进行迭代，测量从那里到补丁外每个像素的距离（使用函数距离），然后寻找到每个其他补丁的最小距离。它有效，但非常耗时。这也非常低效，因为我对每个距离都测量了两次，而且因为我也在测量不需要的距离（如果从补丁 A 到补丁 C 的唯一方法穿过补丁 B，我不需要补丁 A 和补丁之间的距离C）。
下面是我正在使用的代码...
感谢任何建议...
卡洛斯·阿尔贝托
gdistance 代码
library(gdistance)
library(raster)
# setting the patches
mF <- raster(nrows=10, ncols=20)
mF[] <- 0; mF[4:8,3] <- 2; mF[9,14:18] <- 3; mF[3,9:12] <- 1
# and the cost function
mg <- mF <= 0
# and the transition matrix
tr1 <- transition(1/mg, transitionFunction=mean, directions=16)
tr1C <- geoCorrection(tr1, type="c")
# getting coordinates of sampling points
dF1 <- as.data.frame(mF,xy=T, na.rm=T)
dF2 <- dF1[!duplicated(dF1[,3]),]
dF3 <- as.matrix(dF2[,1:2])
rownames(dF3) <- dF2[,3]
# and measuring the cost distance
cbind(dF3, as.matrix(costDistance(tr1C,dF3)))
给这个：
x y 0 1 2 3
0 -171 81 0 2192567 2079216.3 2705664.0
1 -27 45 2192567 0 2727389.7 3353837.4
2 -135 27 2079216 2727390 0.0 626447.7
3 63 -63 2705664 3353837 626447.7 0.0
重点关注点：1.这个值是什么意思？如何将它们与km联系起来？ 2.为什么到0class的距离会随着纬度的增加而增加？如果像素的面积在靠近两极的地方减小，那么到每个图外的距离也应该减小。
光栅编码
region <- (mF > 0) + 0 # the landscape map.
# this is just to reduce the creation/destruction of the variables
patch <- region
biome <- region
biomes <- unique(region)
areas <- area(region)
map.distances <-function (i) {
dA <- data.frame(biome = integer(0),
patch = integer(0),
area = numeric(0))
dD <- data.frame(biome = integer(0),
from = integer(0),
to = integer(0),
dist = numeric(0))
biome[] <- NA_integer_
# creating the patches
biome[region == i] <- i
biomeC <- clump(biome, directions=8)
dA <- rbind(dA, cbind(biome = i,
zonal(areas, biomeC, 'sum')))
patches <- as.integer(unique(biomeC))
# in each patch...
for (j in patches[-1]) {
patch[] <- NA_integer_
patch[biomeC == j] <- 1L
# get the distances from the patch
dists <- distance(patch)
d <- zonal(dists, biomeC, "min")
f <- j > d[,1]
# and combine the info
dD <- rbind(dD, data.frame(from = j,
to = d[f,1],
dist = d[f,2], biome = i))
}
return(list(edges=dD, vertices=dA))
}
# applying to the same map as before gives:
mpd <- map.distances(i=1)
rownames(mpd$edges) <- NULL
mpd$edges
from to dist biome
1 2 1 9210860 1
2 3 1 12438366 1
3 3 2 5671413 1
看到距离不是线性相关的。
这是另一种 raster 的方法， 可能 更有效（但对于真实（大型）数据集可能有问题）。我正在使用更简单的（平面）光栅来更好地理解结果：
library(raster)
# setting the patches
mF <- raster(nrows=10, ncols=20, xmn=0, xmx=20, ymn=0, ymx=10, crs='+proj=utm +zone=10 +datum=WGS84')
mF[] <- 0; mF[4:8,3] <- 2; mF[9,14:18] <- 3; mF[3,9:12] <- 1
dF <- as.data.frame(mF, xy=TRUE, na.rm=TRUE)
dF <- dF[dF[,3] > 0,]
pd <- pointDistance(dF[,1:2], lonlat=FALSE)
pd <- as.matrix(as.dist(pd))
diag(pd) <- NA
a <- aggregate(pd, dF[,3,drop=FALSE], min, na.rm=TRUE)
a <- t(a[,-1])
a <- aggregate(a, dF[,3,drop=FALSE], min, na.rm=TRUE)[, -1]
diag(a) <- 0
a
# V1 V2 V3
#1 0.000000 6.082763 6.324555
#2 6.082763 0.000000 11.045361
#3 6.324555 11.045361 0.000000
这里有 gdistance:
# I think the cost is the same everywhere:
mg <- setValues(mF, 1)
# and the transition matrix
tr1 <- transition(1/mg, transitionFunction=mean, directions=16)
tr1C <- geoCorrection(tr1, type="c")
cd <- as.matrix(costDistance(tr1C, as.matrix(dF[,1:2])))
b <- aggregate(cd, dF[,3,drop=FALSE], min, na.rm=TRUE)
b <- t(b[,-1])
b <- aggregate(b, dF[,3,drop=FALSE], min, na.rm=TRUE)[, -1]
diag(b) <- 0
b
# V1 V2 V3
#1 0.000000 6.236068 6.472136
#2 6.236068 0.000000 11.236068
#3 6.472136 11.236068 0.000000
您可以通过仅使用补丁的边界来减少要处理的点数。考虑：
mF[] <- NA; mF[4:8,3] <- 2; mF[9,14:18] <- 3; mF[3,9:12] <- 1
mF[1:5, 1:5] = 2
plot(boundaries(mF, type='inner'))
您也可以先创建多边形
mF[] <- NA; mF[4:8,3] <- 2; mF[9,14:18] <- 3; mF[3,9:12] <- 1
p <- rasterToPolygons(mF, dissolve=TRUE)
gDistance(p, byid=T)
# 1 2 3
#1 0.00000 5 5.09902
#2 5.00000 0 10.00000
#3 5.09902 10 0.00000
我一直在尝试找到一个快速函数来使用 R 同时测量栅格中多个补丁之间的距离。特别是，我想测量到每个方向上最近的补丁的距离（不仅是最近的一）。由于在每个方向上识别最近的一个可能很耗时，因此到每个补丁的距离也可以解决这个问题。
首先我使用了 gDistance，但结果并不直观（请参见下面的示例）。特别是，很难将电导与光栅中正方形之间的实际距离联系起来。
然后我尝试使用栅格包对每个补丁进行迭代，测量从那里到补丁外每个像素的距离（使用函数距离），然后寻找到每个其他补丁的最小距离。它有效，但非常耗时。这也非常低效，因为我对每个距离都测量了两次，而且因为我也在测量不需要的距离（如果从补丁 A 到补丁 C 的唯一方法穿过补丁 B，我不需要补丁 A 和补丁之间的距离C）。 下面是我正在使用的代码...
感谢任何建议...
卡洛斯·阿尔贝托
gdistance 代码
library(gdistance)
library(raster)
# setting the patches
mF <- raster(nrows=10, ncols=20)
mF[] <- 0; mF[4:8,3] <- 2; mF[9,14:18] <- 3; mF[3,9:12] <- 1
# and the cost function
mg <- mF <= 0
# and the transition matrix
tr1 <- transition(1/mg, transitionFunction=mean, directions=16)
tr1C <- geoCorrection(tr1, type="c")
# getting coordinates of sampling points
dF1 <- as.data.frame(mF,xy=T, na.rm=T)
dF2 <- dF1[!duplicated(dF1[,3]),]
dF3 <- as.matrix(dF2[,1:2])
rownames(dF3) <- dF2[,3]
# and measuring the cost distance
cbind(dF3, as.matrix(costDistance(tr1C,dF3)))
给这个：
x y 0 1 2 3
0 -171 81 0 2192567 2079216.3 2705664.0
1 -27 45 2192567 0 2727389.7 3353837.4
2 -135 27 2079216 2727390 0.0 626447.7
3 63 -63 2705664 3353837 626447.7 0.0
重点关注点：1.这个值是什么意思？如何将它们与km联系起来？ 2.为什么到0class的距离会随着纬度的增加而增加？如果像素的面积在靠近两极的地方减小，那么到每个图外的距离也应该减小。
光栅编码
region <- (mF > 0) + 0 # the landscape map.
# this is just to reduce the creation/destruction of the variables
patch <- region
biome <- region
biomes <- unique(region)
areas <- area(region)
map.distances <-function (i) {
dA <- data.frame(biome = integer(0),
patch = integer(0),
area = numeric(0))
dD <- data.frame(biome = integer(0),
from = integer(0),
to = integer(0),
dist = numeric(0))
biome[] <- NA_integer_
# creating the patches
biome[region == i] <- i
biomeC <- clump(biome, directions=8)
dA <- rbind(dA, cbind(biome = i,
zonal(areas, biomeC, 'sum')))
patches <- as.integer(unique(biomeC))
# in each patch...
for (j in patches[-1]) {
patch[] <- NA_integer_
patch[biomeC == j] <- 1L
# get the distances from the patch
dists <- distance(patch)
d <- zonal(dists, biomeC, "min")
f <- j > d[,1]
# and combine the info
dD <- rbind(dD, data.frame(from = j,
to = d[f,1],
dist = d[f,2], biome = i))
}
return(list(edges=dD, vertices=dA))
}
# applying to the same map as before gives:
mpd <- map.distances(i=1)
rownames(mpd$edges) <- NULL
mpd$edges
from to dist biome
1 2 1 9210860 1
2 3 1 12438366 1
3 3 2 5671413 1
看到距离不是线性相关的。
这是另一种 raster 的方法， 可能 更有效（但对于真实（大型）数据集可能有问题）。我正在使用更简单的（平面）光栅来更好地理解结果：
library(raster)
# setting the patches
mF <- raster(nrows=10, ncols=20, xmn=0, xmx=20, ymn=0, ymx=10, crs='+proj=utm +zone=10 +datum=WGS84')
mF[] <- 0; mF[4:8,3] <- 2; mF[9,14:18] <- 3; mF[3,9:12] <- 1
dF <- as.data.frame(mF, xy=TRUE, na.rm=TRUE)
dF <- dF[dF[,3] > 0,]
pd <- pointDistance(dF[,1:2], lonlat=FALSE)
pd <- as.matrix(as.dist(pd))
diag(pd) <- NA
a <- aggregate(pd, dF[,3,drop=FALSE], min, na.rm=TRUE)
a <- t(a[,-1])
a <- aggregate(a, dF[,3,drop=FALSE], min, na.rm=TRUE)[, -1]
diag(a) <- 0
a
# V1 V2 V3
#1 0.000000 6.082763 6.324555
#2 6.082763 0.000000 11.045361
#3 6.324555 11.045361 0.000000
这里有 gdistance:
# I think the cost is the same everywhere:
mg <- setValues(mF, 1)
# and the transition matrix
tr1 <- transition(1/mg, transitionFunction=mean, directions=16)
tr1C <- geoCorrection(tr1, type="c")
cd <- as.matrix(costDistance(tr1C, as.matrix(dF[,1:2])))
b <- aggregate(cd, dF[,3,drop=FALSE], min, na.rm=TRUE)
b <- t(b[,-1])
b <- aggregate(b, dF[,3,drop=FALSE], min, na.rm=TRUE)[, -1]
diag(b) <- 0
b
# V1 V2 V3
#1 0.000000 6.236068 6.472136
#2 6.236068 0.000000 11.236068
#3 6.472136 11.236068 0.000000
您可以通过仅使用补丁的边界来减少要处理的点数。考虑：
mF[] <- NA; mF[4:8,3] <- 2; mF[9,14:18] <- 3; mF[3,9:12] <- 1
mF[1:5, 1:5] = 2
plot(boundaries(mF, type='inner'))
您也可以先创建多边形
mF[] <- NA; mF[4:8,3] <- 2; mF[9,14:18] <- 3; mF[3,9:12] <- 1
p <- rasterToPolygons(mF, dissolve=TRUE)
gDistance(p, byid=T)
# 1 2 3
#1 0.00000 5 5.09902
#2 5.00000 0 10.00000
#3 5.09902 10 0.00000