#MODIFIED ON: 12/19/2015            

#1) Make this is a script/function callable from the shell/bas

#2) clean up temporary files, it builds currently on the disk

#3) fix output folder for some of output files: create a mosaic output folder if doesn't exist?

#4) create a helper function for inputs/arguments to automate...?? Could also be in the assessment stage

#reg1   : North America

#reg23" : Europe + Asia

#reg4"  : South America

#reg5   : Africa

#reg6   : Oceania+ South East Asia

#

function_mosaicing <-"global_run_scalingup_mosaicing_function_12192015.R"

in_dir_script <-"/home/parmentier/Data/IPLANT_project/env_layers_scripts" #NCEAS UCSB

#in_dir_script <- "/nobackupp8/bparmen1/env_layers_scripts" #NASA NEX

#Data is on ATLAS or NASA NEX

#PARAM 1

in_dir <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_pred_1992_12072015" #NCEAS

#in_dir <- "/nobackupp8/bparmen1/output_run10_1500x4500_global_analyses_pred_1992_12072015" #NEX

in_dir_tiles <- file.path(in_dir,"tiles") #this is valid both for Atlas and NEX

#infile_mask <- "/nobackupp8/bparmen1/regions_input_files/r_mask_reg4.tif"

infile_mask <- "/data/project/layers/commons/NEX_data/regions_input_files/r_mask_reg4.tif"

#in_dir can be on NEX or Atlas

#mosaic_python <- "/nobackupp6/aguzman4/climateLayers/sharedCode/"

#python_bin <- "/nobackupp6/aguzman4/climateLayers/sharedModules2/bin"

mosaic_python <- "/data/project/layers/commons/NEX_data/sharedCode" #PARAM 26

python_bin <- "/usr/bin" #PARAM 27

#for residuals...

list_models <- NULL #PARAM 30

#list_models <- paste(var_pred,"~","1",sep=" ") #if null then this is the default...

#This is is assuming a list of file for a region!! 

#this is where the main function for mosaicing region starts!!

#use reg4 to test the code for now, redo later for any regions!!!

k<-2

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

  region_selected <- region_names[k]

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

  #### First generate rmse images for each date and tile for the region

  lf_mosaic <- reg_lf_mosaic[[k]] #list of files to mosaic by regions

                                #There a 28 files for reg4, South America

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

  ################### PART I: Accuracy layers by tiles #############

  #first generate accuracy layers using tiles definitions and output from the accuracy assessment

  #tb <- list_param$tb #fitting or validation table with all days

  #metric_name <- "rmse" #RMSE, MAE etc.

  #pred_mod_name <- "mod1"

  #y_var_name 

  #interpolation_method #c("gam_CAI") #PARAM3

  days_to_process <- day_to_mosaic

  #NA_flag_val <- list_param$NA_flag_val

  #file_format <- list_param$file_format

  out_dir_str <- out_dir

  out_suffix_str <- out_suffix

  lf <- lf_mosaic

  #Improved by adding multicores option

  num_cores_tmp <- num_cores

  list_param_accuracy_metric_raster <- list(lf,tb,metric_name,pred_mod_name,y_var_name,interpolation_method,

  names(list_param_accuracy_metric_raster) <- c("lf","tb","metric_name","pred_mod_name","y_var_name","interpolation_method",

  list_raster_created_obj <- lapply(1:length(day_to_mosaic),FUN=create_accuracy_metric_raster,

  #debug(create_accuracy_metric_raster)

  #list_raster_created_obj <- lapply(1:1,FUN=create_accuracy_metric_raster,

  #                                  list_param=list_param_accuracy_metric_raster)

  #raster_created_obj <- create_accuracy_metric_raster(1, list_param_accuracy_metric_raster)

  #Extract list of files for rmse and date 1 (19920101), there should be 28 raster images

  lf_accuracy_raster <- lapply(1:length(list_raster_created_obj),FUN=function(i){unlist(list_raster_created_obj[[i]]$list_raster_name)}) 

  #Plot as quick check

  #r1 <- raster(lf_mosaic[[1]][1]) 

  #plot(r1)

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

  ### Now create accuracy surfaces from residuals...

  ## Create accuracy surface by kriging

  num_cores_tmp <-num_cores

  lf_day_tiles  <- lf_mosaic #list of raster files by dates

  data_df <- data_day_v # data.frame table/spdf containing stations with residuals and variable

  #df_tile_processed  #tiles processed during assessment usually by region

  #var_pred  #variable being modeled

  #if not list of models is provided generate one

  if(is.null(list_models)){

    list_models <- paste(var_pred,"~","1",sep=" ")

  }

  #use_autokrige #if TRUE use automap/gstat package

  #y_var_name  #"dailyTmax" #PARAM2

  #interpolation_method #c("gam_CAI") #PARAM3, need to select reg!!

  #date_processed #can be a monthly layer

  #num_cores #number of cores used

  #NA_flag_val 

  #file_format 

  out_dir_str <- out_dir

  out_suffix_str <- out_suffix 

  days_to_process <- day_to_mosaic

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

  df_tile_processed$reg <- basename(dirname(df_tile_processed$path_NEX))

  ##By regions, selected earlier

  #for(k in 1:length(region_names)){

  df_tile_processed_reg <- subset(df_tile_processed,reg==region_selected)#use reg4

  list_create_accuracy_residuals_raster_obj <- lapply(1:length(day_to_mosaic),FUN=create_accuracy_residuals_raster,

                                  list_param=list_param_create_accuracy_residuals_raster)

  #list_create_accuracy_residuals_raster_obj <- lapply(1:1,FUN=create_accuracy_residuals_raster,

  #                                list_param=list_param_create_accuracy_residuals_raster)

  #note that three tiles did not produce a residuals surface!!! find out more later, join the output

  #to df_raste_tile to keep track of which one did not work...

  #lf_accuracy_residuals_raster <- as.character(unlist(lapply(1:length(list_create_accuracy_residuals_raster_obj),FUN=function(i,x){unlist(extract_from_list_obj(x[[i]]$list_pred_res_obj,"raster_name"))},x=list_create_accuracy_residuals_raster_obj))) 

  lf_accuracy_residuals_raster <- lapply(1:length(list_create_accuracy_residuals_raster_obj),FUN=function(i,x){as.character(unlist(extract_from_list_obj(x[[i]]$list_pred_res_obj,"raster_name")))},x=list_create_accuracy_residuals_raster_obj)

  #Plot as quick check

  #r1 <- raster(lf_mosaic[[1]][1]) 

  #list_create_accuracy_residuals_raster_obj

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

  #### PART 2: GENETATE MOSAIC FOR LIST OF FILES #####

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

  #### Mosaic tiles for the variable predicted and accuracy metric

  #methods availbable:use_sine_weights,use_edge,use_linear_weights

  #only use edge method for now

  #loop to dates..., make this a function...

  list_mosaic_obj <- vector("list",length=length(day_to_mosaic))

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

    #

    mosaic_method <- "use_edge_weights" #this is distance from edge

    out_suffix_tmp <- paste(interpolation_method,y_var_name,day_to_mosaic[i],out_suffix,sep="_")

    #debug(mosaicFiles)

    #can also loop through methods!!!

    #python_bin <- "/usr/bin/" #python gdal bin, on Atlas NCEAS

    #python_bin <- "/nobackupp6/aguzman4/climateLayers/sharedModules/bin" #on NEX

    #gdal_merge_sum_noDataTest.py

    mosaic_edge_obj_prediction <- mosaicFiles(lf_mosaic[[i]],

    mosaic_method <- "use_edge_weights" #this is distance from edge

    out_suffix_tmp <- paste(interpolation_method,metric_name,day_to_mosaic[i],out_suffix,sep="_")

    #debug(mosaicFiles)

    #can also loop through methods!!!

    mosaic_edge_obj_accuracy <- mosaicFiles(lf_accuracy_raster[[i]],

    list_mosaic_obj[[i]] <- list(prediction=mosaic_edge_obj_prediction,accuracy=mosaic_edge_obj_accuracy)

    ### produce residuals mosaics

    mosaic_method <- "use_edge_weights" #this is distance from edge

    out_suffix_tmp <- paste(interpolation_method,"kriged_residuals",day_to_mosaic[i],out_suffix,sep="_")

    #lf_tmp<-list.files(pattern="*kriged_residuals.*.tif",full.names=T)

    lf_tmp <- lf_accuracy_residuals_raster[[i]]

    #lf_accuracy_residuals_raster[[i]]

    #debug(mosaicFiles)

    mosaic_edge_obj_residuals <- mosaicFiles(lf_tmp,

                                        match_extent=match_extent,

    list_mosaic_obj[[i]] <- list(prediction=mosaic_edge_obj_prediction,accuracy=mosaic_edge_obj_accuracy,mosaic_edge_obj_residuals)

  }

  ##End of mosaicing function for region predictions