R: calculate Euclidean distance between two raster layers pixels

Question

I have two raster layers for same area. I need to find the Euclidean distance between cell of coarse resolution raster and cell of fine resolution raster that fall within each cell of the pixels from my coarse resolution raster. For example:

The red square is the pixel of coarse resolution raster while the blue squares are the pixels of fine resolution raster. The black dot is the centroid of coarse resolution raster and the blue dots are the centroids of fine resolution raster.

There are similar questions posted, but the difference with my question is that I don't want to compute the nearest distances between raster cells.

My coarse resolution raster has a pixel size of 460m and my fine resolution raster of 100m. What I have done so far is to create point symbols from the centroids of the raster cells for both rasters. How can I compute the Euclidean distance between each coarse pixel and its corresponding fine pixels?

library(terra)

fr = rast("path/fine_image.tif") # fine resolution raster
cr = rast("path/coarse_image.tif") # coarse resolution raster

fr_p = as.points(fr, 
                 values = T, 
                 na.rm = T, 
                 na.all = F) # fine resolution points

cr_p = as.points(cr, 
                 values = T, 
                 na.rm = T, 
                 na.all = F) # coarse resolution points

I am not sure how to proceed from here. Any recommendations?

Here are my rasters:

fr = rast(ncols=108, nrows=203, nlyrs=1, xmin=583400, xmax=594200, ymin=1005700, ymax=1026000, names=c('B10_median'), crs='EPSG:7767')

cr = rast(ncols=23, nrows=43, nlyrs=1, xmin=583280, xmax=593860, ymin=1006020, ymax=1025800, names=c('coarse_image'), crs='EPSG:7767')

The solution came from the @michael answer and the output raster (after cropping and masking with a polygon shp) looks like this:

where the yellow squares are the cells from the coarse raster and the raster underneath it's the output from the code in the answer section.

score 3 · Accepted Answer · answered Nov 02 '22 at 05:58

3

This is a bit hacky but I think it might do what you want...

# Raster at fine resolution where values are cell indices
fr_cells <- fr
values(fr_cells) <- 1:ncell(fr)

# Second raster at fine resolution where values are indices of
# the surrounding coarse res cell (if there is one)
fr_cr <- fr
fr_xy <- xyFromCell(fr, 1:ncell(fr))
values(fr_cr) <- extract(cr, fr_xy, cells = TRUE)[, "cell"]

# Function to calculate distance given a pair of cell indices
fn <- function(x) {
  fr_xy <- xyFromCell(fr, x[1])
  cr_xy <- xyFromCell(cr, x[2])
  
  sqrt( sum( (fr_xy - cr_xy)^2 ) )
}

fr_dist <- app(c(fr_cells, fr_cr), fun = fn)

answered Nov 02 '22 at 05:58

michael

762
7
13

I don't know how to post a picture in the comments, so I posted the output on my question. It seems that your solution is the right answer so I will accept it. – Nikos Nov 02 '22 at 12:23
Very clever (even if perhaps more complicated than need be) – Robert Hijmans Nov 02 '22 at 17:22

Robert Hijmans · Answer 2 · 2022-11-03T02:50:52.770

3

You can use terra::distance for that

Example data

library(terra)    
fr <- rast(ncols=108, nrows=203, nlyrs=1, xmin=583400, xmax=594200, ymin=1005700, ymax=1026000, names='B10_median', crs='EPSG:7767')
cr <- rast(ncols=23, nrows=43, nlyrs=1, xmin=583280, xmax=593860, ymin=1006020, ymax=1025800, names='coarse_image', crs='EPSG:7767')

Solution

pts <- as.points(cr, values=FALSE, na.rm=F)
crs(pts) <- crs(cr)    
d <- distance(fr, pts)

Illustration

plot(d)
zoom(d, col=gray((1:255)/255))
lines(cr, col="red", lwd=2)

Note that this approach also computes the distance to the center of the nearest cell in cr for cells that are not covered by cr. You could remove those values with

dm <- mask(d, as.polygons(ext(cr)))
dm
#class       : SpatRaster 
#dimensions  : 203, 108, 1  (nrow, ncol, nlyr)
#resolution  : 100, 100  (x, y)
#extent      : 583400, 594200, 1005700, 1026000  (xmin, xmax, ymin, ymax)
#coord. ref. : WGS 84 / Maharashtra (EPSG:7767) 
#source(s)   : memory
#name        : B10_median 
#min value   :      0.000 
#max value   :    311.127

edited Nov 03 '22 at 02:50

answered Nov 02 '22 at 17:12

Robert Hijmans

40,301
4
55
63

I don't know how I can post an image in the comments, but using your method I get a pixel value range (0-474.130), while using the method from the other comment I am getting range values (0-311.12). I am not sure who is right though – Nikos Nov 02 '22 at 17:22
See my expanded answer. – Robert Hijmans Nov 02 '22 at 17:58
I don't know how you got that range, but I am still getting 0-474.130 even now that I used `dm <- mask(d, as.polygons(ext(cr)))`. The only difference is that I am not using the code under Illustration, but I plot the results in `QGIS`. Does this make any difference? – Nikos Nov 02 '22 at 19:57
I think there must be a problem with how you are looking at the distances Nikos. When I ran the neat solution from @RobertHijmans with your rasters it gave exactly the same distances as my long-winded method. – michael Nov 03 '22 at 02:16

R: calculate Euclidean distance between two raster layers pixels

2 Answers2