如何在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是一个栅格文件,其中Bathymetry和从GEBCO数据集中提取和转换的区域的高度,而TempMinPatNow是从worldclim中提取的1月相同区域的最低温度.下面是数据集的图:

这个问题的目标

为了从过去的冰川最大值下调过去的数据,我需要模拟当前的气候如果海平面与过去相同的情况.为了做到这一点,我使用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乘25个单元,总共产生625个单元,但是它们表明如果移动窗口的数据少于170个单元,则不应执行回归,并且应该最多有624个单元格,以确保我们只对靠近海岸的区域进行建模.具有移动窗口的这种多元回归的结果应该是具有局部截距的叠加,以及上面示出的等式中的每个Betas的局部估计.我发现如何使用以下代码使用该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层,这些层可以基于填充间隙的其他层.

关于如何在这方面采取行动的任何想法？

谢谢

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