使用 gstat 和 automap 包的克里金法 - 复制教程时出现的问题
Kriging using the gstat and automap packages - problems when copying a tutorial
对于我研究地点的几个位置，我有水深 (m) 数据，我正在尝试使用克里金法将深度插入到我没有数据的位置。我发现了一个有用的博客 post 由 Dr. Wilke here 撰写，并尝试将他的代码应用到我的数据集。但是，当我 运行 所有代码并绘制结果时，所有绘图都返回为空（灰色框；深度样本返回时带有灰点，而不是在连续尺度上反映深度的颜色）。谁能告诉我我做错了什么？下面是我试过的代码。我怀疑问题出在第 2 步，因为这里的代码输出开始偏离博客代码输出。
设置
depth <- structure(list(X = c(668875.407301466, 675555.822548817, 677165.984257575,
671427.139500822, 670565.010067336, 676196.888461476, 669435.428473297,
677288.114994538, 679355.228144992, 677003.998083734, 677585.341252976,
679549.867927876, 677659.023812075, 679730.677304853, 677749.442051316,
669617.624346464, 671569.53262659, 673306.072271015, 670280.828937677,
677894.966356866, 677192.79791361, 672207.85937462, 674600.067203517,
674244.297764696, 673763.964359285, 669723.003482714, 677106.449499354,
671172.528780568, 673040.137152061, 672415.374390262, 675286.888306849,
674272.583595863, 672978.2232579, 672207.731105759, 673057.984061267,
673098.008775404, 670469.846198017, 670461.213900251, 671449.443237703,
671476.810705895, 672429.201320869, 672426.444205806, 673397.966944431,
673397.001177047, 676506.198002199, 679013.454950302, 668360.513631175,
669309.125013407, 675540.056558172, 673946.851010242, 668736.385947233,
669581.113779484, 669531.260804133, 668572.915079444, 668562.271762613
), Y = c(2843306.80511914, 2846996.53521373, 2847483.31959936,
2843061.483504, 2842817.70752641, 2853858.66887037, 2843862.35002759,
2844357.28754673, 2844326.38872235, 2843760.75211182, 2845096.92928085,
2844999.32573033, 2845739.40590226, 2845654.35288888, 2846405.37514488,
2846251.04899388, 2856292.70424562, 2855857.84577301, 2848289.20711897,
2848868.04137537, 2851322.64880942, 2853330.72832414, 2854696.08607969,
2850951.06030575, 2849409.17353187, 2840685.32365573, 2837738.38794808,
2850715.9222808, 2846526.89253892, 2842945.78513155, 2844114.61778685,
2844613.08986426, 2844813.27348601, 2845405.92798845, 2846082.85831243,
2848010.02032868, 2841631.05374873, 2840583.99941216, 2841647.00399675,
2840617.02709824, 2841650.84197688, 2840627.11766348, 2841420.84288763,
2840419.29474805, 2843303.19107449, 2843694.67766026, 2838715.91143071,
2838736.83548694, 2839850.31020878, 2842589.10652076, 2840663.87712246,
2842046.19194071, 2839662.52414273, 2842072.11298559, 2839670.12536644
), Z = c(27, 2, 1, 1.5, 1.5, 4, 6, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5,
1.5, 1.5, 14, 26, 4, 12, 3, 3, 7, 4, 1, 1, 14.4, 2, 6, 1, 1,
1, 1, 1, 1, 1, 1, 5, 4, 3, 3, 2, 4, 2, 3, 2, 2, 27.8, 10, 2,
1.5, 28, 8, 16.6, 28, 28.2)), row.names = c(NA, -55L), class = c("tbl_df",
"tbl", "data.frame"))
X 和 Y 是包含空间坐标的列，Z 包含这些位置的深度值。
# Convert to {sf} because that is the best way to store spatial points
depth_sf <- st_as_sf(depth, coords = c("X", "Y"), crs = 32617) %>%
cbind(st_coordinates(.))
创建变异函数
# We will discuss later, what Z~1 does actually mean in this context
v_emp_OK <- gstat::variogram(
Z~1,
as(depth_sf, "Spatial") # switch from {sf} to {sp}
)
# automap's autofitVariogram actually produces more info than we need. I will only keep the var_model part.
v_mod_OK <- automap::autofitVariogram(Z~1, as(depth_sf, "Spatial"))$var_model
# To inspect the automatic fit that was chosen for us we can use automap's excellent build in methods for base::plot
plot(automap::autofitVariogram(Z~1, as(depth_sf, "Spatial")))
定义目标网格
# technically we already have a grid from the initial dataset, but as we are still working under the pretense that our only available data are the simulated observation wells, we will construct our grid from that object.'
# Step 1: define a grid based on the bounding box of our observations
grd_100_sf <- depth_sf %>%
st_bbox() %>%
st_as_sfc() %>%
st_make_grid(
cellsize = c(100, 100), # 100m pixel size
what = "centers"
) %>%
st_as_sf() %>%
cbind(., st_coordinates(.))
#Step 2: making our grid work for gstat
grd_100_sp <- as(grd_100_sf, "Spatial") # converting to {sp} format
gridded(grd_100_sp) <- TRUE # informing the object that it is a grid
grd_100_sp <- as(grd_100_sp, "SpatialPixels") # specifying what kind of grid
克里金法
# Ordinary Kriging
OK <- krige(
Z~1, # Z is our variable and "~1" means "depends on mean"
as(depth_sf, "Spatial"), # input data in {sp} format
grd_100_sp, # locations to interpolate at
model = v_mod_OK # the variogram model fitted above
)
## [using ordinary kriging]
# Simple Kriging
SK <- krige(
Z~1, # Z still depends on mean
beta = mean(depth$Z), # but this time we know the mean's value
as(depth_sf, "Spatial"), # input data in {sp} format
grd_100_sp, # locations to interpolate at
model = v_mod_OK # the variogram model fitted above
)
## [using simple kriging]
# Universal Kriging
# Implementing this method is somewhat different.
# we no longer assume that Z is essentially depending on a single mean but
# rather on the position of the interpolation location within our target grid
UK <- krige(
Z~coords.x1+coords.x2, # Think "Z~X+Y" but {sp} conversion alters variable naming
as(depth_sf, "Spatial"), # input data in {sp} format (`X` --> `coords.x1`)
grd_100_sp, # locations to interpolate at
model = autofitVariogram( # we need an appropriate variogram fit
Z~X+Y, # here we can keep "X+Y" - it's just how it is
as(depth_sf, "Spatial")
)$var_model
)
## [using universal kriging]
# I'll also add an inverse distance weighted model to provide a baseline
# for model evaluation
# Note how the only difference to Ordinary Kriging is the absence of a
# fitted variogram model
idwres <- idw(
Z~1, # idw also depends on mean
as(depth_sf, "Spatial"), # input data in {sp} format
grd_100_sp, # locations to interpolate at
)
检查结果
# A function to plot rasters
plot_my_gstat_output <- function(raster_object, object_name){
df <- rasterToPoints(raster_object) %>% as_tibble()
colnames(df) <- c("X", "Y", "Z")
ggplot(df, aes(x = X, y = Y, fill = Z)) +
geom_raster() +
ggtitle(label = object_name) +
scale_fill_viridis(option = "B", limits = c(50, 100)) +
theme_void() +
theme(
plot.title = element_text(hjust = 0.5)
)
}
p_idw <- plot_my_gstat_output(raster(idwres), "IDW")
p_SK <- plot_my_gstat_output(raster(SK), "Simple Kriging")
p_OK <- plot_my_gstat_output(raster(OK), "Ordinary Kriging")
p_UK <- plot_my_gstat_output(raster(UK), "Universal Kriging")
# I also want to display sampling locations
p_depth <- ggplot(
data = depth,
mapping = aes(x = X, y = Y, color = Z)
) +
geom_point(size = 3) +
scale_color_viridis(option = "B", limits = c(50, 100)) +
ggtitle(label = "Depths Sampled") +
theme_void() +
theme(
plot.title = element_text(hjust = 0.5)
)
# This works because of library(patchwork)
(p_depth + p_idw) /
(p_SK + p_OK + p_UK) +
plot_layout(guides = 'collect')
看起来罪魁祸首在代码的 ggplot 部分。使用 scale_color_viridis() 时，您需要更改 limits 以适合您的数据。由于数据的深度为 0 到 28，因此您应该将限制从 c(50, 100) 更改为 c(0, 30).
# A function to plot rasters
plot_my_gstat_output <- function(raster_object, object_name){
df <- rasterToPoints(raster_object) %>% as_tibble()
colnames(df) <- c("X", "Y", "Z")
ggplot(df, aes(x = X, y = Y, fill = Z)) +
geom_raster() +
ggtitle(label = object_name) +
scale_fill_viridis(option = "B", limits = c(0, 30)) +
theme_void() +
theme(
plot.title = element_text(hjust = 0.5)
)
}
p_idw <- plot_my_gstat_output(raster(idwres), "IDW")
p_SK <- plot_my_gstat_output(raster(SK), "Simple Kriging")
p_OK <- plot_my_gstat_output(raster(OK), "Ordinary Kriging")
p_UK <- plot_my_gstat_output(raster(UK), "Universal Kriging")
# I also want to display sampling locations
p_depth <- ggplot(
data = depth,
mapping = aes(x = X, y = Y, color = Z)
) +
geom_point(size = 3) +
scale_color_viridis(option = "B", limits = c(0, 30)) +
ggtitle(label = "Depths Sampled") +
theme_void() +
theme(
plot.title = element_text(hjust = 0.5)
)
# This works because of library(patchwork)
(p_depth + p_idw) /
(p_SK + p_OK + p_UK) +
plot_layout(guides = 'collect')
输出
对于我研究地点的几个位置，我有水深 (m) 数据，我正在尝试使用克里金法将深度插入到我没有数据的位置。我发现了一个有用的博客 post 由 Dr. Wilke here 撰写，并尝试将他的代码应用到我的数据集。但是，当我 运行 所有代码并绘制结果时，所有绘图都返回为空（灰色框；深度样本返回时带有灰点，而不是在连续尺度上反映深度的颜色）。谁能告诉我我做错了什么？下面是我试过的代码。我怀疑问题出在第 2 步，因为这里的代码输出开始偏离博客代码输出。
设置
depth <- structure(list(X = c(668875.407301466, 675555.822548817, 677165.984257575,
671427.139500822, 670565.010067336, 676196.888461476, 669435.428473297,
677288.114994538, 679355.228144992, 677003.998083734, 677585.341252976,
679549.867927876, 677659.023812075, 679730.677304853, 677749.442051316,
669617.624346464, 671569.53262659, 673306.072271015, 670280.828937677,
677894.966356866, 677192.79791361, 672207.85937462, 674600.067203517,
674244.297764696, 673763.964359285, 669723.003482714, 677106.449499354,
671172.528780568, 673040.137152061, 672415.374390262, 675286.888306849,
674272.583595863, 672978.2232579, 672207.731105759, 673057.984061267,
673098.008775404, 670469.846198017, 670461.213900251, 671449.443237703,
671476.810705895, 672429.201320869, 672426.444205806, 673397.966944431,
673397.001177047, 676506.198002199, 679013.454950302, 668360.513631175,
669309.125013407, 675540.056558172, 673946.851010242, 668736.385947233,
669581.113779484, 669531.260804133, 668572.915079444, 668562.271762613
), Y = c(2843306.80511914, 2846996.53521373, 2847483.31959936,
2843061.483504, 2842817.70752641, 2853858.66887037, 2843862.35002759,
2844357.28754673, 2844326.38872235, 2843760.75211182, 2845096.92928085,
2844999.32573033, 2845739.40590226, 2845654.35288888, 2846405.37514488,
2846251.04899388, 2856292.70424562, 2855857.84577301, 2848289.20711897,
2848868.04137537, 2851322.64880942, 2853330.72832414, 2854696.08607969,
2850951.06030575, 2849409.17353187, 2840685.32365573, 2837738.38794808,
2850715.9222808, 2846526.89253892, 2842945.78513155, 2844114.61778685,
2844613.08986426, 2844813.27348601, 2845405.92798845, 2846082.85831243,
2848010.02032868, 2841631.05374873, 2840583.99941216, 2841647.00399675,
2840617.02709824, 2841650.84197688, 2840627.11766348, 2841420.84288763,
2840419.29474805, 2843303.19107449, 2843694.67766026, 2838715.91143071,
2838736.83548694, 2839850.31020878, 2842589.10652076, 2840663.87712246,
2842046.19194071, 2839662.52414273, 2842072.11298559, 2839670.12536644
), Z = c(27, 2, 1, 1.5, 1.5, 4, 6, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5,
1.5, 1.5, 14, 26, 4, 12, 3, 3, 7, 4, 1, 1, 14.4, 2, 6, 1, 1,
1, 1, 1, 1, 1, 1, 5, 4, 3, 3, 2, 4, 2, 3, 2, 2, 27.8, 10, 2,
1.5, 28, 8, 16.6, 28, 28.2)), row.names = c(NA, -55L), class = c("tbl_df",
"tbl", "data.frame"))
X 和 Y 是包含空间坐标的列，Z 包含这些位置的深度值。
# Convert to {sf} because that is the best way to store spatial points
depth_sf <- st_as_sf(depth, coords = c("X", "Y"), crs = 32617) %>%
cbind(st_coordinates(.))
创建变异函数
# We will discuss later, what Z~1 does actually mean in this context
v_emp_OK <- gstat::variogram(
Z~1,
as(depth_sf, "Spatial") # switch from {sf} to {sp}
)
# automap's autofitVariogram actually produces more info than we need. I will only keep the var_model part.
v_mod_OK <- automap::autofitVariogram(Z~1, as(depth_sf, "Spatial"))$var_model
# To inspect the automatic fit that was chosen for us we can use automap's excellent build in methods for base::plot
plot(automap::autofitVariogram(Z~1, as(depth_sf, "Spatial")))
定义目标网格
# technically we already have a grid from the initial dataset, but as we are still working under the pretense that our only available data are the simulated observation wells, we will construct our grid from that object.'
# Step 1: define a grid based on the bounding box of our observations
grd_100_sf <- depth_sf %>%
st_bbox() %>%
st_as_sfc() %>%
st_make_grid(
cellsize = c(100, 100), # 100m pixel size
what = "centers"
) %>%
st_as_sf() %>%
cbind(., st_coordinates(.))
#Step 2: making our grid work for gstat
grd_100_sp <- as(grd_100_sf, "Spatial") # converting to {sp} format
gridded(grd_100_sp) <- TRUE # informing the object that it is a grid
grd_100_sp <- as(grd_100_sp, "SpatialPixels") # specifying what kind of grid
克里金法
# Ordinary Kriging
OK <- krige(
Z~1, # Z is our variable and "~1" means "depends on mean"
as(depth_sf, "Spatial"), # input data in {sp} format
grd_100_sp, # locations to interpolate at
model = v_mod_OK # the variogram model fitted above
)
## [using ordinary kriging]
# Simple Kriging
SK <- krige(
Z~1, # Z still depends on mean
beta = mean(depth$Z), # but this time we know the mean's value
as(depth_sf, "Spatial"), # input data in {sp} format
grd_100_sp, # locations to interpolate at
model = v_mod_OK # the variogram model fitted above
)
## [using simple kriging]
# Universal Kriging
# Implementing this method is somewhat different.
# we no longer assume that Z is essentially depending on a single mean but
# rather on the position of the interpolation location within our target grid
UK <- krige(
Z~coords.x1+coords.x2, # Think "Z~X+Y" but {sp} conversion alters variable naming
as(depth_sf, "Spatial"), # input data in {sp} format (`X` --> `coords.x1`)
grd_100_sp, # locations to interpolate at
model = autofitVariogram( # we need an appropriate variogram fit
Z~X+Y, # here we can keep "X+Y" - it's just how it is
as(depth_sf, "Spatial")
)$var_model
)
## [using universal kriging]
# I'll also add an inverse distance weighted model to provide a baseline
# for model evaluation
# Note how the only difference to Ordinary Kriging is the absence of a
# fitted variogram model
idwres <- idw(
Z~1, # idw also depends on mean
as(depth_sf, "Spatial"), # input data in {sp} format
grd_100_sp, # locations to interpolate at
)
检查结果
# A function to plot rasters
plot_my_gstat_output <- function(raster_object, object_name){
df <- rasterToPoints(raster_object) %>% as_tibble()
colnames(df) <- c("X", "Y", "Z")
ggplot(df, aes(x = X, y = Y, fill = Z)) +
geom_raster() +
ggtitle(label = object_name) +
scale_fill_viridis(option = "B", limits = c(50, 100)) +
theme_void() +
theme(
plot.title = element_text(hjust = 0.5)
)
}
p_idw <- plot_my_gstat_output(raster(idwres), "IDW")
p_SK <- plot_my_gstat_output(raster(SK), "Simple Kriging")
p_OK <- plot_my_gstat_output(raster(OK), "Ordinary Kriging")
p_UK <- plot_my_gstat_output(raster(UK), "Universal Kriging")
# I also want to display sampling locations
p_depth <- ggplot(
data = depth,
mapping = aes(x = X, y = Y, color = Z)
) +
geom_point(size = 3) +
scale_color_viridis(option = "B", limits = c(50, 100)) +
ggtitle(label = "Depths Sampled") +
theme_void() +
theme(
plot.title = element_text(hjust = 0.5)
)
# This works because of library(patchwork)
(p_depth + p_idw) /
(p_SK + p_OK + p_UK) +
plot_layout(guides = 'collect')
看起来罪魁祸首在代码的 ggplot 部分。使用 scale_color_viridis() 时，您需要更改 limits 以适合您的数据。由于数据的深度为 0 到 28，因此您应该将限制从 c(50, 100) 更改为 c(0, 30).
# A function to plot rasters
plot_my_gstat_output <- function(raster_object, object_name){
df <- rasterToPoints(raster_object) %>% as_tibble()
colnames(df) <- c("X", "Y", "Z")
ggplot(df, aes(x = X, y = Y, fill = Z)) +
geom_raster() +
ggtitle(label = object_name) +
scale_fill_viridis(option = "B", limits = c(0, 30)) +
theme_void() +
theme(
plot.title = element_text(hjust = 0.5)
)
}
p_idw <- plot_my_gstat_output(raster(idwres), "IDW")
p_SK <- plot_my_gstat_output(raster(SK), "Simple Kriging")
p_OK <- plot_my_gstat_output(raster(OK), "Ordinary Kriging")
p_UK <- plot_my_gstat_output(raster(UK), "Universal Kriging")
# I also want to display sampling locations
p_depth <- ggplot(
data = depth,
mapping = aes(x = X, y = Y, color = Z)
) +
geom_point(size = 3) +
scale_color_viridis(option = "B", limits = c(0, 30)) +
ggtitle(label = "Depths Sampled") +
theme_void() +
theme(
plot.title = element_text(hjust = 0.5)
)
# This works because of library(patchwork)
(p_depth + p_idw) /
(p_SK + p_OK + p_UK) +
plot_layout(guides = 'collect')
输出