#MODIFIED ON: 01/01/2016            

in_dir <- "" #PARAM 0

#day_to_mosaic <- c("19920101","19920102","19920103") #PARAM7

#CRS_WGS84 <-    CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84 #CONSTANT1

#tile_size <- "1500x4500" #PARAM 11

#region_name <- "world" #PARAM 13

countries_shp <-"/data/project/layers/commons/NEX_data/countries.shp" #PARAM 13, copy this on NEX too

reg_modified <- TRUE #PARAM 15

region_name <- c("reg4") #reference region to merge if necessary, if world all the regions are together #PARAM 16

#use previous files produced in step 1a and stored in a data.frame

df_assessment_files <- "df_assessment_files_reg4_1984_run_global_analyses_pred_12282015.txt" #PARAM 17

threshold_missing_day <- c(367,365,300,200) #PARM18

#i <- year_predicted

run_assessment_plotting_prediction_fun <-function(list_param_run_assessment_plotting){

  in_dir <- list_param_run_assessment_plotting$in_dir #PARAM 0

  y_var_name <- list_param_run_assessment_plotting$y_var_name #PARAM1

  interpolation_method <- list_param_run_assessment_plotting$interpolation_method #c("gam_CAI") #PARAM2

  out_suffix <- list_param_run_assessment_plotting$out_suffix #PARAM3

  out_dir <- list_param_run_assessment_plotting$out_dir #

  create_out_dir_param <- list_param_run_assessment_plotting$create_out_dir_param # FALSE #PARAM 5

  mosaic_plot <- list_param_run_assessment_plotting$mosaic_plot #FALSE #PARAM6

  proj_str<- list_param_run_assessment_plotting$proj_str #CRS_WGS84 #PARAM 8 #check this parameter

  file_format <- list_param_run_assessment_plotting$file_format #".rst" #PARAM 9

  NA_value <- list_param_run_assessment_plotting$NA_value #-9999 #PARAM10

  multiple_region <- list_param_run_assessment_plotting$multiple_region # <- TRUE #PARAM 12

  countries_shp <- list_param_run_assessment_plotting$countries_shp #<- "world" #PARAM 13

  plot_region <- list_param_run_assessment_plotting$plot_region #<- TRUE

  num_cores <- list_param_run_assessment_plotting$num_cores # 6 #PARAM 14

  reg_modified <- list_param_run_assessment_plotting$reg_modified #<- TRUE #PARAM 15

  region <- list_param_run_assessment_plotting$region #<- c("reg4") #reference region to merge if necessary #PARAM 16

  region_name <- list_param_run_assessment_plotting$region_name #<- "world" #PARAM 13

  df_assessment_files <- list_param_run_assessment_plotting$df_assessment_files #PARAM 16

  threshold_missing_day <- list_param_run_assessment_plotting$threshold_missing_day #PARM18

  NA_flag_val <- NA_value

  #tile_size <- "1500x4500" #PARAM 11

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

  ###Table 1: Average accuracy metrics

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

  #df_assessment_files, note that in_dir indicate the path of the textfiles

  summary_metrics_v <- file.path(in_dir,basename(df_assessment_files$files[1]))

  tb <- file.path(in_dir, basename(df_assessment_files$files[2]))

  tb_s <-file.path(in_dir, basename(df_assessment_files$files[4]))

  tb_month_s <- file.path(in_dir,basename(df_assessment_files$files[3]))

  pred_data_month_info <- file.path(in_dir, basename(df_assessment_files$files[10]))

  pred_data_day_info <- file.path(in_dir, basename(df_assessment_files$files[11]))

  df_tile_processed <- file.path(in_dir, basename(df_assessment_files$files[12]))

  #add column for tile size later on!!!

  tb$pred_mod <- as.character(tb$pred_mod)

  summary_metrics_v$pred_mod <- as.character(summary_metrics_v$pred_mod)

  summary_metrics_v$tile_id <- as.character(summary_metrics_v$tile_id)

  df_tile_processed$tile_id <- as.character(df_tile_processed$tile_id)

  tb_month_s$pred_mod <- as.character(tb_month_s$pred_mod)

  tb_month_s$tile_id<- as.character(tb_month_s$tile_id)

  tb_s$pred_mod <- as.character(tb_s$pred_mod)

  tb_s$tile_id <- as.character(tb_s$tile_id)

  #multiple regions?

  if(multiple_region==TRUE){

    df_tile_processed$reg <- basename(dirname(as.character(df_tile_processed$path_NEX)))

    tb <- merge(tb,df_tile_processed,by="tile_id")

    tb_s <- merge(tb_s,df_tile_processed,by="tile_id")

    tb_month_s<- merge(tb_month_s,df_tile_processed,by="tile_id")

    summary_metrics_v <- merge(summary_metrics_v,df_tile_processed,by="tile_id")

  }

  tb_all <- tb

  summary_metrics_v_all <- summary_metrics_v 

  #table(summary_metrics_v_all$reg)

  #table(summary_metrics_v$reg)

  #table(tb_all$reg)

  #table(tb$reg)

  ############ PART 2: PRODUCE FIGURES ################

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

  ### Figure 1: plot location of the study area with tiles processed

  #df_tiled_processed <- na.omit(df_tile_processed) #remove other list of folders irrelevant

  #list_shp_reg_files <- df_tiled_processed$shp_files

  list_shp_reg_files<- as.character(df_tile_processed$shp_files)

  #list_shp_reg_files <- file.path("/data/project/layers/commons/NEX_data/",out_dir,

  #          as.character(df_tile_processed$tile_coord),"shapefiles",basename(list_shp_reg_files))

  #list_shp_reg_files <- file.path("/data/project/layers/commons/NEX_data/",out_dir,

                                  #"shapefiles",basename(list_shp_reg_files))

  #table(summary_metrics_v$reg)

  #table(summary_metrics_v$reg)

  ### Potential function starts here:

  #function(in_dir,out_dir,list_shp_reg_files,title_str,region_name,num_cores,out_suffix,out_suffix)

  ### First get background map to display where study area is located

  #can make this more general later on..should have this already in a local directory on Atlas or NEX!!!!

  centroids_pts <- vector("list",length(list_shp_reg_files))

  shps_tiles <- vector("list",length(list_shp_reg_files))

  #collect info: read in all shapfiles

  #This is slow...make a function and use mclapply??

  #/data/project/layers/commons/NEX_data/output_run6_global_analyses_09162014/shapefiles

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

    #path_to_shp <- dirname(list_shp_reg_files[[i]])

    path_to_shp <- file.path(out_dir,"/shapefiles")

    layer_name <- sub(".shp","",basename(list_shp_reg_files[[i]]))

    shp1 <- try(readOGR(path_to_shp, layer_name)) #use try to resolve error below

    #shp_61.0_-160.0

    #Geographical CRS given to non-conformant data: -186.331747678

    #shp1<-readOGR(dirname(list_shp_reg_files[[i]]),sub(".shp","",basename(list_shp_reg_files[[i]])))

    if (!inherits(shp1,"try-error")) {

    }else{

      pt <- shp1

      centroids_pts[[i]] <- pt

    }

    shps_tiles[[i]] <- shp1

    #centroids_pts[[i]] <- centroids

  #fun <- function(i,list_shp_files)

  #coord_names <- c("lon","lat")

  #l_ras#t <- rasterize_df_fun(test,coord_names,proj_str,out_suffix=out_suffix,out_dir=".",file_format,NA_flag_val,tolerance_val=0.000120005)

  #remove try-error polygons...we loose three tiles because they extend beyond -180 deg

  tmp <- shps_tiles

  shps_tiles <- remove_errors_list(shps_tiles) #[[!inherits(shps_tiles,"try-error")]]

  #shps_tiles <- tmp

  length(tmp)-length(shps_tiles) #number of tiles with error message

  tmp_pts <- centroids_pts 

  centroids_pts <- remove_errors_list(centroids_pts) #[[!inherits(shps_tiles,"try-error")]]

  #centroids_pts <- tmp_pts 

  #plot info: with labels

  res_pix <-1200

  col_mfrow <- 1 

  png(filename=paste("Figure1_tile_processed_region_",region_name,"_",out_suffix,".png",sep=""),

      width=col_mfrow*res_pix,height=row_mfrow*res_pix)

  plot(reg_layer)

  #Add polygon tiles...

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

    shp1 <- shps_tiles[[i]]

    pt <- centroids_pts[[i]]

    if(!inherits(shp1,"try-error")){

      plot(shp1,add=T,border="blue")

      #plot(pt,add=T,cex=2,pch=5)

      label_id <- df_tile_processed$tile_id[i]

      text(coordinates(pt)[1],coordinates(pt)[2],labels=i,cex=1.3,font=2,col=c("red"))

    }

  }

  #title(paste("Tiles ", tile_size,region_name,sep=""))

  #unique(summaty_metrics$tile_id)

  #text(lat-shp,)

  #union(list_shp_reg_files[[1]],list_shp_reg_files[[2]])

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

  ### Figure 2: boxplot of average accuracy by model and by tiles

  tb$tile_id <- factor(tb$tile_id, levels=unique(tb$tile_id))

  model_name <- as.character(unique(tb$pred_mod))

  ## Figure 2a

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

    res_pix <- 480

    col_mfrow <- 1

    row_mfrow <- 1

    png(filename=paste("Figure2a_boxplot_with_oultiers_by_tiles_",model_name[i],"_",out_suffix,".png",sep=""),

        width=col_mfrow*res_pix,height=row_mfrow*res_pix)

    boxplot(rmse~tile_id,data=subset(tb,tb$pred_mod==model_name[i]))

    title(paste("RMSE per ",model_name[i]))

    dev.off()

  }

  ## Figure 2b

  #wtih ylim and removing trailing...

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

    res_pix <- 480

    col_mfrow <- 1

    row_mfrow <- 1

    png(filename=paste("Figure2b_boxplot_scaling_by_tiles","_",model_name[i],"_",out_suffix,".png",sep=""),

        width=col_mfrow*res_pix,height=row_mfrow*res_pix)

    model_name <- unique(tb$pred_mod)

    boxplot(rmse~tile_id,data=subset(tb,tb$pred_mod==model_name[i])

            ,ylim=c(0,4),outline=FALSE)

    title(paste("RMSE per ",model_name[i]))

    dev.off()

  }

  #bwplot(rmse~tile_id, data=subset(tb,tb$pred_mod=="mod1"))

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

  ### Figure 3: boxplot of average RMSE by model acrosss all tiles

  ## Figure 3a

  png(filename=paste("Figure3a_boxplot_overall_region_with_oultiers_",model_name[i],"_",out_suffix,".png",sep=""),

      width=col_mfrow*res_pix,height=row_mfrow*res_pix)

  boxplot(rmse~pred_mod,data=tb)#,names=tb$pred_mod)

  title("RMSE per model over all tiles")

  ## Figure 3b

  png(filename=paste("Figure3b_boxplot_overall_region_scaling_",model_name[i],"_",out_suffix,".png",sep=""),

      width=col_mfrow*res_pix,height=row_mfrow*res_pix)

  boxplot(rmse~pred_mod,data=tb,ylim=c(0,5),outline=FALSE)#,names=tb$pred_mod)

  title("RMSE per model over all tiles")

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

  ### Figure 4: plot predicted tiff for specific date per model

  #y_var_name <-"dailyTmax"

  #index <-244 #index corresponding to Sept 1

  # if (mosaic_plot==TRUE){

  #   index  <- 1 #index corresponding to Jan 1

  #   date_selected <- "20100901"

  #   name_method_var <- paste(interpolation_method,"_",y_var_name,"_",sep="")

  # 

  #   pattern_str <- paste("mosaiced","_",name_method_var,"predicted",".*.",date_selected,".*.tif",sep="")

  #   lf_pred_list <- list.files(pattern=pattern_str)

  # 

  #   for(i in 1:length(lf_pred_list)){

  #     

  #   

  #     r_pred <- raster(lf_pred_list[i])

  #   

  #     res_pix <- 480

  #     col_mfrow <- 1

  #     row_mfrow <- 1

  #   

  #     png(filename=paste("Figure4_models_predicted_surfaces_",model_name[i],"_",name_method_var,"_",data_selected,"_",out_suffix,".png",sep=""),

  #        width=col_mfrow*res_pix,height=row_mfrow*res_pix)

  #   

  #     plot(r_pred)

  #     title(paste("Mosaiced",model_name[i],name_method_var,date_selected,sep=" "))

  #     dev.off()

  #   }

  #   #Plot Delta and clim...

  # 

  #    ## plotting of delta and clim for later scripts...

  # 

  # }

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

  ### Figure 5: plot accuracy ranked 

  #Turn summary table to a point shp

  list_df_ac_mod <- vector("list",length=length(model_name))

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

    ac_mod <- summary_metrics_v[summary_metrics_v$pred_mod==model_name[i],]

    ### Ranking by tile...

    df_ac_mod <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("pred_mod","rmse","mae","tile_id")]

    list_df_ac_mod[[i]] <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("rmse","mae","tile_id")]

    res_pix <- 480

    col_mfrow <- 1

    row_mfrow <- 1

    png(filename=paste("Figure5_ac_metrics_ranked_",model_name[i],"_",out_suffix,".png",sep=""),

        width=col_mfrow*res_pix,height=row_mfrow*res_pix)

    x<- as.character(df_ac_mod$tile_id)

    barplot(df_ac_mod$rmse, names.arg=x)

    #plot(ac_mod1,cex=sqrt(ac_mod1$rmse),pch=1,add=T)

    #plot(ac_mod1,cex=(ac_mod1$rmse1)*2,pch=1,add=T)

    title(paste("RMSE ranked by tile for ",model_name[i],sep=""))

    dev.off()

  }

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

  ### Figure 6: plot map of average RMSE per tile at centroids

  ### Without 

  #list_df_ac_mod <- vector("list",length=length(lf_pred_list))

  list_df_ac_mod <- vector("list",length=2)

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

    ac_mod <- summary_metrics_v[summary_metrics_v$pred_mod==model_name[i],]

    #r_pred <- raster(lf_list[i])

    res_pix <- 1200

    #res_pix <- 480

    col_mfrow <- 1

    row_mfrow <- 1

    png(filename=paste("Figure6_ac_metrics_map_centroids_tile_",model_name[i],"_",out_suffix,".png",sep=""),

        width=col_mfrow*res_pix,height=row_mfrow*res_pix)

    #plot(r_pred)  

    #plot(reg_layer)

    #plot(ac_mod1,cex=sqrt(ac_mod1$rmse),pch=1,add=T)

    #plot(ac_mod,cex=(ac_mod$rmse^2)/10,pch=1,col="red",add=T)

    #coordinates(ac_mod) <- ac_mod[,c("lon","lat")] 

    coordinates(ac_mod) <- ac_mod[,c("lon.x","lat.x")] #solve this later

    p_shp <- layer(sp.polygons(reg_layer, lwd=1, col='black'))

    #title("(a) Mean for 1 January")

    p <- bubble(ac_mod,"rmse",main=paste("Averrage RMSE per tile and by ",model_name[i]))

    p1 <- p+p_shp

    print(p1)

    #plot(ac_mod1,cex=(ac_mod1$rmse1)*2,pch=1,add=T)

    #title(paste("Averrage RMSE per tile and by ",model_name[i]))

    dev.off()

    ### Ranking by tile...

    #df_ac_mod <- 

    list_df_ac_mod[[i]] <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("rmse","mae","tile_id")]

  }

  ## Number of tiles with information:

  sum(df_tile_processed$metrics_v) #26,number of tiles with raster object

  length(df_tile_processed$metrics_v) #26,number of tiles in the region

  sum(df_tile_processed$metrics_v)/length(df_tile_processed$metrics_v) #80 of tiles with info

  #coordinates

  #coordinates(summary_metrics_v) <- c("lon","lat")

  coordinates(summary_metrics_v) <- c("lon.y","lat.y")

  #threshold_missing_day <- c(367,365,300,200)

  nb<-nrow(subset(summary_metrics_v,model_name=="mod1"))  

  sum(subset(summary_metrics_v,model_name=="mod1")$n_missing)/nb #33/35

  ## Make this a figure...

  #plot(summary_metrics_v)

  #Make this a function later so that we can explore many thresholds...

  j<-1 #for model name 1

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

    #summary_metrics_v$n_missing <- summary_metrics_v$n == 365

    #summary_metrics_v$n_missing <- summary_metrics_v$n < 365

    summary_metrics_v$n_missing <- summary_metrics_v$n < threshold_missing_day[i]

    summary_metrics_v_subset <- subset(summary_metrics_v,model_name=="mod1")

    #res_pix <- 1200

    res_pix <- 960

    col_mfrow <- 1

    row_mfrow <- 1

    png(filename=paste("Figure7a_ac_metrics_map_centroids_tile_",model_name[j],"_","missing_day_",threshold_missing_day[i],

                       "_",out_suffix,".png",sep=""),

        width=col_mfrow*res_pix,height=row_mfrow*res_pix)

    model_name[j]

    p_shp <- layer(sp.polygons(reg_layer, lwd=1, col='black'))

    #title("(a) Mean for 1 January")

    p <- bubble(summary_metrics_v_subset,"n_missing",main=paste("Missing per tile and by ",model_name[j]," for ",

                                                                threshold_missing_day[i]))

    p1 <- p+p_shp

    try(print(p1)) #error raised if number of missing values below a threshold does not exist

    #plot(ac_mod1,cex=(ac_mod1$rmse1)*2,pch=1,add=T)

    #title(paste("Averrage RMSE per tile and by ",model_name[i]))

    dev.off()

  }

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

  ### Figure 7: Number of predictions: daily and monthly

  #xyplot(rmse~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v),

  #                                           pred_mod!="mod_kr"),type="b")

  #xyplot(n~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v),

  #                                           pred_mod!="mod_kr"),type="h")

  #cor

  # 

  ## Figure 7a

  png(filename=paste("Figure7a_number_daily_predictions_per_models","_",out_suffix,".png",sep=""),

      width=col_mfrow*res_pix,height=row_mfrow*res_pix)

  xyplot(n~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v),

                                          pred_mod!="mod_kr"),type="h")

  table(tb$pred_mod)

  table(tb$index_d)

  table(subset(tb,pred_mod!="mod_kr"))

  table(subset(tb,pred_mod=="mod1")$index_d)

  #aggregate()

  tb$predicted <- 1

  test <- aggregate(predicted~pred_mod+tile_id,data=tb,sum)

  xyplot(predicted~pred_mod | tile_id,data=subset(as.data.frame(test),

                                                  pred_mod!="mod_kr"),type="h")

  as.character(unique(test$tile_id)) #141 tiles

  dim(subset(test,test$predicted==365 & test$pred_mod=="mod1"))

  histogram(subset(test, test$pred_mod=="mod1")$predicted)

  unique(subset(test, test$pred_mod=="mod1")$predicted)

  table((subset(test, test$pred_mod=="mod1")$predicted))

  #LST_avgm_min <- aggregate(LST~month,data=data_month_all,min)

  histogram(test$predicted~test$tile_id)

  #table(tb)

  ## Figure 7b

  #png(filename=paste("Figure7b_number_daily_predictions_per_models","_",out_suffix,".png",sep=""),

  #    width=col_mfrow*res_pix,height=row_mfrow*res_pix)

  #xyplot(n~month | tile_id + pred_mod,data=subset(as.data.frame(tb_month_s),

  #                                           pred_mod!="mod_kr"),type="h")

  #xyplot(n~month | tile_id,data=subset(as.data.frame(tb_month_s),

  #                                           pred_mod="mod_1"),type="h")

  #test=subset(as.data.frame(tb_month_s),pred_mod="mod_1")

  #table(tb_month_s$month)

  #dev.off()

  #

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

  ##### Figure 8: Breaking down accuracy by regions!! #####

  #summary_metrics_v <- merge(summary_metrics_v,df_tile_processed,by="tile_id")

  ## Figure 8a

  res_pix <- 480

  png(filename=paste("Figure8a_boxplot_overall_separated_by_region_with_oultiers_",model_name[i],"_",out_suffix,".png",sep=""),

      width=col_mfrow*res_pix,height=row_mfrow*res_pix)

  p<- bwplot(rmse~pred_mod | reg, data=tb,

             main="RMSE per model and region over all tiles")

  print(p)

  ## Figure 8b

  png(filename=paste("Figure8b_boxplot_overall_separated_by_region_scaling_",model_name[i],"_",out_suffix,".png",sep=""),

      width=col_mfrow*res_pix,height=row_mfrow*res_pix)

  boxplot(rmse~pred_mod,data=tb,ylim=c(0,5),outline=FALSE)#,names=tb$pred_mod)

  title("RMSE per model over all tiles")

  p<- bwplot(rmse~pred_mod | reg, data=tb,ylim=c(0,5),

             main="RMSE per model and region over all tiles")

  print(p)

  dev.off()

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

  #### Figure 9: plotting mosaics of regions ###########

  ## plot mosaics...

  #First collect file names

  #names(lf_mosaics_reg) <- l_reg_name

  #This part should be automated...

  #plot Australia

  #lf_m <- lf_mosaics_reg[[2]]

  #out_dir_str <- out_dir

  #reg_name <- "reg6_1000x3000"

  #list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str,out_suffix=out_suffix)

  #list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str,out_suffix=out_suffix,l_dates=day_to_mosaic)

  #debug(plot_daily_mosaics)

  #lf_m_mask_reg6_1000x3000 <- plot_daily_mosaics(1,list_param=list_param_plot_daily_mosaics)

  #lf_m_mask_reg6_1000x3000 <- mclapply(1:length(lf_m),FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 10)

  ################## WORLD MOSAICS NEEDS MAJOR CLEAN UP OF CODE HERE

  ##make functions!!

  ###Combine mosaics with modified code from Alberto

  #use list from above!!

  # test_list <-list.files(path=file.path(out_dir,"mosaics"),    

  #            pattern=paste("^world_mosaics.*.tif$",sep=""), 

  # )

  # #world_mosaics_mod1_output1500x4500_km_20101105_run10_1500x4500_global_analyses_03112015.tif

  # 

  # #test_list<-lapply(1:30,FUN=function(i){lapply(1:x[[i]]},x=lf_mosaics_mask_reg)

  # #test_list<-unlist(test_list)

  # #mosaic_list_mean <- vector("list",length=1)

  # mosaic_list_mean <- test_list 

  # out_rastnames <- "world_test_mosaic_20100101"

  # out_path <- out_dir

  # 

  # list_param_mosaic<-list(mosaic_list_mean,out_path,out_rastnames,file_format,NA_flag_val,out_suffix)

  # names(list_param_mosaic)<-c("mosaic_list","out_path","out_rastnames","file_format","NA_flag_val","out_suffix")

  # #mean_m_list <-mclapply(1:12, list_param=list_param_mosaic, mosaic_m_raster_list,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement

  # 

  # lf <- mosaic_m_raster_list(1,list_param_mosaic)

  # 

  # debug(mosaic_m_raster_list)

  # mosaic_list<-list_param$mosaic_list

  # out_path<-list_param$out_path

  # out_names<-list_param$out_rastnames

  # file_format <- list_param$file_format

  # NA_flag_val <- list_param$NA_flag_val

  # out_suffix <- list_param$out_suffix

  ##Now mosaic for mean: should reorder files!!

  #out_rastnames_mean<-paste("_",lst_pat,"_","mean",out_suffix,sep="")

  #list_date_names<-c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")

  #mosaic_list<-split(tmp_str1,list_date_names)

  #new_list<-vector("list",length(mosaic_list))

  # for (i in 1:length(list_date_names)){

  #    j<-grep(list_date_names[i],mosaic_list,value=FALSE)

  #    names(mosaic_list)[j]<-list_date_names[i]

  #    new_list[i]<-mosaic_list[j]

  #  }

  #  mosaic_list_mean <-new_list #list ready for mosaicing

  #  out_rastnames_mean<-paste(list_date_names,out_rastnames_mean,sep="")

  ### Now mosaic tiles...Note that function could be improved to use less memory

  ################## PLOTTING WORLD MOSAICS ################

  #lf_world_pred <-list.files(path=file.path(out_dir,"mosaics"),    

  #           pattern=paste("^world_mosaics.*.tif$",sep=""),full.names=T) 

  lf_world_pred <-list.files(path=file.path(out_dir,"mosaics"),    

                             pattern=paste("^reg5.*.",".tif$",sep=""),full.names=T) 

  l_reg_name <- unique(df_tile_processed$reg)

  lf_world_pred <-list.files(path=file.path(out_dir,l_reg_name[[i]]),    

                             pattern=paste(".tif$",sep=""),full.names=T) 

  prefix_str <- "Figure10_clim_world_mosaics_day_"

  l_dates <-day_to_mosaic

  screenRast=FALSE

  debug(plot_screen_raster_val)

  world_m_list1<- plot_screen_raster_val(1,list_param_plot_screen_raster)

  s_pred <- stack(lf_raster_fname)

  res_pix <- 1500

  col_mfrow <- 3 

  row_mfrow <- 2

  png(filename=paste("Figure10_levelplot_combined_",region_name,"_",out_suffix,".png",sep=""),

      width=col_mfrow*res_pix,height=row_mfrow*res_pix)

  levelplot(s_pred,layers=1:6,col.regions=rev(terrain.colors(255)),cex=4)

  dev.off()

  #   blues<- designer.colors(6, c( "blue", "white") )

  # reds <- designer.colors(6, c( "white","red")  )

  # colorTable<- c( blues[-6], reds)

  # breaks with a gap of 10 to 17 assigned the white color

  # brks<- c(seq( 1, 10,,6), seq( 17, 25,,6)) 

  # image.plot( x,y,z,breaks=brks, col=colorTable)

  #

  #lf_world_mask_reg <- vector("list",length=length(lf_world_pred))

  #for(i in 1:length(lf_world_pred)){

  #

  #  lf_m <- lf_mosaics_reg[i]

  #  out_dir_str <- out_dir

  #reg_name <- paste(l_reg_name[i],"_",tile_size,sep="") #make this automatic

  #lapply()

  #  list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=tile_size,out_dir_str=out_dir_str,out_suffix=out_suffix,l_dates=day_to_mosaic)

  #lf_m_mask_reg4_1500x4500 <- mclapply(1:2,FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 6)

  #lf_world_mask_reg[[i]] <- mclapply(1:length(lf_m),FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 10)

  }

  ############# NEW MASK AND DATA

  ## Plot areas and day predicted as first check

  l_reg_name <- unique(df_tile_processed$reg)

  #(subset(df_tile_processed$reg == l_reg_name[i],date)

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

    lf_world_pred<-list.files(path=file.path(out_dir,l_reg_name[[i]]),    

                            pattern=paste(".tif$",sep=""),full.names=T) 

    #mosaic_list_mean <- test_list 

    #out_rastnames <- "world_test_mosaic_20100101"

    #out_path <- out_dir

    #lf_world_pred <- list.files(pattern="world.*2010090.*.tif$")

    #lf_raster_fname <- list.files(pattern="world.*2010*.*02162015.tif$",full.names=T)

    lf_raster_fname <- lf_world_pred

    prefix_str <- paste("Figure10_",l_reg_name[i],sep="")

    l_dates <- basename(lf_raster_fname)

    tmp_name <- gsub(".tif","",l_dates)

    tmp_name <- gsub("gam_CAI_dailyTmax_predicted_mod1_0_1_","",tmp_name)

    #l_dates <- tmp_name

    l_dates <- paste(l_reg_name[i],"_",tmp_name,sep="")

    screenRast=TRUE

    list_param_plot_screen_raster <- list(lf_raster_fname,screenRast,l_dates,out_dir,prefix_str,out_suffix)

    names(list_param_plot_screen_raster) <- c("lf_raster_fname","screenRast","l_dates","out_dir_str","prefix_str","out_suffix_str")

    #undebug(plot_screen_raster_val)

    #world_m_list1<- plot_screen_raster_val(1,list_param_plot_screen_raster)

    #world_m_list <- mclapply(11:30, list_param=list_param_plot_screen_raster, plot_screen_raster_val,mc.preschedule=FALSE,mc.cores = num_cores) #This is the end bracket from mclapply(...) statement

    world_m_list <- mclapply(1:length(l_dates), list_param=list_param_plot_screen_raster, plot_screen_raster_val,mc.preschedule=FALSE,mc.cores = num_cores) #This is the end bracket from mclapply(...) statement

    #s_pred <- stack(lf_raster_fname)

    #res_pix <- 1500

    #col_mfrow <- 3 

    #row_mfrow <- 2

    #png(filename=paste("Figure10_levelplot_combined_",region_name,"_",out_suffix,".png",sep=""),

    #  width=col_mfrow*res_pix,height=row_mfrow*res_pix)

    #levelplot(s_pred,layers=1:6,col.regions=rev(terrain.colors(255)),cex=4)

    #dev.off()