如何在 R 中使用外推栅格

Question

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

需要的包和数据

对于本示例，我将一些数据上传到 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)

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

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

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

我发现了几个 r 包来填补空白，例如 Gapfill 包，甚至找到了一个 review of methods 来填补空白，但它们中的大多数用于插值，并且主要用于 NDVI 层，这些层可以基于其他层，其中填补了空白。

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

谢谢

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))