#MODIFIED ON: 12/16/2015            

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

library(automap)                             # Kriging automatic fitting of variogram using gstat

library(rgeos)                               # Geometric, topologic library of functions

function_mosaicing <-"global_run_scalingup_mosaicing_function_12162015b.R"

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

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

load_obj <- function(f)

{

  env <- new.env()

  nm <- load(f, env)[1]

  env[[nm]]

}

create_dir_fun <- function(out_dir,out_suffix){

  if(!is.null(out_suffix)){

    out_name <- paste("output_",out_suffix,sep="")

    out_dir <- file.path(out_dir,out_name)

  }

  #create if does not exists

  if(!file.exists(out_dir)){

    dir.create(out_dir)

  }

  return(out_dir)

}

out_suffix_str <- "run10_1500x4500_global_analyses_pred_1992_12072015" #PARAM 7

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

pred_mod_name <- "mod1" #PARAM 9

var_pred <- "res_mod1" #used in residuals mapping #PARAM 10

out_dir <- in_dir #PARAM 11

create_out_dir_param <- FALSE #PARAM 12

#if daily mosaics NULL then mosaicas all days of the year #PARAM 13

file_format <- ".tif" #PARAM 14

NA_value <- -9999 #PARAM 15

NA_flag_val <- NA_value #PARAM 16

num_cores <- 6 #PARAM 17                 

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

use_autokrige <- F #PARAM 19

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

tb_accuracy_name <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_pred_1992_12072015/tb_diagnostic_v_NA_run10_1500x4500_global_analyses_pred_1992_12072015.txt" #PARAM 21

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" #PARAM 22

data_day_v_name <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_pred_1992_12072015/data_day_v_NAM_run10_1500x4500_global_analyses_pred_1992_12072015.txt" #PARAM 23

data_day_s_name <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_pred_1992_12072015/data_day_s_NAM_run10_1500x4500_global_analyses_pred_1992_12072015.txt" ##PARAM 24

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" ##PARAM 25

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

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

algorithm <- "python" #PARAM 28 #if R use mosaic function for R, if python use modified gdalmerge script from Alberto Guzmann

match_extent <- "FALSE" #PARAM 29 #try without matching!!!

tb <- read.table(tb_accuracy_name,sep=",")

data_day_s <- read.table(file.path(data_day_s_name),sep=",") #daily testing/validation stations by dates and tiles

data_day_v <- read.table(file.path(data_day_v_name),sep=",") #daily training stations by dates and tiles

#Improved by adding multicores option

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

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

## Create accuracy surface by kriging

num_cores_tmp <- 6

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

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

#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

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

#var_pred <- "res_mod1"

#list_models<- c("res_mod1 ~ 1")

##By regions

i<-1 #loop by days/date to process!!

list_param_create_accuracy_residuals_raster <- list(lf_day_tiles,data_df,df_tile_processed,

                    var_pred,list_models,use_autokrige,y_var_name,interpolation_method,

                    days_to_process,num_cores_tmp,NA_flag_val,file_format,out_dir_str,

                    out_suffix_str) 

names(list_param_create_accuracy_residuals_raster) <- c("lf_day_tiles","data_df","df_tile_processed",

                    "var_pred","list_models","use_autokrige","y_var_name","interpolation_method",

                    "days_to_process","num_cores_tmp","NA_flag_val","file_format","out_dir_str",

                    "out_suffix_str") 

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

#                                  list_param=list_param_create_accuracy_residuals_raster_obj)

debug(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)

#create_accuracy_residuals_raster_obj <- create_accuracy_metric_raster(1, list_param_create_accuracy_residuals_raster_obj)