#Part 1: Script produces plots for every selected date

#Part 2: Examine 

#AUTHOR: Benoit Parmentier                                                                       

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

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

### Parameters and arguments

##Paths to inputs and output

#Select relevant dates and load R objects created during the interpolation step

#names(list_param_results_analyses)<-c("in_path","out_path","script_path","raster_prediction_obj", "interpolation_method",

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

plots_assessment_by_date<-function(j,list_param){

  interpolation_method <- list_param$interpolation_method

  raster_prediction_obj<-list_param$raster_prediction_obj

  method_mod_obj<-raster_prediction_obj$method_mod_obj

    y_var_month <-"TMin"

  }

  names(s_raster)<-covar_names               #Assigning names to the raster layers: making sure it is included in the extraction

  ## Prepare study area  mask: based on LC12 (water)

  #Get raster stack of interpolated surfaces

  index<-i_dates[[j]]

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

  sampling_dat<-method_mod_obj[[index]]$sampling_dat

  #The names of covariates can be changed...

  LST_month<-paste("mm_",month,sep="") # name of LST month to be matched

  s_raster<-dropLayer(s_raster,pos)      # If it exists drop layer

  r1<-raster(s_raster,layer=pos)             #Select layer from stack

  #Get mask image!!

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

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

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

  plot(data_month[[y_var_month]],data_month$LST,xlab=paste("Station mo ",y_var_month,sep=""),ylab=paste("LST mo ",y_var_month,sep=""))

  title(paste("LST vs ", y_var_month,"for",datelabel,sep=" "))

  abline(0,1)

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

  plot(data_s[[y_var_name]]~data_s[[y_var_month]],xlab=paste("Month") ,ylab=paste("Daily for",datelabel),main="across stations in VE")

  nb_point<-paste("ns=",length(data_s[[y_var_month]]),sep="")

  nb_point2<-paste("ns_obs=",length(data_s[[y_var_month]])-sum(is.na(data_s[[y_var_name]])),sep="")

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

  model_name<-"mod_kr" #can be looped through models later on...

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

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

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

  dev.off()

  ## Figure 4a: prediction raster images

  #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=":")

  #plot(rast_pred_temp)

  dev.off()

  ## Figure 4b: prediction raster images

  #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=":")

  plot(rast_pred_temp)

  dev.off()

  ## Figure 5: training and testing stations used

  plot(raster(rast_pred_temp,layer=5))

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

  plot(data_v,col="red",cex=1.2,pch=1,add=TRUE)

  ## Figure 6: monthly stations used

  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")

  ## Figure 7: delta surface and bias

  if (interpolation_method=="gam_fus"){

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

    delta_rast<-raster(method_mod_obj[[index]]$delta) #only one delta image!!!

    names(delta_rast)<-"delta"

    rast_temp_date<-stack(bias_rast,delta_rast)

    #bias_d_rast<-raster("fusion_bias_LST_20100103_30_1_10d_GAM_fus5_all_lstd_02082013.rst")

    plot(rast_temp_date)

    dev.off()

  }

  if (interpolation_method=="gam_CAI"){

    clim_rast<-stack(clim_method_mod_obj[[index]]$clim)

    rast_temp_date<-stack(clim_rast,delta_rast)

    #bias_d_rast<-raster("fusion_bias_LST_20100103_30_1_10d_GAM_fus5_all_lstd_02082013.rst")

    plot(rast_temp_date)

  #histogram(rast_pred_temp)

  list_output_analyses<-list(metrics_s,metrics_v)

  return(list_output_analyses)

}