我正在尝试使用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中提取的同一区域1月份的最低温度。下面是这些数据集的图表：
这个问题的目标
为了缩小过去冰川时期的数据，我需要模拟当前气候如果海平面与过去相同会是什么样子。为了做到这一点，我使用GEBCO数据，并找出海岸线的大致位置。根据上述文章中的方法，以下是要遵循的前三个步骤：
- 创建一个高于海平面20米的陆地DEM
- 在移动窗口中计算多元线性回归
- 将系数外推到海洋
第3点是我一直在努力解决的问题，我将展示我如何完成前两个步骤，并展示我一直在寻找的解决第3点的方法。
1. 创建一个高于海平面20米的陆地DEM
为了做到这一点，我使用了栅格图像，并使用以下代码将所有低于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")
下图显示了生成的堆栈：
如论文所述，移动窗口应为25x25个单元格，共计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 栈：
我卡在这里了
将系数外推到海洋
现在的目标是将系数堆栈的前4个栅格（截距，高程估计，经度估计和纬度估计）外推到上次冰河期时海岸线应该在的位置，当时海平面下降了120米。
MaxGlacier <- BatPatagonia
values(MaxGlacier) <- ifelse(values(MaxGlacier) < -120, NA,1)
以下地图显示了预计海岸线：
poisson_grid_fill函数通过求解Poisson方程来填补间隙。但我想保持简单，尝试在同一种语言中完成所有操作。我还在Python中找到了类似的函数。我会满意任何有效的外推过程，我不限制搜索算法。
我发现了几个R包，如Gapfill包，并找到了一个review of methods来填补间隙，但它们大多用于插值，主要用于基于其他图层的NDVI图层，其中填充了间隙。
我该如何继续？
Elev20 = mask(BatPatagonia, BatPatagonia <= 20, maskvalue = 1)、Longitud = mask(init(BatPatagonia, 'x'), Elev20)和Latitud = mask(init(BatPatagonia, 'y'), Elev20)。;-) - ztl