Кригинг в одной пространственной точке? - PullRequest
Новая
       12

Кригинг в одной пространственной точке?

0 голосов
/ 23 мая 2018

Этот вопрос может быть связан с моим плохим знанием кригинга - возможно ли рассчитать значение кригеда в одном пространственном месте?Насколько я понимаю, типичный метод кригинга использует пространственную корреляцию, встроенную в разреженные данные, для интерполяции значений во всех точках регулярной сетки, которая имеет такой же пространственный экстент, что и данные.Я хотел бы знать, могу ли я выполнить эту интерполяцию в определенной точке (которая может не попасть в сетку).В качестве примера, приведенный ниже код влияет на концентрацию меди в данных о мезе и накладывает значения на карту мез.Я хотел бы знать, как вычислить значение в криге в "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]

enter image description here

1 Ответ

0 голосов
/ 23 мая 2018

Вы можете предсказать одну точку, и точка может быть за пределами тренировочной сетки, но статистически это не правильно, и если вы уйдете слишком далеко, вы получите только среднее значение.

Пример прогнозированияТочка вне обучения (новые данные): 

> 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"
...