---
title: "R Scripts related to Analyzing spatio-temporal data with R: 
Everything you always wanted to know -- but were afraid to ask"
author: Ghislain Geniaux and Davide Martinetti and Edith Gabriel and Eric Parent and  Nicolas Desassis and Denis Allard and Thomas Opitz and Liliane Bel and Thomas Romary
date: "`r format(Sys.time(), '%d %B %Y')`"
output:
  pdf_document:
    highlight: tango
    toc: yes
  html_document:
    toc: yes
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(cache = TRUE, tidy = FALSE, dev = 'pdf')
par(mar = c(5, 4, 1, 1) + .1)
rep_base = "Z:/dallard/RESSTE/atelier/SPAT2017_Workshop/Tutorial"
setwd(rep_base)
options(warn=-1)
```

1) Installing R packages
====
```{r install_packages, eval=FALSE }
install.packages('raster')
install.packages('shiny')
install.packages('spacetime')
install.packages('leaflet')
install.packages('foreign')
install.packages('maptools')
install.packages('gstat')
install.packages('plm')
install.packages('sp')
install.packages('xts')
install.packages('rgeos')
install.packages('diveMove')
install.packages('maps')
install.packages('googleVis')
install.packages('rgdal')
install.packages('ggmap')
install.packages('data.table')
install.packages('mapdata')
install.packages('cshapes')
install.packages('rworldmap')
install.packages('dplyr')
install.packages('RColorBrewer')
install.packages('ncdf4')
install.packages('rworldmap')
install.packages('plotKML')
install.packages('tseries')
install.packages('ggfortify')
install.packages('ggplot2')
install.packages('CompRandFld')
install.packages('fields')
install.packages('spam')
install.packages('scatterplot3d')
install.packages('RandomFields')
```


2) Importing data and converting to STSDF and SDIDF spacetime formats.
====

2.a) Importing stations
---
```{r import_data_station, eval=FALSE,tidy=FALSE}

### AirBase stations database (cleaned and purged) and their plot
stations <- read.csv("data/AirBase_v8_stations.csv",header=TRUE,sep="\t")
names(stations)

coord <- cbind(stations$station_longitude_deg,stations$station_latitude_deg)
stations <-stations[which(coord[,2]>20),]
coord <- cbind(stations$station_longitude_deg,stations$station_latitude_deg)

not_duplicated <- which(!duplicated(coord))
stations <- stations[not_duplicated,]
coord <- cbind(stations$station_longitude_deg,stations$station_latitude_deg)

plot(coord)
save(stations,file='data/stations.Rdata')
```

2.b) Importing hourly data
---
```{r import_data_station_obs_hourly, eval=FALSE}
library(data.table)

file=list.files('data/OBS')
length(file)
head(file)

### Select only file containing hourly data (file names with more than 15-17 characters)
filePM10 <- file[grep('PM10',file)]
filePM10 <- filePM10[nchar(filePM10)>17]
filePM25 <- file[grep('PM25',file)]
filePM25 <- filePM25[nchar(filePM25)>17]
fileO3   <- file[grep('O3',file)]
fileO3   <- fileO3[nchar(fileO3)>15]
fileNO2  <- file[grep('NO2',file)]
fileNO2  <- fileNO2[nchar(fileNO2)>16]
filePM10 <- paste('data/OBS/',filePM10,sep="")
filePM25 <- paste('data/OBS/',filePM25,sep="")
fileO3   <- paste('data/OBS/',fileO3,sep="")
fileNO2  <- paste('data/OBS/',fileNO2,sep="")

date_station <- gsub(".csv","",filePM10)
date_station <- gsub("PM10_","",date_station)

### Checking if all files have the same lenght and order
{
  length(filePM10)
  length(filePM25)
  length(fileO3)
  length(fileNO2)
  filePM10h <- gsub(".csv","",filePM10)
  filePM10h <- gsub("PM10_","",filePM10h)
  filePM25h <- gsub(".csv","",filePM25)
  filePM25h <- gsub("PM25_","",filePM25h)
  fileO3h   <- gsub(".csv","",fileO3)
  fileO3h   <- gsub("O3_","",fileO3h)
  fileNO2h  <- gsub(".csv","",fileNO2)
  fileNO2h  <- gsub("NO2_","",fileNO2h)
  all(filePM10h==filePM25h)
  all(filePM10h==fileO3h)
  all(filePM10h==fileNO2h)
}


### Preparing list for rbindList
PM10l <- vector(mode = "list", length = length(filePM10))
PM25l <- vector(mode = "list", length = length(filePM25))
O3l   <- vector(mode = "list", length = length(fileO3))
NO2l  <- vector(mode = "list", length = length(fileNO2))

for(i in 1:length(PM10l)){
  date  <- strptime(substr(date_station[i],10,19),format= "%Y%m%d%H")
  PM10t <- try(read.table(filePM10[i],header = FALSE,sep=" ",na.strings="-999",
                          col.names = c("ID","long","lat","PM10"),
                          colClasses = c("factor","numeric","numeric","numeric")))
  PM10t <- cbind(PM10t,rep(date,nrow(PM10t)))
  PM10l[[i]] <- PM10t
  PM25t <- try(read.table(filePM25[i],header = FALSE,sep=" ",na.strings="-999",
                          col.names = c("ID","long","lat","PM25"),
                          colClasses=c("factor","numeric","numeric","numeric")))
  PM25t <- cbind(PM25t,rep(date,nrow(PM25t)))
  PM25l[[i]] <- PM25t
  NO2t  <- try(read.table(fileNO2[i],header = FALSE,sep=" ",na.strings="-999",
                          col.names = c("ID","long","lat","NO2"),
                          colClasses=c("factor","numeric","numeric","numeric")))
  NO2t <- cbind(NO2t,rep(date,nrow(NO2t)))
  NO2l[[i]] <- NO2t
  O3t <- try(read.table(fileO3[i],header = FALSE,sep=" ",na.strings="-999",
                        col.names = c("ID","long","lat","O3"),
                        colClasses=c("factor","numeric","numeric","numeric")))
  O3t <- cbind(O3t,rep(date,nrow(O3t)))
  O3l[[i]] <- O3t
  cat(i,' ')
}

### Using rbindlist for faster row binding of all dataframes
PM10 <- rbindlist(PM10l)
names(PM10)[5] <- 'date'
PM10$ID2 <- paste(PM10$ID,PM10$date,sep='_') 
PM25 <- rbindlist(PM25l)
names(PM25)[2:3] <- c("long1","lat1")
names(PM25)[5]   <- 'date1' 
PM25$ID2 <- paste(PM25$ID,PM25$date1,sep='_') 
NO2 <- rbindlist(NO2l)
names(NO2)[2:3] <- c("long2","lat2")
names(NO2)[5]   <- 'date2'  
NO2$ID2 <- paste(NO2$ID,NO2$date2,sep='_') 
O3 <- rbindlist(O3l)
names(O3)[2:3] <- c("long3","lat3")
names(O3)[5]   <- 'date3'  
O3$ID2 <- paste(O3$ID,O3$date3,sep='_') 


### Using data.table for faster merge
PM10 <- data.table(PM10)
PM25 <- data.table(PM25)
NO2  <- data.table(NO2)
O3   <- data.table(O3)
res  <- merge(PM10[,.(PM10,ID2,date,long,lat)],PM25[,.(PM25,ID2,date1,long1,lat1)],
              by="ID2",all=TRUE)
res  <- merge(res,NO2[,.(NO2,ID2,date2,long2,lat2)] ,by="ID2",all=TRUE)
res  <- merge(res, O3[,.(O3,ID2,date3,long3,lat3)] ,by="ID2",all=TRUE)
res$ID <- sapply(res$ID2, function(x) strsplit(x,"_")[[1]][1])
res$ID <- as.factor(res$ID)
res    <- data.frame(res)

### Finding locations and dates avoiding NAs
nodate <- is.na(res$date)
res$date[nodate] <- res$date1[nodate]
nodate <- is.na(res$date)
res$date[nodate] <- res$date2[nodate]
nodate <- is.na(res$date)
res$date[nodate] <- res$date3[nodate]

nolong <- is.na(res$long)
res$long[nolong] <- res$long1[nolong]
nolong <- is.na(res$long)
res$long[nolong] <- res$long2[nolong]
nolong <- is.na(res$long)
res$long[nolong] <- res$long3[nolong]

nolat <- is.na(res$lat)
res$lat[nolat] <- res$lat1[nolat]
nolat <- is.na(res$lat)
res$lat[nolat] <- res$lat2[nolat]
nolat <- is.na(res$lat)
res$lat[nolat] <- res$lat3[nolat]

### Building the final dataframe
OBS_hourly <- res[!is.na(res$date),c('ID','long','lat','date','PM10','PM25','NO2','O3')]
summary(OBS_hourly)

save(OBS_hourly,file='data/OBS_hourly.Rdata')
```

2.c) Importing daily data
---
```{r import_data_station_obs, eval=FALSE}
file <- list.files('data/OBS')
length(file)
head(file)

filePM10 <- file[grep('PM10',file)]
filePM10 <- filePM10[nchar(filePM10)<=17]
filePM25 <- file[grep('PM25',file)]
filePM25 <- filePM25[nchar(filePM25)<=17]
fileO3   <- file[grep('O3',file)]
fileO3   <- fileO3[nchar(fileO3)<=15]
fileNO2  <- file[grep('NO2',file)]
fileNO2  <- fileNO2[nchar(fileNO2)<=16]

date_station <- gsub(".csv","",filePM10)
date_station <- gsub("PM10_","",date_station)

### Checking if all files have the same length and order
{
  length(filePM10)
  length(filePM25)
  length(fileO3)
  length(fileNO2)
  filePM10h <- gsub(".csv","",filePM10)
  filePM10h <- gsub("PM10_","",filePM10h)
  filePM25h <- gsub(".csv","",filePM25)
  filePM25h <- gsub("PM25_","",filePM25h)
  fileO3h   <- gsub(".csv","",fileO3)
  fileO3h   <- gsub("O3_","",fileO3h)
  fileNO2h  <- gsub(".csv","",fileNO2)
  fileNO2h  <- gsub("NO2_","",fileNO2h)
  all(filePM10h==filePM25h)
  all(filePM10h==fileO3h)
  all(filePM10h==fileNO2h)
}

PM10l <- vector(mode = "list", length = length(filePM10))
PM25l <- vector(mode = "list", length = length(filePM25))
O3l   <- vector(mode = "list", length = length(fileO3))
NO2l  <- vector(mode = "list", length = length(fileNO2))

filePM10 <- paste('data/OBS/',filePM10,sep="")
filePM25 <- paste('data/OBS/',filePM25,sep="")
fileO3   <- paste('data/OBS/',fileO3,sep="")
fileNO2  <- paste('data/OBS/',fileNO2,sep="")

for(i in 1:length(PM10l)){
      date  <- strptime(date_station[i],format= "%Y%m%d")
      PM10t <- try(read.table(filePM10[i],header = FALSE,sep=" ",na.strings="-999",
      col.names = c("ID","long","lat","PM10"),colClasses = c("factor","numeric","numeric","numeric")))
      PM10t <- cbind(PM10t,rep(date,nrow(PM10t)))
      PM10l[[i]] <- PM10t
      PM25t <- try(read.table(filePM25[i],header = FALSE,sep=" ",na.strings="-999",
      col.names = c("ID","long","lat","PM25"),colClasses = c("factor","numeric","numeric","numeric")))
      PM25t <- cbind(PM25t,rep(date,nrow(PM25t)))
      PM25l[[i]] <- PM25t
      NO2t <- try(read.table(fileNO2[i],header = FALSE,sep=" ",na.strings="-999",
      col.names = c("ID","long","lat","NO2"),colClasses = c("factor","numeric","numeric","numeric")))
      NO2t <- cbind(NO2t,rep(date,nrow(NO2t)))
      NO2l[[i]] <- NO2t
      O3t <- try(read.table(fileO3[i],header = FALSE,sep=" ",na.strings="-999",
      col.names = c("ID","long","lat","O3"),colClasses = c("factor","numeric","numeric","numeric")))
      O3t <- cbind(O3t,rep(date,nrow(O3t)))
      O3l[[i]] <- O3t
      cat(i,'')
}

### Using rbindlist for faster row binding of all loaded dataframes
PM10 <- rbindlist(PM10l)
names(PM10)[5]='date'
PM10$ID2 <- paste(PM10$ID,PM10$date,sep='_') 
PM25 <- rbindlist(PM25l)
names(PM25)[2:3] <- c("long1","lat1")
names(PM25)[5]   <- 'date1' 
PM25$ID2 <- paste(PM25$ID,PM25$date1,sep='_') 
NO2 <- rbindlist(NO2l)
names(NO2)[2:3] <- c("long2","lat2")
names(NO2)[5]   <- 'date2'  
NO2$ID2 <- paste(NO2$ID,NO2$date2,sep='_') 
O3 <- rbindlist(O3l)
names(O3)[2:3] <- c("long3","lat3")
names(O3)[5]   <- 'date3'  
O3$ID2 <- paste(O3$ID,O3$date3,sep='_') 

### Using data.table for faster merge
PM10 <- data.table(PM10)
PM25 <- data.table(PM25)
NO2  <- data.table(NO2)
O3   <- data.table(O3)
res <- merge(PM10[,.(PM10,ID2,date,long,lat)],PM25[,.(PM25,ID2,date1,long1,lat1)],
           by="ID2",all=TRUE)
res <- merge(res,NO2[,.(NO2,ID2,date2,long2,lat2)] ,by="ID2",all=TRUE)
res <- merge(res, O3[,.(O3,ID2,date3,long3,lat3)] ,by="ID2",all=TRUE)

res$ID<-sapply(res$ID2, function(x) strsplit(x,'_')[[1]][1])
res$ID<-as.factor(res$ID)

res=data.frame(res)

### Selecting the filled dates and locations
nodate <- is.na(res$date)
res$date[nodate] <- res$date1[nodate]
nodate <- is.na(res$date)
res$date[nodate] <- res$date2[nodate]
nodate=is.na(res$date)
res$date[nodate] <- res$date3[nodate]

nolong <- is.na(res$long)
res$long[nolong] <- res$long1[nolong]
nolong <- is.na(res$long)
res$long[nolong] <- res$long2[nolong]
nolong <- is.na(res$long)
res$long[nolong] <- res$long3[nolong]

nolat <- is.na(res$lat)
res$lat[nolat] <- res$lat1[nolat]
nolat <- is.na(res$lat)
res$lat[nolat] <- res$lat2[nolat]
nolat <- is.na(res$lat)
res$lat[nolat] <- res$lat3[nolat]

OBS_daily <- res[!is.na(res$date),c('ID','long','lat','date','PM10','PM25','NO2','O3')]
summary(OBS_daily)
save(OBS_daily,file='data/OBS_daily.Rdata')
```

2.d) Importing CHIMERE grid data
---
```{r import_chimere_data,eval=FALSE}
library(ncdf4)

file <- list.files('data/CHM')
length(file)
head(file)

file <- file[grep('chim_',file,fixed=TRUE)]
file <- file[grep('.nc',file)]
file <- paste('data/CHM/',file,sep='')

a <- vector(mode = "list", length = length(file))

for(i in 1:length(file)){
  ncfile <- nc_open(file[i])
  print(file[i])
  lat  <- as.vector(ncvar_get(ncfile,"lat"))
  lon  <- as.vector(ncvar_get(ncfile,"lon"))
  PM10 <- as.vector(ncvar_get(ncfile,"PM10"))
  PM25 <- as.vector(ncvar_get(ncfile,"PM25"))
  NO2  <- as.vector(ncvar_get(ncfile,"NO2"))
  O3   <- as.vector(ncvar_get(ncfile,"O3"))
  time <- strptime(strsplit(file[i],'.',fixed = TRUE)[[1]][2],format = "%Y%m%d")
  dfr  <- data.frame(lon=lon,lat=lat,PM10=PM10,PM25=PM25,NO2=NO2,O3=O3,
                     time=rep(time,each=length(lon)))
  a[[i]]<-dfr
}
CHM <- rbindlist(a)
CHM <- data.frame(CHM)

save(CHM,file='data/CHM.Rdata')
```

2.e) Converting dataframes to spacetime objects
---

```{r coerce_data_to_spacetime_class,eval=TRUE}
rm(list=ls())
library(spacetime)
load('data/OBS_hourly.Rdata'); OBS_hourly <- OBS_horaire
load('data/OBS_daily.Rdata'); OBS_daily <- OBS_jour
load('data/CHM.Rdata')
load('data/stations.Rdata')

### WARNING: daily data contain 2 stations with only 214 recorded days.
### This can cause problems for the construction of a ST object

STIDF_daily <- stConstruct(OBS_daily,space=c('long','lat'),time='date',
                          SpatialObj=SpatialPoints(OBS_daily[,c('long','lat')]))
STFDF_daily    <- as(STIDF_daily, "STFDF") 
STFDF_30d  <- STFDF_daily[,"2014-01-01::2014-01-31"]


### The entire dataset is too large and difficult to handle
### STIDF_obs=stConstruct(OBS_hourly,space=c('long','lat'),time='date')
### with long format and STIDF dimension s and t are wrong.
### For monitoring station data, spacetime::stplot is particularly useful if
### it knows that these are repeated measurements for fixed stations.
### Provided with an STIDF, stplot assumes they are not. 
### Pebesma https://stat.ethz.ch/pipermail/r-sig-geo/2015-July/023101.html)
#STIDF_obs =as(STIDF_obs, "STFDF") ### too long
  
### Considering two subsamples of hourly data of 30 or 7 days

STIDF_hourly_30d <- stConstruct(OBS_hourly[OBS_hourly$date>= "2013-12-31" &
                    OBS_hourly$date <= "2014-02-01",],space=c('long','lat'),time='date')
STIDF_hourly_30d <- as(STIDF_hourly_30d, "STSDF") 

STIDF_hourly_7d <- stConstruct(OBS_hourly[OBS_hourly$date>= "2013-12-31" & 
                    OBS_hourly$date <= "2014-01-08",],space=c('long','lat'),time='date')
STIDF_hourly_7d <- as(STIDF_hourly_7d, "STSDF")

CHM_daily <- stConstruct(CHM[CHM$time>= "2013-12-31" & CHM$time <= "2014-02-01",],
                     space=c('lon','lat'),time='time')
CHM_daily <- as(CHM_daily, "STSDF")


add_ID<-function(STFDF){
  ID  <- paste('ID_',1:(dim(STFDF)[1]),sep='')
  IDs <- rep(ID,dim(STFDF)[2])
  attr(STFDF,'data') <- cbind(IDs,attr(STFDF,'data'))
  return(STFDF)
}
CHM_daily=add_ID(CHM_daily)

### Adding projections
library(sp)
proj4string(CHM_daily)   <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
proj4string(STFDF_daily) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
proj4string(STIDF_hourly_30d) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
proj4string(STIDF_hourly_7d)  <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

save(CHM_daily,file='data/CHM_daily.Rdata')
save(STFDF_daily,file='data/STFDF_daily.Rdata')
save(STIDF_hourly_30d,file='data/STIDF_hourly_30d.Rdata')
save(STIDF_hourly_7d,file='data/STIDF_hourly_7d.Rdata')

```

```{r spacetime_features}
rm(list=ls())
load('data/STFDF_daily.Rdata')


### Understanding STFDF object
dim(STFDF_daily)
str(STFDF_daily)
names(STFDF_daily@data)
dim(STFDF_daily@data)
as.numeric(dim(STFDF_daily)[1]*dim(STFDF_daily)[2])

### Select and subset of STFDF object: STFDF[location_subset,time_subset,var_selection]
### Some examples:
#### 1) Select first location, first day, second variable
STFDF_daily[1,1,2]
#### object of class numeric
class(STFDF_daily[1,1,2])

#### 2) Select first two locations, first day, second variable
STFDF_daily[1:2,1,2]
####  object of class SpatialPointDataFrame
class(STFDF_daily[1:2,1,2])

#### 3) Select first location, first two days, second variable
STFDF_daily[1,1:2,2]
####  object of class xts (time serie)
class(STFDF_daily[1,1:2,2])

#### 4) Select first ten locations, first two days and ID + PM10 variables
STFDF_daily[1:10,1:2,1:2]
####  object of class STFDF
class(STFDF_daily[1:10,1:2,1:2])

#### 5) select 10 days time interval (using date format) and two variables (PM10 and 03)
STFDF_daily[1,"2014-01-01::2014-01-10",c('PM10','O3')]

#### 6) select one station by ID (not as natural as expected), 10 days by date and 2 variables
IDs=levels(STFDF_daily$ID)
STFDF_daily[which(IDs %in% c('FR01011')),"2014-01-01::2014-01-10",c('PM10','O3')]

#### 7) select 2 stations by ID (not as natural as expected), 10 days by date and 2 variables
temp<-STFDF_daily[which(IDs %in% c('FR01011','FR03063')),"2014-01-01::2014-01-10",c('PM10','O3')]

### Get dataframe from spacetime STFDF/STIDF
mydata<-temp@data
summary(mydata)
```

3) Plots of STFDF/STIDF
===

3.a) Static plots 
---
```{r plot_map_statique}
#####################################
rm(list=ls())
load('data/STFDF_daily.Rdata')
library(RColorBrewer)
library(sp)

###  Some examples of spplot
#### 1) Plot all locations, 1 day (by date), 1 variable (O3)
spplot(STFDF_daily[,"2014-01-01",'O3'],col.regions=brewer.pal(10,"Spectral"),cuts=10,
       cex=ifelse(is.na(STFDF_daily[['O3']]),0,1))

#### 2) The same, but point size is proportional to O3 value for that station
cx03 = attr(STFDF_daily[,"2014-01-01",'O3'],'data')/max(attr(STFDF_daily[,"2014-01-01",'O3'],
            'data'),na.rm=T) * 5
spplot(STFDF_daily[,"2014-01-01",'O3'],cex=ifelse(is.na(STFDF_daily[['O3']]),0, 0.4*as.numeric(cx03$O3)),scales=list(draw=TRUE))

#### 3) Plot all locations, 1 day (by date), 2 variables (O3 and NO2)
#### point size is proportional to NO2 value, color to O3 value
cxNO2 = attr(STFDF_daily[,"2014-01-01",'NO2'],'data')/max(attr(STFDF_daily[,"2014-01-01",'NO2']
            ,'data'),na.rm=T) * 5
cxNO2[is.na(cxNO2)]<-0
spplot(STFDF_daily[,"2014-01-01",'O3'],cex=ifelse(is.na(STFDF_daily[['O3']]),0, 0.4*as.numeric(cxNO2$NO2)),scales=list(draw=TRUE))

### Download and add map of countries 
par(mfrow=c(1,1))
library(rworldmap)
library(maptools)
library(ggmap)
world <- getMap(resolution="low")
proj4string(world)
europe.limits <- geocode(c("Spain","Ireland","Germany","Italy"))
plot(world,xlim=range(europe.limits$lon),ylim=range(europe.limits$lat),asp=1,
     main='NO2 colour with O3 circle size \n 2014-01-01')
colours <- brewer.pal(10,"Spectral")
O3      <- attr(STFDF_daily[,"2014-01-01",'O3'],'data')
plot(STFDF_daily[,"2014-01-01",'O3'],col=colours[findInterval(O3$O3,
     quantile(O3$O3,seq(0,1,length.out =10),na.rm=T), all.inside=TRUE)],
     cex= 0.4*as.numeric(cxNO2$NO2),scales=list(draw=TRUE),add=TRUE,pch=19)



### Same as before but using Lambert 93 (L93) projection
splimit <- SpatialPoints(europe.limits)
proj4string(splimit) <- proj4string(world)
splimit_l93 <- spTransform(splimit,CRSobj=CRS("+init=epsg:2154"))
worldL93    <- spTransform(world,CRSobj=CRS("+init=epsg:2154"))
plot(worldL93,xlim=range(splimit_l93$lon),ylim=range(splimit_l93$lat),asp=1,
     main='NO2 colour and O3 circle size \n 2014-01-01 \n proj LAMBERT 93')
STFDF_dayL93 <- STFDF_daily
STFDF_dayL93 <- spTransform(STFDF_dayL93,CRSobj=CRS("+init=epsg:2154"))
plot(STFDF_dayL93[,"2014-01-01",'O3'],col=colours[findInterval(O3$O3, 
    quantile(O3$O3,seq(0,1,length.out =10),na.rm=T), all.inside=TRUE)],
    cex= 0.4*as.numeric(cxNO2$NO2),scales=list(draw=TRUE),add=TRUE,pch=19)


### The same, using spplot (easier and centered around our data)
spplot(STFDF_daily[,"2014-01-01",'O3'],cex=0.4*as.numeric(cxNO2$NO2),
       scales=list(draw=TRUE),pch=19,sp.layout =list(world,first=TRUE))

##################################

### Some examples of stplot
#### 1) Plot all locations, 4 day (by date), 1 variable (O3)
stplot(STFDF_daily[,"2014-01-01::2014-01-04",'O3'],col.regions=
       brewer.pal(10,"Spectral"),cuts=10,cex=rep(1,length(STFDF_daily[,"2014-01-01::2014-01-04",'O3'])))

#### 2) Plot GRID data, 1 day (by date), 1 variable (O3)
load('data/CHM_daily.Rdata')
chim_O3_1d<-CHM_daily[,"2014-01-01",'O3']
spplot(chim_O3_1d,col.regions=brewer.pal(10,"Spectral"),cuts=10)

#### 3) with grid data as real gridded object
spplot(chim_O3_1d,col.regions=brewer.pal(10,"Spectral"),cuts=10)

#### 4) with a STFDF we use gridded(STFDF@sp), 7 days
chim_O3_7j <- CHM_daily[,"2014-01-01::2014-01-07",'O3']
stplot(chim_O3_7j,col.regions=brewer.pal(10,"Spectral"),
       cuts=10,cex=rep(1,length(chim_O3_7j[['O3']])),
       sp.layout=list(world,first=FALSE),number=7)
#### 5) add and plot new variable (ratio between NO2 and PM10) 
NO2   <- attr(STFDF_daily[,"2014-01-01",'NO2'],'data')
PM10  <- attr(STFDF_daily[,"2014-01-01",'PM10'],'data')
ratio <- NO2$NO2/PM10$PM10
ratio[is.na(ratio)]<-0
spplot(STFDF_daily[,"2014-01-01",'PM10'],cex=ifelse(is.na(CHM_daily[['PM10']]),0,0.4*as.numeric(ratio)),
       scales=list(draw=TRUE))
```

3.b) Static plots using klm (plotklm)
---
```{r plotklm, eval=FALSE}

### May require an update of Google Earth
library(plotKML)
O3df    <- attr(STFDF_daily[,"2014-01-01","O3"],'data')
NO2df   <- attr(STFDF_daily[,"2014-01-01","NO2"],'data')
PM10df  <- attr(STFDF_daily[,"2014-01-01","PM10"],'data')
colours <- SAGA_pal[[1]]

### plotKML with colours for O3
plotKML(STFDF_daily[,"2014-01-01",'O3'],colours=colours[findInterval(O3df$O3, 
        quantile(O3df$O3,seq(0,1,length.out =10),na.rm=T), all.inside=TRUE)])

### same for NO2 with colour_scale
plotKML(STFDF_daily[,"2014-01-01",'NO2'])

```

3.c) Time series plots
---
```{r time_by_station, eval=TRUE}
rm(list=ls())
library(tseries)
library(ggfortify)
library(ggplot2)
load('data/STFDF_daily.Rdata')
load('data/OBS_daily.Rdata'); OBS_daily <- OBS_jour

### One time serie plot can represent the data corresponding to one or
### multiple locations as well as one or multiple variables. 

#### 1) single location, multiple variables
IDs <- levels(OBS_daily$ID)
ts  <- STFDF_daily[which(IDs %in% c('FR01011')),,c("PM10","NO2","O3")]
colours <- brewer.pal(3,"Set1")
irr     <- irts(index(ts),as.matrix(ts[,1:3]))

p1 <- autoplot(irr, facets = FALSE)
p1 <- p1+ggtitle('FR01011')
p1

#### 2) Two locations, multiple variables
ts2  <- STFDF_daily[which(IDs %in% c('FR03063')),,c("PM10","NO2","O3")]
irr2 <- irts(index(ts2),as.matrix(ts2[,1:3]))
p2   <- autoplot(irr2, facets = FALSE)
p2   <- p2+ggtitle('FR03063')
new("ggmultiplot",plots=list(p1,p2),ncol=1)


#### 3)  Multiple locations >3, single variable
##### PM10
init <- STFDF_daily[which(IDs==IDs[1]),,c("PM10")]
   a <- matrix(0,ncol=10,nrow=length(init[,1]))
for(i in 1:10 ){
  ID    <- IDs[i]
  a[,i] <- STFDF_daily[which(IDs==ID),,c("PM10")][,1]
}
colnames(a) <- IDs[1:10]
irr <- irts(index(init),a)
p3  <- autoplot(irr, facets = FALSE)+ggtitle('PM10')
p3

####  4) Multiple locations >3, multiples variables (grouped by variable)

#####  O3
init <- STFDF_daily[which(IDs==IDs[1]),,c("O3")]
a    <- matrix(0,ncol=10,nrow=length(init[,1]))
for(i in 1:10 ){
  ID    <- IDs[i]
  a[,i] <- STFDF_daily[which(IDs==ID),,c("O3")][,1]
}
colnames(a) <- IDs[1:10]
irr <- irts(index(init),a)
p4  <- autoplot(irr, facets = FALSE)+ggtitle('O3')

##### NO2
init <- STFDF_daily[which(IDs==IDs[1]),,c("NO2")]
a    <- matrix(0,ncol=10,nrow=length(init[,1]))
for(i in 1:10 ){
  ID    <- IDs[i]
  a[,i] <- STFDF_daily[which(IDs==ID),,c("NO2")][,1]
}
colnames(a) <- IDs[1:10]
irr <- irts(index(init),a)
p5  <- autoplot(irr, facets = FALSE)+ggtitle('NO2')

new("ggmultiplot",plots=list(p3,p4,p5),ncol=1)

####  5) As before, but grouping by location
IDs <- levels(OBS_daily$ID)
ts  <- STFDF_daily[which(IDs %in% c('FR03080')),"2014-01-01::2014-04-01",c("PM10","NO2","O3")]
colours <- brewer.pal(3,"Set1")
irr <- irts(index(ts),as.matrix(ts[,1:3]))
p1  <- autoplot(irr, facets = FALSE)
p1  <- p1+ggtitle('FR03080')

ts2  <- STFDF_daily[which(IDs %in% c('FR08209')),"2014-01-01::2015-04-01",c("PM10","NO2","O3")]
irr2 <- irts(index(ts2),as.matrix(ts2[,1:3]))
p2   <- autoplot(irr2, facets = FALSE)
p2   <- p2+ggtitle('FR08209')

ts2  <- STFDF_daily[which(IDs %in% c('FR03084')),"2014-01-01::2014-04-01",c("PM10","NO2","O3")]
irr2 <- irts(index(ts2),as.matrix(ts2[,1:3]))
p3   <- autoplot(irr2, facets = FALSE)
p3   <- p3+ggtitle('FR03084')

new("ggmultiplot",plots=list(p1,p2,p3),ncol=1)

```

3.d) Animated plot using stplot() -- RESULTS NOT DISPLAYED --
---
```{r animated_stplot, eval=FALSE}
###  Animated plot -- NOT DISPLAYED --
### 10 days, two variables, with sp.layout, but refresh too slow
cxNO2_10d <- attr(STFDF_daily[,"2014-01-01::2014-01-10",'NO2'],'data')/
                  max(attr(STFDF_daily[,"2014-01-01::2014-01-10",'NO2'],'data'),na.rm=T) * 5
stplot(STFDF_daily[,"2014-01-01::2014-01-10",'O3'],cex=0.4* as.numeric(cxNO2_10d$NO2),
       col.regions=brewer.pal(10,"Spectral"),cuts=10, 
       animate=0.2,do.repeat=F,sp.layout =list(world,first=TRUE))

### 10 days, two columns, without sp.layout: better refresh
cxNO2_4d <- attr(STFDF_daily[,"2014-01-01::2014-01-04",'NO2'],'data')
                 /max(attr(STFDF_daily[,"2014-01-01::2014-01-04",'NO2'],'data'),na.rm=T) * 5
stplot(STFDF_daily[,"2014-01-01::2014-01-04",'O3'],cex=0.4* as.numeric(cxNO2_4d$NO2),
       col.regions=brewer.pal(10,"Spectral"),cuts=10, animate=0.2,do.repeat=F)

### 10 days, one variable for color one for size, with HTML output
### The animation is saved in mpeg format with high quality 
### (requires installation of command line tool imagemagick)
stplot(STFDF_daily[,"2014-01-01::2014-01-10",'O3'],cex=0.4* as.numeric(cxNO2_10d$NO2),
       col.regions=brewer.pal(10,"Spectral"),cuts=10,sp.layout =world, animate=0.2,do.repeat=F)

for(i in c('01','02','03','04','05','06','07','08','09','10')){
 par(mar = c(3, 3, 2, 0.5), mgp = c(2, 0.5, 0), tcl = -0.3, cex.axis = 0.8 
     ,cex.lab = 0.8, cex.main = 1)
 cxNO2_1d <- attr(STFDF_daily[,paste("2014-01-",i,sep=''),'NO2'],'data')
                 /max(attr(STFDF_daily[,paste("2014-01-",i,sep=''),'NO2'],'data'),na.rm=T) * 5
 jpeg(paste('images/se',i,'.jpg',sep=''),width=5, height=5, units="in",res=600)
 cat(paste('se',i,'.jpg',sep=''),'\n')
 print(spplot(STFDF_daily[,paste("2014-01-",i,sep=''),'O3'],cex=0.4*as.numeric(cxNO2_1d$NO2),
             col.regions=brewer.pal(10,"Spectral"),cuts=10,sp.layout =world))
dev.off()
}   
system('convert  -delay 25 -quality 100  images/*.jpg  animation10days.mpeg')
system('open animation10days.mpeg')


### Animated plotKML with colours NO2 
library(plotKML)
ID <- STFDF_daily[,"2014-01-01",'ID']
plotKML(STFDF_daily[,"2014-01-01",'ID'],colour_scale=SAGA_pal[[1]],
labels='')

```

3.e) Plot of spatio-temporal dependencies
---
```{r spatio_temporal_dep}
rm(list=ls())
library(gstat)
load('data/STFDF_daily.Rdata')


### Temporal auto-correlation 
### choose station with no NA for O3 for the period 2014-01-01::2014-01-30
rn <- row.names(STFDF_daily@sp)[c(26:28,36)] 

## Multiple stations, one variable 
## lags are expressed in secondes, not days
par(mfrow=c(2,2))
for(i in rn) acf(STFDF_daily[i,"2014-01-01::2014-01-30","O3"]$O3,main=i)

### Cross correlation

## Cross correlation between multiple stations and one variable
par(mfrow=c(1,1))

acf(as(STFDF_daily[rn,"2014-01-01::2014-01-30","O3"],"xts"))

## Cross correlation between multiple stations and three variables

acf(STFDF_daily[5,"2014-01-01:::2014-01-30",c("O3","NO2","PM10")][,1:3])



### Spatial auto-correlation

## Spatial auto-correlation for fixed date and one variable

library(gstat)
v <- variogram(O3~lat+long,STFDF_daily[!is.na(STFDF_daily[,"2014-01-01","O3"]$O3),"2014-01-01","O3"])
plot(v,main="O3")


## Spatial auto-correlation for fixed date and multiple variables
g <- gstat(NULL, "O3", O3 ~ 1, data = STFDF_daily[!is.na(STFDF_daily[,"2014-01-01","O3"]$O3)
           ,"2014-01-01","O3"])
g <- gstat(g, "NO2", NO2 ~ 1, data = STFDF_daily[!is.na(STFDF_daily[,"2014-01-01","NO2"]$NO2)
            ,"2014-01-01","NO2"])
g <- gstat(g, "PM10", PM10 ~ 1, data = STFDF_daily[!is.na(STFDF_daily[,"2014-01-01","PM10"]$PM10)
            ,"2014-01-01","PM10"])
g <- gstat(g, "PM25", PM25 ~ 1, data = STFDF_daily[!is.na(STFDF_daily[,"2014-01-01","PM25"]$PM25)
            ,"2014-01-01","PM25"])

vm <- variogram(g)
plot(vm)

## mean variogram on multiple dates

lst <- lapply(1:5,function(i){x=STFDF_daily[,i,"O3"]; x$ti=i ; rownames(x@coords)=NULL;x})
pts <- do.call(rbind,lst)

vv <- variogram(O3~ti,pts[!is.na(pts$O3),])
plot(vv, main="O3")

## Spatio-temporal auto-correlation

vst <- variogram(O3~lat+long,STFDF_daily[!is.na(STFDF_daily[,"2014-01-01::2014-01-31","O3"]$O3),
                                     "2014-01-01::2014-01-31","O3"])
plot(vst,main="O3")

vvst <- variogram(O3~1,STFDF_daily[!is.na(STFDF_daily[,"2014-01-01::2014-12-31","O3"]$O3),
                               "2014-01-01::2014-12-31","O3"],
                  width=5,cutoff=200,tlags=0:5)
plot(vvst)
plot(vvst,map=FALSE)
plot(vvst,wireframe=T)

```

3.f) Interactive plots with Shiny
---
```{r shiny_leaflet,eval=FALSE}
setwd(paste(rep_base,'/data',sep=""))
library(raster)
library(tseries)
library(ggfortify)
library(gstat)
library(ggplot2)
library(shiny)
library(spacetime)
library(leaflet)
library(RColorBrewer)
library(sp)
library(rgeos)
library(rworldmap)
library(maptools)
library(ggmap)

load('data/STFDF_daily.Rdata')
ST2 <- STFDF_daily[,"2014-01-01::2014-01-30",c("ID","NO2","O3","PM10")]

### This Shiny application creates an interface that the user can navigate and that
### contains different objects,such as, maps, plots, summaries, etc, used to explore our data set. 

### The user has the possibility to interact with the inteface by choosing, e.g., 
### the station to consider, the time interval or the variables to explore, and the 
### resulting outputs reacts according to the choice 

### The basic ingredients to construct a Shiny application are: 
### 1) the UI, where the layout of the interface is defined,
### 2) the SERVER, which contains the "instructions" on how the UI should work.

### In the UI (user interface) we control:
### 1) the VISUAL ASPECT of the application (color, disposition of panels, etc.)
### 2) the INPUT variables: the objects that the user can handle and modify (db, dates, etc.)
### 3) the OUTPUT: what the app will show to the user (plots, maps, graphs, etc.)

### In the SERVER we:
### 1) explain what to do with input variables and how to generate the output
### 2) control the reactivity of the app (more complicated, but highly efficient)


runApp(list(
 ui = fluidPage(
### ui will contain the UI of the app.
### UI is of type "FluidPage": layout consisting of rows which in turn include columns.
   sidebarLayout(
### Create a layout with a sidebar and main area.      
     sidebarPanel(
### Create a sidebar panel containing input controls that can in turn be passed to sidebarLayout
### sidebarPanel() contains the objects that the user can interact with, i.e. the inputs.
### To each input correspond a unique name, that can be called in the SERVER part by input$name
### Every time that an input is modified by the user, all functions in SERVER that make use of
### that input will react to the change (unless differently specified).
### Objects of different class can become inputs, such as numerical values, dataframes, files, etc.
### Different control widgets, listed below, are used to interact with the user:
###actionButton	-> Action Button
###checkboxGroupInput	-> A group of check boxes
###checkboxInput	-> A single check box
###dateInput	-> A calendar to aid date selection
###dateRangeInput	-> A pair of calendars for selecting a date range
###fileInput	-> A file upload control wizard
###numericInput	-> A field to enter numbers
###radioButtons	-> A set of radio buttons
###selectInput	-> A box with choices to select from
###sliderInput	-> A slider bar
###submitButton	-> A submit button
###textInput	-> A field to enter text
       
         selectInput('dataset', 'Choose Dataset', ls()[sapply(mget(ls(), .GlobalEnv), 
                        function(x) 'STFDF' %in% is(x))]),
### Create a select list with all active STFDF objects. Choose one to start working.
         selectInput('vars', 'Variables', "",multiple=TRUE),
### Choose the variables of the chosen STFDF.
         dateRangeInput(inputId="tim_int",label="Choose date interval"),
### Set the time interval using two calendars.
         radioButtons("choiceID","Choice of locations using...",c(IDs = "byIDS",map = "bymap"),
                      selected=NULL),conditionalPanel(condition = "input.choiceID == 'byIDS'",
                          selectInput('IDs', 'Which station', "",multiple=TRUE)),
         conditionalPanel(condition = "input.choiceID == 'bymap'", leafletOutput("mymap")),
### Create a radio button to select the type of choice for the station:
### If Ids, it will visualize the list of all IDs amongst which the user can choose
### If map, it will show a leaflet map with marked dots in correspondence of the stations
### radioButtons() works in combination with conditionalPanel(), for the conditional choices
         actionButton("buttonUpdate","Update")
### creates a button that updates the choices made by the user and start the computations
### for generating the output. This functionality can be avoided in case that output generation
### is very fast, allowing hence an interactive and immediate control of the user on the outputs.
### We prefered to add it because our outputs are sometimes slow to compute.
### It is a way of controlling the reactivity of the application.
         ),
     mainPanel(tabsetPanel(
       tabPanel("Summary",textOutput("choices")),
       tabPanel("Cartography",plotOutput("plot_carto"),plotOutput("plot_carto2")),
       tabPanel("Plots",plotOutput("plot1")),
       tabPanel("Dependence Structure",plotOutput("dep1"),plotOutput("dep2"))
### mainPanel is the place where the outputs will be visualized. We divided in 4 tabs (tabPanel())
### Outputs are generated by the functions in the SERVER part.
### To each output correspond a unique name, that can be called in the SERVER part by output$name
### Outputs can take different forms, such as text, images, plots, maps, etc.
### htmlOutput	
### imageOutput	
### plotOutput	
### tableOutput	
### textOutput	
### uiOutput	
### verbatimTextOutput
         )
       ))),
 server = function(session,input, output){
### the SERVER function contains all the instruction on what to do once the inputs have been chosen   

   triggerd.by <<- NULL
   data1=reactive({get(input$dataset)})
### The chosen dataset "input$dataset" is set as data1 (class STFDF)
   observeEvent(input$dataset,{
     df1<-attr(data1(),'data')
     dt1<-attr(data1(),'time')
     updateSelectInput(session, "vars",choices = names(df1))
     updateDateRangeInput(session,"tim_int",start =min(index(dt1)),end=max(index(dt1)),
                          min=min(index(dt1)),max=max(index(dt1)))
     updateSelectInput(session, "IDs",choices = levels(df1$ID))
### Once the dataset has been chosen, the widgets of "Which variable","time interval" and " available 
### stations" are automatically updated to reflect the choice of the database
   })

   observeEvent({"input.choiceID == 'bymap'"},{
     ds1=as(data1(),"Spatial")
     output$mymap <- renderLeaflet({leaflet(ds1) %>% 
         addTiles(urlTemplate = "http://{s}.tile.opentopomap.org/{z}/{x}/{y}.png") %>%
         setView(lng =mean(coordinates(ds1)[,1]),
                 lat = mean(coordinates(ds1)[,2]), zoom = 8)%>%
         addCircles(color="red")})
### If the choice of the stations is made "by map", a leaflet map will appear in the sidebar
   })

   observeEvent(input$buttonUpdate,{
### All function contained in observeEvent(input$buttonUpdate{}) will execute only after the 
### "Update" button has been clicked. It is an easy way of controlling the reactivity.
### Obviously, the outputs depends on the chosen inputs:
     nvars <- length(input$vars)
     time_interval <- input$tim_int
     int_in <- input$tim_int[1]
     int_en <- input$tim_int[2]
     int_date<-paste(int_in,int_en,sep="::")
     ndays <- as.numeric(int_en-int_in)


     if(input$choiceID == 'bymap'){
       ds1=as(data1()[,1:2],"Spatial")
       a<-matrix(unlist(input$mymap_bounds)[c(4,3,2,1)],2,2)
       b_poly <- as(extent(a), "SpatialPolygons")
       proj4string(b_poly)<-proj4string(ds1)
       inters_sp<-!is.na(over(ds1,b_poly))
       npoints <- sum((inters_sp==1)*1)
       d<-data1()[inters_sp,int_date,input$vars]
       d2<-data1()[inters_sp,int_date,c("ID",input$vars)]
     }
     if(input$choiceID == 'byIDS'){
       d<-data1()[attr(data1(),"data")$ID %in% input$IDs,int_date,input$vars]
       d2<-data1()[attr(data1(),"data")$ID %in% input$IDs,int_date,c("ID",input$vars)]
       npoints <- length(input$IDs)
       chosenID <- attr(data1(),"data")$ID %in% input$IDs
     }

     output$choices = renderText({
       paste("You have chosen ",nvars," variable(s).\n",
             "The time interval is of ",as.numeric(int_en-int_in)," days.\n",
             "You are considering ",npoints," points/stations\n",sep="")
     })
### The summary of the chosen subset is print using "renderText()".

     world<-getMap(resolution="low")

### output$plot_carto and output$plot_carto2 are cartographies plotted using
### the functions presented in Sections 3.a and 3.d. They are plotted using "renderPlot()".
### output$plot1 is a static plot as presented in Section 3.c for time series.
### output$dep1 and output$dep2 are dependency plots as explained in Section 3.e.

###Map: n stations, 1 day, 1 variable
     if (npoints>1 && ndays==0 && nvars==1 ){
       output$plot_carto <- renderPlot({spplot(d,col.regions=brewer.pal(10,"Spectral"),
                                               cuts=10,sp.layout=world)})
     }

###Map: n stations, 1 day, 2 variables
     if (npoints>1 && ndays==0 && nvars==2 ){
       output$plot_carto <- renderPlot({
         cx = attr(d,"data")[,input$vars[2]]/max(attr(d,"data")[,input$vars[2]],na.rm=T) * 5
         spplot(d[,input$vars[1]],cex=0.4*as.numeric(cx),scales=list(draw=TRUE),
                pch=19,sp.layout=world)
       })
       output$plot_carto2 <- renderPlot({
         cx = attr(d,"data")[,input$vars[1]]/max(attr(d,"data")[,input$vars[1]],na.rm=T) * 5
         spplot(d[,input$vars[2]],cex=0.4*as.numeric(cx),scales=list(draw=TRUE),
                pch=19,sp.layout=world)
       })
     }

###Animated Map: n stations, n days, 1 variable
###Time seriesPlot: n stations, n days, 1 variable
     if (npoints>1 && ndays>0 && nvars==1 ){
       output$plot_carto <- renderPlot({stplot(d,col.regions=brewer.pal(10,"Spectral"),
                                               cuts=10,sp.layout=world) })
       output$plot_carto2 <- renderPlot({stplot(d,col.regions=brewer.pal(10,"Spectral"),
                                      cuts=10,sp.layout=world, animate=0.2,do.repeat=F)})
       output$plot1 <- renderPlot({
         tim<-attr(d2,"time")
         dt<-attr(d2,"data")
         n<-dim(dt)[1]/length(tim)
         a<-split(dt,dt$ID)
         m<-matrix(0,ncol=n,nrow=length(tim))
         for (i in 1:n){
           if (!all(is.na(a[[i]][,2]))) m[,i]<-as.numeric(a[[i]][,2])
         }
         colnames(m)<-levels(d2$ID)
         ir<-irts(index(tim),m)
         p=autoplot(ir, facets = FALSE)+ggtitle(as.character(input$vars))
         p
       })

     }

###Animated Map: n stations, 1 day, 2 vars
###Time series Plot: n stations, 1 day, 2 vars
     if (npoints>1 && ndays>0 && nvars==2 ){
       output$plot_carto <- renderPlot({
         cx = attr(d,'data')[,input$vars[2]]/max(attr(d,'data')[,input$vars[2]],na.rm=T) * 5
         stplot(d[,,input$vars[1]],cex=0.4* as.numeric(cx),col.regions=brewer.pal(10,"Spectral"),
                cuts=10,sp.layout=world, animate=0.2,do.repeat=F)
       })

       output$plot_carto2 <- renderPlot({
         cx = attr(d,'data')[,input$vars[1]]/max(attr(d,'data')[,input$vars[1]],na.rm=T) * 5
         stplot(d[,,input$vars[2]],cex=0.4* as.numeric(cx),col.regions=brewer.pal(10,"Spectral"),
                cuts=10,sp.layout=world, animate=0.2,do.repeat=F)
       })

     }

### 1 station, multiple variables
     if(npoints==1 && ndays>0 && nvars>=2 ){
       output$plot1 <- renderPlot({
         colours=brewer.pal(3,"Set1")
         irr<-irts(index(d),as.matrix(d[,1:nvars]))
         p1<-autoplot(irr, facets = FALSE)
         p1})
     }

     if(npoints>1 && ndays>0 && nvars>=2 ){
       output$plot1 <- renderPlot({
         tim<-attr(d2,"time")
         dt<-attr(d2,"data")
         n<-dim(dt)[1]/length(tim)
         k<-length(input$vars)
         a<-split(dt,dt$ID)
         p=vector(mode="list",length = k)
         for (j in 1:k){
           m<-matrix(0,ncol=n,nrow=length(tim))
           for (i in 1:n){
             if (!all(is.na(a[[i]][,j+1]))) m[,i]<-as.numeric(a[[i]][,j+1])}
           colnames(m)<-levels(dt$ID)
           ir<-irts(index(tim),m)
           p[[j]]=autoplot(ir, facets = FALSE)+ggtitle(as.character(input$vars[j]))
         }
         new("ggmultiplot",plots=p,ncol=1)})
     }

     if (npoints>1 && ndays>3 && nvars==1){
       output$dep1 <- renderPlot({
         acf(as(d,"xts"),na.action = na.pass)  
       })
     }

     if (npoints==1 && ndays>3 && nvars==1){
       output$dep1 <- renderPlot({
         acf(d[,input$vars],na.action = na.pass)  
       })
     }

     if (ndays>3 && nvars>1 && npoints==1){
       output$dep1 <- renderPlot({
         acf(as(d,"xts"),na.action = na.pass)  
       })
     }

   })
 }))
```

4.) Estimating a space-time covariance function
===
4.1) Data preparation
---

```{r data preparation,eval=TRUE,echo=TRUE}
rm(list=ls())
library(rgdal)
library(sp)
library(spacetime)
######### reload saved data objects (if necessary)
load('data/OBS_daily.Rdata'); OBS_daily <- OBS_jour
load('data/CHM.Rdata')
load('data/stations.Rdata')
######### extract French background stations
IDs_bg   <- stations$station_european_code[stations$type_of_station=="Background"]
IDs_fr   <- stations$station_european_code[stations$country_name=="France"]
IDs_bgfr <- intersect(IDs_bg,IDs_fr)
idx_stations <- which(stations$station_european_code %in% IDs_bgfr)
coord    <- cbind(stations$station_longitude_deg,stations$station_latitude_deg)[idx_stations,]
dim(coord) #around 700 French background stations
######### select data for fitting according to background stations and time window
######### we choose a three month period (but not longer) to keep the data size reasonably and to avoid problems with handling NA data for CompRandFld
tstart   <- "2014-01-01"
tend     <- "2014-03-31"
OBS_daily <- OBS_daily[OBS_daily$ID %in% IDs_bgfr,]
OBS_sel  <- stConstruct(OBS_daily,space=c('long','lat'),time='date',
                        SpatialObj=SpatialPoints(OBS_daily[,c('long','lat')]))
OBS_sel  <- as(OBS_sel,"STFDF") 
proj4string(OBS_sel)="+init=epsg:4326"
OBS_sel  <- OBS_sel[,paste0(tstart,"::",tend)]
######### transform to Lambert 93 coordinates system
OBS_sel@sp <- spTransform(OBS_sel@sp, CRS("+init=epsg:2154")) 

######### interpolate gridded CHIMERE values to selected space-time observation points
library(fields)
#first, get CHIMERE predictions for selected time period
CHM_sel  <- CHM[as.character(CHM$time)>=tstart&as.character(CHM$time)<=tend,]
grid_CHM <- unique(cbind(CHM_sel$lon,CHM_sel$lat))
CHM_mat  <- matrix(CHM_sel$PM10, nrow=nrow(grid_CHM))
dim(CHM_mat)
lon_grid <- sort(unique(CHM_sel$lon))
lat_grid <- sort(unique(CHM_sel$lat))
idx_sites_obs <- match(unique(OBS_sel@data$ID),IDs_bgfr)
#function to interpolate gridded CHIMERE data (for a fixed day i):
fun2interp <- function(i){
  tmp <- list(x=lon_grid,y=lat_grid,
              z=matrix(CHM_mat[,i],length(lon_grid),length(lat_grid)))
  interp.surface(tmp, coord[idx_sites_obs,])
}
CHM_sites <- sapply(1:ncol(CHM_mat),fun2interp)
######### save as variable CHM into the STFDF object:
OBS_sel@data$CHM <- as.numeric(CHM_sites)

######### save PM10 residuals with respect to CHM as variable PM10res into the STFDF object:
OBS_sel@data$PM10res <- OBS_sel@data$PM10-as.numeric(CHM_sites)
```

4.2) Estimation with *CompRandFld*
---

```{r Estimation with CompRandFld,eval=TRUE,echo=TRUE, tidy=TRUE}
library(CompRandFld)
######### data must follow the format t x d (with t observation times and d sites; currently no missing data allowed)
######### to obtain a set of stations and time points without missing data, we proceed through an iterative removal of stations and time points with many missing data
names(OBS_sel@data)
data_mat <- matrix(OBS_sel@data$PM10res,ncol=dim(OBS_sel)[1],nrow=dim(OBS_sel)[2],byrow=T)
dim(data_mat)
mean(is.na(data_mat)) #around 44 per cent of missing data
######### define a list data2keep, to be modified subsequently, to retain a subset of non NA data
data2keep=list(data_mat=data_mat,times=1:nrow(data_mat),stations=1:ncol(data_mat))
######### define a function cleandata that removes all rows / columns with NA proportion above propNA
cleandata=function(data2keep,propNA){
  #remove sites and days with more a proportion of more than propNA NA values
  idx2keep=apply(is.na(data2keep$data_mat),1,mean) < propNA
  data2keep$times=data2keep$times[idx2keep]
  data2keep$data_mat  <- data2keep$data_mat[idx2keep,]
  idx2keep=apply(is.na(data2keep$data_mat), 2, mean) < propNA
  data2keep$stations=data2keep$stations[idx2keep]
  data2keep$data_mat  <- data2keep$data_mat[,idx2keep]
  data2keep
}
data2keep=cleandata(data2keep,.9)
data2keep=cleandata(data2keep,.2)
data2keep=cleandata(data2keep,.1)
data2keep=cleandata(data2keep,.05)
data2keep=cleandata(data2keep,.025)
data2keep=cleandata(data2keep,.01)
data2keep=cleandata(data2keep,.10^{-10})
data_mat=data2keep$data_mat
dim(data_mat) #we keep 70 stations and 103 days
######### now update the coordinate and stations object for fitting
coord_fit <- (OBS_sel@sp@coords/1000)[data2keep$stations,] #Lambert 93 coordinates in  km
stations_fit <- stations[idx_stations[data2keep$stations],]
plot(coord_fit,asp=1,xlab="(km)",ylab="(km)",pch=4) 
times_fit <- (1:dim(OBS_sel)[2])[data2keep$times]
######### save selected data for later use
save(stations_fit, file="data/stations_fit.Rdata")
save(times_fit,file="data/times_fit.Rdata")
save(data_mat,file="data/data_mat.Rdata")

######### we assume that the PM10res data are stationary over space and time
######### we start by estimating the empirical mean and the empirical variance from data
mean_est <- mean(data_mat)
mean_est
var_est  <- var(as.numeric(data_mat))
var_est
sqrt(var_est) #estimated standard deviation
######### now calculate an empirical space-time variogram (see example(EVariogram) for examples of plotting the result)
vgm_emp  <- EVariogram(data=data_mat, coordx=coord_fit, coordt=times_fit, cloud=F, maxdist=1000, maxtime=12)
######### full likelihood estimation would be really hard and slow:
prod(dim(data_mat))
######### now we define some parametric models to estimate (with initial values and some of the parameters fixed)
#space-time separable exponential model without nugget
cormod1 <- "exp_exp"
fixed1  <- list(mean=mean_est,nugget=0,sill=var_est)
start1  <- list(scale_s=200,scale_t=2)
#space-time separable exponential model with nugget
cormod2 <- "exp_exp"
fixed2  <- list(mean=mean_est,sill=var_est)
start2  <- list(scale_s=200,scale_t=2,nugget=0)
#Gneiting model with powers fixed to 1, without nugget
cormod3 <- "gneiting"
fixed3  <- list(sill=var_est,mean=mean_est,nugget=0,power_s=1,power_t=1)
start3  <- list(scale_s=200,scale_t=2,sep=.5)
#Gneiting model with powers fixed to 1, with nugget
cormod4 <- "gneiting"
fixed4  <- list(sill=var_est,mean=mean_est,power_s=1,power_t=1)
start4  <- list(scale_s=200,scale_t=2,sep=.5)
#Gneiting model with powers fixed to .5, without nugget
cormod5 <- "gneiting"
fixed5  <- list(sill=var_est,mean=mean_est,nugget=0,power_s=.5,power_t=.5)
start5  <- list(scale_s=200,scale_t=2,sep=.5)
#Gneiting model with powers fixed to .5, with nugget
cormod6 <- "gneiting"
fixed6  <- list(sill=var_est,mean=mean_est,power_s=.5,power_t=.5)
start6  <- list(scale_s=200,scale_t=2,sep=.5)
#Gneiting model where powers are estimated, without nugget
#difficult to propose good starting values, but we will take estimated parameters from preceding models 
cormod7 <- "gneiting"
fixed7  <- list(sill=var_est,mean=mean_est,nugget=0)
#Gneiting model where powers are estimated, with nugget
#difficult to propose good starting values, but we will take estimated parameters from preceding models 
cormod8 <- "gneiting"
fixed8  <- list(sill=var_est,mean=mean_est)

######### pairwise likelihood fitting of the first model with space and time cut-off distances (should take less than a minute)
fit1PL <- FitComposite(data=data_mat,coordx=coord_fit,coordt=times_fit,maxdist=300,maxtime=4,corrmodel=cormod1,likelihood="Marginal",type="Pairwise",fixed=fixed1,start=start1,varest=T)
fit1PL
######### visual comparison of empirical and fitted variogram
library(scatterplot3d)
Covariogram(fit1PL,vario=vgm_emp,show.vario=T,pch=20)
######### weighted least squares fitting
fit1WLS <- WLeastSquare(data=data_mat,coordx=coord_fit,coordt=times_fit,maxdist=300,maxtime=4,corrmodel=cormod1,fixed=fixed1,start=start1,weighted=T)
fit1WLS
######### comparison of relative difference in PL and WLS fits: less than 25%
(fit1PL$param-fit1WLS$param)/pmax(abs(fit1PL$param),abs(fit1WLS$param))

######### now we fit all models with weighted least squares or pairwise likelihood
######### smaller values for maxtime, maxdist to go faster (100,2, instead of 300,4)
######### Uncomment following lines for runningthe fitting procedures
# fitsWLS <- fitsPL <- vector("list", 8)
# for(i in 1:8){
#   cat("fitting model number ", i, "\n")
#   #need to fix sensible initial values for Gneiting models (7,8) with many parameters to estimate
#   if(i==7){
#     startWLS=c(as.list(fitsWLS[[5]]$param),power_s=.5,power_t=.5)
#   }else if(i==8){
#     startWLS=c(as.list(fitsWLS[[6]]$param),power_s=.5,power_t=.5)
#   }else{
#     startWLS=get(paste0("start",i))
#   }
#   fitsWLS[[i]] <- WLeastSquare(data=data_mat,coordx=coord_fit,coordt=times_fit,
#                                maxdist=300,maxtime=4,corrmodel=get(paste0("cormod",i)),
#                                fixed=get(paste0("fixed",i)),start=startWLS,weighted=T)
#   if(i==7){
#     startPL=c(as.list(fitsPL[[5]]$param),power_s=.5,power_t=.5)
#   }else if(i==8){
#     startPL=c(as.list(fitsPL[[6]]$param),power_s=.5,power_t=.5)
#   }else{
#     startPL=as.list(fitsWLS[[i]]$param)
#   }
#   fitsPL[[i]] <- FitComposite(data=data_mat,coordx=coord_fit,coordt=times_fit,
#                               maxdist=300,maxtime=4,corrmodel=get(paste0("cormod",i)),
#                               likelihood="Marginal",type="Pairwise",fixed=get(paste0("fixed",i)),
#                               start=startPL,varest=T) 
# }
# ######### save lists of fit objects
# save(fitsWLS,file="data/fitsWLS.Rdata")
# save(fitsPL,file="data/fitsPL.Rdata")
load(file="data/fitsWLS.Rdata")
load(file="data/fitsPL.Rdata")

######### which pairwise likelihood model is best in terms of the CLIC?
unlist(sapply(fitsPL,getElement,"clic"))
which.min(unlist(sapply(fitsPL,getElement,"clic")))

######### test some hypotheses for embedded Gneiting models fitted with PL
fit8 <- fitsPL[[8]]
fit7 <- fitsPL[[7]]
fit5 <- fitsPL[[5]]
HypoTest(fit8,fit7,fit5,statistic="Wald")
#result: submodel 5 is rejected, while the absence of a nugget in model 7 (in comparison to model 8 ) cannot be rejected
```

4.3) Estimation with *gstat*
---

```{r Estimation with gstat,eval=TRUE,echo=TRUE,tidy=FALSE}
library(gstat)
######### define a function variot for plotting the marginal temporal 
######### variogram at distance 0
variot <- function(model,vario=vario){
  tseq <- seq(0,400,20)
  grid <- expand.grid(seq(0,50000,10),tseq)
  dg   <- data.frame(spacelag=grid[,1],timelag=grid[,2])
  table<- variogramSurface(model,dg)
  nx   <- length(seq(0,50000,10))
  ny   <- length(seq(0,15,1))
  plot(vario[vario$dist==0,5],vario[vario$dist==0,3],ylim=range(c(0,vario[vario$dist==0,3],
       table[table[,1]==0,3]),na.rm=T))
  lines(tseq,table[table[,1]==0,3])
}

######### we calculate the empirical space-time variogram
vario <- variogramST(PM10res~1,data=OBS_sel)

######### separable model, least squares fit
separableModel <- vgmST("separable",space=vgm(0.9,"Exp",10000,0.1),
                        time=vgm(0.9,"Exp",3.5,0.1),sill=40)
separable_fit  <- fit.StVariogram(model=separableModel,object=vario)
plot(vario,separable_fit,all=T,map=F) 
#(this does not look very convincing)
#(maybe some stability issues in the search for optimal parameters)
######### separable model, "manual" fit ("eye fit")
sill <- 250
separable_man <- vgmST("separable",space=vgm(.18,"Exp",1.5e5,0.01,
                       add.to=vgm(.12,"Exp",8e3,.01)),time=vgm(80/sill,"Exp",1670,19/sill,
                       add.to=vgm(66/sill,"Exp",85,0)),sill=sill)
plot(vario,separable_man,all=T,map=F)
variot(separable_man,vario)
######### product sum model, least squares fit
ProductSum     <- vgmST("productSum",space =vgm(9,"Exp",8e3,1),time=vgm(8,"Exp",106,2),k=2)
ProductSum_fit <- fit.StVariogram(model=ProductSum,object=vario)
plot(vario,ProductSum_fit,all=T,map=F)
variot(ProductSum_fit,vario)
######### product sum model, manual fit
ProductSum_man <- vgmST("productSum",space=vgm(12.5,"Exp",2e5,0,add.to=vgm(10,"Exp",9e3,0)),
                        time=vgm(35,"Exp",800,9.5,add.to=vgm(36,"Exp",85,0)),k=0.035) 
plot(vario,ProductSum_man,all=T,map=F)
variot(ProductSum_man,vario)
######### metric model, "manual" fit
Metric <- vgmST("metric",joint=vgm(60,"Exp",8e3,10,add.to=vgm(70,"Exp",1e5,0,stAni=1000)),
                stAni=150);model=Metric
plot(vario,Metric,map=F,all=T) 
variot(Metric,vario)
######### sum-metric model, manual fit
SumMetric <- vgmST("sumMetric",space=vgm(30,"Exp",20000),time=vgm(40,"Exp",70,15),
                   joint=vgm(50,"Exp",1.6e5,0),stAni=1000)
plot(vario,SumMetric,map=F,all=T)
plot(vario,SumMetric,map=T,diff=T)
variot(SumMetric,vario)
######### create list of all models fitted with gstat
models <- list(separable_fit,separable_man,ProductSum_fit,ProductSum_man,Metric,SumMetric)
######### save list of all models
save(models,file="data/fits_gstat.Rdata")
```


5) Prediction
===
5.1) Data Preparation 
---
```{r pred data preparation, eval=TRUE,echo=TRUE, tidy=TRUE}
rm(list=ls())
library(CompRandFld)
library(fields)
library(gstat)
library(rgdal)
library(sp)
library(spam)
library(spacetime)
library(RColorBrewer)
#library(raster)
library(rworldmap)
library(maptools)
library(ggmap)

# data preparation
#-----------------
DATA = 'data/'
load('data/data_mat.Rdata')
load('data/stations_fit.Rdata')
load('data/times_fit.Rdata')
load('data/CHM.Rdata')
load("DATA/STFDF_daily.Rdata")

### Getting the CHIMERE grid coordinates
coord.chm    <- expand.grid(unique(CHM$lon),unique(CHM$lat))
coord.spchm <- SpatialPoints(coord.chm, proj4string=CRS("+init=epsg:4326"))
# transform to Lambert 93 coordinate projection for approximately correct spatial geodesic distances
coord.spchm <- spTransform(coord.spchm, CRSobj=CRS("+init=epsg:2154"))
loc_to_pred <- coord.spchm@coords/1000   # prediction
CHM_mat <- matrix(CHM$PM10, nrow=nrow(loc_to_pred))
dim(CHM_mat)

### Getting the stations coordinates
coord.sta <- cbind(stations_fit$station_longitude_deg,
                   stations_fit$station_latitude_deg)
coord.spsta  <- SpatialPoints(coord.sta, proj4string=CRS("+init=epsg:4326"))
# transform to Lambert 93 coordinate projection for approximately correct spatial geodesic distances
coord.spsta <- spTransform(coord.spsta, CRSobj=CRS("+init=epsg:2154"))
coord.obs   <- coord.spsta@coords/1000   # coordinates of the selected stations

world <- getMap(resolution="low")
proj4string(world)
worldL93        <- spTransform(world,CRSobj=CRS("+init=epsg:2154"))
europe.limits <- geocode(c("Spain","Greece","Ireland","Norway","Ukraine"))

ns <- ncol(data_mat)    # number of stations
nT <- nrow(data_mat)    # number of times

times_to_pred <- c(88,90:92)   # prediction times
times_to_pred.st <- STFDF_daily@time[times_to_pred]
tobs = 88:90 # observation times
times.st <- STFDF_daily@time[tobs]
ind.obs = which(is.element(times_fit,tobs))
dobs = data_mat[ind.obs,]

mean.est <-   mean(data_mat)
mean.est
var.est  <-   var(as.numeric(data_mat))
var.est

data2 <- as.data.frame(as.vector(t(dobs)))
names(data2) <- 'PM10res'
REScor <- STFDF(sp=coord.spsta,time=times.st,data=data2)
stplot(REScor[,"2014-03-29::2014-03-31",'PM10res'],
       col.regions=tim.colors(100),colorkey=TRUE, 
       sp.layout=list(worldL93,first=FALSE))

```

5.2) A) Kriging with *CompRandFld*
---
```{r Kriging with CompRandFld,eval=TRUE,echo=TRUE, tidy=TRUE} 
load("data/fitsPL.Rdata") # load the fitted models
nummod    <- 8
corrmodel <- fitsPL[[nummod]]$corrmodel
param <- c(as.list(fitsPL[[nummod]]$param),as.list(fitsPL[[nummod]]$fixed))

pr <- Kri(loc=loc_to_pred,time=times_to_pred,coordx=coord.obs,coordt=tobs,
          corrmodel=corrmodel, param=param, data=dobs)

preddata <- as.data.frame(cbind(as.double(t(pr$pred)),
                                as.double(t(sqrt(pr$varpred))),
                                as.double(CHM_mat[,times_to_pred]),
                                as.double(CHM_mat[,times_to_pred]+t(pr$pred))))

names(preddata) <- c("Predicted residuals","Stdev predicted residuals",
                     "CHIMERE","CHIMERE corrected")

CHMcor = STFDF(sp=coord.spchm,time=times_to_pred.st,data=preddata)

ivar = 1
stplot(CHMcor[,1:4,ivar],col.regions=tim.colors(100),colorkey=TRUE,
    sp.layout=list(coord.spsta,worldL93,cex=.5,pch=19,first=FALSE,col=1))
ivar = 2
stplot(CHMcor[,1:4,ivar],col.regions=tim.colors(100),colorkey=TRUE,
    sp.layout=list(coord.spsta,worldL93,cex=.5,pch=19,first=FALSE,col=1))

cuts = c(seq(0,100,5),seq(110,150,10))
ivar = 3
stplot(CHMcor[,1:4,ivar],col.regions=tim.colors(100),colorkey=TRUE,
       cuts=cuts,main ="CHIMERE",sp.layout=list(coord.spsta,worldL93,
       cex=.5,pch=19,first=FALSE,col=1))

ivar = 4
stplot(CHMcor[,1:4,ivar],col.regions=tim.colors(100),colorkey=TRUE,
       cuts=cuts, main = "CHIMERE corrected",
       sp.layout=list(coord.spsta,worldL93,   
       cex=.5,pch=19,first=FALSE,col=1))

```
5.2) B) Cross-validation with *CompRandFld*
```{r Cross validation,eval=TRUE,echo=TRUE, tidy=TRUE}


err  <- matrix(NA,ns,8)  # LOO all stations
time <- times_to_pred[2] # predict 31 March 2014 
for(nummod in 1:8){
  corrmodel <- fitsPL[[nummod]]$corrmodel
  param <- c(as.list(fitsPL[[nummod]]$param),as.list(fitsPL[[nummod]]$fixed))
  for ( ista in 1:ns){
    loc    <- coord.obs[ista,]  
    coordx <- coord.obs[-ista,]
    data   <- dobs[,-ista]
    prvc   <- Kri(loc=loc,time=time,coordx=coordx,coordt=tobs,
                corrmodel=corrmodel, param=param, data=data)
    err[ista,nummod] = data_mat[ind.obs[3],ista]-prvc$pred
  }
}
RMSE <- sqrt(apply(err^2,2,mean))
RMSE

```

5.2) C) External validation with *CompRandFld*

```{r CaompRndFld external validation, eval=TRUE,echo=TRUE, tidy=TRUE}

## ----CompRndFld external validation, eval=TRUE,echo=TRUE, tidy=TRUE-----
load('data/OBS_daily.Rdata'); OBS_daily <- OBS_jour
load('data/stations.Rdata')

# extract French background stations
IDs.bg   <- stations$station_european_code[stations$type_of_station=="Background"]
IDs.fr   <- stations$station_european_code[stations$country_name=="France"]
IDs.bgfr <- intersect(IDs.bg,IDs.fr)
idx.stations <- which(stations$station_european_code %in% IDs.bgfr)
coord.bgfr   <- cbind(stations$station_longitude_deg,stations$station_latitude_deg)[idx.stations,]

OBS_bgfr <- OBS_daily[OBS_daily$ID %in% IDs.bgfr,]
OBS_sel2 <- stConstruct(OBS_bgfr,space=c('long','lat'),time='date',
                        SpatialObj=SpatialPoints(OBS_bgfr[,c('long','lat')]))
OBS_sel2 <- as(OBS_sel2,"STFDF")
proj4string(OBS_sel2) <- "+init=epsg:4326"
# transforming to Lambert 93 coordinates system
OBS_sel2@sp <- spTransform(OBS_sel2@sp, CRS("+init=epsg:2154"))

# defining time window for data used for validation
tstart <- "2014-04-01"
tend   <- "2014-04-02" # at least two dates !

OBS_sel2 <- OBS_sel2[,paste0(tstart,"::",tend)]

### interpolating gridded CHIMERE values to space-time observation points

CHM_sel  <- CHM[as.character(CHM$time)>=tstart&as.character(CHM$time)<=tend,]
CHM_mat2 <- matrix(CHM_sel$PM10, nrow=nrow(coord.chm))
dim(CHM_mat2)
lon.grid <- sort(unique(CHM_sel$lon))
lat.grid <- sort(unique(CHM_sel$lat))
ind.sites.obs <- match(unique(OBS_sel2@data$ID),IDs.bgfr)
fun <- function(i){
  tmp <- list(x=lon.grid,y=lat.grid,z=matrix(CHM_mat2[,i],length(lon.grid),length(lat.grid)))
  interp.surface(tmp, coord.bgfr[ind.sites.obs,])
}
CHM_sites <- sapply(1:ncol(CHM_mat2),fun)
OBS_sel2@data$CHM <- as.numeric(CHM_sites) # saving as variable CHM

time   <- times_to_pred[3]
loc    <- OBS_sel2@sp@coords/1000

RMSE_ext <- NULL
val <- (OBS_sel2[,1,'PM10'])$PM10
indexnona <- which(!is.na(val))
CHMval    <- (OBS_sel2[,1,'CHM'])$CHM[indexnona]
val <- val[indexnona]
for(nummod in 1:8){
  corrmodel <- fitsPL[[nummod]]$corrmodel
  param <- c(as.list(fitsPL[[nummod]]$param),as.list(fitsPL[[nummod]]$fixed))
  prvc <- Kri(loc=loc,time=time,coordx=coord.obs,coordt=tobs,
              corrmodel=corrmodel, param=param, data=dobs)
  RMSE_ext <- c(RMSE_ext,sqrt(mean((val-prvc$pred[indexnona]-CHMval)^2)))
}
RMSE_ext
### error made by CHIMERE
sqrt(mean((val-CHMval)^2))
```


5.3) A) Kriging with *gstat*
---
```{r kriging with gstat,eval=TRUE,echo=TRUE, tidy=TRUE}
load('data/fits_gstat.Rdata')

### Spatio-Temporal format 
predGridTime <- STF(coord.spchm,times_to_pred.st)
nummod <- 4 # Product sum manual
krig <- krigeST(PM10res~1,
                data = REScor[,seq(as.Date("2014-03-29"),as.Date("2014-03-31"),by="day")],
                newdata = predGridTime, modelList = models[[nummod]],
                computeVar=TRUE)

preddata_gstat <- as.data.frame(cbind(as.double(krig@data[[1]]),
                  sqrt(as.double(krig@data[[2]])), 
                  as.double(CHM_mat[,times_to_pred]+krig@data[[1]])))
names(preddata_gstat) <- c("Predicted residuals",
                           "Stdev predicted residuals",
                           "CHIMERE corrected")

CHMcor_gstat = STFDF(sp=coord.spchm,time=times_to_pred.st,data=preddata_gstat)

stplot(CHMcor_gstat[,,1],col.regions=tim.colors(100),colorkey=TRUE,
       sp.layout=list(coord.spsta,worldL93,cex=.5,pch=19,first=FALSE,col=1))

stplot(CHMcor_gstat[,,2],col.regions=tim.colors(100),colorkey=TRUE,
       sp.layout=list(coord.spsta,worldL93,cex=.5,pch=19,first=FALSE,col=1))

cuts = c(seq(0,100,5),seq(110,150,10))
stplot(CHMcor_gstat[,,3],col.regions=tim.colors(100),colorkey=TRUE,
       cuts=cuts, main = "CHIMERE corrected",
       sp.layout=list(coord.spsta,worldL93,cex=.5,pch=19,first=FALSE,col=1))
stplot(krig,col.regions=tim.colors(100),colorkey=TRUE,
       main="residuals prediction",
       sp.layout=list(coord.spsta,worldL93,cex=.5,pch=19,first=FALSE,col=1))

```
5.3) B) Cross Validation with *gstat*
```{r cross validation, eval=TRUE,echo=TRUE, tidy=TRUE}
dates = STFDF_daily@time[times_to_pred[1:2]]
err   <- matrix(NA,ns,6)
suppressMessages(for(nummod in 1:6){
  corrmodel <- models[[nummod]]
  for ( ista in 1:ns){
    newdata <- STF(coord.spsta[ista],dates)
    data = REScor[-ista,seq(as.Date("2014-03-29"),as.Date("2014-03-31"),by="day")]
    krig = krigeST(PM10res~1, data = data, newdata = newdata, 
                   modelList = corrmodel)
    err[ista,nummod] = data_mat[nT,ista]-krig@data[2,]
  }
})
RMSE_gs <- sqrt(apply(err^2,2,mean))
RMSE_gs

```
5.3) C) External validation with *gstat*
```{r  gstat external validation,eval=TRUE,echo=TRUE, tidy=TRUE}
dates <- STFDF_daily@time[times_to_pred[3:4]]
RMSE_extgs <- NULL
suppressMessages(newdata <- STF(OBS_sel2@sp,dates))
for(nummod in 1:6){
  corrmodel <- models[[nummod]]
  krig <- krigeST(PM10res~1, data = REScor[,seq(as.Date("2014-03-29"),as.Date("2014-03-31"),by="day")],
                  newdata <- newdata, modelList = corrmodel)
  RMSE_extgs <- c(RMSE_extgs,sqrt(mean((val-(krig[,1,1])$var1.pred[indexnona]-CHMval)^2)))
}
RMSE_extgs

```

