匹配R中数值的像素的栅格位置/坐标

Question

有没有一种方法可以提取与特定值匹配的像素的网格位置或（最好为具有显式范围的栅格）点/质心坐标？我几乎有一个非常低效的工作流程，将其转换为矩阵并使用which(mtrx == max(mtrx), arr.ind = TRUE)来获取矩阵位置，但这会（a）丢失地理空间信息和（b）在矩阵转换过程中导致数据旋转90度，两者都需要额外的代码使其工作并显著减慢计算速度。是否有等效的栅格工作流程？

- geotheory

---

你能澄清一下你的数据最初是以哪种形式给出的吗？ - Georg Schnabel

也许我漏掉了什么，但是问题不应该普遍适用于所有栅格对象吗，无论其起源如何？假设rastr = raster(volcano); rastr@extent = extent(-76,-74,42,44)，我们想要定位rastr@data@max的值。 - geotheory

Answer 1

示例数据：

library(raster)
set.seed(0)
r <- raster(ncols=10, nrows=10)
r[] <- sample(50, 100, replace=T)

现在执行以下步骤：

p <- rasterToPoints(r, function(x) x == 11)

为了获取：

       x   y layer
[1,]   18  81    11
[2,] -126  63    11
[3,]  -90  45    11
[4,]   54 -63    11

如果您想要具有最大值的单元格

vmax = maxValue(r)
p <- rasterToPoints(r, function(x) all.equal(x, vmax)

（请注意，不要使用@data@max）

- Robert Hijmans

Answer 2

我不明白为什么你要强制将其转换为矩阵？也许我没有理解你的问题，但如果我正确理解了你的意思，你可以直接查询栅格值，然后强制转换为点来获取地理位置。

require(raster)
r <- raster(ncols=100, nrows=100)
  r[] <- runif(ncell(r), 0,1)

# Coerce < max to NA and coerce result to points 
rMax <- r
m = maxValue(r) 
rMax[rMax != m] <- NA 
( r.pts <- rasterToPoints (rMax) )  

# You could also use the raster specific Which or which.max functions. 
i <- which.max(r)
  xy.max <- xyFromCell(r, i)
    plot(r)
      points(xy.max, pch=19, col="black")

# Or for a more general application of Which
i <- Which(r >= 0.85, cells=TRUE)
  xy.max <- xyFromCell(r, i)
    plot(r)
      points(xy.max, pch=19, col="black")

# If you prefer a raster object set cells=FALSE
i <- Which(r >= 0.85, cells=FALSE)
  plot(i)

有多个栅格函数可以让您向其中传递自定义或基本函数。您可能想查看“focal”（局部算子）或“calc”。您可能还需要阅读与栅格相关的帮助文件。

- Jeffrey Evans

---

谢谢，这比我的当前工作流程要好：pts = data.frame(rasterToPoints(rastr)); pts[pts[,3] == rastr@data@max,]。我仍然想知道是否有一个栅格等效的 which 函数来返回像素组。如果没有，我会接受这个答案。 - geotheory

你想要从哪个最大值获取arr.ind？你所说的像素组是什么意思？ - mdsumner

Answer 3

扩展Jeffrey的答案，您可以使用以下代码选择最低栅格值的最后一个实例：

r <- raster(ncols=12, nrows=12)
set.seed(0)
r[] <- round(runif(ncell(r))*0.7 )
rc <- clump(r)
rc[12,8]<-1
plot(rc)

xy.min<-data.frame(xyFromCell(rc,max(which.min(rc))))
xy.min$dat<-1
coordinates(xy.min)<-~x+y
points(xy.min,lwd=2)

- jsta