Calculating distance between each pixel of rasterstack image?

Question

I am new on R and I have an image of class "rasterstack" of only 3 colors: Red, Green and blue as illustrated:

I use this function: image = stack(img.red, img.green, img.blue)

I would like to find the distance between each pixel of one color to each others pixel of another color. By example the distance between each red pixel, to each blue pixel:

• Distance between red pixel 1 to blue pixel 1 ,pixel 2 ,pixel 3...
• Distance between red pixel 2 to blue pixel 1 ,pixel 2 ,pixel 3...
• etc

And create a table with all these distances. Then I could do some statistic as to know the distribution of the distance, the average min distance between red and blue, etc

I have no idea how to do it?

Answer 1

If you have a raster red, in this case 35 cells in a 5x7 raster:

> red
class       : RasterLayer 
dimensions  : 5, 7, 35  (nrow, ncol, ncell)
resolution  : 0.1428571, 0.2  (x, y)
extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : in memory
names       : layer 
values      : 0, 255  (min, max)

You can get a vector of the values:

> red[]
 [1]   0   0 255   0   0 255   0 255   0   0   0   0   0   0   0   0 255   0   0
[20]   0   0   0   0 255   0   0   0   0   0   0   0   0 255   0 255

and you can convert all 35 cell centres to SpatialPoints:

> as(red,"SpatialPoints")
class       : SpatialPoints 
features    : 35 
extent      : 0.07142857, 0.9285714, 0.1, 0.9  (xmin, xmax, ymin, ymax)
coord. ref. : NA 

To then get only the red centres, subset the above:

> as(red,"SpatialPoints")[red[]==255]
class       : SpatialPoints 
features    : 7 
extent      : 0.07142857, 0.9285714, 0.1, 0.9  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
> 

So that's the coordinates of the 7 red pixels in my red layer.

Repeat for green and blue to get three SpatialPoints objects and then compute distance using rgeos::gDistance:

> redpts = as(red,"SpatialPoints")[red[]==255]
> greenpts = as(green,"SpatialPoints")[green[]==255]
> library(rgeos)
> gDistance(redpts, greenpts, byid=TRUE)
           1         2         3         4         5         6         7
1  0.1428571 0.5714286 0.2457807 0.4247448 0.6167724 0.9075646 1.0724757
2  0.1428571 0.2857143 0.4729413 0.4247448 0.6167724 0.8126550 0.9075646
3  0.2857143 0.1428571 0.6054177 0.4915614 0.6645545 0.8000000 0.8494896
[ 12 rows for the 12 green pixels ]

Its very useful when you are working on problems to create small example, both so that it is easy for you to see what is going on and also to share with others. For example, I created the 5x7 rasters with this code:

Make a 5x7 matrix with random 1,2,3 values:

rgb = matrix(sample(1:3,35,TRUE),5,7)

Turn into a raster object:

rgb = raster(rgb)

Extract the 1s, 2s, and 3s into three rasters:

red = (rgb==1)*255
green = (rgb==2)*255
blue = (rgb==3)*255

The scaling by 255 is so that plotRGB(stack(red,green,blue)) produces a truly coloured plot. Please try and make sample data in your questions to save us time answering them.

Answer 2

It is unclear if you are working with an RGB multiband composite or a single raster with 3 classes that you are calling red, green, blue. Please edit your question to clarify. Here I address a single band multi-class raster. The workflow would be somewhat different if multivariate and I would ask what are you after using distance across bands? In this case, something like a non-metric MDS seems like it would be in order.

I ran into the distance function in raster and terra being somewhat untenable for reasonably sized problems. As such, I added a rasterDistance function to the spatialEco package that, is functionally the same as raster::distance but, uses the RANN package to speed up the distance calculation. Here is the example from help.

Create example data

library(raster)
library(spatialEco)
r <- raster(ncol=100,nrow=100)
  r[] <- sample(c(0,1), ncell(r), replace = TRUE)
    majority <- function(x){
     m <- table(x)
     names(m)[which.max(m)][1]
    }

r <- focal(r, matrix(1,11,11, byrow=TRUE), majority) 

Perform class distance analysis

 pts <- rasterToPoints(r, spatial=TRUE)
   cls <- pts[pts$layer == "1",] 
 d <- rasterDistance(pts, cls, reference = r, scale=TRUE)
   dev.new(height=8,width=11)
     plot(d)
       points(cls,pch=19,cex=0.5)  

My solution does not actually directly answer your question because the result is a raster representing a given class distance and not a pairwise distance matrix. @Spacedman provides the correct answer to create an actual distance matrix. However, in reality, I question how one would leverage an "all combinations" set of distance matrices representing raster data. Say, you have a 500x500 raster (n=250,000 is fairly small) with 3 equal classes (~n=83,333). Each of your 3 resulting distance matrices will be 83,333 x 250,000. How are you planning on analyzing this?