使用R測量每個修補程序之間的距離

Question

我一直在試圖找到一個快速函數來使用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.爲什麼到0級的距離隨着緯度增加？如果像素面積減小到接近極點，則每個圖形外部的距離也應該減小。

光柵碼

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 

看到距離不是線性相關。

Answer 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