library(sp)
library(spacetime)
library(gstat)
library(plyr)
library(CAST)
library(doParallel)
library(ranger)
# preparing data
data(dtempc)
data(stations)
data(regdata) # covariates, made by mete2STFDF function
regdata@sp@proj4string <- CRS('+proj=longlat +datum=WGS84')
data(tvgms) # ST variogram models
data(tregcoef) # MLR coefficients
lonmin=18 ;lonmax=22.5 ; latmin=40 ;latmax=46
serbia = point.in.polygon(stations$lon, stations$lat, c(lonmin,lonmax,lonmax,lonmin),
c(latmin,latmin,latmax,latmax))
st = stations[ serbia!=0, ] # stations in Serbia approx.
crs = CRS('+proj=longlat +datum=WGS84')
# create STFDF
stfdf <- meteo2STFDF(obs = dtempc,
stations = st,
crs = crs)
# Cross-validation for mean temperature for days "2011-07-05" and "2011-07-06"
# global model is used for regression and variogram
# Overlay observations with covariates
time <- index(stfdf@time)
covariates.df <- as.data.frame(regdata)
names_covar <- names(tregcoef[[1]])[-1]
for (covar in names_covar){
nrowsp <- length(stfdf@sp)
regdata@sp=as(regdata@sp,'SpatialPixelsDataFrame')
ov <- sapply(time, function(i)
if (covar %in% names(regdata@data)) {
if (as.Date(i) %in% as.Date(index(regdata@time))) {
over(stfdf@sp, as(regdata[, i, covar], 'SpatialPixelsDataFrame'))[, covar]
} else {
rep(NA, length(stfdf@sp))
}
} else {
over(stfdf@sp, as(regdata@sp[covar], 'SpatialPixelsDataFrame'))[, covar]
}
)
ov <- as.vector(ov)
if (all(is.na(ov))) {
stop(paste('There is no overlay of data with covariates!', sep = ""))
}
stfdf@data[covar] <- ov
}
# Remove stations out of covariates
for (covar in names_covar){
# count NAs per stations
numNA <- apply(matrix(stfdf@data[,covar],
nrow=nrowsp,byrow= FALSE), MARGIN=1,
FUN=function(x) sum(is.na(x)))
rem <- numNA != length(time)
stfdf <- stfdf[rem,drop= FALSE]
}
# Remove dates out of covariates
rm.days <- c()
for (t in 1:length(time)) {
if(sum(complete.cases(stfdf[, t]@data)) == 0) {
rm.days <- c(rm.days, t)
}
}
if(!is.null(rm.days)){
stfdf <- stfdf[,-rm.days]
}
### Example with STFDF and without parallel processing and without refitting of variogram
results <- cv.strk(data = stfdf,
obs.col = 1, # "tempc"
data.staid.x.y.z = c(1,NA,NA,NA),
reg.coef = tregcoef[[1]],
vgm.model = tvgms[[1]],
sp.nmax = 20,
time.nmax = 2,
type = "LLO",
k = 5,
seed = 42,
refit = FALSE,
progress = TRUE
)
# stplot(results[,,"pred"])
summary(results)
# accuracy
acc.metric.fun(results@data$obs, results@data$pred, "R2")
acc.metric.fun(results@data$obs, results@data$pred, "RMSE")
acc.metric.fun(results@data$obs, results@data$pred, "MAE")
acc.metric.fun(results@data$obs, results@data$pred, "CCC")
Run the code above in your browser using DataLab