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)

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 cont

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 an

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 ju

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 valida

threshold.res

critical threshold for removing potential outliers

progress

if FALSE remove progress bar

Value

An list object containing:
predan object of STFDF-class with column contains prediction or variance
cvcross validation information for points used in prediction, as object of STFDF-class
outpotential outliers, returned as vector of row names of x$cv@sp, only returned if out.remove=FALSE
remstremoved locations as an object of Spatial-class, if out.remove=TRUE
remobsremoved locations with observations as an object of STFDF-class, if out.remove=TRUE

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