I have been trying to find a quick function to measure the distance between several patches in a raster simultaneously using R. In particular, I want to measure the distance to the closest patch in every direction (not only to the closest one). Since identifying the closest one in each direction can be time consuming, the distance to every single patch could also solve the problem.
First I used gDistance, but the results were not intuitive (see example below). In particular, it is hard to relate the conductance with the actual distance between squares in the raster.
Then I tried with the raster packages iterating by every patch, measuring the distance from there to every single pixel out of the patch (using the function distance) and then looking for the minimum distance to every single other patch. It works, but it's extremely time consuming. It is also extremely inefficient because I'm measuring every distance twice and because I'm measuring unneeded distances also (if the only way to go from patch A to patch C crosses patch B, I don't need the distance between patch A and C). Below is the code I'm using...
Thanks for any advice...
Carlos Alberto
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)))
give this:
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
Key concerns: 1. what does the value mean? how to relate them to km? 2. why the distance to the 0 class increase with the latitude? if the area of the pixels reduce closer to the poles, the distance to outside each plot should decrease also.
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
see that the distances are not linearly correlated.
Here is another approach with raster, perhaps more efficient (but probably problematic for a real (large) data set). I am using a simpler (planar) raster to better understand the result:
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
And here with 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
You could reduce the number of points to work with by only using the boundaries of the patches. Consider:
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'))
You could also create polygons first
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