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

Question

我正在尝试使用 this 中的方法降低气候条件文章使用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 数据并确定或多或少是海岸。根据上面引用的文章中的方法，这是要遵循的前三个步骤:

为海拔 20 米以上的土地创建 DEM

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

将系数外推到海洋

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

1. 为海拔 20 米的陆地创建 DEM

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

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 个栅格(截距、海拔Est、经度Est 和纬度Est)外推到根据最后一次冰期最大值(较浅 120 米)的海岸位置

MaxGlacier <- BatPatagonia

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

投影海岸线如下图所示:



作者将系数投影到海岸的方式是通过使用 poisson_grid_fill 求解泊松方程来填补空白。 NCL 语言来自 NCAR .但我想保持简单，并尝试用同一种语言来做所有事情。我也在 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))