如何在 R 中使用外推栅格
How to extrapolate a raster using in R
我正在尝试使用 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 中提取的同一区域 1 月份的最低温度。数据集的图表如下所示：
本题的目标
为了缩减上一次冰川盛期的过去数据，我需要模拟如果海平面与过去相同，当前气候会是什么样子。为了做到这一点，我使用 GEBCO 数据并或多或少地弄清楚海岸是什么。根据上面引用的文章中的方法，这是要遵循的前三个步骤：
- 为海拔 20 米以上的土地创建 DEM
- 计算移动中的多元线性回归window
- 将系数外推到海洋
第 3 点是我一直在努力做的事情，我将展示我是如何做到前 2 点的，并展示我一直在寻找什么来解决第 3 点
1。为海拔 20 米的陆地创建 DEM
为了做到这一点，我采用了 BatPatagonia 栅格，并使用以下代码将 20 米以下的所有值替换为 NA 值：
Elev20 <- BatPatagonia
values(Elev20) <- ifelse(values(Elev20) <= 20, NA, values(Elev20))
生成的栅格如下图所示
2。在移动 window
中计算多元线性回归
根据第 2591 页的手稿，下一步是在移动 window 中使用以下公式对超过 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")
得到的栈如下图所示：
如论文所述，移动 window 应该是 25 x 25 个单元格，总共有 625 个单元格，但是他们指出，如果移动 window 的单元格少于 170 个对于数据，不应执行回归，并且它最多应具有 624 个单元格，以确保我们仅对靠近海岸的区域进行建模。这个带有移动 window 的多元回归的结果应该是一个带有局部截距的堆栈，以及上面所示等式中每个贝塔的局部估计。我发现如何使用 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 个栅格（截距、elevationEst、longitudeEst 和 latitudeEst）外推到海岸应该根据最后一个冰川最大值浅 120 米的位置
MaxGlacier <- BatPatagonia
values(MaxGlacier) <- ifelse(values(MaxGlacier) < -120, NA,1)
预计的海岸线如下图所示：
作者将系数投影到海岸的方法是通过使用来自 NCAR. But I would like to keep it simple and try to do all in the same language. I also found a similar function in python.
的 NCL 语言的 poisson_grid_fill 求解泊松方程来填补空白
我对任何运行良好的外推过程都很满意，我不会将我的搜索限制在该算法上。
我发现了几个填充间隙的 r 包，例如 Gapfill package and even found a review of methods 来填充间隙，但它们中的大多数用于插值，并且主要用于 NDVI 图层，这些图层可以基于填充间隙的其他图层。
关于如何在这方面取得进展有什么想法吗？
谢谢
回想几十年前我的物理本科时代，我们使用拉普拉斯松弛来解决这些类型的泊松方程问题。我不确定，但我想这也可能是 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))
我正在尝试使用 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 中提取的同一区域 1 月份的最低温度。数据集的图表如下所示：
本题的目标
为了缩减上一次冰川盛期的过去数据，我需要模拟如果海平面与过去相同，当前气候会是什么样子。为了做到这一点，我使用 GEBCO 数据并或多或少地弄清楚海岸是什么。根据上面引用的文章中的方法，这是要遵循的前三个步骤：
- 为海拔 20 米以上的土地创建 DEM
- 计算移动中的多元线性回归window
- 将系数外推到海洋
第 3 点是我一直在努力做的事情，我将展示我是如何做到前 2 点的，并展示我一直在寻找什么来解决第 3 点
1。为海拔 20 米的陆地创建 DEM
为了做到这一点，我采用了 BatPatagonia 栅格，并使用以下代码将 20 米以下的所有值替换为 NA 值：
Elev20 <- BatPatagonia
values(Elev20) <- ifelse(values(Elev20) <= 20, NA, values(Elev20))
生成的栅格如下图所示
2。在移动 window
中计算多元线性回归根据第 2591 页的手稿，下一步是在移动 window 中使用以下公式对超过 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")
得到的栈如下图所示：
如论文所述，移动 window 应该是 25 x 25 个单元格，总共有 625 个单元格，但是他们指出，如果移动 window 的单元格少于 170 个对于数据，不应执行回归，并且它最多应具有 624 个单元格，以确保我们仅对靠近海岸的区域进行建模。这个带有移动 window 的多元回归的结果应该是一个带有局部截距的堆栈，以及上面所示等式中每个贝塔的局部估计。我发现如何使用 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 个栅格（截距、elevationEst、longitudeEst 和 latitudeEst）外推到海岸应该根据最后一个冰川最大值浅 120 米的位置
MaxGlacier <- BatPatagonia
values(MaxGlacier) <- ifelse(values(MaxGlacier) < -120, NA,1)
预计的海岸线如下图所示：
作者将系数投影到海岸的方法是通过使用来自 NCAR. But I would like to keep it simple and try to do all in the same language. I also found a similar function in python.
的 NCL 语言的poisson_grid_fill 求解泊松方程来填补空白
我对任何运行良好的外推过程都很满意，我不会将我的搜索限制在该算法上。
我发现了几个填充间隙的 r 包，例如 Gapfill package and even found a review of methods 来填充间隙，但它们中的大多数用于插值，并且主要用于 NDVI 图层，这些图层可以基于填充间隙的其他图层。
关于如何在这方面取得进展有什么想法吗？
谢谢
回想几十年前我的物理本科时代，我们使用拉普拉斯松弛来解决这些类型的泊松方程问题。我不确定，但我想这也可能是 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))