Remote sensing rainfall products: part two: SM2RAIN

SM2RAIN is satellite rainfall etimation algorithm developed by Luca Brocca and his research team. The idea is that it is possible to invert the soil moisture imagary to reanalyse the rainfall. References about it can be found here, here, here and references therein. I am using this rainfall product to estimate spatial rainfall over upper Ble Nile basin. Lets have qiuck look to the data.
As usual necessary library

library(raster)
## Loading required package: sp
library(rasterVis)
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
## Loading required package: hexbin
library(raster)
library(rgdal)
## rgdal: version: 0.9-1, (SVN revision 518)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/proj
install.packages("R.matlab")
## Error in contrib.url(repos, "source"): trying to use CRAN without setting a mirror
library(R.matlab)
## R.matlab v3.2.0 (2015-02-24) successfully loaded. See ?R.matlab for help.
## 
## Attaching package: 'R.matlab'
## 
## The following objects are masked from 'package:base':
## 
##     getOption, isOpen
library(lattice)

Data import: I download the dataset, which was originally in matlab format, to my pc and import it to R

SM2RF<-readMat("/~/P_SM2RAIN_ECVSM_BluNile_DISTR.mat")
proj <- CRS('+proj=longlat +ellps=WGS84')
SM2RF2 <- cbind(SM2RF$lonubn, SM2RF$latubn, t(SM2RF$Psim)) 
SM2RF3 <- rasterFromXYZ(SM2RF2)

The dataset is raster layer for 8766 days. To indix each layer by the time-steps, we can use time index as follows:

Timeseries <- seq(as.Date('1990-01-01'), as.Date('2013-12-31'), 'day')
SM2RF4 <- setZ(SM2RF3, Timeseries)

To extract point data based on location of ground measurments, we can create spatial points and extract the data from the layer using these spatial points
To create the spatial points of measurment stations:

ADET<-cbind(37.47,11.27); ADET<-SpatialPoints(ADET) 
BAHRDAR<-cbind(37.42,11.6); BAHRDAR<-SpatialPoints(BAHRDAR) 
CHAGNI<-cbind(36.00,10.95); CHAGNI<-SpatialPoints(CHAGNI) 
DANGILA<-cbind(36.42,11.12); DANGILA<-SpatialPoints(DANGILA) 
DEBREMARKOS<-cbind(37.67,10.33);DEBREMARKOS<-SpatialPoints(DEBREMARKOS) 
MOTTA<-cbind(37.87,11.08); MOTTA<-SpatialPoints(MOTTA) 
PAWE<-cbind(36.05,11.15);  PAWE<-SpatialPoints(PAWE) 
AYKEL<-cbind(37.05,12.53); AYKEL<-SpatialPoints(AYKEL) 
GONDER<-cbind(37.42,12.55); GONDER<-SpatialPoints(GONDER)
DEBREBRIHAN<-cbind(39.58,9.63);DEBREBRIHAN<-SpatialPoints(DEBREBRIHAN)

To extract the time series rainfall (1990-2013 daily rainfall) and create the dataframe

SM2RF_station<-data.frame(
time=seq(as.Date('1990-01-01'), as.Date('2013-12-31'), 'day'),
ADET=as.vector(extract(SM2RF3, ADET)),
BAHRDAR=as.vector(extract(SM2RF3, BAHRDAR)),
CHAGNI=as.vector(extract(SM2RF3, CHAGNI)),
DANGILA=as.vector(extract(SM2RF3, DANGILA)),
DEBREMARKOS=as.vector(extract(SM2RF3,DEBREMARKOS)),
MOTTA=as.vector(extract(SM2RF3, MOTTA)),
PAWE=as.vector(extract(SM2RF3, PAWE)),
AYKEL=as.vector(extract(SM2RF3, AYKEL)),
GONDER=as.vector(extract(SM2RF3, GONDER)),
DEBREBRIHAN=as.vector(extract(SM2RF3, DEBREBRIHAN)))

Visulisation: to visualize sample of time series and stations,subset of the datset is used. Here is the first 11 years, 1990-2000, from top to down, each plot for one year.

SM2RF_station2<-subset(SM2RF_station[0:4018, ])
xyplot(ADET+PAWE+MOTTA+GONDER+CHAGNI+DEBREMARKOS~time|equal.count(as.numeric(time), 11, overlap = 0.1),
SM2RF_station2,  
type = c("a","g"),   , strip = FALSE,
xlab="time", ylab="SM2RAIN rainfall",
as.table = TRUE, layout = c(1, 11),
auto.key = list(points = FALSE, lines = TRUE, columns = 3),
scales = list(x = list(relation = "sliced", axs = "i"),
              y = list(alternating = FALSE)))