The special assumptions of homogeneity in kriging interpolation have led to the use of several types of suitable models to inspect any particular trend before attempting a variogram stimation. While extending the median-polish capacity of table variance analysis in the study of space-time processes with irregular sampling, we have proposed to implement it in the extraction of space-time trends of climatic variables, in this case monthly precipitation. This technique, introduced in cressie 1991, blends an iterative algorithm for removing trends by computing medians throught the spatial and temporal domain and the estimation of the variogram to apply kriging (william 2017). Therefore, this notes attempts to give a brief introduction of STMedianPolish package, giving insight of how every function works in the inspection of trends and how this is integrated into the gstat workflow for spatio temporal kriging.

Examination of monthly precipitation records in a southwestern basin of Colombia.

This data analysis is based on the monthly precipitation records from January 2000 to January 2010 taken from 102 pluviometric stations operated by the Institute of Hydrology, Meteorology and Environmental Studies of Colombia (IDEAM, for its acronym in Spanish). The monitoring network covering this region is irregular due to the physical and meteorological conditions of the different zones, and potential needs of the network users

library(forecast)
library(spacetime)
library(sp)
library(STMedianPolish)
data(DemMeta)
data(HZRMeta)
proj4string(HZRMeta)<-CRS(proj4string(DemMeta))
polygon0 = polygons(HZRMeta)
data(Metadb)
xy0 = SpatialPoints(Metadb[,2:4],CRS(proj4string(DemMeta)))

Rotation of the study area

For practical purposes, the area of study is rotated 40° at northest direction in order to make the sweeps in the east and north spatial direction equivalet to the tables in which the method is based. The following code is published in RPuds in the following web site “https://rpubs.com/geospacedman/rotatespatial”

rotateProj = function(spobj, angle) {
  # get bounding box as spatial points object
  boxpts = SpatialPoints(t(bbox(spobj)), proj4string = CRS(proj4string(spobj)))
  # convert to lat-long
  boxLL = bbox(spTransform(boxpts, CRS("+init=epsg:4326")))
  # find the centre
  llc = apply(boxLL, 1, mean)
  # construct the proj4 string
  prj = paste0("+proj=omerc +lat_0=", llc[2], " +lonc=", llc[1], " +alpha=", 
               angle, " +gamma=0.0 +k=1.000000 +x_0=0.000 +y_0=0.000 +ellps=WGS84 +units=m ")
  # return as a CRS:
  CRS(prj)
}
rtp<-rotateProj(HZRMeta,40)
polygon1 = polygons(spTransform(HZRMeta,rtp))
xy1 = spTransform(xy0,rtp)

par(mfrow=c(1,2))
plot(polygon0,border="blue",axes=TRUE,las=1)
points(xy0, pch = 3, cex = 0.5, col = "red")
plot(polygon1,border="blue",axes=TRUE,las=1)
points(xy1, pch = 3, cex = 0.5, col = "red")

Metadb[,2:4]<-xy1@coords

Grid prediction

time

library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
x<-matrix(0,1,121)
for(i in 1:121){
  x[,i] <- 2000 + (seq(0, 120)/12)[i]
}
x<-as.Date (as.yearmon(x), frac = 1)
time = as.POSIXct(x, tz = "GMT")
length(time)
## [1] 121

space

Gridxy<- spsample(polygon1, cellsize=4000, n=5000,"regular")
plot(polygon1)
points(Gridxy, pch = 3, cex = .5, col = "red")

Intersection with the DEM

Gridxyz<-data.frame(Gridxy,over(spTransform(Gridxy,proj4string(DemMeta)),DemMeta))
colnames(Gridxyz)<-c("East", "North","height")
Grid_pred <- STF(sp=SpatialPoints(Gridxyz,CRS(proj4string(xy1))), time=time[c(102,103)])

Create an spatio - temporal object with regular data

MPST<-ConstructMPst(sqrt(0.5+Metadb[,-c(1:4)]),time,pts=Metadb[,2:4],Delta=c(7,6,5))
Mpplot(MPST)

Visualization of spatial traces

We have defined simply 3 iterarions here.

MpSTData<-MedianPolishM(MPST,eps=0, maxiter=3, na.rm=TRUE)
plot(MpSTData)

Temporal prediction effect

library(forecast)
a=ts(MpSTData$effects[[4]], frequency = 12, start = 2000,  end = c(2010,1))
fit <- auto.arima(a)
plot(forecast(fit,h=12,level=95),main="",
     xlab="Years",ylab="Estimates of Effect Temporal"
     ,axes = TRUE,panel.first = grid(10, 10),pch = 16)

Trend Surface

TendenciaGrilla<-splineMPST(Gridxyz,Ef_t=time[c(102,103)],MPST,eps=0.01, maxiter=2)

this study keeps in a simple form the decomposition of four factors that exert a separate and additive influence on the variable.

IDs = paste("ID",1:nrow(TendenciaGrilla))
mydata = data.frame(values = TendenciaGrilla[,5], ID=IDs)
Trend.ST = STFDF(SpatialPixels(Gridxy),time[c(102,103)],mydata)
stplot(Trend.ST,col.regions=bpy.colors(40),par.strip.text = list(cex=0.7),main="Trend surface")

Remove trend, get residuals

residuals<-removetrendMPst(MPST,eps=0.01, maxiter=2)
rain.loc<-Metadb[,c("Station","East","North","Height")]
coordinates(rain.loc) = ~East+North+Height
proj4string(rain.loc) = CRS(proj4string(xy1))
rain_residual = stConstruct(data.frame(Res=residuals[,7]), space = list(values = 1),
                            time, SpatialObj = rain.loc,interval=TRUE)

Variography

Empirical variogram

#VRes = variogram(values ~ 1, rain_residual, cutoff=20000,tlags=0:4,width=4000,
#               assumeRegular=TRUE, na.omit=TRUE)

data(VRes)
print(plot(VRes, ylab = "time lag (days)"), split = c(1,1,1,2), more = TRUE)
library(lattice)
print(plot(VRes,wireframe=T),split = c(1,2,1,2))

Calculation of the Variogram, sum-product model generalized

FitPar_st = function(p, gfn, v, trace = FALSE, ...) {
  mod = gfn(v$spacelag, v$timelag,p,  ...)
  resid = v$gamma - mod
  if (trace)
    print(c(p, MSE = mean(resid^2)))
  mean(resid^2)
}

#Spatial Model
ModSpatial = function(h,p){p[2]*(1-exp(-h/p[3]))}

#Temporal Model
ModTemporal = function(u,p){p[4]*(1-exp(-u/p[5]))+ p[6]*(1-cos(pi*u/180))+
    p[7]*abs(sin(pi*u/180))}

VariogST<-function(h,u,p){ModTemporal(u,p)+ModSpatial(h,p)+p[8]*ModTemporal(u,p)*ModSpatial(h,p)}

#Initial parameters
p1<-c(2,14.5,13900,5.9,29,1.55,3.7,-0.07)
pars.st = optim(p1, FitPar_st, method =  "BFGS",
                gfn = VariogST, v = VRes, hessian=TRUE)
fit_Variog_ST<-VRes
fit_Variog_ST$gamma<-VariogST(VRes$spacelag, VRes$timelag, pars.st$par)
plot(fit_Variog_ST,wireframe=T)

Covariogram

CS = function(h,p){p[2]*exp(-h/p[3])}
CT = function(u,p){p[4]*exp(-u/p[5])+ p[6]*cos(pi*u/180)+p[7]*(1-abs(sin(pi*u/180)))}

Ct_0<-max(ModTemporal(VRes$timelag,pars.st$par))
Cs_0<-max(ModSpatial(VRes$spacelag,pars.st$par))
C00<-max(VariogST(VRes$spacelag, VRes$timelag,pars.st$par))

k1<--pars.st$par[8]
k2<-(C00-Ct_0)/Cs_0
k3<-(C00-Cs_0)/Ct_0 

k2+k1*Ct_0==k3+k1*Cs_0
## [1] TRUE
#Spatio - temporal model of cavariance.
CST<-function(h,u,p){k3*CT(u,p)+k2*CS(h,p)+k1*CT(u,p)*CS(h,p)}

fit_Covariog_ST<-VRes
fit_Covariog_ST$gamma<-CST(VRes$spacelag, VRes$timelag,pars.st$par)
plot(fit_Covariog_ST,wireframe=T)

Kriging space time

Estimation of the spatio-temporal anisotropy, used only for the selection of the k-nearest neighbours and asumming a metric spatio-temporal space Based on the Gstat function estiStAni

library(gstat)
stAni<-estiStAni(VRes, interval=c(10, 150))
PredictValue<-krige0STlocalMP(data=rain_residual,newdata=Grid_pred,p=pars.st$par,model=CST,k=15,stAni)

Residuals Interpolation on June and july 2008

IDs = paste("ID",1)
mydata = data.frame(PredictValue[,5], ID=IDs)
wind.ST1 = STFDF(SpatialPixels(Gridxy),time[c(102,103)],mydata)
stplot(wind.ST1,col.regions=bpy.colors(40),par.strip.text = list(cex=0.7)
       ,main="Kriging ordinary residuals: Prediction surface")

#Original Values

TendenciaGrilla$forecast<-(mydata[,1]+TendenciaGrilla[,5]-0.5)^2

IDs = paste("ID",1:length(TendenciaGrilla[,c("forecast")]))
mydata1 = data.frame(values =TendenciaGrilla[,c("forecast")] , ID=IDs)
wind.ST1 = STFDF(SpatialPixels(Gridxy),time[c(102,103)],mydata1)
stplot(wind.ST1,col.regions=bpy.colors(40),par.strip.text = list(cex=0.7)
       ,main="Forecast: Prediction surface Monthly Precipitation")