#Part 1: Script produces plots for every selected date

#Part 2: Examine 

#AUTHOR: Benoit Parmentier                                                                       

#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#???--   

out_path<- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/output_data/"

infile_covar<-"covariates__venezuela_region__VE_01292013.tif" #this is an output from covariate script

date_selected<-c("20000101") ##This is for year 2000!!!

#out_prefix<-"_365d_GAM_fus5_all_lstd_03132013"

#out_prefix<-"_365d_GAM_fus5_all_lstd_03142013"                #User defined output prefix

#gam_fus_mod<-load_obj("gam_fus_mod_365d_GAM_fus5_all_lstd_02202013.RData")

#validation_obj<-load_obj("gam_fus_validation_mod_365d_GAM_fus5_all_lstd_02202013.RData")

#clim_obj<-load_obj("gamclim_fus_mod_365d_GAM_fus5_all_lstd_02202013.RData")

             "nobs_09","nobs_10","nobs_11","nobs_12")

covar_names<-c(rnames,lc_names,lst_names)

                     "infile_covar","covar_names","date_selected","var","out_prefix")

setwd(in_path)

### PLOTTING RESULTS FROM VENEZUELA INTERPOLATION FOR ANALYSIS

#j=1

plots_assessment_by_date<-function(j,list_param){

  date_selected<-list_param$date_selected

  #validation_obj<-load_obj("gam_fus_validation_mod_365d_GAM_fus5_all_lstd_02202013.RData")

  #clim_obj<-load_obj("gamclim_fus_mod_365d_GAM_fus5_all_lstd_02202013.RData")

  method_mod_obj<-raster_prediction_obj$gam_fus_mod #change later for any model type

  #validation_obj<-raster_prediction_obj$validation_obj

  validation_obj<-raster_prediction_obj$gam_fus_validation_mod #change later for any model type

  clim_obj<-raster_prediction_obj$gamclim_fus_mod #change later for any model type

  if (var=="TMAX"){

  }

  if (var=="TMIN"){

  }

  ## Read covariate stack...

  LC_mask_rec[is.na(LC_mask_rec)]<-0

  #determine index position matching date selected

  for (j in 1:length(date_selected)){

    for (i in 1:length(method_mod_obj)){

      if(method_mod_obj[[i]]$sampling_dat$date==date_selected[j]){  

  #Get raster stack of interpolated surfaces

  index<-i_dates[[j]]

  rast_pred_temp<-stack(pred_temp) #stack of temperature predictions from models 

  #Get validation metrics, daily spdf training and testing stations, monthly spdf station input

  rmse<-metrics_v$rmse[nrow(metrics_v)]

  rmse_f<-metrics_s$rmse[nrow(metrics_s)]  

            out_prefix,".png", sep=""))

  abline(0,1)

  mean_bias<-paste("Mean LST bias= ",format(mean(data_month$LSTD_bias,na.rm=TRUE),digits=3),sep="")

  #Add the number of data points on the plot

  legend("topleft",legend=c(mean_bias,nb_point),bty="n")

  ## Figure 2: Daily_tmax_monthly_TMax_scatterplot, modify for TMin!!

            out_prefix,".png", sep=""))

  #Add the number of data points on the plot

  legend("topleft",legend=c(nb_point,nb_point2,nb_point3),bty="n",cex=0.8)

  dev.off()

  ## Figure 3: Predicted_tmax_versus_observed_scatterplot 

  #This is for mod_kr!! add other models later...

  x_range<-range(c(data_s[[y_var_name]],data_v[[y_var_name]]),na.rm=T)

  col_t<- c("black","red")

  pch_t<- c(1,2)

       xlab=paste("Actual daily for",datelabel),ylab="Pred daily", 

       ylim=y_range,xlim=x_range,col=col_t[1],pch=pch_t[1])

  grid(lwd=0.5, col="black")

  abline(0,1)

  legend("topleft",legend=c("training","testing"),pch=pch_t,col=col_t,bty="n",cex=0.8)

  rmse_str1<-paste("RMSE= ",format(rmse,digits=3),sep="")

  rmse_str2<-paste("RMSE_f= ",format(rmse_f,digits=3),sep="")

  #Add the number of data points on the plot

  dev.off()

            out_prefix,".png", sep=""))

  #paste(metrics_v$pred_mod,format(metrics_v$rmse,digits=3),sep=":")

  names(rast_pred_temp)<-paste(metrics_v$pred_mod,format(metrics_v$rmse,digits=3),sep=":")

  levelplot(rast_pred_temp)

  dev.off()

  png(paste("Raster_prediction_plot",sampling_dat$date,"_",sampling_dat$prop,"_",sampling_dat$run_samp,

            out_prefix,".png", sep=""))

  #paste(metrics_v$pred_mod,format(metrics_v$rmse,digits=3),sep=":")

  plot(rast_pred_temp)

  dev.off()

            out_prefix,".png", sep=""))

  plot(raster(rast_pred_temp,layer=5))

  plot(data_s,col="black",cex=1.2,pch=2,add=TRUE)

         pch=c(2,1),bty="n")

  dev.off()

            out_prefix,".png", sep=""))

  plot(raster(rast_pred_temp,layer=5))

  plot(data_month,col="black",cex=1.2,pch=4,add=TRUE)

  title("Monthly ghcn station in Venezuela for January")

  dev.off()

            "_",sampling_dat$run_samp[i],out_prefix,".png", sep=""))

  bias_rast<-stack(clim_obj[[index]]$bias)