Finding minimum distance between two raster layer pixels in R

Question

I have two thematic raster layers r1 and r2 for same area each following same classification scheme and has 16 classes. I need to find minimum distance between cell of r1 and cell of r2 but with same value. E.g. nth cell in r1 has value 10 and coordinates x1,y1. And in r2, there are 2 cells with value 10 and coordinates x1+2,y1+2 and x1-0.5,y1-0.5. Thus the value that I need for this cell would be 0.5,0.5.

I tried distance from raster package but it gives distance, for all cells that are NA, to the nearest cell that is not NA. I am confused as to how can I include second raster layer into this.

Answer 1

You can use knn from class package so that for each cell of r1 find index of nearest cell of r2 with the same category:

library(class)
library(raster)
#example of two rasters
r1 <- raster(ncol =  600, nrow = 300)
r2 <- raster(ncol =  600, nrow = 300)
#fill each with categories that rabge from 1 to 16
r1[] <- sample(1:16, ncell(r1), T)
r2[] <- sample(1:16, ncell(r2), T)
# coordinates of cells extracted
xy = xyFromCell(r1, 1:ncell(r1))
#multiply values of  raster with a relatively large number so  cells thet belong
#to each category have smaller distance with reagrd to other categories.
v1 = values(r1) *  1000000
v2 = values(r2) *  1000000
# the function returns indices of nearest cells
out = knn(cbind(v2, xy) ,cbind(v1, xy) ,1:ncell(r1), k=1)

Answer 2

So, use rasterToPoints to extract SpatialPoints object for unique thematic class. Then use the sp::spDists function to find the distance between your points.

library(raster)


r1 <- raster( nrow=10,ncol=10)
r2 <- raster( nrow=10,ncol=10)

set.seed(1)
r1[] <- ceiling(runif(100,0,10))
r2[] <- ceiling(runif(100,0,10))

dist.class <- NULL
for(i in unique(values(r1))){
p1 <- rasterToPoints(r1, fun=function(xx) xx==i, spatial=T)
p2 <- rasterToPoints(r2, fun=function(xx) xx==i, spatial=T)
dist.class[i] <- min(spDists(p1,p2))
}
cbind(class = unique(values(r1)),dist.class)

The loop may not be efficient for you. If it's a problem, wrap it into a function and lapply it. Also, be carefull with your class, if they aren't 1:10, my loop won't work. If your projection is in degree, you will probably need the geosphere package to get accurate results. But the best in that case I think is to use a projection in meters.

Answer 3

A memory safe approach using the raster-package would be to use the layerize() function to split up your raster value into a stack of binary rasters (16 in your case) and then use the distance() function to compute distances in the layers of r2, masking them with the respective layers of r1. Something like this:

layers1 <- layerize(r1, falseNA=TRUE)
layers2 <- layerize(r2, falseNA=TRUE)

# now you can loop over the layers (use foreach loop if you want 
# to speed things up using parallel processing)

dist.stack <- layers1

for (i in 1:nlayers(r1)) {
    dist.i <- distance(layers2[[i]])
    dist.mask.i <- mask(dist, layers1[[i]])
    dist.stack[[i]] <- dist.mask.i
}

# if you want pairwise distances for all classes in one layer, simply
# combine them using sum()

dist.combine <- sum(dist.stack, na.rm=TRUE)