0

我有两个完全重叠的栅格（相同的范围和像元大小）。对于一个栅格中的每个像元（即对于每个 XY），我想确定在栅格之间给定阈值差异内到最近像元的欧几里得地理距离。

换句话说：raster1 和 raster2 测量一些变量 Z。我有一个 Z 值 (t) 的阈值差异，它构成了 raster1 和 raster2 之间的“匹配”值（或“足够接近”）。对于 raster1 中的每个参考像元，我需要 1) 找到 raster2 中 Z 值为 abs(Z2-Z1) 的所有像元

每个栅格有约 2600 万个像元，其中约 1000 万个具有非 NA 值。对于这个问题，我想出了一个基于非栅格的解决方法，但只能通过将栅格转换为 XYZ 表/向量并为每个参考单元格执行循环功能。对于我正在处理的数据大小来说，这在计算上过于密集（需要大约 10 天的时间来处理！）。但是，为了帮助理解我的问题，该代码如下：

library(SDMTools)
c.in <- asc2dataframe("reference.asc"); names(c.in) <- c("X","Y","Z")
f.in <- asc2dataframe("destination.asc"); names(f.in) <- c("X","Y","Z")

x=c.in$X
y=c.in$Y
c=c.in$Z
f=f.in$Z
dist=vector(length=length(c))
threshold <- 0.01

id <- 1:length(c)
for (i in length(id)) {
  # First, find all rows within the threshold
  t <- id[abs(f-c[i])<threshold]
  # Second, find the distance to the closest row
  dist[i] <- round(sqrt(min((x[t]-x[i])^2+(y[t]-y[i])^2)))
}

library(raster)
dist.rast <- rasterFromXYZ(x,y,dist)
4

1 回答 1

1

您可以将超过阈值的值设置为 NA，然后使用包中的distance()函数旁边的direction()函数计算“as-the-crow-flies”的最短距离raster。在结合栅格分辨率进行小型毕达哥拉斯计算之后，您将拥有两个栅格图层或矩阵，指定每个像元的确切位置（距离和方向）。为简单起见，您可能需要事先从栅格中移除空间投影，以移除距离和方向计算的椭球分量。它们稍后很容易添加回来。如果这一切都太慢了，我建议sparseMatrix()在 IDL 或 MATLAB 中尝试或编码。如果您使用带有 MKL 优化的 RRO，那应该会提高矩阵运算的性能。

顺便说一句，你做了很好的研究:) 替我向安德烈亚斯问好。

于 2016-01-19T20:59:34.457 回答