library(magrittr)
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:magrittr':
##
## extract
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(FRK)
##
## Attaching package: 'FRK'
## The following object is masked from 'package:raster':
##
## distance
library(spacetime)
library(sp)
library(rgdal)
## rgdal: version: 1.2-8, (SVN revision 663)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
## Path to GDAL shared files: /usr/share/gdal/2.2
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.2-5
library(mapview)
## Loading required package: leaflet
library(fst)
# devtools::install_github("andrewzm/dggrids")
setwd("/home/michael/Dropbox/BGU/Itai_Kloog/p_46_FRK_Kriging_test/")
# Read data
dat = read.fst("01_raw/modX2009.fst")
# Rename
dat$lon = dat$Longitude
dat$lat = dat$Latitude
dat$Longitude = NULL
dat$Latitude = NULL
# Subset days
dat = dat[dat$day < (min(dat$day) + 11), ]
# Standard deviation of measurements (cannot be estimated from data in spatio-temporal kriging)
dat$std = 1
# Observations to 'spacetime' object
STObj = stConstruct(
x = dat, # dataset
space = c("lon","lat"), # spatial fields
time = "day", # time field
crs = CRS("+proj=longlat +ellps=sphere"), # CRS
interval = TRUE # time reflects an interval
)
# Time to POSIXct
STObj@time = xts::xts(zoo::coredata(STObj@time), as.POSIXct(as.character(time(STObj@time))))
###################################
# Plot - locations
sp = as(STObj, "Spatial")
plot(sp)
###################################
# Prediction (BAU) grid
grid_BAUs = auto_BAUs(
manifold = STsphere(), # space-time field on sphere
data = STObj,
cellsize = c(0.1, 0.1, 1),
# cellsize = c(1, 1, 1),
type = "grid",
tunit = "days" # time spacing in days
)
grid_BAUs$fs = 1
###################################
# Plot
sp = as(grid_BAUs, "Spatial")
plot(sp)
###################################
# Basis functions
G_spatial = auto_basis(
manifold = sphere(), # functions on sphere
data = as(STObj,"Spatial"), # collapse time out
nres = 3, # use three DGGRID resolutions
prune = 15, # prune basis functions
type = "bisquare", # bisquare basis functions
subsamp = 2000, # use only 2000 data points for pruning
isea3h_lo = 2
)
G_temporal = local_basis(
manifold = real_line(), # functions on real line
loc = matrix(c(2,7,12)), # location parameter
scale = rep(3,3), # scale parameter
type = "Gaussian"
)
G_spacetime = TensorP(G_spatial,G_temporal)
###############################################################################################
# Cross-validation
result = NULL
for(i in unique(STObj[["fold"]])) {
# Model
f = error ~ 1 # formula for fixed effects
S = SRE(
f = f, # formula
data = list(STObj[STObj[["fold"]] != i, ]), # spatio-temporal object
basis = G_spacetime, # space-time basis functions
BAUs = grid_BAUs, # space-time BAUs
est_error = FALSE, # do not estimate measurement error
average_in_BAU = TRUE
) # average data that fall inside BAUs
# Fit
S = SRE.fit(
SRE_model = S, # SRE model
n_EM = 2, # max. EM iterations
tol = 0.01
) # convergence criteria
# Predict
grid_BAUs = predict(S, obs_fs = TRUE) # predict only at select days
# Extract
x = as.data.frame(STObj[STObj[["fold"]] == i, ])
y = over(STObj[STObj[["fold"]] == i, ], grid_BAUs)
y$error_pred = y$mu
z = cbind(x, y[, "error_pred", drop = FALSE])
result = rbind(result, z)
}
## Normalising basis function evaluations at BAU level ...
## Binning data ...
## Binned data in 0.491 seconds
## Binned data in 0.481 seconds
## Binned data in 0.441 seconds
## Binned data in 0.445 seconds
## Binned data in 0.433 seconds
## Binned data in 0.415 seconds
## Binned data in 0.48 seconds
## Binned data in 0.412 seconds
## Binned data in 0.439 seconds
## Binned data in 0.446 seconds
## Binned data in 0.439 seconds
##
|
| | 0%
|
|================================ | 50%
|
|=================================================================| 100%
## Maximum EM iterations reached
## Normalising basis function evaluations at BAU level ...
## Binning data ...
## Binned data in 0.49 seconds
## Binned data in 0.444 seconds
## Binned data in 0.416 seconds
## Binned data in 0.516 seconds
## Binned data in 0.408 seconds
## Binned data in 0.454 seconds
## Binned data in 0.456 seconds
## Binned data in 0.449 seconds
## Binned data in 0.439 seconds
## Binned data in 0.486 seconds
## Binned data in 0.461 seconds
##
|
| | 0%
|
|================================ | 50%
|
|=================================================================| 100%
## Maximum EM iterations reached
## Normalising basis function evaluations at BAU level ...
## Binning data ...
## Binned data in 0.534 seconds
## Binned data in 0.482 seconds
## Binned data in 0.453 seconds
## Binned data in 0.55 seconds
## Binned data in 0.481 seconds
## Binned data in 0.507 seconds
## Binned data in 0.461 seconds
## Binned data in 0.442 seconds
## Binned data in 0.42 seconds
## Binned data in 0.489 seconds
## Binned data in 0.431 seconds
##
|
| | 0%
|
|================================ | 50%
|
|=================================================================| 100%
## Maximum EM iterations reached
## Normalising basis function evaluations at BAU level ...
## Binning data ...
## Binned data in 0.511 seconds
## Binned data in 0.463 seconds
## Binned data in 0.432 seconds
## Binned data in 0.498 seconds
## Binned data in 0.424 seconds
## Binned data in 0.449 seconds
## Binned data in 0.52 seconds
## Binned data in 0.431 seconds
## Binned data in 0.459 seconds
## Binned data in 0.485 seconds
## Binned data in 0.458 seconds
##
|
| | 0%
|
|================================ | 50%
|
|=================================================================| 100%
## Maximum EM iterations reached
## Normalising basis function evaluations at BAU level ...
## Binning data ...
## Binned data in 0.505 seconds
## Binned data in 0.427 seconds
## Binned data in 0.451 seconds
## Binned data in 0.473 seconds
## Binned data in 0.464 seconds
## Binned data in 0.438 seconds
## Binned data in 0.509 seconds
## Binned data in 0.471 seconds
## Binned data in 0.426 seconds
## Binned data in 0.501 seconds
## Binned data in 0.431 seconds
##
|
| | 0%
|
|================================ | 50%
|
|=================================================================| 100%
## Maximum EM iterations reached
## Normalising basis function evaluations at BAU level ...
## Binning data ...
## Binned data in 0.509 seconds
## Binned data in 0.474 seconds
## Binned data in 0.424 seconds
## Binned data in 0.499 seconds
## Binned data in 0.423 seconds
## Binned data in 0.451 seconds
## Binned data in 0.488 seconds
## Binned data in 0.459 seconds
## Binned data in 0.453 seconds
## Binned data in 0.499 seconds
## Binned data in 0.459 seconds
##
|
| | 0%
|
|================================ | 50%
|
|=================================================================| 100%
## Maximum EM iterations reached
## Normalising basis function evaluations at BAU level ...
## Binning data ...
## Binned data in 0.561 seconds
## Binned data in 0.462 seconds
## Binned data in 0.456 seconds
## Binned data in 0.498 seconds
## Binned data in 0.482 seconds
## Binned data in 0.427 seconds
## Binned data in 0.548 seconds
## Binned data in 0.532 seconds
## Binned data in 0.453 seconds
## Binned data in 0.469 seconds
## Binned data in 0.485 seconds
##
|
| | 0%
|
|================================ | 50%
|
|=================================================================| 100%
## Maximum EM iterations reached
## Normalising basis function evaluations at BAU level ...
## Binning data ...
## Binned data in 0.527 seconds
## Binned data in 0.504 seconds
## Binned data in 0.501 seconds
## Binned data in 0.461 seconds
## Binned data in 0.449 seconds
## Binned data in 0.422 seconds
## Binned data in 0.495 seconds
## Binned data in 0.427 seconds
## Binned data in 0.481 seconds
## Binned data in 0.662 seconds
## Binned data in 0.522 seconds
##
|
| | 0%
|
|================================ | 50%
|
|=================================================================| 100%
## Maximum EM iterations reached
## Normalising basis function evaluations at BAU level ...
## Binning data ...
## Binned data in 0.533 seconds
## Binned data in 0.483 seconds
## Binned data in 0.45 seconds
## Binned data in 0.522 seconds
## Binned data in 0.482 seconds
## Binned data in 0.451 seconds
## Binned data in 0.538 seconds
## Binned data in 0.463 seconds
## Binned data in 0.494 seconds
## Binned data in 0.505 seconds
## Binned data in 0.489 seconds
##
|
| | 0%
|
|================================ | 50%
|
|=================================================================| 100%
## Maximum EM iterations reached
## Normalising basis function evaluations at BAU level ...
## Binning data ...
## Binned data in 0.527 seconds
## Binned data in 0.454 seconds
## Binned data in 0.484 seconds
## Binned data in 0.489 seconds
## Binned data in 0.531 seconds
## Binned data in 0.446 seconds
## Binned data in 0.527 seconds
## Binned data in 0.447 seconds
## Binned data in 0.474 seconds
## Binned data in 0.534 seconds
## Binned data in 0.454 seconds
##
|
| | 0%
|
|================================ | 50%
|
|=================================================================| 100%
## Maximum EM iterations reached
###############################################################################################
# Cross-validation Plot
plot(error_pred ~ error, result)
fit = lm(error_pred ~ error, result)
summary(fit)
##
## Call:
## lm(formula = error_pred ~ error, data = result)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0958 -0.2779 -0.0589 0.2880 1.7954
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.016554 0.014157 1.169 0.242
## error 0.005507 0.004585 1.201 0.230
##
## Residual standard error: 0.5453 on 1482 degrees of freedom
## Multiple R-squared: 0.0009725, Adjusted R-squared: 0.0002984
## F-statistic: 1.443 on 1 and 1482 DF, p-value: 0.2299
abline(fit, col = "red")
###################################
# Variogram - one day
sp = as(STObj[STObj@data$dayint == "14245"], "Spatial")
bubble(sp, "error")
v = automap::autofitVariogram(error ~ 1, sp)
plot(v)
###############################################################################################
# Ordinary Kriging
library(gstat)
library(automap)
sp = as(STObj[STObj@data$dayint == "14245"], "Spatial")
model_data = sp[sp$fold != 3, ]
test_data = sp[sp$fold == 3, ]
v = autofitVariogram(error ~ 1, model_data)
g = gstat(formula = error ~ 1, model = v$var_model, data = model_data)
pred = predict(g, test_data)
## [using ordinary kriging]
plot(pred$var1.pred ~ test_data$error)
fit = lm(error_pred ~ error, result)
summary(fit)
##
## Call:
## lm(formula = error_pred ~ error, data = result)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0958 -0.2779 -0.0589 0.2880 1.7954
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.016554 0.014157 1.169 0.242
## error 0.005507 0.004585 1.201 0.230
##
## Residual standard error: 0.5453 on 1482 degrees of freedom
## Multiple R-squared: 0.0009725, Adjusted R-squared: 0.0002984
## F-statistic: 1.443 on 1 and 1482 DF, p-value: 0.2299
abline(fit, col = "red")
# ###################################
# # Plot - Space-time
# stplot(grid_BAUs[, , "mu"])
# ###################################
# # Plot - Space only
# grid_BAUs[, 3, "mu"] %>% spplot
# ###################################