####################################  Cimate Layers Predctions  #######################################

############################  Script for assessment of scaling up on NEX: part 1a ##############################

#MODIFIED ON: 12/29/2015            

# - Separate call in a master script for assessment

# - add second stage in the master script for assessment

# - add mosaicing in the master script for assessment

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

run_assessment_prediction_fun <-function(i,list_param_run_assessment_prediction){

  ##Function to predict temperature interpolation with 12 input parameters

  #12 parameters used in the data preparation stage and input in the current script

  #

  #1) in_dir1 : llcoactation of rroot directory containging tiles

  #2) region_names : region_names #e.g. c("reg23","reg4") 

  #3) y_var_name : list_param_run_assessment_prediction$y_var_name # e.g. dailyTmax" #PARAM3

  #4) interpolation_method :  e.g. #c("gam_CAI") #PARAM4

  #5) out_prefix : e.g. "run_global_analyses_pred_12282015" #PARAM5

  #6) out_dir <- list_param_run_assessment_prediction$out_dir #<- "/nobackupp8/bparmen1/" #PARAM6

  #7) create_out_dir_param: if true a new dirrectory  is craeted for the outputs 

  #8) CRS_locs_WGS84: coordinates  CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84, #PARAM8

  #9) list_year_predicted : list of years predicted e.g. 1984:2004

  #10) file_format : format for mosaiced files #PARAM10

  #11) NA_flag_val : No data value, #PARAM11

  #12) num_cores : number of cores used #PARAM13

  ###Loading R library and packages     

  ### Loading R library and packages        

  #library used in the workflow production:

  library(gtools)                              # loading some useful tools 

  library(mgcv)                                # GAM package by Simon Wood

  library(sp)                                  # Spatial pacakge with class definition by Bivand et al.

  library(spdep)                               # Spatial pacakge with methods and spatial stat. by Bivand et al.

  library(rgdal)                               # GDAL wrapper for R, spatial utilities

  library(gstat)                               # Kriging and co-kriging by Pebesma et al.

  library(fields)                              # NCAR Spatial Interpolation methods such as kriging, splines

  library(raster)                              # Hijmans et al. package for raster processing

  library(gdata)                               # various tools with xls reading, cbindX

  library(rasterVis)                           # Raster plotting functions

  library(parallel)                            # Parallelization of processes with multiple cores

  library(maptools)                            # Tools and functions for sp and other spatial objects e.g. spCbind

  library(maps)                                # Tools and data for spatial/geographic objects

  library(reshape)                             # Change shape of object, summarize results 

  library(plotrix)                             # Additional plotting functions

  library(plyr)                                # Various tools including rbind.fill

  library(spgwr)                               # GWR method

  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

  library(gridExtra)

  #Additional libraries not used in workflow

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

  library(colorRamps)

  #### FUNCTION USED IN SCRIPT

  #This is read in in the master script

  #function_analyses_paper1 <- "global_run_scalingup_assessment_part1_functions_02112015.R" #PARAM12

  #script_path <- "/nobackupp8/bparmen1/env_layers_scripts" #path to script

  #source(file.path(script_path,function_analyses_paper1)) #source all functions used in this script 

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

  #### Parameters and constants  

  #Make this a function

  #reg1 (North Am), reg2(Europe),reg3(Asia), reg4 (South Am), reg5 (Africa), reg6 (Australia-Asia)

  #master directory containing the definition of tile size and tiles predicted

  in_dir1 <- list_param_run_assessment_prediction$in_dir1 

  region_names <- list_param_run_assessment_prediction$region_names #e.g. c("reg23","reg4") 

  y_var_name <- list_param_run_assessment_prediction$y_var_name # e.g. dailyTmax" #PARAM3

  interpolation_method <- list_param_run_assessment_prediction$interpolation_method #c("gam_CAI") #PARAM4

  out_prefix <- list_param_run_assessment_prediction$out_prefix #"run_global_analyses_pred_12282015" #PARAM5

  out_dir <- list_param_run_assessment_prediction$out_dir #<- "/nobackupp8/bparmen1/" #PARAM6

  create_out_dir_param <-list_param_run_assessment_prediction$create_out_dir_param #<- TRUE #PARAM7

  CRS_locs_WGS84 <- list_param_run_assessment_prediction$CRS_locs_WGS84 #<- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84, #PARAM8

  list_year_predicted <- list_param_run_assessment_prediction$list_year_predicted #<- 1984:2004

  #list_param_run_assessment_prediction$year_predicted <- list_year_predicted[1]

  file_format <- list_param_run_assessment_prediction$file_format #<- ".tif" #format for mosaiced files #PARAM10

  NA_flag_val <- list_param_run_assessment_prediction$NA_flag_val #<- -9999  #No data value, #PARAM11

  num_cores <- list_param_run_assessment_prediction$num_cores #<- 6 #number of cores used #PARAM13

  ########################## START SCRIPT #########################################

  #Need to make this a function to run as a job...

  ######################## PART0: Read content of predictions first.... #####

  #function looped over i, correspoding to year predicted

  year_predicted <- list_param_run_assessment_prediction$list_year_predicted[i] 

  list_outfiles <- vector("list", length=6) #collect names of output files

  in_dir_list <- list.dirs(path=in_dir1,recursive=FALSE) #get the list regions processed for this run

  #basename(in_dir_list)

#                       y=in_dir_list) 

  in_dir_list_all  <- unlist(lapply(in_dir_list,function(x){list.dirs(path=x,recursive=F)}))

  in_dir_list <- in_dir_list_all

  #in_dir_list <- in_dir_list[grep("bak",basename(basename(in_dir_list)),invert=TRUE)] #the first one is the in_dir1

  #this was changed on 10052015 because the shapefiles were not matching!!!

  #in_dir_list_tmp <- list.dirs(path=in_dir1,recursive=FALSE) #get the list regions processed for this run

  #in_dir_subset <- in_dir_list_tmp[grep("subset",basename(in_dir_list_tmp),invert=FALSE)] #select directory with shapefiles...

  #in_dir_shp <- file.path(in_dir_subset,"shapefiles")

  #in_dir_list_tmp <- list.dirs(path=in_dir1,recursive=FALSE) #get the list regions processed for this run

  in_dir_subset <- in_dir_list_all[grep("subset",basename(in_dir_list_all),invert=FALSE)] #select directory with shapefiles...

  in_dir_shp <- file.path(in_dir_subset,"shapefiles")

  #select only directories used for predictions

  in_dir_reg <- in_dir_list[grep(".*._.*.",basename(in_dir_list),invert=FALSE)] #select directory with shapefiles...

  #in_dir_reg <- in_dir_list[grep("july_tiffs",basename(in_dir_reg),invert=TRUE)] #select directory with shapefiles...

  in_dir_list <- in_dir_reg

  in_dir_list <- in_dir_list[grep("bak",basename(basename(in_dir_list)),invert=TRUE)] #the first one is the in_dir1

  #list of shapefiles used to define tiles

  in_dir_shp_list <- list.files(in_dir_shp,".shp",full.names=T)

  ## load problematic tiles or additional runs

  #modify later...

  #system("ls /nobackup/bparmen1")

  if(create_out_dir_param==TRUE){

    out_dir <- create_dir_fun(out_dir,out_prefix)

    setwd(out_dir)

  }else{

    setwd(out_dir) #use previoulsy defined directory

  }

  ##raster_prediction object : contains testing and training stations with RMSE and model object

  in_dir_list_tmp <- file.path(in_dir_list,year_predicted)

  list_raster_obj_files <- lapply(in_dir_list_tmp,FUN=function(x){list.files(path=x,pattern="^raster_prediction_obj.*.RData",full.names=T)})

  basename(dirname(list_raster_obj_files[[1]]))

  list_names_tile_coord <- lapply(list_raster_obj_files,FUN=function(x){basename(dirname(x))})

  list_names_tile_id <- paste("tile",1:length(list_raster_obj_files),sep="_")

  names(list_raster_obj_files)<- list_names_tile_id

  #one level up

  lf_covar_obj <- lapply(in_dir_list,FUN=function(x){list.files(path=x,pattern="covar_obj.*.RData",full.names=T)})

  lf_covar_tif <- lapply(in_dir_list,FUN=function(x){list.files(path=x,pattern="covar.*.tif",full.names=T)})

  #sub_sampling_obj_daily_gam_CAI_10.0_-75.0.RData

  #sub_sampling_obj_gam_CAI_10.0_-75.0.RData

  lf_sub_sampling_obj_files <- lapply(in_dir_list,FUN=function(x){list.files(path=x,pattern=paste("^sub_sampling_obj_",interpolation_method,".*.RData",sep=""),full.names=T)})

  lf_sub_sampling_obj_daily_files <- lapply(in_dir_list_tmp,FUN=function(x){list.files(path=x,pattern="^sub_sampling_obj_daily.*.RData",full.names=T)})

  ## This will be part of the raster_obj function

  #debug(create_raster_prediction_obj)

  #out_prefix_str <- paste(basename(in_dir_list),out_prefix,sep="_") 

  #lf_raster_obj <- create_raster_prediction_obj(in_dir_list,interpolation_method, y_var_name,out_prefix_str,out_path_list=NULL)

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

  ######## PART 1: Generate tables to collect information:

  ######## over all tiles in North America 

  ##Function to collect all the tables from tiles into a table

  ###Table 1: Average accuracy metrics

  ###Table 2: daily accuracy metrics for all tiles

  #First create table of tiles under analysis and their coord

  df_tile_processed <- data.frame(tile_coord=basename(in_dir_list))

  df_tile_processed$tile_id <- unlist(list_names_tile_id) #Arbitrary tiling number!!

  df_tile_processed$path_NEX <- in_dir_list

  df_tile_processed$year_predicted <- year_predicted

  df_tile_processed$sub_sampling_clim  <- lf_sub_sampling_obj_files

  df_tile_processed$sub_sampling_daily  <- lf_sub_sampling_obj_daily_files

  #lf_sub_sampling_obj_files

  ##Quick exploration of raster object

  #Should be commented out to make this a function

  #robj1 <- try(load_obj(list_raster_obj_files[[3]])) #This is an example tile

  #robj1 <- load_obj(lf_raster_obj[4]) #This is tile tile

  #names(robj1)

  #names(robj1$method_mod_obj[[2]]) #for January 1, 2010

  #names(robj1$method_mod_obj[[2]]$dailyTmax) #for January

  #names(robj1$method_mod_obj[[11]]) #for January 1, 2010

  #names(robj1$method_mod_obj[[11]]$dailyTmax) #for January

  #names(robj1$clim_method_mod_obj[[1]]$data_month) #for January

  #names(robj1$validation_mod_month_obj[[1]]$data_s) #for January with predictions

  #Get the number of models predicted

  #nb_mod <- length(unique(robj1$tb_diagnostic_v$pred_mod))#

  #list_formulas <- (robj1$clim_method_mod_obj[[1]]$formulas)

  #dates_predicted <- (unique(robj1$tb_diagnostic_v$date))

  #list_tb_diagnostic_v <- mclapply(lf_validation_obj,FUN=function(x){try( x<- load_obj(x)); try(extract_from_list_obj(x,"metrics_v"))},mc.preschedule=FALSE,mc.cores = 6)                           

  #names(list_tb_diagnostic_v) <- list_names_tile_id

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

  #### Table 1: Average accuracy metrics per tile and predictions

  #can use a maximum of 6 cores on the NEX Bridge

  #For 28 tiles and  28 RData boject it takes 15-16 min

  #summary_metrics_v_list <- mclapply(list_raster_obj_files[5:6],FUN=function(x){try( x<- load_obj(x)); try(x[["summary_metrics_v"]]$avg)},mc.preschedule=FALSE,mc.cores = 2)                           

  summary_metrics_v_list <- mclapply(list_raster_obj_files,FUN=function(x){try( x<- load_obj(x)); try(x[["summary_metrics_v"]]$avg)},mc.preschedule=FALSE,mc.cores = num_cores)                         

  #summary_metrics_v_list <- lapply(summary_metrics_v_list,FUN=function(x){try(x$avg)})

  names(summary_metrics_v_list) <- list_names_tile_id

  summary_metrics_v_tmp <- remove_from_list_fun(summary_metrics_v_list)$list

  df_tile_processed$metrics_v <- as.integer(remove_from_list_fun(summary_metrics_v_list)$valid)

  #Now remove "try-error" from list of accuracy)

  summary_metrics_v_NA <- do.call(rbind.fill,summary_metrics_v_tmp) #create a df for NA tiles with all accuracy metrics

  #tile_coord <- lapply(1:length(summary_metrics_v_list),FUN=function(i,x){rep(names(x)[i],nrow(x[[i]]))},x=summary_metrics_v_list)

  #add the tile id identifier

  tile_id_tmp <- lapply(1:length(summary_metrics_v_tmp),

                        FUN=function(i,x,y){rep(y[i],nrow(x[[i]]))},x=summary_metrics_v_tmp,y=names(summary_metrics_v_tmp))

  #adding tile id summary data.frame

  summary_metrics_v_NA$tile_id <-unlist(tile_id_tmp)

  summary_metrics_v_NA$n <- as.integer(summary_metrics_v_NA$n)

  summary_metrics_v_NA <- merge(summary_metrics_v_NA,df_tile_processed[,1:2],by="tile_id")

  tx<-strsplit(as.character(summary_metrics_v_NA$tile_coord),"_")

  lat<- as.numeric(lapply(1:length(tx),function(i,x){x[[i]][1]},x=tx))

  long<- as.numeric(lapply(1:length(tx),function(i,x){x[[i]][2]},x=tx))

  summary_metrics_v_NA$lat <- lat

  summary_metrics_v_NA$lon <- long

  list_out_files

  write.table(as.data.frame(summary_metrics_v_NA),

              file=file.path(out_dir,paste("summary_metrics_v2_NA_",year_predicted,"_",out_prefix,".txt",sep="")),sep=",")

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

  ###Table 2: daily validation/testing accuracy metrics for all tiles

  #this takes about 15min for 28 tiles (reg4)

  #tb_diagnostic_v_list <- lapply(list_raster_obj_files,FUN=function(x){x<-load_obj(x);x[["tb_diagnostic_v"]]})                           

  tb_diagnostic_v_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(x[["tb_diagnostic_v"]])},mc.preschedule=FALSE,mc.cores = num_cores)                           

  names(tb_diagnostic_v_list) <- list_names_tile_id

  tb_diagnostic_v_tmp <- remove_from_list_fun(tb_diagnostic_v_list)$list

  #df_tile_processed$tb_diag <- remove_from_list_fun(tb_diagnostic_v_list)$valid

  tb_diagnostic_v_NA <- do.call(rbind.fill,tb_diagnostic_v_tmp) #create a df for NA tiles with all accuracy metrics

  tile_id_tmp <- lapply(1:length(tb_diagnostic_v_tmp),

                        FUN=function(i,x,y){rep(y[i],nrow(x[[i]]))},x=tb_diagnostic_v_tmp,y=names(tb_diagnostic_v_tmp))

  tb_diagnostic_v_NA$tile_id <- unlist(tile_id_tmp) #adding identifier for tile

  tb_diagnostic_v_NA <- merge(tb_diagnostic_v_NA,df_tile_processed[,1:2],by="tile_id")

  write.table((tb_diagnostic_v_NA),

              file=file.path(out_dir,paste("tb_diagnostic_v_NA_",year_predicted,"_",out_prefix,".txt",sep="")),sep=",")

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

  ###Table 3: monthly fit/training accuracy information for all tiles

  ## Monthly fitting information

  tb_month_diagnostic_s_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(x[["tb_month_diagnostic_s"]])},mc.preschedule=FALSE,mc.cores = num_cores)                           

  names(tb_month_diagnostic_s_list) <- list_names_tile_id

  tb_month_diagnostic_s_tmp <- remove_from_list_fun(tb_month_diagnostic_s_list)$list

  #df_tile_processed$tb_diag <- remove_from_list_fun(tb_diagnostic_v_list)$valid

  tb_month_diagnostic_s_NA <- do.call(rbind.fill,tb_month_diagnostic_s_tmp) #create a df for NA tiles with all accuracy metrics

  tile_id_tmp <- lapply(1:length(tb_month_diagnostic_s_tmp),

                        FUN=function(i,x,y){rep(y[i],nrow(x[[i]]))},x=tb_month_diagnostic_s_tmp,y=names(tb_month_diagnostic_s_tmp))

  tb_month_diagnostic_s_NA$tile_id <- unlist(tile_id_tmp) #adding identifier for tile

  tb_month_diagnostic_s_NA <- merge(tb_month_diagnostic_s_NA,df_tile_processed[,1:2],by="tile_id")

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

  tb_month_diagnostic_s_NA$month<-strftime(date_f, "%m")          # current month of the date being processed

  write.table((tb_month_diagnostic_s_NA),

              file=file.path(out_dir,paste("tb_month_diagnostic_s_NA_",year_predicted,"_",out_prefix,".txt",sep="")),sep=",")

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

  ###Table 4: daily fit/training accuracy information with predictions for all tiles

  ## daily fit info:

  tb_diagnostic_s_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(x[["tb_diagnostic_s"]])},mc.preschedule=FALSE,mc.cores = num_cores)                           

  names(tb_diagnostic_s_list) <- list_names_tile_id

  tb_diagnostic_s_tmp <- remove_from_list_fun(tb_diagnostic_s_list)$list

  #df_tile_processed$tb_diag <- remove_from_list_fun(tb_diagnostic_v_list)$valid

  tb_diagnostic_s_NA <- do.call(rbind.fill,tb_diagnostic_s_tmp) #create a df for NA tiles with all accuracy metrics

  tile_id_tmp <- lapply(1:length(tb_diagnostic_s_tmp),

                        FUN=function(i,x,y){rep(y[i],nrow(x[[i]]))},x=tb_diagnostic_s_tmp,y=names(tb_diagnostic_s_tmp))

  tb_diagnostic_s_NA$tile_id <- unlist(tile_id_tmp) #adding identifier for tile

  tb_diagnostic_s_NA <- merge(tb_diagnostic_s_NA,df_tile_processed[,1:2],by="tile_id")

  write.table((tb_diagnostic_s_NA),

              file=file.path(out_dir,paste("tb_diagnostic_s_NA_",year_predicted,"_",out_prefix,".txt",sep="")),sep=",")

  ##### Table 5: Add later on: daily info

  ### with also data_s and data_v saved!!!

  #Insert here...compute input and predicted ranges to spot potential errors?

  ### Make this part a function...this is repetitive

  ##### SPDF of Monhtly Station info

  #load data_month for specific tiles

  #10.45pm

  #data_month <- extract_from_list_obj(robj1$clim_method_mod_obj,"data_month")

  #names(data_month) #this contains LST means (mm_1, mm_2 etc.) as well as TMax and other info

  #data_month_s_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(x$validation_mod_month_obj[["data_s"]])},mc.preschedule=FALSE,mc.cores = 6)                           

  data_month_s_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(extract_from_list_obj(x$validation_mod_month_obj,"data_s"))},mc.preschedule=FALSE,mc.cores = 6)                           

  #test <- mclapply(list_raster_obj_files[1:6],FUN=function(x){try(x<-load_obj(x));try(extract_from_list_obj(x$validation_mod_month_obj,"data_s"))},mc.preschedule=FALSE,mc.cores = 6)                           

  names(data_month_s_list) <- list_names_tile_id

  data_month_tmp <- remove_from_list_fun(data_month_s_list)$list

  #df_tile_processed$metrics_v <- remove_from_list_fun(data_month_s_list)$valid

  tile_id <- lapply(1:length(data_month_tmp),

                    FUN=function(i,x){rep(names(x)[i],nrow(x[[i]]))},x=data_month_tmp)

  data_month_NAM <- do.call(rbind.fill,data_month_tmp) #combined data_month for "NAM" North America

  data_month_NAM$tile_id <- unlist(tile_id)

  write.table((data_month_NAM),

              file=file.path(out_dir,paste("data_month_s_NAM_",year_predicted,"_",out_prefix,".txt",sep="")),sep=",")

  #Get validation data?? Find other object from within the dir 

  #Som region don't have validation data at monthly time scale.

  #### SPDF of daily Station info

  #load data_month for specific tiles

  #data_month <- extract_from_list_obj(robj1$clim_method_mod_obj,"data_month")

  #names(data_month) #this contains LST means (mm_1, mm_2 etc.) as well as TMax and other info

  data_day_s_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(extract_from_list_obj(x$validation_mod_obj,"data_s"))},mc.preschedule=FALSE,mc.cores = num_cores)    

  data_day_v_list <- mclapply(list_raster_obj_files,FUN=function(x){try(x<-load_obj(x));try(extract_from_list_obj(x$validation_mod_obj,"data_v"))},mc.preschedule=FALSE,mc.cores = num_cores)    

  names(data_day_s_list) <- list_names_tile_id

  names(data_day_v_list) <- list_names_tile_id

  data_day_s_tmp <- remove_from_list_fun(data_day_s_list)$list

  data_day_v_tmp <- remove_from_list_fun(data_day_v_list)$list

  #df_tile_processed$metrics_v <- remove_from_list_fun(data_month_s_list)$valid

  tile_id <- lapply(1:length(data_day_s_tmp),

                    FUN=function(i,x){rep(names(x)[i],nrow(x[[i]]))},x=data_day_s_tmp)

  data_day_s_NAM <- do.call(rbind.fill,data_day_s_tmp) #combined data_month for "NAM" North America

  data_day_s_NAM$tile_id <- unlist(tile_id)

  tile_id <- lapply(1:length(data_day_v_tmp),

                    FUN=function(i,x){rep(names(x)[i],nrow(x[[i]]))},x=data_day_v_tmp)

  data_day_v_NAM <- do.call(rbind.fill,data_day_v_tmp) #combined data_month for "NAM" North America

  data_day_v_NAM$tile_id <- unlist(tile_id)

  write.table((data_day_s_NAM),

              file=file.path(out_dir,paste("data_day_s_NAM_",year_predicted,"_",out_prefix,".txt",sep="")),sep=",")

  write.table((data_day_v_NAM),

              file=file.path(out_dir,paste("data_day_v_NAM_",year_predicted,"_",out_prefix,".txt",sep="")),sep=",")

  #### Recover subsampling data

  #For tiles with many stations, there is a subsampling done in terms of distance (spatial pruning) and 

  #in terms of station numbers if there are still too many stations to consider. This is done at the 

  #daily and monthly stages.

  #lf_sub_sampling_obj_files <- lapply(in_dir_list,FUN=function(x){list.files(path=x,pattern=paste("^sub_sampling_obj_",interpolation_method,".*.RData",sep=""),full.names=T)})

  #lf_sub_sampling_obj_daily_files <- lapply(in_dir_list,FUN=function(x){list.files(path=x,pattern="^sub_sampling_obj_daily.*.RData",full.names=T)})

  #sub_sampling_obj <- try(load_obj(lf_sub_sampling_obj_files[[3]])) #This is an example tile

  #data_removed contains the validation data...

  #this data can be used for validation of the product. Note that it may be missing for some tiles

  #as no stations are removed if there are too enough stations in the tile

  #this will need to be checked later on...

  data_month_v_subsampling_list <- mclapply(lf_sub_sampling_obj_files,FUN=function(x){try(x<-load_obj(x));try(extract_from_list_obj(x$validation_mod_month_obj,"data_removed"))},mc.preschedule=FALSE,mc.cores = 6)                           

  #test <- mclapply(list_raster_obj_files[1:6],FUN=function(x){try(x<-load_obj(x));try(extract_from_list_obj(x$validation_mod_month_obj,"data_s"))},mc.preschedule=FALSE,mc.cores = 6)                           

  names(data_month_v_subsampling_list) <- list_names_tile_id

  data_month_v_subsampling_tmp <- remove_from_list_fun(data_month_v_subsampling_list)$list

  #df_tile_processed$metrics_v <- remove_from_list_fun(data_month_s_list)$valid

  #if no stations have been removed then there are no validation stations !!!

  if(length(data_month_v_subsampling_tmp)!=0){

    tile_id <- lapply(1:length(data_month_v_subsampling_tmp),

                      FUN=function(i,x){try(rep(names(x)[i],nrow(x[[i]])))},x=data_month_v_subsampling_tmp)

    data_month_v_subsmapling_NAM <- do.call(rbind.fill,ddata_month_v_subsampling_tmp) #combined data_month for "NAM" North America

    data_month_v_subsampling_NAM$tile_id <- unlist(tile_id)

    write.table((data_month_v_subsampling_NAM),

                file=file.path(out_dir,paste("data_month_v_subsampling_NAM_",year_predicted,"_",out_prefix,".txt",sep="")),sep=",")

  }

  ## Do the same for daily...

  ## End of potential function started in line 317...this section will be cut down for simplicity

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

  ####### PART 3: EXAMINE STATIONS AND MODEL FITTING ###

  ### Stations and model fitting ###

  #summarize location and number of training and testing used by tiles

  #names(robj1$clim_method_mod_obj[[1]]$data_month) # monthly data for January

  #names(robj1$validation_mod_month_obj[[1]]$data_s) # daily for January with predictions

  #note that there is no holdout in the current run at the monthly time scale:

  #robj1$clim_method_mod_obj[[1]]$data_month_v #zero rows for testing stations at monthly timescale

  #load data_month for specific tiles

  data_month <- extract_from_list_obj(robj1$clim_method_mod_obj,"data_month")

  #names(data_month) #this contains LST means (mm_1, mm_2 etc.) as well as TMax and other info

  use_day=TRUE

  use_month=TRUE

  #list_raster_obj_files <- c("/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/output10Deg/reg1//30.0_-100.0/raster_prediction_obj_gam_CAI_dailyTmax30.0_-100.0.RData",

  #                    "/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/output10Deg/reg1//30.0_-105.0/raster_prediction_obj_gam_CAI_dailyTmax30.0_-105.0.RData")

  list_names_tile_id <- df_tile_processed$tile_id

  list_raster_obj_files[list_names_tile_id]

  #list_names_tile_id <- c("tile_1","tile_2")

  list_param_training_testing_info <- list(list_raster_obj_files[list_names_tile_id],use_month,use_day,list_names_tile_id)

  names(list_param_training_testing_info) <- c("list_raster_obj_files","use_month","use_day","list_names_tile_id")

  list_param <- list_param_training_testing_info

  #debug(extract_daily_training_testing_info)

  #pred_data_info <- extract_daily_training_testing_info(1,list_param=list_param_training_testing_info)

  pred_data_info <- mclapply(1:length(list_raster_obj_files[list_names_tile_id]),FUN=extract_daily_training_testing_info,list_param=list_param_training_testing_info,mc.preschedule=FALSE,mc.cores = num_cores)

  #pred_data_info <- mclapply(1:length(list_raster_obj_files[list_names_tile_id][1:6]),FUN=extract_daily_training_testing_info,list_param=list_param_training_testing_info,mc.preschedule=FALSE,mc.cores = 6)

  #pred_data_info <- lapply(1:length(list_raster_obj_files),FUN=extract_daily_training_testing_info,list_param=list_param_training_testing_info)

  #pred_data_info <- lapply(1:length(list_raster_obj_files[1]),FUN=extract_daily_training_testing_info,list_param=list_param_training_testing_info)

  pred_data_info_tmp <- remove_from_list_fun(pred_data_info)$list #remove data not predicted

  ##Add tile nanmes?? it is alreaready there

  #names(pred_data_info)<-list_names_tile_id

  pred_data_month_info <- do.call(rbind,lapply(pred_data_info_tmp,function(x){x$pred_data_month_info}))

  pred_data_day_info <- do.call(rbind,lapply(pred_data_info_tmp,function(x){x$pred_data_day_info}))

  #putput inforamtion in csv !!

  write.table(pred_data_month_info,

              file=file.path(out_dir,paste("pred_data_month_info_",year_predicted,"_",out_prefix,".txt",sep="")),sep=",")

  write.table(pred_data_day_info,

              file=file.path(out_dir,paste("pred_data_day_info_",year_predicted,"_",out_prefix,".txt",sep="")),sep=",")

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

  ####### PART 4: Get shapefile tiling with centroids ###

  #get shape files for the region being assessed:

  list_shp_world <- list.files(path=in_dir_shp,pattern=".*.shp",full.names=T)

  l_shp <- gsub(".shp","",basename(list_shp_world))

  l_shp <- sub("shp_","",l_shp)

  #l_shp <- unlist(lapply(1:length(list_shp_world),

  #                       FUN=function(i){paste(strsplit(list_shp_world[i],"_")[[1]][3:4],collapse="_")}))

  l_shp <- unlist(lapply(1:length(l_shp),

                         FUN=function(i){paste(strsplit(l_shp[i],"_")[[1]][1:2],collapse="_")}))

  df_tiles_all <- as.data.frame(as.character(unlist(list_shp_world)))

  df_tiles_all$tile_coord <- l_shp

  #names(df_tiles_all) <- "list_shp_world"

  names(df_tiles_all) <- c("shp_files","tile_coord")

  matching_index <- match(basename(in_dir_list),l_shp)

  list_shp_reg_files <- list_shp_world[matching_index]

  df_tile_processed$shp_files <-list_shp_reg_files

  #df_tile_processed$shp_files <- ""

  #df_tile_processed$tile_coord <- as.character(df_tile_processed$tile_coord)

  #test <- df_tile_processed

  #test$shp_files <- NULL

  #test3 <- merge(test,df_tiles_all,by=c("tile_coord"))

  #test3 <- merge(df_tiles_all,test,by=c("tile_coord"))

  #merge(df_tile_processed,df_tiles_all,by="shp_files")

  tx<-strsplit(as.character(df_tile_processed$tile_coord),"_")

  lat<- as.numeric(lapply(1:length(tx),function(i,x){x[[i]][1]},x=tx))

  long<- as.numeric(lapply(1:length(tx),function(i,x){x[[i]][2]},x=tx))

  df_tile_processed$lat <- lat

  df_tile_processed$lon <- long

  #put that list in the df_processed and also the centroids!!

  write.table(df_tile_processed,

              file=file.path(out_dir,paste("df_tile_processed_",year_predicted,"_",out_prefix,".txt",sep="")),sep=",")

  write.table(df_tiles_all,

              file=file.path(out_dir,paste("df_tiles_all_",year_predicted,"_",out_prefix,".txt",sep="")),sep=",")

  #Copy to local home directory on NAS-NEX

  #

  dir.create(file.path(out_dir,"shapefiles"))

  file.copy(list_shp_world,file.path(out_dir,"shapefiles"))

  #save a list of all files...

  write.table(df_tiles_all,

              file=file.path(out_dir,"shapefiles",paste("df_tiles_all_",year_predicted,"_",out_prefix,".txt",sep="")),sep=",")

  ###Prepare files for copying back?

  ## Prepare list of files to return...

  return(1)