In MariekeDirk/RSOI: Reduce optimal interpolation

knitr::opts_chunk$set(echo = TRUE)

RSOI

Reduced space optimal interpolation workshop 2019

For an interpolation area the following files need to be prepared:

empty grid layers for each month of the reconstruction period
observations for the 30 year reference period
observations for the reconstruction period

An overview of the observations

Indonesia

library(mapview)
data("sp.RR_stations")
mapview(sp.RR_stations)

The Netherlands

data("sp.TN_stations")
mapview(sp.TN_stations)

Datasets in this project

This project contains meteorological data from two countries: the Netherlands and Indonesia. From the Netherlands temperature observations since 1905 are included. For Indonesia precipitation observations are stored in data("RR_obs") starting in the year 1930.

#load Temperature data for the Netherlands
data("TN_observations") #observations
data("TN_reference")    #reference observations
data("sp_grid_nl")      #empty array
#data("ok_NL")           #ordinary kriging climatology (not used in RSOI)

#load Precipitation data for Indonesia
data("RR_obs")          #observations
data("RR_ref")          #reference observations
data("grid_matrix")     #with the dimensions from the raster stack [nrow,ncol,nlayers,fill=values]
data("RR_ref_grid")     #reference grid interpolated raster(not used in RSOI)

RSOI prepare: data availability

We are going to start with the Indonesia data for the analysis.

library(gridrsoi)
library(ggplot2)

RR_obs_count<-subset(RR_obs,select=-IT_DATETIME)
out<-apply(RR_obs_count,2,FUN=function(x){ as.numeric(!is.na(x))})
sums_out<-colSums(out)

df.out<-data.frame(RR_obs$IT_DATETIME,out)
df.out.long<-df.out %>% tidyr::gather(key=stations,value=value,-RR_obs.IT_DATETIME)
# ggplot(df.out.long,aes(RR_obs.IT_DATETIME,value,fill=value,group=stations))+geom_tile()+
#   #scale_x_discrete(expand = c(0, 0)) +
#   scale_y_discrete(expand = c(0, 1))

RSOI prepare: PCA analysis

For the PCA analysis you need a reconstruction on a grid for the reference period. The reconstruction could for example be calculated using ordinary kriging or IDW.

trafo <- function(x) {sqrt(x)}
back.trafo <- function( x) { x[x<0] <- 0; x^2 }


RR_pca<-gridrsoi::gridrsoi.prepare(grid=grid_matrix,transf = trafo,back.transf = back.trafo,maxL = 80 )

scree.plot(RR_pca$pca,num.pcs = 10)

Reconstruction using RSOI

Note: reconstruction of Indonesia does not work here since there are missing values in the recontruction period!

RR_obs_val<-as.matrix(RR_obs[,2:ncol(RR_obs)])
attr(RR_obs_val,"dimnames")<-list(as.character(RR_obs$IT_DATETIME),
                                               names(RR_obs[,2:ncol(RR_obs)]))
attr(RR_obs_val,"x")<-as.numeric(attr(RR_obs,"x")$lat)
attr(RR_obs_val,"y")<-as.numeric(attr(RR_obs,"y")$lon)

# RR_ref<-RR_ref[complete.cases(RR_ref),]
RR_ref_val<-as.matrix(RR_ref[,2:ncol(RR_ref)])
attr(RR_ref_val,"dimnames")<-list(as.character(RR_ref$IT_DATETIME),names(RR_ref[,2:ncol(RR_ref)]))
attr(RR_ref_val,"x")<-as.numeric(attr(RR_ref,"x")$lat)
attr(RR_ref_val,"y")<-as.numeric(attr(RR_ref,"y")$lon)

res<-gridrsoi.reconstruct(stdat.recon = RR_obs_val, 
                          stdat.calib = RR_ref_val,
                          preparation = RR_pca,
                          L = 6,
                          do.what = c("grid","xval"))