pred.strk: Spatio-temporal regression kriging

Description

Function for spatio-temporal regression kriging prediction based on krigeST. The prediction is made for raster objects, i.e. for STFDF-class objects.

Usage

pred.strk(temp,zcol = 1, newdata, pred.id = "tempPred", zero.tol = 0, 
  dynamic.cov = c(1, 2), static.cov = c(1, 2),
reg.coef=list(
tmean=c(-0.126504415,0.4051734447,0.4943247727,0.0001837527,-0.0189207588),  
tmin = c(-0.9825601517,0.5672140021,0.3344561638, 0.0003119777,-0.0243629638),
tmax = c(1.7873573081,0.350228076, 0.5569091092, 0.0002571338,-0.0012988123)
)[['tmean']],
vgm.model=list( tmean=vgmST("sumMetric",
                                space=vgm( 14.13, "Sph", 5903, 1.933),
                                time =vgm(0, "Sph",  0.1, 0),
                                joint=vgm(9.06, "Sph", 2054, 0.474),
                                 stAni=497.9),
                   tmin = vgmST("sumMetric",
                                space=vgm( 22.682, "Sph", 5725, 3.695),
                                time =vgm(0, "Sph",  0.1, 0),
                                joint=vgm(9.457, "Sph",1888, 1.67),
                                 stAni=485),
                    tmax = vgmST("sumMetric",
                                 space=vgm( 8.31, "Sph", 4930, 2.872),
                                 time =vgm(0, "Sph",  0.1, 0),
                                 joint=vgm(11.175, "Sph", 2117, 1.75),
                                 stAni=527) ) [['tmean']] ,
tiling = FALSE, ntiles = 64, parallel.processing = FALSE, cpus = 3, 
sp.nmax = 18, time.nmax = 2,fast = FALSE, computeVar = FALSE, 
do.cv = FALSE, only.cv = FALSE, out.remove = FALSE, threshold.res = 15,progress=TRUE)

Value

An list object containing:

pred: an object of STFDF-class with column contains prediction or variance
cv: cross validation information for points used in prediction, as object of STFDF-class
out: potential outliers, returned as vector of row names of x$cv@sp, only returned if out.remove=FALSE
remst: removed locations as an object of Spatial-class, if out.remove=TRUE
remobs: removed locations with observations as an object of STFDF-class, if out.remove=TRUE

Arguments

temp: object of STFDF-class containing dependent variable (observations) in space and time.
zcol: variable column name or number showing position of dependent variable in temp@data
newdata: dynamic and static covariates as STFDF-class object; spatial and temporal overlay with temp object must be possible
pred.id: identifier of new variable
zero.tol: distance values less than or equal to this threshold value locations are considered as duplicates, see rm.dupl, duplicates are removed to avoid singular covariance matrices in kriging.
dynamic.cov: vector of variable column names or numbers showing position of dynamic covariates in newdat@data; dynamic covariates are spatio-temporal explanatory variables, changing in space and time domain
static.cov: vector of variable column names or numbers showing position of static covariates in newdata@data@sp; static covariates are spatial explanatory variables changing just in space; static in time dimension
reg.coef: linear regression coefficients; order is assumed as intercept, dynamic.cov, static.cov. Coefficients can be specified by user; depending on type, number and order of dynamic and static covariates. At the moment the function contains regression coefficient for mean, minimum and maximum temperature calculated globally for GSOD and ECA&D data set on geometrical temperature trend, MODIS LST-8 day, elevation and TWI, see regdata. Coefficients for mean temperature are defined by default.
vgm.model: spatio-temporal variogram of regression residuals, see vgmST. At the moment the function contains spatio-temporal variogram model on residuals for mean, minimum and maximum temperature calculated globally for GSOD and ECA&D data set. Regression residuals on geometrical temperature trend, MODIS LST-8 day, elevation and TWI, see regdata. Ranges are in km. Spatio-temporal variogram for mean temperatures is defined by default. User can specified own variogram model as vgmST object.
tiling: for simplified local kriging. Area is divided in tiles and kriging calculation is done for each tile separately, number of observation used per tile is defined with sp.nmax and time.nmax. Default is TRUE. If FALSE just temporal local kriging will be applied defined with time.nmax, sp.nmax will be ignored.
ntiles: number of tiles. Default is 64. Each tile at minimum should contain less observations than sp.nmax, ideally each tile should contain observations falling in neighboring tiles.
parallel.processing: if TRUE parallel processing is performed via sfLapply
cpus: number of processing units
sp.nmax: number of nearest spatial observations that should be used for a kriging prediction for each tile
time.nmax: number of nearest time observations that should be used for a kriging prediction
fast: if TRUE tiling, tiling is done twice to avoid edge effect
computeVar: if TRUE, just variance is computed
do.cv: if TRUE, cross validation leave-one-station-out is performed
only.cv: if TRUE, only cross validation leave-one-station-out is performed without prediction
out.remove: if TRUE, potential outliers are removed. Removing procedure is iterative, all location with residual higher than defined threshold (treshold.res) are selected. Only location with highest cross validation residual is removed, than cross validation is done again, the procedure removing one by one location run until all locations have residuals under defined threshold.
threshold.res: critical threshold for removing potential outliers
progress: if FALSE remove progress bar

Author

Milan Kilibarda kili@grf.bg.ac.rs

References

Kilibarda, M., T. Hengl, G. B. M. Heuvelink, B. Graeler, E. Pebesma, M. Percec Tadic, and B. Bajat (2014), Spatio-temporal interpolation of daily temperatures for global land areas at 1 km resolution, J. Geophys. Res. Atmos., 119, 2294-2313, doi:10.1002/2013JD020803.

Examples

Run this code

# prepare data
# load observation - data.frame of mean temperatures
data(dtempc) 
# str(dtempc)
data(stations)
library(sp)
library(spacetime)
library(gstat)

# str(stations)
## lonmin,lonmax,lonmax, lonmin   latmin, latnmin,latmax,latmax
serbia= point.in.polygon(stations$lon, stations$lat, c(18,22.5,22.5,18), c(40,40,46,46))
st= stations[ serbia!=0, ]
# create STFDF
temp <- meteo2STFDF(dtempc,st)
rm(dtempc)
# str(temp)
# Adding CRS
temp@sp@proj4string <- CRS('+proj=longlat +datum=WGS84')

# load covariates for mean temperatures
data(regdata)
# str(regdata)
regdata@sp@proj4string <- CRS('+proj=longlat +datum=WGS84')

# Calculate prediction of mean temperatures for "2011-07-05" 
# global model is used for regression and variogram
# load precalculated variograms
data(tvgms)
data(tregcoef)
res= pred.strk(temp,zcol=1, newdata= regdata[,1,drop=FALSE], 
              reg.coef=tregcoef[[1]] ,vgm.model=tvgms[[1]], progress=FALSE )

## plot prediction
# stplot(res$pred, col.regions=bpy.colors())


# t1= temp[regdata@sp,]
# # create fake observations
# t1@data$tempc[seq(1,120,by=8)] =35
# 
# 
# res= pred.strk(t1,zcol=1, newdata= regdata[,1:2], 
#                reg.coef=tregcoef[[1]], vgm.model=tvgms[[1]] , 
#                threshold.res=5, do.cv=T, out.remove = T)
# # plot cross validation residuals 
# stplot(res$cv[,,'resid.cv'] , col.regions=bpy.colors())
# 
# # plot locations of removed stations
# spplot(res$remst, zcol='station_name' , col.regions=bpy.colors())
# #plot removed stations as time-series
# row.names(res$remobs@sp) = res$remst$station_name
# res$remobs[,1:2,c('tempc','pred.cv')]
# stplot(res$remobs[,1:2,c('tempc','pred.cv')], mode='tp')

## Calculate prediction of mean temperature for "2011-07-05" "2011-07-06"
## only MODIS is used as covariate

# modisVGM =vgmST("sumMetric",space=vgm( 18.27, "Sph", 6000, 3.22),
#                                           time =vgm(0, "Sph",  0.1, 0),
#                                           joint=vgm(8.34, "Sph", 2349, 1.80),
#                                           stAni=583)
# attr(modisVGM,"temporal unit") = "days"                                           

# rkmod <-  pred.strk(temp,zcol=1, newdata= STFDF(regdata@sp,
#                time=as.POSIXct("2011-07-05"), endTime=as.POSIXct("2011-07-06"), 
#                data=regdata[,1]@data) , threshold.res=10, 
#              dynamic.cov='modis', static.cov=NULL,
#               reg.coef= c(-0.23,0.7303284),
#               vgm.model= modisVGM  )
                                           
## coefficients and variogram is calculated globally for GSOD and ECA&D obs. for 2011 year 

# stplot(rkmod$pred, col.regions=bpy.colors())

## parallel processing
# library(snowfall)
# rkmod <-  pred.strk(temp,zcol=1, 
#                   newdata= STFDF(regdata@sp,
#                    time=as.POSIXct("2011-07-05"), endTime=as.POSIXct("2011-07-06"), 
#                    data=regdata[,1]@data) ,
#                     threshold.res=10, 
#                     dynamic.cov='modis', static.cov=NULL,
#                     reg.coef= c(-0.23,0.7303284),
#                     vgm.model= modisVGM, parallel.processing=TRUE)

Run the code above in your browser using DataLab