计算两个栅格图层像素之间的欧几里得距离

Question

我有两个栅格图层，表示同一区域。我需要找到粗分辨率栅格的每个像素单元和细分辨率栅格的每个像素单元之间的欧几里得距离，这些单元都在粗分辨率栅格的像素单元内。例如：

红色正方形是粗分辨率栅格的像素，而蓝色正方形是细分辨率栅格的像素。黑点是粗分辨率栅格的质心，蓝点是细分辨率栅格的质心。

已经有类似的问题发布了，但我的问题与他们的不同之处在于我不想计算栅格单元之间的最近距离。

我的粗分辨率栅格具有460米的像素大小，而我的细分辨率栅格具有100米的像素大小。到目前为止，我已经从两个栅格单元的质心创建了点符号。如何计算每个粗像素与其相应的细像素之间的欧几里得距离？

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

我不确定该如何继续。有什么建议吗？

这是我的光栅图像：

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')

解决方案来自@michael的答案，输出的栅格图像（经过裁剪并使用多边形shp进行遮罩）如下所示：

“黄色正方形是粗栅格中的单元格，下面的栅格是答案部分代码的输出结果。”

- Nikos

Answer 1

这个有点儿技巧性，但我认为它可能做到你想要的...

# 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)

- michael

---

我不知道如何在评论中发布图片，所以我在我的问题上发布了输出。看起来你的解决方案是正确答案，所以我会接受它。 - Nikos

非常聪明（即使可能比必要的更复杂） - Robert Hijmans

Answer 2

你可以使用terra::distance进行计算。

示例数据。

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')

解决方案

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

插图

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

请注意，这种方法还会计算未被cr覆盖的单元格到最近单元格中心的距离。您可以使用以下方法删除这些值。

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 

- Robert Hijmans

---

我不知道如何在评论中发布图片，但是使用您的方法，我获得了像素值范围（0-474.130），而使用另一条评论的方法则获得了范围值（0-311.12）。我不确定哪个方法是正确的。 - Nikos

请查看我扩展的答案。 - Robert Hijmans

我不知道你是怎么得到那个范围的，但是即使我使用了 dm <- mask(d, as.polygons(ext(cr)))，我仍然得到的是0-474.130。唯一的区别是我没有使用Illustration下面的代码，而是在QGIS中绘制了结果。这会有什么影响吗？ - Nikos

我认为你在计算距离时可能存在问题，Nikos。当我使用@RobertHijmans的简洁解决方案并输入你的栅格数据时，得到的距离与我冗长的方法完全相同。 - michael