Dailymean.Rmd

---
title: Daily Mean Radiation Products from Ground-based observations and Satellite
  data
author: "Marieke Dirksen"
date: "September, 2016"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,message=FALSE,warning=FALSE)
#options(defaultPackages = "/usr/people/dirksen/R-20160725/x86_64-redhat-linux-gnu-library/3.3/")

```

#Loading Packages
The data from the hdf5 files is loaded with the R package rhdf5. From the .h5 files a raster is created, filped and croped using the R's package raster. The data from the ground based observations is read with the data.table function fread. For the geostatistical analysis we need the automap package. The machine learning algorithms from the caret package are used. 
```{r}
library(adehabitat)
library(automap)
library(caret)
# library(caretEnsemble)
library(data.table)
# library(doParallel)
# library(foreach)
library(GSIF)
library(kernlab)
library(maptools)
library(reshape)
library(raster)
library(rgdal)
library(rhdf5)
library(SDMTools)
data(wrld_simpl)
```

#Color settings for plotting routine
For all the plotting routines we use the same color settings. 
```{r}
kleur.breaks<-seq(0,500,by=10)
kleur.cols<-terrain.colors(length(kleur.breaks-1))
```

#Functions
Defining some basic functions. 
```{r,message=FALSE}
rmse<-function(sim,obs){
  rmse=sqrt(mean((sim-obs)^2,na.rm=TRUE))
}

r2<-function(residual,observed){
  teller <- sum(residual^2)
  noemer <- sum((observed-mean(observed))^2)
  r2 <- 1 - teller/noemer
  return(r2)
}
```

## Surface irradiance
Loading the attributes "direct irradiance" and "diffuse irradiance" to calculate the total irradiance. Also the satellite grid coordinates "lat" and "lon" are loaded.
```{r,message=FALSE}
time<-as.POSIXct("2014-07-04")
file<-"/net/pc150398/nobackup_1/users/meirink/siccs_wouter/SICCS/daymean/harmonie_proj/daymean_reproj_20140704.h5"
h5ls(file)
data.direct.irradiance<-h5read(file,"direct irradiance")
data.diffuse.irradiance<-h5read(file,"diffuse irradiance")
data.total.irradiance<-data.direct.irradiance+data.diffuse.irradiance
data.total.irradiance<-t(data.total.irradiance)

data.lat<-h5read(file,"/lat") #Latitude
data.lon<-h5read(file,"/lon") #Longitude
```

## Satellite grid
Now we have all the satellite data we create a raster object. 
```{r,message=FALSE}
data.lat[which(data.lat==-999)]<-NA # replace no data value with NA (this case -999)
data.lon[which(data.lon==-999)]<-NA # replace no data value with NA (this case -999)

r<-raster(data.total.irradiance,crs=CRS("+init=epsg:4326"),
      xmn=min(data.lon,na.rm=T),
      xmx=max(data.lon,na.rm=T),
      ymn=min(data.lat,na.rm=T),
      ymx=max(data.lat,na.rm=T))
rr<-flip(r,direction='y')
plot(rr,main=paste("Datum=",time),col=kleur.cols,breaks=kleur.breaks,legend=F)
plot(wrld_simpl,add=TRUE)
```

#Reprojecting raster and cropping
```{r,message=FALSE}
pro=CRS("+init=epsg:28992")
#Natural Earth dataset: unprojected shape files
mymap.unpro=readOGR(dsn='Rdata/NaturalEarthData/ne_10m_admin_0_countries',layer="ne_10m_admin_0_countries") # Read in (unprojected) map data
mymap.pro=spTransform(mymap.unpro, pro) # Reproject the map

mymap.unpro_lakes=readOGR(dsn='Rdata/NaturalEarthData/ne_10m_lakes',layer="ne_10m_lakes") # Read in (unprojected) map data
mymap.pro_lakes=spTransform(mymap.unpro_lakes, pro) # Reproject the map

rr<-projectRaster(rr,crs=pro)
r.NED<-crop(rr,extent(12621.630033977,278621.630033977,305583.0457758,620583.0457758))

plot(r.NED,col=kleur.cols,breaks=kleur.breaks,main=paste("Datum=",time),legend=F)
plot(mymap.pro,add=TRUE)
plot(mymap.pro_lakes,add=TRUE)
```

<!-- From the [KNMI daily record](http://projects.knmi.nl/klimatologie/daggegevens/index.cgi) the following data is available from this day: -->
<!-- ![Figure](/nobackup/users/dirksen/radiation/figures/20140707.png) -->

#Data from observations
From the observations a subset is made (2014/07/04), coordinates are stored in a seperate file and merged with the dataset. 
```{r,message=FALSE}
# coords<-fread("Rdata/radiation_KNMI_day.csv")
# coords<-subset(coords,select=c("DS_LAT","DS_LON","DS_CODE"))
# coords<-unique(coords)
#saveRDS(coords,file="Rdata/coordsKNMI.rda")
coords<-readRDS("Rdata/coordsKNMI.rda")

obs<-fread("Rdata/radiation_KNMI_day_v2.csv")
obs$IT_DATETIME<-as.POSIXct(obs$IT_DATETIME,format="%Y%m%d_240000_000000")

obs.subset<-obs[which(IT_DATETIME==time),]
obs.subset$Q<-obs.subset$REH1.Q24*10000/(24*3600)
obs.subset<-merge(obs.subset,coords,by="DS_CODE")
obs.subset<-na.omit(obs.subset,cols=c("DS_LAT","DS_LON","Q"))

coordinates(obs.subset)<-~DS_LON+DS_LAT
proj4string(obs.subset)<-CRS("+init=epsg:4326")
obs.subset<-spTransform(obs.subset,pro)

kleur.data<-cut(obs.subset$Q,breaks=kleur.breaks,include.lowest=TRUE,labels=FALSE)
plot(obs.subset,bg=(kleur.cols[kleur.data]),
     col="darkgrey",pch=21)
text(obs.subset$DS_LON,obs.subset$DS_LAT,round(obs.subset$Q,0),pos=3,cex=0.7)
plot(mymap.pro,add=TRUE)
plot(mymap.pro_lakes,add=TRUE)
```

##Test and Train set
For the splitting of the train and test set we want all the data from 1 day, as we only have 33 observations. A post on splitting with timeslices can be found [here](http://stackoverflow.com/questions/24758218/time-series-data-spliting-and-model-evaluation). Also on Github some more background information can be [found](http://topepo.github.io/caret/data-splitting.html). 

> Note that these methods can also be included in the trainControl function (see stackoverflow link). The trainControl would look something like: trainControl(method='cv', index=createFolds(obs,list = TRUE)). Here we explore how the functions work for subsetting the data. Another important discussion on cross-validation and data splitting can be found [here](http://stats.stackexchange.com/questions/188955/is-it-necessary-to-split-dataset-for-cross-validation). 

```{r}
obs<-obs[complete.cases(obs$REH1.Q24),]
#timestamps<-unique(obs$IT_DATETIME)

#Create Folds
folds<-createFolds(obs,2)

trainFolds<-folds[[1]]
testFolds<-folds[[2]]

trainSet<-obs[trainFolds]
testSet<-obs[testFolds]

print(head(trainSet))
print(head(testSet))

#Time slices: slow
# I.Sub.TimeSlices<-createTimeSlices(1:nrow(obs),initialWindow = 32,horizon=32,fixedWindow = TRUE)
# 
# str(I.Sub.TimeSlices,max.level=1)
# 
# 
# trainSlices<-I.Sub.TimeSlices[[1]]
# testSlices<-I.Sub.TimeSlices[[2]]
# 
# trainSet<-obs[trainSlices[[1]]]
# testSet<-obs[testSlices[[1]]]
# 
# print(trainSet)
# print(testSet)

#create Data Partition
I.Sub<-createDataPartition(obs$REH1.Q24,p=0.05,list=FALSE)


```

#Compare the two products
```{r,message=FALSE}
rASC<-asc.from.raster(r.NED)
spdf<-asc2spixdf(rASC)
proj4string(spdf)<-pro

var.Q<-subset(obs.subset,select=Q)
sat.var<-over(var.Q,spdf)

n<-names(sat.var)
diff<-sat.var[n]-var.Q$Q

plot(r.NED,col=kleur.cols,breaks=kleur.breaks,main=paste("Datum=",time),legend=F)
plot(obs.subset,bg=(kleur.cols[kleur.data]),
     col="darkgrey",pch=21,add=TRUE)
text(obs.subset$DS_LON,obs.subset$DS_LAT,round(obs.subset$Q,0),pos=3,cex=0.7)
text(obs.subset$DS_LON,obs.subset$DS_LAT,round(diff[[n]],0),pos=1,cex=0.7,col="red")
plot(mymap.pro,add=TRUE)
plot(mymap.pro_lakes,add=TRUE)
```

#Geostatistical Approach: Kriging interpolation with trend
```{r,message=FALSE}
# Kriging
mxdkrige=Inf # maxdist Krige
# over functions
gridded(spdf)=FALSE;gridded(spdf)=TRUE;fullgrid(spdf) = TRUE
slot(slot(spdf, "grid"), "cellsize") <-rep(mean(slot(slot(spdf, "grid"), "cellsize")), 2)
# over Distshore on Var
distshore.ov=over(obs.subset,spdf)
# Copy the values to Var )
  var = obs.subset
  var$var=distshore.ov$var

  #Prepare input
  field = spdf
  field@data = cbind(field@data, coordinates(field))
  names(field@data) = c("s","x","y")
  var$x = over(var,field)$x
  var$y = over(var,field)$y
  var$s = over(var,field)$s

  # Remove nodata from dataframe based on missing distshore
  var = var[!is.na(var$var),]


ked_exp <- autoKrige(Q~var, var, spdf,maxdist=mxdkrige, model = c("Exp"), na.action=na.pass, fix.values=c(NA,NA,NA), miscFitOptions = list(merge.small.bins = FALSE))  #log(distshore)

# Krige Cross validation
ked_exp.cv <- autoKrige.cv(Q~var, var, model = c("Exp"),maxdist=mxdkrige,fix.values=c(NA,NA,NA), miscFitOptions = list(merge.small.bins = FALSE),verbose=c(FALSE,FALSE))
teller <- sum(ked_exp.cv$krige.cv_output$residual^2)
noemer <- sum((var$var-mean(var$var))^2)
ked_exp.r2 <- 1 - teller/noemer
ked.zscoremean <- mean(ked_exp.cv$krige.cv_output$zscore)
ked.zscore.var <- var(ked_exp.cv$krige.cv_output$zscore)
ked_exp.rmse<-rmse(ked_exp.cv$krige.cv_output$var1.pred,ked_exp.cv$krige.cv_output$observed)
plot(ked_exp,col=kleur.cols,breaks=kleur.breaks,sp.layout=list(pts=list("sp.points",obs.subset,pch=21),mymap.pro,mymap.pro_lakes))

print(paste("R2=",round(ked_exp.r2,2),"RMSE=",round(ked_exp.rmse,2)))
print(head(ked_exp$krige_output))
```

# Caret: Machine Learning Algorithms
For grid predictions with the caret package a combination of the caret::train and raster::predict functions is used. The input for the raster::predict function is either a raster Stack or raster Brick.

* Methods [here](http://topepo.github.io/caret/train-models-by-tag.html#Two_Class_Only.html).
* Some methods of the caret package run in parallel, background information can be found [here](https://cran.r-project.org/web/packages/doParallel/vignettes/gettingstartedParallel.pdf).
* Bug in Caret "Error in e$fun(obj, substitute(ex), parent.frame(), e$data) :
  worker initialization failed: there is no package called ‘caret’" can be fixed like [this](http://stackoverflow.com/questions/21029019/parallel-execution-of-train-in-caret-fails-with-function-not-found)

The machine learning algorithms (MLA) used in this example are all regression based, as we have only 2 variables (ground-based measurements vs. satellite derived product). We compare: 2 linear models, a support vector machine, treebag and cubist. Below a short description of the models and their tuning methods.

### Linear Models
Two linear models are compared: the linear model (lm) and gaussprLinear.

### Support Vector Machines
Sigma is depending on the predictor. The function "sigest" provides the value for sigma.

### Treebag
The treebag model also has no tuneGrid.

### Cubist
Cubist is a rule-based model but differs from other tree models. The final model combines models using a linear combination of two models. The different models are weighted based on their RMSE. The final model makes up the initial set of rules. The model tunes itself using committees and neighbors.

## Data preperation
```{r,message=FALSE,warning=FALSE}
#cl<-makeCluster(6)
#registerDoParallel(cl)

obs.subset.Q<-subset(obs.subset,select="Q")
grid<-as(r.NED,"SpatialPixelsDataFrame")

ov<-over(obs.subset.Q,grid)
ov<-cbind(data.frame(obs.subset.Q["Q"]),ov)
ov<-rename(ov,c("DS_LON"="x"))
ov<-rename(ov,c("DS_LAT"="y"))
```

## Train control settings and tuning parameters
* The github post on [model training and tuning](http://topepo.github.io/caret/model-training-and-tuning.html).
* A list with all the available regression models can be found [online](http://topepo.github.io/caret/available-models.html), search for regression.
```{r}
control<-trainControl(method="cv",number=10,repeats=3) #setting a 10 fold-cross validation (best performing)
length<-10 #for the tuneLength of the models

#uncertainty measurements (calibration): 10W/m2, uncertainty sat: 30-50W/m2.
sigmaRangeReduced<-sigest(as.matrix(ov$layer))[1]
svmRadialRGridReduced<-expand.grid(.sigma=sigmaRangeReduced,.C=2^(seq(-4,4)))
svmLinearRGridReduced<-expand.grid(C=2^(seq(-4,4)))
ctreeGridReduced<-expand.grid(mincriterion=seq(from=0.1,to=0.5,by=0.1))
knnGridReduced<-expand.grid(.k=3:15)
gbmGridReduced<-expand.grid(.shrinkage=c(0.1),.n.trees=1:5,.interaction.depth=1:6,.n.minobsinnode=1:5)
earthGridReduced <- data.frame(.degree = 1, .nprune = (2:4)*2)
```

```{r}
#Linear Models
set.seed(50)
m1.lm<-caret::train(Q~layer,data=ov,method="lm",preProcess=c("center","scale","BoxCox"),tuneLength=length,trControl=control)

set.seed(50)
m2.glm<-caret::train(Q~layer,data=ov,method="glm",preProcess=c("center","scale","BoxCox"),tuneLength=length,trControl=control)

set.seed(50)
m3.gaussprLinear<-caret::train(Q~layer,data=ov,method="gaussprLinear",preProcess=c("center","scale","BoxCox"),verbose=FALSE,tuneLength=length,trControl=control)

#Support Vector Models
set.seed(50)
m4.svmRadial<-caret::train(Q~layer,data=ov,method="svmRadial",preProcess=c("center","scale","BoxCox"),verbose=FALSE,tuneLength=length,trControl=control,tuneGrid=svmRadialRGridReduced)

set.seed(50)
m5.svmLinear<-caret::train(Q~layer,data=ov,method="svmLinear",preProcess=c("center","scale","BoxCox"),verbose=FALSE,tuneLength=length,trControl=control,tuneGrid=svmLinearRGridReduced)

#Tree models
set.seed(50)
m6.treebag<-caret::train(Q~layer,data=ov,method="treebag",preProcess=c("center","scale","BoxCox"),verbose=FALSE,tuneLength=length,trControl=control)

set.seed(50)
m7.cubist<-caret::train(Q~layer,data=ov,method="cubist",preProcess=c("center","scale","BoxCox"),verbose=FALSE,tuneLength=length,trControl=control)

set.seed(50)
m8.ctree<-caret::train(Q~layer,data=ov,method="ctree",preProcess=c("center","scale","BoxCox"),tuneLength=length,trControl=control,tuneGrid=ctreeGridReduced)

#K-nearest neighbors
set.seed(50)
m9.knn<-caret::train(Q~layer,data=ov,method="knn",preProcess=c("center","scale","BoxCox"),verbose=FALSE,tuneLength=length,trControl=control,tuneGrid=knnGridReduced)

# Others
set.seed(50)
m10.earth<-caret::train(Q~layer,data=ov,method="earth",preProcess=c("center","scale","BoxCox"),tuneLength=length,trControl=control,tuneGrid=earthGridReduced)

set.seed(50)
m11.gbm<-caret::train(Q~layer,data=ov,method="gbm",preProcess=c("center","scale","BoxCox"),tuneLength=length,trControl=control,tuneGrid=gbmGridReduced,verbose=FALSE)
```

```{r}
results<-resamples(list(lm=m1.lm,
                        glm=m2.glm,
                        gaussprLinear=m3.gaussprLinear,
                        svmRadial=m4.svmRadial,
                        svmLinear=m5.svmLinear,
                        treebag=m6.treebag,
                        cubist=m7.cubist,
                        ctreer=m8.ctree,
                        knn=m9.knn,
                        earth=m10.earth,
                        gbm=m11.gbm))
summary(results)
bwplot(results,scales=list(relation="free"),xlim=list(c(0,70),c(0,1)))
modelDifferences<-diff(results)
print(modelDifferences) #see str(modelDifferences) for values. Large p-values: the models fail to show any diference in performance (Kuhn&Johnson 2013, p.101 )
```

## Model selection
In the previous section we found that differences between the models are too small to select the best performing model. What we don't want is an under- or over-predicting model. Comparing predictions, observations and residuals helps to find the best fit.

Comparing the models showed the following results:

* The model fits all plot around the 1:1 line
* The range of the predictions is generally smaller than the observations (especially for gaussprLinear and rlm)
* Residuals are generally equally distributed (linear models have a slight tendency towards an under-fit)
* The cubist and treebag explain most of the variance in the observations
* Also after tuning (add tuneGrid) the svmRadial has a good performance

```{r}
final.model<-m5.svmLinear

st<-stack(grid)
p1.caret<-raster::predict(model=final.model,object=st)

p1.ASC<-asc.from.raster(p1.caret)
p1.spdf<-asc2spixdf(p1.ASC)
proj4string(p1.spdf)<-pro

p1.var<-over(var.Q,p1.spdf)

n<-names(p1.var)
p1.diff<-p1.var[n]-var.Q$Q

observed<-data.frame(var.Q$Q)
predicted<-p1.var[n]

axisRange<-extendrange(c(observed,predicted))
obs.vs.pred<-data.frame(observed,predicted)
pred.vs.res<-data.frame(predicted,p1.diff)

plot(obs.vs.pred,xlab="observed",ylab="predicted",ylim=axisRange,xlim=axisRange,asp=1)
abline(0, 1, col = "darkgrey", lty = 2)

plot(pred.vs.res,xlab="predicted",ylab="residual",asp=1)
abline(h = 0, col = "darkgrey", lty = 2)

plot(p1.caret,col=kleur.cols,breaks=kleur.breaks,main=paste("Datum =",time,"\n","method = ",final.model$method),legend=F)
plot(obs.subset,bg=(kleur.cols[kleur.data]),
     col="darkgrey",pch=21,add=TRUE)
text(obs.subset$DS_LON,obs.subset$DS_LAT,round(obs.subset$Q,0),pos=3,cex=0.7)
text(obs.subset$DS_LON,obs.subset$DS_LAT,round(p1.diff[[n]],0),pos=1,cex=0.7,col="red")
plot(mymap.pro,add=TRUE)
plot(mymap.pro_lakes,add=TRUE)

obs.vs.pred<-data.frame(observed,predicted)
pred.vs.res<-data.frame(predicted,p1.diff)
axisRange <- extendrange(c(observed, predicted))
```

## Comparing models
```{r}
st$lm<-raster::predict(model=m1.lm,object=st)
st$glm<-raster::predict(model=m2.glm,object=st)
st$gaussprLinear<-raster::predict(model=m3.gaussprLinear,object=st)
st$svmRadial<-raster::predict(model=m4.svmRadial,object=st)
st$svmLinear<-raster::predict(model=m5.svmLinear,object=st)
st$treebag<-raster::predict(model=m6.treebag,object=st)
st$cubist<-raster::predict(model=m7.cubist,object=st)
st$ctree<-raster::predict(model=m8.ctree,object=st)
st$knn<-raster::predict(model=m9.knn,object=st)
st$earth<-raster::predict(model=m10.earth,object=st)
st$gbm<-raster::predict(model=m11.gbm,object=st)
st$kriging<-ked_exp$krige_output$var1.pred
print(st)

spplot(st,col.regions=kleur.cols,at=kleur.breaks,
       sp.layout=list(pts=list("sp.points",obs.subset,pch=21),mymap.pro,mymap.pro_lakes))

```

<!-- ## Statistics from the first test run (2014) -->
<!-- A first test run shows the Machine learnin algorithms have statistically mutch better correlations that the kriging model. Although the RMSE is low R2 values are dramatically low. -->

<!-- ```{r} -->
<!-- statistics_first_run<-fread("/nobackup/users/dirksen/Radiation_Obs_Satellite/output/statistical_summary/RMSE_Rsquared.csv",header=T) -->
<!-- statistics_first_run<-subset(statistics_first_run,select=-c(1:3)) -->

<!-- I.RMSE<-statistics_first_run[,grep(".rmse",colnames(statistics_first_run),value=TRUE)] -->
<!-- I.Rsquared<-statistics_first_run[,grep(".r2",colnames(statistics_first_run),value=TRUE)] -->

<!-- RMSE<-subset(statistics_first_run,select=I.RMSE) -->
<!-- Rsquared<-subset(statistics_first_run,select=I.Rsquared) -->

<!-- summary.RMSE<-as.data.frame( t(sapply(RMSE, function(cl) list(mean=mean(cl,na.rm=TRUE), -->
<!--                                                                         sd=sd(cl,na.rm=TRUE), -->
<!--                                                                         min=min(cl,na.rm = TRUE), -->
<!--                                                                         max=max(cl,na.rm=TRUE))))) -->

<!-- summary.Rsquared<-as.data.frame( t(sapply(Rsquared, function(cl) list(mean=mean(cl,na.rm=TRUE), -->
<!--                                                                         sd=sd(cl,na.rm=TRUE), -->
<!--                                                                         min=min(cl,na.rm = TRUE), -->
<!--                                                                         max=max(cl,na.rm=TRUE))))) -->
<!-- print(summary.RMSE) -->
<!-- print(summary.Rsquared) -->
<!-- ``` -->


<!-- ```{r} -->
<!-- SICCS_vs_obs<-fread("/nobackup/users/dirksen/Radiation_Obs_Satellite/output/summary_differences/SICCS.csv",header=T) -->
<!-- SICCS_vs_obs<-subset(SICCS_vs_obs,select=-c(1:3)) -->

<!-- difference<-SICCS_vs_obs[,grep(".rmse",colnames(statistics_first_run),value=TRUE)] -->
<!-- I.Rsquared<-statistics_first_run[,grep(".r2",colnames(statistics_first_run),value=TRUE)] -->

<!-- difference<-subset(SICCS_vs_obs,select="difference") -->

<!-- summary.diff<-as.data.frame( t(sapply(difference, function(cl) list(mean=mean(cl,na.rm=TRUE), -->
<!--                                                                         sd=sd(cl,na.rm=TRUE), -->
<!--                                                                         min=min(cl,na.rm = TRUE), -->
<!--                                                                         max=max(cl,na.rm=TRUE))))) -->

<!-- print(summary.diff) -->
<!-- ``` -->


<!-- ## Differences for the grids on a clear sky day and others -->
<!-- A spatial comparisson between the grids and observations is made. 4 clear sky days 2008-05-07 until 2008-05-10 are selected. Also 3 other days from the report of van Tiggelen 2014 are selected. -->
<!-- ```{r} -->
<!-- fun<-function() { -->
<!--   plot(mymap.pro,add=TRUE) -->
<!--   plot(mymap.pro_lakes,add=TRUE) -->
<!-- } -->
<!-- # For the 4 clearsky reference days there are very small differences with observations -->
<!-- clearsky_ref1<-readRDS("/nobackup/users/dirksen/Radiation_Obs_Satellite/output/predictions/2008-05-07_st.rds") -->
<!-- clearsky_ref2<-readRDS("/nobackup/users/dirksen/Radiation_Obs_Satellite/output/predictions/2008-05-08_st.rds") -->
<!-- clearsky_ref3<-readRDS("/nobackup/users/dirksen/Radiation_Obs_Satellite/output/predictions/2008-05-09_st.rds") -->
<!-- clearsky_ref4<-readRDS("/nobackup/users/dirksen/Radiation_Obs_Satellite/output/predictions/2008-05-10_st.rds") -->

<!-- #Jurgen van Tiggelen used from 2010 the following days in his report: May12, June 21 and July 31 -->
<!-- #From all these days and the surrounding days differences are huge -->
<!-- #Is something wrong with the dataset? -->
<!-- May12<-readRDS("/nobackup/users/dirksen/Radiation_Obs_Satellite/output/predictions/2010-05-12_st.rds") -->
<!-- June21<-readRDS("/nobackup/users/dirksen/Radiation_Obs_Satellite/output/predictions/2010-06-21_st.rds") -->
<!-- July31<-readRDS("/nobackup/users/dirksen/Radiation_Obs_Satellite/output/predictions/2010-07-31_st.rds") -->
<!-- July30<-readRDS("/nobackup/users/dirksen/Radiation_Obs_Satellite/output/predictions/2010-07-30_st.rds") -->
<!-- Aug1<-readRDS("/nobackup/users/dirksen/Radiation_Obs_Satellite/output/predictions/2010-08-01_st.rds") -->

<!-- stacklayer<-June21 -->
<!-- kleur.min<-min(minValue(stacklayer)) -->
<!-- kleur.max<-max(maxValue(stacklayer))+10 -->
<!-- kleur.breaks<-seq(kleur.min,kleur.max,by=10) -->
<!-- kleur.cols<-colorRampPalette(c("green","yellow","red"))(length(kleur.breaks-1)) -->

<!-- arg<-list(at=range(kleur.breaks),labels=round(range(kleur.breaks),0)) -->
<!-- plot(stacklayer,col=kleur.cols,breaks=kleur.breaks,axis.arg=arg,addfun=fun) -->


<!-- ``` -->
<!-- ## Ensemble predictions with Caret -->
<!-- As we saw, the previous predictions from the caret package are already nice. But, can we further improve the model by combining several models? Here we explore the caret Ensemble package ([online example](https://cran.r-project.org/web/packages/caretEnsemble/vignettes/caretEnsemble-intro.html)). As the differences between the models are small an ensemble prediction is not expected to improve the prediction. Though, lets try and make a code and see how the ensemble model performs: -->
<!-- ```{r} -->
<!-- set.seed(50) -->
<!-- model_list<-caretList(REH1.Q~layer, -->
<!--                       data=ov, -->
<!--                       trControl=control, -->
<!--                       methodList=c("lm","cubist","treebag"), -->
<!--                       preProcess=c("center","scale","BoxCox"), -->
<!--                       tuneLength=length) -->
<!--                       #tuneList=list(svm=caretModelSpec(method="svmRadial", -->
<!--                       #                                 tuneGrid=svmRGridReduced))) -->
<!-- #NOTE: trControl not fine-tuned for caretEnsemble and caretStack! -->
<!-- model_ensemble<-caretEnsemble(model_list) -->

<!-- #not specified a rf ensemble is generated -->
<!-- lm_ensemble<-caretStack(model_list,trControl=trainControl(method='cv'),tuneGrid = expand.grid(.mtry=2:5),tuneLength=length) -->

<!-- print(lm_ensemble) -->
<!-- #ens_pred<-predict(model_ensemble,newdata=grid) -->
<!-- ens_pred<-predict(lm_ensemble,newdata=grid) -->
<!-- # predictors<-names(st)[names(st) != "layer"] -->
<!-- # final_ensemble<-caret::train(subset(st,predictors),subset(st,"layer"),method='treebag',trControl=control) -->
<!-- ``` -->

<!-- ## Creating a raster from the prediction -->
<!-- ```{r} -->
<!-- ens_raster<-raster(matrix(ens_pred,nrow=nrow(r.NED),ncol=ncol(r.NED))) -->
<!-- ens_raster@extent<-extent(r.NED) -->
<!-- proj4string(ens_raster)<-pro -->
<!-- plot(ens_raster,col=kleur.cols,breaks=kleur.breaks,main=paste("Datum =",time,"\n","Ensemble"),legend=F) -->
<!-- plot(mymap.pro,add=TRUE) -->
<!-- plot(mymap.pro_lakes,add=TRUE) -->
<!-- ``` -->