#MODIFIED ON: 12/12/2015            

function_mosaicing <-"global_run_scalingup_mosaicing_function_12122015.R"

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

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

out_suffix <- paste(region_name,"_","run10_1500x4500_global_analyses_pred_1992_12072015",sep="") #PARAM 6

out_suffix_str <- "run10_1500x4500_global_analyses_pred_1992_12072015"

region_names <- c("reg23","reg4") #selected region names, #PARAM2 

###Make a separate folder for masks by regions...

#tb_accuracy_name <- "/nobackupp8/bparmen1/output_run10_1500x4500_global_analyses_pred_1992_10052015/tb_diagnostic_v_NA_run10_1500x4500_global_analyses_pred_1992_10052015.txt"

tb_accuracy_name <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_pred_1992_12072015"

data_month_s_name <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_pred_1992_12072015/data_month_s_NAM_run10_1500x4500_global_analyses_pred_1992_12072015.txt"

df_tile_processed_name <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_pred_1992_12072015/df_tile_processed_run10_1500x4500_global_analyses_pred_1992_12072015.txt"

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

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

mosaic_python <- "/data/project/layers/commons/NEX_data/sharedCode"

python_bin <- "/usr/bin"

# accuracy table by tiles

# textfiles of stations by month

data_month_s <- read.table(file.path(data_month_s_name),sep=",")

df_tile_processed <- read.table( df_tile_processed_name,sep=",")

#Take into account multiple region in some cases!!!      

           pattern=paste(".*.",day_to_mosaic[i],".*.tif$",sep=""),full.names=T,recursive=T)})

#interpolation_method_day_function_multisampling.R

## Create accuracy surface by kriging

data_tmp <- subset(data_month_s,month==1)

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

var_pred <- "res_mod1"

list_models<- c("res_mod1 ~ 1")

df_tile_processed

#list_models<-c("y_var ~ lat*lon + elev_s + N_w*E_w",

#               "y_var ~ lat*lon + elev_s + DISTOC",

#               "y_var ~ lat*lon + elev_s + LST",

#               "y_var ~ lat*lon + elev_s + LST + I(LST*LC1)")

##By regions

create_accuracy_metric_raster_kriging <- function(data_df,df_tile_processed,var_pred,lf_ref_rast,out_dir,out_suffix){

  #This function generates kriged values for a giving set of formula and data input.

  #Inputs:

  #lf: list of raster files

  #tb: data.frame table #fitting or validation table with all days

  #metric_name: accuracy metric selected to be mapped, RMSE, MAE etc.

  #pred_mod_name: model selected, such as mod1, mod_kr etc.

  #y_var_name: variable being modeled e.g."dailyTmax", dailyTmin, precip  

  #interpolation_method: names of the interpolation/modeling method

  #date_processed: day being processed , e.g. 19920101

  #num_cores : number of cores used in the parallelization

  #NA_flag_val: value used as flag in the raster 

  #file_format: e.g. tif., .rst

  #out_dir_str: output directory

  #out_suffix_str: output suffix

  #Outputs:

  #raster list of weights and product of wegihts and inuts

  #TODO: 

  # - improve efficiency

  #

  #promote to spdf

  list_formulas <- lapply(list_models,as.formula,env=.GlobalEnv) #mulitple arguments passed to lapply!!

  #data_training <- data_df 

  coordinates(data_df)<-cbind(data_df$x,data_df$y)

  #proj4string(data_s)<-proj_str

  #coordinates(data_v)<-cbind(data_v$x,data_v$y)

  #proj4string(data_v)<-proj_str

  lf_day_tiles <- lf_mosaic[[i]] #i index for time...

  #Now match the correct tiles with data used in kriging...

  #match the correct tile!!! df_tile_processed

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

  list_tile_coord <- as.character(df_tile_processed$tile_coord)

  pattern_str <- glob2rx(paste("*",list_tile_coord,"*","*.tif",sep=""))

  keywords_str <- pattern_str

  tmp_str2 <-unlist(lapply(keywords_str,grep,lf_day_tiles,value=TRUE))

  df_raster_pred_tiles_tmp <- data.frame(files =tmp_str2, tile_coord=list_tile_coord)

  df_raster_pred_tiles <- merge(df_raster_pred_tiles_tmp,df_tile_processed,by="tile_coord")

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

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

  df_raster_pred_tiles$files <- as.character(df_raster_pred_tiles$files)

  j <- 1 #loops across tiles from set of files

  r_stack <- df_raster_pred_tiles$files[j]

  tile_selected <- as.character(df_raster_pred_tiles$tile_id[j])

  data_training <- subset(data_df,tile_id==tile_selected) #df_raster_pred_tiles$files

  #r_stack <- stack(lf_day[j])

  out_filename <- "test.tif"

  #debug(predict_auto_krige_raster_model)

  test <- predict_auto_krige_raster_model(list_formulas,r_stack,data_training,out_filename)

  #identify residuals

  #call kriging function

  #write output

}