如何在R中基于栅格逐网格RMSE条件绘制点？

问题分析 #

你需要对两个栅格逐网格计算RMSE，当结果超过阈值a时，在对应坐标绘制标记点。现有嵌套循环仅收集RMSE值，未处理坐标记录与绘图逻辑，且循环效率较低。

高效解决方案（推荐） #

使用raster包的矢量化操作替代嵌套循环，简洁且高效，步骤如下：

1. 确保raster1和raster2为同投影、同分辨率、同行列数的多层栅格（单个网格需包含多组对比值，否则单网格RMSE无实际意义）。
2. 定义逐网格RMSE计算函数
3. 生成全栅格RMSE结果
4. 提取超标点坐标并绘图

library(raster)

# 自定义逐网格RMSE计算函数
grid_rmse <- function(x, y) {
  sqrt(mean((x - y)^2, na.rm = TRUE))
}

# 计算全栅格的RMSE结果
rmse_raster <- overlay(raster1, raster2, fun = grid_rmse)

# 设置阈值a（替换为你的实际阈值）
a <- 5

# 提取RMSE大于阈值的网格中心坐标
high_rmse_points <- rasterToPoints(rmse_raster, fun = function(x) x > a)

# 绘图：先绘制底图（示例用raster1的均值），再叠加超标点
plot(mean(raster1), main = "RMSE超过阈值的网格点")
points(high_rmse_points[, 1:2], col = "red", pch = 16, cex = 0.8)

基于原有循环的修改方案（不推荐，仅作参考） #

如果坚持使用嵌套循环，需在循环中记录超标网格的坐标，最后统一绘图：

library(raster)

nrows <- nrow(raster1)
ncols <- ncol(raster1)
rmse <- c()
points_coords <- data.frame(x = numeric(), y = numeric())
a <- 5  # 替换为你的阈值

for (i in 1:nrows) {
  for (j in 1:ncols) {
    # 获取当前网格的所有对比值（多层栅格）
    ras1_vals <- values(raster1[i, j])
    ras2_vals <- values(raster2[i, j])
    
    # 计算当前网格的RMSE
    current_rmse <- sqrt(mean((ras1_vals - ras2_vals)^2, na.rm = TRUE))
    rmse <- c(rmse, current_rmse)
    
    # 若RMSE超过阈值，记录网格中心坐标
    if (current_rmse > a) {
      cell_idx <- cellFromRowCol(raster1, i, j)
      coords <- xyFromCell(raster1, cell_idx)
      points_coords <- rbind(points_coords, coords)
    }
  }
}

# 绘图展示
plot(mean(raster1), main = "RMSE超标网格点")
points(points_coords, col = "blue", pch = 16, cex = 0.8)

注意事项 #

• 若使用单波段栅格，单个网格的两个值计算RMSE等价于绝对值（sqrt((x-y)^2) = |x-y|），建议确认需求是否为多栅格/时间序列的差异计算。
• 矢量化操作（如overlay）比嵌套循环效率高数十倍，处理大型栅格时优势更明显。

内容的提问来源于stack exchange，提问作者CovetTachi