I have a rasterStack obtained from a stack of 3 raster layers by:

    image = stack(A, B, C)

I would like to calculate the distance between the pixels of layer A and pixels of layer C. After obtaining SpatialPoints of these 2 layers I used:

    A_pts = as(red,"SpatialPoints")[red[]==1]
    C_pts = as(blue,"SpatialPoints")[blue[]==1]
    Distance_AC <- gDistance(A_pts, C_pts, byid=TRUE)

But as the number of pixels are huge the calculating process cannot finish.

Therefore I would like to calulate the distance between the layer C and only the pixels of the layer A that are in contact with layer C = the distance between pixels of layer A and the pixels of the border between A/C.

How could get these pixels between layer A and C as illustrated in the following image ?

enter image description here

  • If the size of your problem is a problem you should state what that size is in your question. As I understand it you want the distance from each red pixel to the nearest blue pixel? Write some functions to create sample data at various sizes and that will make it easier for us to work on it. Have you looked at specialised nearest-neighbour packages like FNN? – Spacedman May 5 '19 at 11:21
  • FNN can compute the nearest neighbour of a million points from another million points in 6 seconds on my laptop. – Spacedman May 5 '19 at 11:30
  • Or do you really want the NxM distances between all the red pixels and all the border pixels, for N red pixels and M border pixels? Use an edge-detection filter to get the edge pixels. – Spacedman May 5 '19 at 11:46

Do your layers A and C just have NA where they don't have a value?

i'd suggest using raster::distance. Just run separately for rasters A and C and merge (if relevant)

Here's a dummy for you, ncells = 1 million, using an extremely crude boundary line;


# make 2 rasters (your A and C), with generic values
r1 <- raster(xmn=0, xmx=1000, ymn=0, ymx=1000, ncol=1000, nrow=1000)
r2 <- raster(xmn=0, xmx=1000, ymn=0, ymx=1000, ncol=1000, nrow=1000)

r1[] <- 1
r2[] <- 2

# make a boundary polygon (crude here)

x_coord <- c(300,350,500, 500, 600, 700, 600, 800, 800, 900, 1000,1000)
y_coord <- c(0,200,400, 400, 500, 550, 450, 500, 600, 650, 700,0)
xym <- cbind(x_coord, y_coord)
p = Polygon(xym)
ps = Polygons(list(p),1)
sps = SpatialPolygons(list(ps))

# mask the two rasters with the boundary polygon just to make the data look like yours
r1.m <- mask(r1, sps, inverse = T)
r2.m <- mask(r2, sps)

# in this plot the rasters have been merged for visual aid, you don't need to

enter image description here

# use raster::distance to establish distance between NA and nearest cell of value. 
# My data was masked by the polygon but yours should already contain the NA data
# output is another raster

dists.to.r1 <- distance(r1.m)
dists.to.r2 <- distance(r2.m)

# you now have 2 rasters; distance of NA in r1 to nearest data in r1, and same for r2

# ADD them together (0 in one layer will be > 0 in other)
# this is, in effect, the distance of one raster to the other raster (i.e. nearest NA)

all.dists <- dists.r1 + dists.r2

enter image description here

the distance code take a few seconds with this dummy data

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.