[R-sig-Geo] Obtaining covariance matrix from space time variogram

[Edzer Pebesma]

You can reduce your global problem to a smaller, local problem and get 
it, or you can modify the code of krigeST to save the covariance matrix 
e.g. to your global environment using <<- or assign().

In local kriging the covariance matrix is a thing that you compute and 
then abandon, so there's little point in the function making it 
available to users.

On 18/05/2022 11:33, Philibert Weenink wrote:
> Hi all,
> 
> I'm trying to obtain a full covariance matrix for a STIDF object
> (spacetime library).
> I use the vgmST function from gstat to specify a functional form based
> on regression residuals of a subset.
> 
> While the krigeST function provides the option to compute such a full
> covariance matrix, this only works for global kriging.
> As this requires a lot of computing power, this is no option for me.
> Hence, I am looking for a method to obtain the covariance matrix
> directly.
> 
> While I found various functions to do so with purely spatial
> variograms, this was not the case for space time variograms.
> Any help will be appreciated.
> 
> Kind regards,
> Philibert
> 
> Example Script:
> 
> #--R-Script--#
> 
> library(sp)
> library(spacetime)
> library(gstat)
> data(air)
> 
> suppressWarnings(proj4string(stations) <- CRS(proj4string(stations)))
> rural = STFDF(stations, dates, data.frame(PM10 = as.vector(air)))
> rr <- rural[,"2005-06-01/2005-09-01"]
> rr <- as(rr,"STSDF")
> 
> vgram <- variogramST(formula = PM10 ~ 1, data = rr, tunit = "days",
> tlags = seq(0,3,1), cutoff = 500, cores = 7)
> 
> sumMetricVgm <- vgmST("sumMetric",
>                        space = vgm( 8.5, "Exp", 196.6,  3),
>                        time  = vgm( 25, "Sph",   1.1,  2),
>                        joint = vgm(5, "Exp", 136.6, 12),
>                        stAni = 51.7)
> 
> VAR <- fit.StVariogram(vgram, sumMetricVgm,
>                         lower = c(0,0,0,0,0,0,0,0,0,0),,
>                         method="L-BFGS-B",
>                         fit.method = 7)
> plot(vgram,VAR, all = TRUE)
> 
> #Thus far I have tried the krige method, but this is not an option as
> the newdata contains a larger set.
> kr <- krigeST(PM10 ~ 1, data = rr, newdata = rr, VAR, computeVar =
> TRUE, fullCovariance = TRUE)
> 
> #While the local kriging option is computably feasable it does not
> provide full covariance matrices.
> lkr <- krigeST(PM10 ~ 1, data = rr, newdata = rural, VAR, nmax = 10,
> computeVar = TRUE, fullCovariance = TRUE)
> 
> #--R-Script--#
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

-- 
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081