Измерение расстояний между каждым патчем с использованием R

Question

Я пытаюсь найти быструю функцию для измерения расстояния между несколькими патчами в растре одновременно с использованием R. В частности, я хочу измерить расстояние до ближайшего патча во всех направлениях (не только к ближайшему). Поскольку определение ближайшего в каждом направлении может занять много времени, расстояние до каждого патча также может решить проблему.

Сначала я использовал gDistance, но результаты не были интуитивными (см. Пример ниже). В частности, трудно связать проводимость с фактическим расстоянием между квадратами в растре.

Затем я попробовал с растровыми пакетами, итерациями каждого патча, измеряя расстояние оттуда до каждого пикселя из патча (используя расстояние между функциями), а затем ищем минимальное расстояние до каждого другого патча. Он работает, но это занимает много времени. Это также крайне неэффективно, потому что я измеряю каждое расстояние дважды, и потому что я измеряю ненужные расстояния (если единственный способ перейти от патча A к патчу C пересекает патч B, мне не нужно расстояние между патчем A и С). Ниже приведен код, я использую ...

Спасибо за любые советы ...

Карлос Альберто

код 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. Что означает значение? как связать их с км? 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