R:在两个栅格之间计算的阈值内找到最近的栅格单元

本文介绍了R:在两个栅格之间计算的阈值内找到最近的栅格单元的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！

问题描述

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

I have two perfectly overlapping rasters (same extents and cell size). For every cell in one raster (i.e. for every XY), I would like to determine the Euclidean geographical distance to the closest cell within a given threshold difference between the rasters.

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

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

Put another way: raster1 and raster2 measure some variable Z. I have a threshold difference for Z values (t) which constitutes a "matching" value (or "close enough") between raster1 and raster2. For each reference cell in raster1, I need to 1) find all the cells in raster2 with a Z value of abs(Z2-Z1)

Each raster has ~26 million cells, ~10 million of which have non-NA values. I have come up with a non-raster-based work-around for this problem, but only by converting the rasters to XYZ tables/vectors and performing a looping function for each reference cell. This is much too computationally intensive for the data size that I'm dealing with (takes ~10 days to process!). To assist comprehension of my question, however, that code is as follows:

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)

推荐答案

您可以将跨越阈值的值设置为 NA，然后使用 distance() 计算as-the-crow-flies"最短距离 函数与 raster 包中的 direction() 函数一起使用.然后，在包含栅格分辨率的小型勾股计算之后，您将有两个栅格层或矩阵指定每个像元的确切位置(距离和方向).为简单起见，您可能希望事先从栅格中移除空间投影，以移除距离和方向计算的椭球分量.它们很容易在以后添加回来.如果这一切都太慢，我建议尝试 sparseMatrix() 或在 IDL 或 MATLAB 中编码.如果您使用 MKL 优化附带的 RRO，那应该会提高矩阵运算的性能.

You could set the values crossing the threshold to NA and then compute the 'as-the-crow-flies' shortest distance using the distance() function alongside the direction() function within the raster package. You'd then have two raster layers or matrices specifying the exact location of each cell (distance and direction), after a small Pythagorean calculation incorporating the raster resolution. For the sake of simplicity, you may want to remove the spatial projection from the rasters beforehand to remove the ellipsoidal component of the distance and direction calculations. They are easily added back later. If all this is too slow, I recommend trying sparseMatrix() or coding it in IDL or MATLAB. If you use RRO, which ships with MKL optimization, that should improve the performance of matrix operations.

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

You do excellent research by the way :) Say hi to Andreas for me.

这篇关于R:在两个栅格之间计算的阈值内找到最近的栅格单元的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！

VALUES