单个空间点的克里金法?
Kriging at a single spatial point?
这个问题可能与我对克里格法的了解不足有关——是否可以在单个空间位置计算克里格值?据我了解,典型的克里金方法使用嵌入在稀疏数据中的空间相关性来插入与数据具有相同空间范围的规则网格上所有点的值。我想知道的是我是否可以在特定点(可能不会落在网格上)上进行插值。例如,下面的代码对默兹数据中的铜浓度进行克里金法处理,并将这些值叠加在默兹地图上。我想知道如何计算地图中心附近 "Maasband" 处的克里格值。
谢谢。
# transform meuse data to SpatialPointsDataFrame
suppressMessages(library(sp))
data(meuse)
coordinates(meuse) <- ~ x + y
proj4string(meuse) <- CRS("+proj=stere
+lat_0=52.15616055555555 +lon_0=5.38763888888889
+k=0.999908 +x_0=155000 +y_0=463000
+ellps=bessel +units=m +no_defs
+towgs84=565.2369,50.0087,465.658,
-0.406857330322398,0.350732676542563,-1.8703473836068, 4.0812")
# define a regular grid for kriging
xrange <- range(as.integer(meuse@coords[, 1])) + c(0,1)
yrange <- range(as.integer(meuse@coords[, 2]))
grid <- expand.grid(x = seq(xrange[1], xrange[2], by = 40),
y = seq(yrange[1], yrange[2], by = 40))
coordinates(grid) <- ~ x + y
gridded(grid) <- T
# do kriging
suppressMessages(library(automap))
krg <- autoKrige(formula = copper ~ 1,
input_data = meuse,
new_data = grid)
# extract kriged data
krg_df <- data.frame(krg$krige_output@coords,
pred = krg$krige_output@data$var1.pred)
# transform to SpatialPointsDF & assign original (meuse) projection
krg_spdf <- krg_df
coordinates(krg_spdf) <- ~ x + y
proj4string(krg_spdf) <- proj4string(meuse)
# transform again to longlat coordinates (for overlaying on google map below)
krg_spdf <- spTransform(krg_spdf, CRS("+init=epsg:4326"))
krg_df <- data.frame(krg_spdf@coords, pred = krg_spdf@data$pred)
# get meuse map and overlay kriged data
suppressMessages(library(ggmap))
suppressMessages(library(RColorBrewer))
lon <- range(krg_df$x)
lat <- range(krg_df$y)
meuse_map <- get_map(location = c(lon = mean(lon), lat = mean(lat)),
zoom = 13)
print(ggmap(meuse_map, extent = "normal", maprange = F) +
stat_summary_2d(aes(x = x, y = y, z = pred),
binwidth = c(0.001,0.001),
alpha = 0.5,
data = krg_df) +
scale_fill_gradientn(name = "Copper",
colours = brewer.pal(6, "YlOrRd")) +
coord_cartesian(xlim = lon, ylim = lat, expand = 0) +
theme(aspect.ratio = 1))
# geocode for Maasband
longlat <- as.numeric(geocode("Maasband"))
x0 <- longlat[1]
y0 <- longlat[2]
您可以预测单个点,该点可以在训练网格之外,但从统计上讲它是不正确的,如果距离太远,您将只能得到一个平均值。
训练外点预测示例(newdata):
> range(coordinates(meuse)[,1])
[1] 178605 181390
> range(coordinates(meuse)[,2])
[1] 329714 333611
> newdata <- data.frame(x = 178600,
+ y = 329710)
> coordinates(newdata) <- ~ x + y
> krg <- autoKrige(formula = copper ~ 1,
+ input_data = meuse,
+ new_data = newdata)
[using ordinary kriging]
> krg
$krige_output
coordinates var1.pred var1.var var1.stdev
1 (178600, 329710) 43.10343 538.8824 23.21384
$exp_var
np dist gamma dir.hor dir.ver id
1 17 59.33470 142.2353 0 0 var1
2 36 86.01449 288.4722 0 0 var1
3 114 131.02870 259.1184 0 0 var1
4 149 176.18845 417.4295 0 0 var1
5 184 226.75652 322.5353 0 0 var1
6 711 337.60359 473.3298 0 0 var1
7 830 502.04773 529.3259 0 0 var1
8 1349 713.21485 594.9974 0 0 var1
9 1314 961.27179 628.1739 0 0 var1
10 1139 1213.41157 612.3446 0 0 var1
11 1355 1506.55052 583.7055 0 0 var1
$var_model
model psill range kappa
1 Nug 29.45239 0.0000 0.0
2 Ste 592.45663 365.6081 0.4
$sserr
[1] 78.06413
attr(,"class")
[1] "autoKrige" "list"
这个问题可能与我对克里格法的了解不足有关——是否可以在单个空间位置计算克里格值?据我了解,典型的克里金方法使用嵌入在稀疏数据中的空间相关性来插入与数据具有相同空间范围的规则网格上所有点的值。我想知道的是我是否可以在特定点(可能不会落在网格上)上进行插值。例如,下面的代码对默兹数据中的铜浓度进行克里金法处理,并将这些值叠加在默兹地图上。我想知道如何计算地图中心附近 "Maasband" 处的克里格值。
谢谢。
# transform meuse data to SpatialPointsDataFrame
suppressMessages(library(sp))
data(meuse)
coordinates(meuse) <- ~ x + y
proj4string(meuse) <- CRS("+proj=stere
+lat_0=52.15616055555555 +lon_0=5.38763888888889
+k=0.999908 +x_0=155000 +y_0=463000
+ellps=bessel +units=m +no_defs
+towgs84=565.2369,50.0087,465.658,
-0.406857330322398,0.350732676542563,-1.8703473836068, 4.0812")
# define a regular grid for kriging
xrange <- range(as.integer(meuse@coords[, 1])) + c(0,1)
yrange <- range(as.integer(meuse@coords[, 2]))
grid <- expand.grid(x = seq(xrange[1], xrange[2], by = 40),
y = seq(yrange[1], yrange[2], by = 40))
coordinates(grid) <- ~ x + y
gridded(grid) <- T
# do kriging
suppressMessages(library(automap))
krg <- autoKrige(formula = copper ~ 1,
input_data = meuse,
new_data = grid)
# extract kriged data
krg_df <- data.frame(krg$krige_output@coords,
pred = krg$krige_output@data$var1.pred)
# transform to SpatialPointsDF & assign original (meuse) projection
krg_spdf <- krg_df
coordinates(krg_spdf) <- ~ x + y
proj4string(krg_spdf) <- proj4string(meuse)
# transform again to longlat coordinates (for overlaying on google map below)
krg_spdf <- spTransform(krg_spdf, CRS("+init=epsg:4326"))
krg_df <- data.frame(krg_spdf@coords, pred = krg_spdf@data$pred)
# get meuse map and overlay kriged data
suppressMessages(library(ggmap))
suppressMessages(library(RColorBrewer))
lon <- range(krg_df$x)
lat <- range(krg_df$y)
meuse_map <- get_map(location = c(lon = mean(lon), lat = mean(lat)),
zoom = 13)
print(ggmap(meuse_map, extent = "normal", maprange = F) +
stat_summary_2d(aes(x = x, y = y, z = pred),
binwidth = c(0.001,0.001),
alpha = 0.5,
data = krg_df) +
scale_fill_gradientn(name = "Copper",
colours = brewer.pal(6, "YlOrRd")) +
coord_cartesian(xlim = lon, ylim = lat, expand = 0) +
theme(aspect.ratio = 1))
# geocode for Maasband
longlat <- as.numeric(geocode("Maasband"))
x0 <- longlat[1]
y0 <- longlat[2]
您可以预测单个点,该点可以在训练网格之外,但从统计上讲它是不正确的,如果距离太远,您将只能得到一个平均值。
训练外点预测示例(newdata):
> range(coordinates(meuse)[,1])
[1] 178605 181390
> range(coordinates(meuse)[,2])
[1] 329714 333611
> newdata <- data.frame(x = 178600,
+ y = 329710)
> coordinates(newdata) <- ~ x + y
> krg <- autoKrige(formula = copper ~ 1,
+ input_data = meuse,
+ new_data = newdata)
[using ordinary kriging]
> krg
$krige_output
coordinates var1.pred var1.var var1.stdev
1 (178600, 329710) 43.10343 538.8824 23.21384
$exp_var
np dist gamma dir.hor dir.ver id
1 17 59.33470 142.2353 0 0 var1
2 36 86.01449 288.4722 0 0 var1
3 114 131.02870 259.1184 0 0 var1
4 149 176.18845 417.4295 0 0 var1
5 184 226.75652 322.5353 0 0 var1
6 711 337.60359 473.3298 0 0 var1
7 830 502.04773 529.3259 0 0 var1
8 1349 713.21485 594.9974 0 0 var1
9 1314 961.27179 628.1739 0 0 var1
10 1139 1213.41157 612.3446 0 0 var1
11 1355 1506.55052 583.7055 0 0 var1
$var_model
model psill range kappa
1 Nug 29.45239 0.0000 0.0
2 Ste 592.45663 365.6081 0.4
$sserr
[1] 78.06413
attr(,"class")
[1] "autoKrige" "list"