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

script_path <- "/home/parmentier/Data/IPLANT_project/env_layers_scripts" #path to script

#interpolation_method<-c("gam_fusion") #other otpions to be added later

interpolation_method<-c("gam_CAI") #param 2

CRS_interp <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0" #param 3

#CRS_interp <-"+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs";

#out_region_name<-""

#list_models<-c("y_var ~ s(lat,lon,k=5) + s(elev_s,k=3) + s(LST,k=3)") #param 4

metric_name <- "var_pred" #use RMSE if accuracy

#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_dir <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg5/assessment"

#in_dir_mosaic <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg5/mosaic/mosaic"

in_dir <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg6/assessment"

in_dir_mosaic <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg6/mosaics/mosaic" #predicted mosaic

day_start <- "1984101" #PARAM 12 arg 12

day_end <- "20141231" #PARAM 13 arg 13

#infile_mask <- "/nobackupp8/bparmen1/NEX_data/regions_input_files/r_mask_LST_reg4.tif"

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

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

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

#run_figure_by_year <- TRUE # param 10, arg 7

#list_year_predicted <- "1984,2014"

scaling <- 0.01 #was scaled on 100 

#if scaling is null then perform no scaling!!

#df_centroids_fname <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg5/mosaic/output_reg5_1999/df_centroids_19990701_reg5_1999.txt"

#df_centroids_fname <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/output_reg4_1999/df_centroids_19990701_reg4_1999.txt"

#df_centroids_fname <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg6/mosaic/output_reg6_1984/df_centroids_19840101_reg6_1984.txt"

#/nobackupp6/aguzman4/climateLayers/out/reg1/assessment//output_reg1_1984/df_assessment_files_reg1_1984_reg1_1984.txt

#dates to plot and analyze

#l_dates <- c("19990101","19990102","19990103","19990701","19990702","19990703")

#l_dates <- c("19990101","19990102","19990103","19990104","19990105") 

#df_points_extracted_fname <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg5/mosaic/int_mosaics/data_points_extracted.txt"

#df_points_extracted_fname <- NULL #if null extract on the fly

#r_mosaic_fname <- "r_mosaic.RData"

#r_mosaic_fname <- NULL #if null create a stack from input dir

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

####### PART 1: Read in data ########

out_dir <- in_dir

if (create_out_dir_param == TRUE) {

  out_dir <- create_dir_fun(out_dir,out_suffix)

  setwd(out_dir)

}else{

  setwd(out_dir) #use previoulsy defined directory

}

#setwd(out_dir)

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

############ Using predicting first ##########

## using predictions

#pattern_str <- ".*.tif"

pattern_str <-"*.tif"

lf_raster <- list.files(path=in_dir_mosaic,pattern=pattern_str,recursive=F,full.names=T)

r_stack <- stack(lf_raster,quick=T) #this is very fast now with the quick option!

#save(r_mosaic,file="r_mosaic.RData")

#### check for missing dates from list of tif

###This should be a function!!

#####  This could be moved in a separate file!!

###############  PART4: Checking for mosaic produced for given region ##############

## From list of mosaic files predicted extract dates

## Check dates predicted against date range for a given date range

## Join file information to centroids of tiles data.frame

#list_dates_produced <- unlist(mclapply(1:length(lf_mosaic_list),FUN=extract_date,x=lf_mosaic_list,item_no=13,mc.preschedule=FALSE,mc.cores = num_cores))                         

#list_dates_produced <-  mclapply(1:2,FUN=extract_date,x=lf_mosaic_list,item_no=13,mc.preschedule=FALSE,mc.cores = 2)                         

item_no <- 13

date_start <- day_start

date_end <- day_end

#day_start <- "1984101" #PARAM 12 arg 12

#day_end <- "20141231" #PARAM 13 arg 13

#Using default values for parameters exectpt for num_cores=11 instead of 1

#debug(check_missing)

#test_missing <- check_missing(lf=lf_raster, 

#                              pattern_str=NULL,

#                              in_dir=".", #this is not used if lf is given

#                              date_start="1984101",

#                              date_end="20141231",

#                              item_no=13,

#                              out_suffix="",

#                              num_cores=num_cores,

#                              out_dir=".")

##Run this on reg4 and reg5 after

#Add report by year in text file?

#Using specified values for parameters

test_missing <- check_missing(lf=lf_raster, 

                              pattern_str=NULL,

                              in_dir=in_dir_mosaic,

                              date_start="1984101",

                              date_end="20141231",

                              item_no=13,

                              out_suffix=out_suffix,

                              num_cores=num_cores,

                              out_dir=".")

df_time_series <- test_missing$df_time_series

head(df_time_series)

table(df_time_series$missing)

table(df_time_series$year)

#combine polygon

#http://gis.stackexchange.com/questions/155328/merging-multiple-spatialpolygondataframes-into-1-spdf-in-r

#http://gis.stackexchange.com/questions/116388/count-overlapping-polygons-in-single-shape-file-with-r

#### Use the predictions directory

#By region

#For each polygon/tile find polygon overlapping with count and ID (like list w)

#for each polygon/tile and date find if there is a prediction using the tif (multiply number)

#for each date of year report data in table.

#go through table and hsow if there are missing data (no prediction) or report min predictions for tile set?

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

#### Parameters and constants  

#on ATLAS

#in_dir1 <- "/data/project/layers/commons/NEX_data/test_run1_03232014/output" #On Atlas

#parent output dir : contains subset of the data produced on NEX

#in_dir1 <- "/data/project/layers/commons/NEX_data/output_run6_global_analyses_09162014/output20Deg2"

# parent output dir for the curent script analyes

#out_dir <-"/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/" #On NCEAS Atlas

# input dir containing shapefiles defining tiles

#in_dir_shp <- "/data/project/layers/commons/NEX_data/output_run5_global_analyses_08252014/output/subset/shapefiles"

#out_dir <- "/nobackup/bparmen1/" #on NEX

#in_dir_shp <- "/nobackupp4/aguzman4/climateLayers/output4/subset/shapefiles/"

#in_dir <- "/data/project/layers/commons/NEX_data/reg4_assessment"

#list_in_dir_run <-

#in_dir_list <-  c("output_run_global_analyses_pred_2009_reg4","output_run_global_analyses_pred_2010_reg4",

#                  "output_run_global_analyses_pred_2011_reg4","output_run_global_analyses_pred_2012_reg4",

#                  "output_run_global_analyses_pred_2013_reg4","output_run_global_analyses_pred_2014_reg4")

#in_dir_list_filename <- "/data/project/layers/commons/NEX_data/reg4_assessment/stage6_reg4_in_dir_list_02072016.txt"

in_dir <- "" #PARAM 0

  in_dir1 <- file.path(in_dir1,region_name)

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

    out_dir <- create_dir_fun(out_dir,out_prefix)