r - 如何在 R 中使用外推栅格

Question

我正在尝试使用本文中使用 R 软件的方法来降低气候条件。我快到了，但我错过了几个步骤

需要的包和数据

对于此示例，我将一些数据上传到 archive.org 网站以加载此示例中使用的所需包和数据，使用以下代码：

library(raster)
library(rgdal)

download.file("https://archive.org/download/Downscaling/BatPatagonia.rds", "Bat.rds")
download.file("https://archive.org/download/Downscaling/TempMinPatNow.rds", "Tmin.rds")

BatPatagonia <- readRDS("Bat.rds")
TempMinPatNow <- readRDS("Tmin.rds")

BatPatagonia 是一个栅格文件，其中包含从 GEBCO 数据集提取和转换的区域的水深和海拔高度，而 TempMinPatNow 是从 worldclim 中提取的同一区域一月份的最低温度。数据集的图如下所示：

这个问题的目标

为了从最后一个冰川最大值缩小过去的数据，我需要模拟如果海平面与过去相同，当前的气候会如何。为了做到这一点，我使用 GEBCO 数据并或多或少地弄清楚海岸是什么。根据上面引用的文章中的方法，这是要遵循的前三个步骤：

1. 为海拔 20 米以上的土地创建 DEM
2. 在移动窗口中计算多元线性回归
3. 将系数外推到海洋

第 3 点是我一直在努力做的事情，我将展示我是如何做到前 2 点的，并展示我一直在寻找解决第 3 点的问题

1.为海拔20米的土地创建一个DEM

为了做到这一点，我使用了BatPatagonia栅格，并使用以下代码将 20 米以下的所有值替换为 NA 值：

Elev20 <- BatPatagonia

values(Elev20) <- ifelse(values(Elev20) <= 20, NA, values(Elev20))

生成的栅格如下图所示

2. 在移动窗口中计算多元线性回归

根据第 2591 页的手稿，下一步是在移动窗口中对超过 20 米的高度使用以下公式进行多重线性回归：

我们已经有了高程数据，但我们还需要纬度和经度的栅格，为此我们使用以下代码，首先创建纬度和经度栅格：

Latitud <- BatPatagonia
Longitud <- BatPatagonia

data_matrix <- raster::xyFromCell(BatPatagonia, 1:ncell(BatPatagonia))

values(Latitud) <- data_matrix[, 2]
values(Longitud) <- data_matrix[, 1]

我们将把它乘以海拔超过 20 米的区域的光栅蒙版，这样我们就只得到我们需要的值：

Elev20Mask <- BatPatagonia

values(Elev20Mask) <- ifelse(values(Elev20Mask) <= 20, NA, 1)

Longitud <- Elev20Mask*Longitud

Latitud <- Elev20Mask*Latitud

现在我将使用响应变量和预测变量构建一个堆栈：

Preds <- stack(Elev20, Longitud, Latitud, TempMinPatNow)

names(Preds) <- c("Elev", "Longitud", "Latitud", "Tmin")

生成的堆栈如下图所示：

如论文中所述，移动窗口应为 25 x 25 个单元格，总共有 625 个单元格，但是他们指出，如果移动窗口中包含数据的单元格少于 170 个，则不应执行回归，并且应该最多有 624 个单元，以确保我们只对靠近海岸的区域进行建模。这个带有移动窗口的多重回归的结果应该是一个带有局部截距的堆栈，以及上面显示的等式中每个 Beta 的局部估计。我发现了如何使用以下代码使用该getValuesFocal函数进行此操作（此循环需要一段时间）：

# First we establish the 25 by 25 window

w <- c(25, 25)

# Then we create the empty layers for the resulting stack

intercept <- Preds[[1]]
intercept[] <- NA

elevationEst <- intercept

latitudeEst <- intercept

longitudeEst <- intercept

现在我们开始代码：

for (rl in 1:nrow(Preds)) {
  v <- getValuesFocal(Preds[[1:4]], row = rl, nrows = 1, ngb = w, array = FALSE)
  int <- rep(NA, nrow(v[[1]]))
  x1 <- rep(NA, nrow(v[[1]]))
  x2 <- rep(NA, nrow(v[[1]]))
  x3 <- rep(NA, nrow(v[[1]]))
  x4 <- rep(NA, nrow(v[[1]]))
  for (i in 1:nrow(v[[1]])) {
    xy <- na.omit(data.frame(x1 = v[[1]][i, ], x2 = v[[2]][i, ], x3 = v[[3]][i, 
                                                                         ], y = v[[4]][i, ]))

    if (nrow(xy) > 170 & nrow(xy) <= 624) {
      coefs <- coefficients(lm(as.numeric(xy$y) ~ as.numeric(xy$x1) + 
                             as.numeric(xy$x2) + as.numeric(xy$x3)))
      int[i] <- coefs[1]
      x1[i] <- coefs[2]
      x2[i] <- coefs[3]
      x3[i] <- coefs[4]
    } else {
      int[i] <- NA
      x1[i] <- NA
      x2[i] <- NA
      x3[i] <- NA
    }
  }

  intercept[rl, ] <- int
  elevationEst[rl, ] <- x1
  longitudeEst[rl, ] <- x2
  latitudeEst[rl, ] <- x3

  message(paste(rl, "of", nrow(Preds), "ready"))
}

Coeffs <- stack(intercept, elevationEst, latitudeEst, longitudeEst, (intercept + Preds$Elev * elevationEst + Preds$Longitud * longitudeEst + Preds$Latitud *latitudeEst), Preds$Tmin)

names(Coeffs) <- c("intercept", "elevationEst", "longitudeEst", "latitudeEst", "fitted", "Observed")

此循环的结果是coeffs堆栈，如下所示：

这就是我卡住的地方：

将系数外推到海洋

现在的目标是将 Coeffs 堆栈的前 4 个栅格（intercept、elevationEst、longitudeEst 和 latitudeEst）外推到海岸应位于的位置，根据最后一次冰川最大值（浅 120 米）

MaxGlacier <- BatPatagonia

values(MaxGlacier) <- ifelse(values(MaxGlacier) < -120, NA,1)

预计的海岸线如下图所示：

作者将系数投影到海岸的方式是通过使用NCARpoisson_grid_fill的 NCL 语言求解泊松方程来填补空白。但我想保持简单，并尝试用同一种语言来做所有事情。我还在python中找到了类似的功能。

我会对任何运行良好的外推过程感到满意，我不会将我的搜索限制在该算法上。

我发现了几个填补空白的 r 包，例如Gapfill包，甚至发现了填补空白的方法的评论，但其中大部分是用于插值的，并且主要用于 NDVI 层，这些层可以基于填补空白的其他层。

关于如何继续前进的任何想法？

谢谢

标签： rr-rasterpde

Answer 1

回想几十年我的物理本科时代，我们使用拉普拉斯松弛来解决这些类型的泊松方程问题。我不确定，但我想这也可能是如何poisson_grid_fill工作的。过程很简单。松弛是一个迭代过程，我们计算每个单元格，除了那些形成边界条件的单元格作为水平或垂直相邻单元格的平均值，然后重复直到结果接近稳定解。

在您的情况下，您已经拥有值的单元格提供了您的边界条件，我们可以迭代其他单元格。像这样的东西（这里展示了截距系数——你可以用同样的方式做其他的）：

gaps = which(is.na(intercept)[])
intercept.ext = intercept
w=matrix(c(0,0.25,0,0.25,0,0.25,0,0.25,0), nc=3, nr=3)
max.it = 1000
for (i in 1:max.it) intercept.ext[gaps] = focal(intercept.ext, w=w, na.rm=TRUE)[gaps]
intercept.ext = mask(intercept.ext, MaxGlacier)

编辑

这是嵌入在函数中的相同过程，以演示如何使用while循环，直到达到所需的容差（或超过最大迭代次数）。请注意，此功能是为了演示原理，并没有针对速度进行优化。

gap.fill = function(r, max.it = 1e4, tol = 1e-2, verbose=FALSE) {
  gaps = which(is.na(r)[])
  r.filled = r
  w = matrix(c(0,0.25,0,0.25,0,0.25,0,0.25,0), nc=3, nr=3)
  i = 0
  while(i < max.it) {
    i = i + 1
    new.vals = focal(r.filled, w=w, na.rm=TRUE)[gaps]
    max.residual = suppressWarnings(max(abs(r.filled[gaps] - new.vals), na.rm = TRUE))
    if (verbose) print(paste('Iteration', i, ': residual = ', max.residual))
    r.filled[gaps] = new.vals
    if (is.finite(max.residual) & max.residual <= tol) break
  }
  return(r.filled)
}

intercept.ext = gap.fill(intercept)
intercept.ext = mask(intercept.ext, MaxGlacier)
plot(stack(intercept, intercept.ext))