Spatio-temporal variograms
Guillaume Blanchet & Steve Vissault
November 21th, 2019
Prepare the R environment
We load the packages needed for this practice.
# Spatial
library(sf)
library(sp) # for some analysis
library(dplyr)
# Analysis
library(gstat)
library(spacetime)
# visualization
library(mapview)Get access to the data
Load data
For this section, we will use the Air quality data obtained from the airBase European air quality data base. Daily averages for rural background stations in Germany, 1998-2009. Using a nationwide network of monitoring sites, EPA has developed ambient air quality trends for particle pollution, also called Particulate Matter (PM). PM10 used in this dataset describes inhalable particles, with diameters that are generally 10 micrometers and smaller.
## station_id pm10 dates lng lat
## 1 DEBB051 8.688 2001-06-08 14.05746 52.54136
## 2 DEBB051 7.354 2001-09-10 14.05746 52.54136
## 3 DEBB051 15.729 2001-10-10 14.05746 52.54136
## 4 DEBB051 26.312 2001-05-03 14.05746 52.54136
## 5 DEBB051 13.958 2001-10-28 14.05746 52.54136
## 6 DEBB051 39.188 2001-10-16 14.05746 52.54136Prepare the data
To reduce the computation time, we focus on the period 2005-2010.
pm10_2005_2010 <- subset(pm10, dates > as.Date("2005-01-01") &
dates < as.Date("2010-12-31"))
STIDF_pm10_2005_2010 <- stConstruct(pm10_2005_2010, c("lng","lat"),"dates", crs = CRS("+init=epsg:4326"))
STFDF_pm10_2005_2010 <- as(STIDF_pm10_2005_2010, "STFDF")
summary(STFDF_pm10_2005_2010)## Object of class STFDF
## with Dimensions (s, t, attr): (53, 1825, 2)
## [[Spatial:]]
## Object of class SpatialPoints
## Coordinates:
## min max
## lng 6.28107 14.78617
## lat 47.80847 54.92497
## Is projected: FALSE
## proj4string :
## [+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
## +towgs84=0,0,0]
## Number of points: 53
## [[Temporal:]]
## Index timeIndex
## Min. :2005-01-02 Min. : 1
## 1st Qu.:2006-04-03 1st Qu.: 457
## Median :2007-07-03 Median : 913
## Mean :2007-07-03 Mean : 913
## 3rd Qu.:2008-10-01 3rd Qu.:1369
## Max. :2009-12-31 Max. :1825
## [[Data attributes:]]
## station_id pm10
## DETH026: 1821 Min. : 0.560
## DENI063: 1815 1st Qu.: 9.275
## DEBY047: 1815 Median : 13.854
## DEMV017: 1814 Mean : 16.261
## DEUB028: 1814 3rd Qu.: 20.333
## (Other):65675 Max. :269.079
## NA's :21971 NA's :21971Build empricial spatio-temporal variogram
##
|
| | 0%
|
|=========== | 17%
|
|====================== | 33%
|
|================================ | 50%
|
|=========================================== | 67%
|
|====================================================== | 83%
|
|=================================================================| 100%
Model variogram
Declare the model
The argument sill that is supplied to vgmST (below) defines the joint spatio-temporal sill.
# Declare the model variogram
separableModel <- vgmST("separable",
method = "Nelder-Mead", # no lower & upper needed
space = vgm(psill = 0.9, model = "Exp", range = 123, nugget = 0.1),
time = vgm(psill = 0.9, model = "Exp", range = 2.9, nugget = 0.1),
sill = 100)The function fit.StVariogram will use the optim function to optimize and find the best set of parameters for the model variogram using the mean squared errors.
# fit the separableModel variogram
separableModel <- fit.StVariogram(pm10_vario, separableModel,
method="BFGS")## Warning in optim(extractPar(model), fitFun, ..., method = method, lower =
## lower, : bounds can only be used with method L-BFGS-B (or Brent)# metric model
metricModel <- vgmST("metric",
method = "Nelder-Mead", # no lower & upper needed
joint=vgm(0.9, "Exp", 123, 0.1),
stAni=100)
metricModel <- fit.StVariogram(pm10_vario, metricModel,
method="BFGS")## Warning in optim(extractPar(model), fitFun, ..., method = method, lower =
## lower, : bounds can only be used with method L-BFGS-B (or Brent)But many other possibilities
# product sum model: spatial and temporal nugget will be ignored and kept
# constant at 0. Only a joint nugget is used.
prodSumModel <- vgmST("productSum",
space=vgm(39, "Sph", 343, 0),
time= vgm(36, "Exp", 3, 0),
k=15)
# sum metric model: spatial, temporal and joint nugget will be estimated
sumMetricModel <- vgmST("sumMetric",
space=vgm( 6.9, "Lin", 200, 3.0),
time =vgm(10.3, "Lin", 15, 3.6),
joint=vgm(37.2, "Exp", 84,11.7),
stAni=77.7)
# simplified sumMetric model, only a overall nugget is fitted. The spatial,
# temporal and jont nuggets are set to 0.
simpleSumMetricModel <- vgmST("simpleSumMetric",
space=vgm(20,"Lin", 150, 0),
time =vgm(20,"Lin", 10, 0),
joint=vgm(20,"Exp", 150, 0),
nugget=1, stAni=15)Assess model quality
You can get the mean squared errors within the output of the function fit.StVariogram() with the following command line and compare them.
## [1] 46.4424## [1] 60.661library(lattice)
plot(pm10_vario, list(separableModel, metricModel), all=T, wireframe=T, zlim=c(0,120),
zlab=NULL,
xlab=list("distance (km)", rot=30),
ylab=list("time lag (days)", rot=-35),
scales=list(arrows=F, z = list(distance = 5)))
Spatio-temporal kriging
We subset 3 days within the STSDF in order to reduce the krigging computation time.
## Formal class 'STFDF' [package "spacetime"] with 4 slots
## ..@ data :'data.frame': 159 obs. of 2 variables:
## .. ..$ station_id: Factor w/ 46 levels "DEBB053","DEBE032",..: 1 2 3 4 5 6 7 8 9 10 ...
## .. ..$ pm10 : num [1:159] 11.9 12.4 15.8 14.8 14.3 ...
## ..@ sp :Formal class 'SpatialPoints' [package "sp"] with 3 slots
## .. .. ..@ coords : num [1:53, 1:2] 14.02 13.23 13.65 9.57 7.76 ...
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:53] "659" "6179" "9814" "11736" ...
## .. .. .. .. ..$ : chr [1:2] "lng" "lat"
## .. .. ..@ bbox : num [1:2, 1:2] 6.28 47.81 14.79 54.92
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "lng" "lat"
## .. .. .. .. ..$ : chr [1:2] "min" "max"
## .. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. .. .. ..@ projargs: chr "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
## ..@ time :An 'xts' object on 2005-06-01/2005-06-03 containing:
## Data: int [1:3, 1] 151 152 153
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr "timeIndex"
## Indexed by objects of class: [Date] TZ: UTC
## xts Attributes:
## NULL
## ..@ endTime: POSIXct[1:3], format: "2005-06-02" ...We create the empty spatial grid.
library(raster)
sp_area <- STFDF_pm10_3d@sp
proj <- proj4string(STFDF_pm10_3d@sp)
rs_area <- raster(extent(sp_area), res = 0.05)
sp_grid_area <- SpatialPoints(rasterToPoints(rs_area), CRS(proj))Perform the krigging
DE_pred <- STF(sp=sp_grid_area, time=STFDF_pm10_3d@time)
DE_kriged <- krigeST(pm10~1, data=STIDF_pm10_3d, newdata=DE_pred,
modelList=separableModel)
stplot(DE_kriged)