#different covariates using two baselines. Accuracy methods are added in the the function script to evaluate results.

#Figures, tables and data for the contribution of covariate paper are also produced in the script.

#AUTHOR: Benoit Parmentier                                                                      

#Version: 5

#PROJECT: Environmental Layers project                                     

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

library(automap)                             # Kriging automatic fitting of variogram using gstat

library(rgeos)                               # Geometric, topologic library of functions

#RPostgreSQL                                 # Interface R and Postgres, not used in this script

#Additional libraries not used in workflow

library(pgirmess)                            # Krusall Wallis test with mulitple options, Kruskalmc {pgirmess}  

library(ncf)

        sub(".shp","",basename(met_stations_obj$monthly_covar_ghcn_data)))

ghcn_dat_WGS84 <-spTransform(ghcn_dat,CRS_locs_WGS84)         # Project from WGS84 to new coord. system

interp_area <- readOGR(dsn=dirname(infile_reg_outline),sub(".shp","",basename(infile_reg_outline)))

interp_area_WGS84 <-spTransform(interp_area,CRS_locs_WGS84)         # Project from WGS84 to new coord. system

lst_md<- lst_md - 273.16

lst_mm_01<-subset(s_raster,"mm_01")

#temp.colors <- colorRampPalette(c('blue', 'white', 'red'))

#temp.colors <- colorRampPalette(c('blue', 'lightgoldenrodyellow', 'red'))

#temp.colors <- matlab.like(no_brks)

png(filename=paste("Figure_3_paper_Comparison_daily_monthly_mean_lst",out_prefix,".png",sep=""),width=960,height=480)

par(mfrow=c(1,2))

plot(interp_area,add=TRUE)

title("Mean for January 1")

plot(interp_area,add=TRUE)

title("Mean for month of January")

dev.off()

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

################### Figure 4. MAE/RMSE and distance to closest fitting station.

#plot.gstatVariogram

#Plot a sample variogram, and possibly a fitted model

#model 1 lat,lon and elev

date_selected <- c("20100101","20100901")

png(paste("Figure9_paper_","_variogram_",date_selected[1],"_",date_selected[2],"_",out_prefix,".png", sep=""),

    #height=480*6,width=480*4)

#p3 <- list_plots_spt[[3]]

#plot(p1)

#plot(p241)

#grid.arrange(p1,p2,p3,ncol=1)

dev.off()

#Combine both plot?     + plot info on sill, nugget and range? and most frequent model selected

tb_variogram["date"] <- dates

tb_variogram$month <- add_month_tag(tb_variogram)

#add dates and month??

#hist(tt1$range)

#hist(tt1$psill)

#hist(tt1$Nug)

boxplot(tb_variogram$Nug~tb_variogram$month,outline=F)

boxplot(tb_variogram$psill~tb_variogram$month,outline=F)

#boxplot(tb_variogram$model~tb_variogram$month,outline=F)

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

### Figure 10: map of residuals...for a specific date...

index <- 244

#date_selected <-  c("20100101")

date_selected <-  c("20100901")

#data_s<-validation_mod_obj[[index]]$data_s

date_proc<-strptime(date_selected, "%Y%m%d")   # interpolation date being processed

mo<-as.integer(strftime(date_proc, "%m"))          # current month of the date being processed

year<-as.integer(strftime(date_proc, "%Y"))

datelabel=format(ISOdate(year,mo,day),"%b %d, %Y")

list_p <- vector("list", length(names_mod))

for (k in 1:length(names_mod)){

  model_name <- names_mod[k]

  #p1 <- levelplot(elev,scales = list(draw = FALSE), colorkey = FALSE,col.regions=rev(terrain.colors(255)),contour=T)

  #add legend..par.settings = GrTheme

  cx <- ((data_v[[res_model_name]])*2)

  #p2 <- spplot(data_v,res_model_name, cex=1 * cx,main=paste("Residuals for ",res_model_name," ",datelabel,sep=""),

  #             col.regions=matlab.like(25))

  p2 <- bubble(data_v,res_model_name, main=paste("Residuals for ",res_model_name," ",datelabel,sep=""),

  p3 <- p2 + p1 + p2 #to force legend...

  #p1 <- spplot(interp_area)

  #p3 <- p1+p2 #to force legend...

  list_p[[k]] <- p3

}

    height=480*layout_m[1],width=480*layout_m[2])

dev.off()   

######## NOW GET A ACCURACY BY STATIONS

#number of observations per day

year_nbv <- sapply(list_data_v,FUN=length)

#Convert sp data.frame and combined them in one unique df, see function define earlier

data_v_combined <-convert_spdf_to_df_from_list(list_data_v) #long rownames

t<-melt(data_v_combined,

        measure=names_var, 

mae_tb<-cast(t,id~variable,mae_fun) #join to station location...

met_obj <-load_obj(file.path(in_dir1,met_obj_file_1))

stat_loc<-readOGR(dsn=in_dir1,layer=sub(".shp","",basename(met_obj$loc_stations)))

coordinates(data_v_mae)<-coords                      #Assign coordinates to the data frame

proj4string(data_v_mae)<- proj4string(stat_loc)                #Assign coordinates reference system in PROJ4 format

data_v_mae<-spTransform(data_v_mae,CRS(CRS_interp))     #Project from WGS84 to new coord. system

infile_reg_outline <- "/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/OR83M_state_outline.shp"  #input region outline defined by polygon: Oregon

reg_outline <- readOGR(dsn=dirname(infile_reg_outline),layer=sub(".shp","",basename(infile_reg_outline)))

p + layer(sp.polygons(reg_outline,lwd=0.9,col="darkgray"))

p4 + layer(sp.polygons(reg_outline,lwd=0.9,col="darkgray"))

col_t <- colorRampPalette(c('blue', 'white', 'red'))

plot(res_mod1~LST,data=test)                    

plot(res_mod5~LST,data=test)                    

brks<-c(0,500,1000,1500,2000,2500,4000)

lab_brks<-1:6

elev_rcstat<-cut(data_v_ag$elev,breaks=brks,labels=lab_brks,right=F)

boxplot(data_v_ag$res_mod2~elev_rcstat,ylim=c(-15,15))           

boxplot(data_v_ag$res_mod1~elev_rcstat,ylim=c(-15,15))                      

boxplot(data_v_ag$res_mod5~elev_rcstat,ylim=c(-15,15))                      

brks<-c(-20,-10,0,10,20,30,40)

lab_brks<-1:6

rcstat<-cut(data_v_ag$LST,breaks=brks,labels=lab_brks,right=F)