#MODIFIED ON: 09/02/2016            

#COMMIT: testing plot for both training and testing and introducing a new function

#source(file.path(script_path,"NASA2016_conference_temperature_predictions_function_05032016b.R"))

function_mosaicing_functions <- "global_run_scalingup_mosaicing_function_08232016.R" #Functions used to mosaic predicted tiles

function_mosaicing <-"global_run_scalingup_mosaicing_08222016.R" #main scripts for mosaicing predicted tiles

function_assessment_part1_functions <- "global_run_scalingup_assessment_part1_functions_12282015.R" #PARAM12

function_assessment_part3 <- "global_run_scalingup_assessment_part3_07292016.R"

l_dates <- c("19990101","19990102","19990103","19990104","19990105") #dates to plot and analze

df_data_v_fname <- data.frame(year_str,region_name,dataset="data_v",file=list_data_v_fname)

##start new function

df_combined_data_v <- do.call(rbind,list_df_v_stations) #reading only the years related to the the dates e.g. 1999

df_combined_data_s <- do.call(rbind,list_df_s_stations)

list_df_points_data_v_dates <- vector("list",length=length(l_dates))

list_df_points_data_s_dates <- vector("list",length=length(l_dates))

  list_df_points_data_v_dates[[i]] <- subset(df_combined_data_v,df_combined_data_v$date==l_dates[i])

  list_df_points_data_s_dates[[i]] <- subset(df_combined_data_s,df_combined_data_s$date==l_dates[i])

df_combined_data_v_dates <- do.call(rbind,list_df_points_data_v_dates)

df_combined_data_s_dates <- do.call(rbind,list_df_points_data_s_dates)

####

proj_str <- CRS_locs_WGS84

list_param_plot_stations_val_by_date <- list(l_dates,df_combined_data_v_dates,countries_shp_tmp,CRS_locs_WGS84,

                                             station_type_name,model_name,

                                             var_selected,y_var_name,out_suffix,out_dir)

names(list_param_plot_stations_val_by_date) <- c("l_dates","df_combined_data_points","countries_shp","proj_str","

                                                 station_type_name","model_name",

                                                 "var_selected","y_var_name","out_suffix","out_dir")

debug(plot_stations_val_by_date)

test<- plot_stations_val_by_date(1,list_param = list_param_plot_stations_val_by_date)

plot_stations_val_by_date <- function(i,list_param){

  #

  #

  #function to plot residuals by date

  #

  #####

  l_dates <- list_param$l_dates

  model_name <- list_param$model_name

  station_type_name <- list_param$station_type_name

  var_selected <- list_param$var_selected

  #list_df_points_dates <- 

  df_combined_data_points <- list_param$df_combined_data_points

  countries_shp <- list_param$countries_shp

  proj_str <- list_param$proj_str

  out_suffix <- list_param$out_suffix

  out_dir <- list_param$out_dir

  ### STARTFUNCTION ####

  df_points <- subset(df_combined_data_points,date==date_processed)

  #df_points <- list_df_points_dates[[i]]

  freq_tile_id <- as.data.frame(table(df_points$tile_id))

  freq_tile_id$date <- date_processed

  names(freq_tile_id) <- c("tile_id","freq","date")

  #proj4string(df_points) <- CRS_locs_WGS84

  proj4string(df_points) <- proj_str

  #dim(df_points_day)

  #dim(df_points)

  if(class(countries_shp)!="SpatialPolygonsDataFrame"){

  }else{

    reg_layer <- countries_shp

  #data_stations_var_pred <- aggregate(id ~ date, data = data_stations, min)

  #data_stations_var_pred <- aggregate(id ~ x + y + date + dailyTmax + mod1 + res_mod1 , data = data_stations, min)

  #change the formula later on to use different y_var_name (tmin and precip)

  data_stations_var_pred <- aggregate(id ~ x + y + date + dailyTmax + res_mod1 + tile_id + reg ,data = df_points, mean ) #+ mod1 + res_mod1 , data = data_stations, min)

  dim(data_stations_var_pred)

  #> dim(data_stations_var_pred)

  data_stations_var_pred$date_str <- data_stations_var_pred$date

  data_stations_var_pred$date <- as.Date(strptime(data_stations_var_pred$date_str,"%Y%m%d"))

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

  #data_stations_var_pred$constant <- c(1000,2000)

  #plot(data_stations_var_pred)

  res_pix <- 800

  #plot(data_stations_var_pred

  #p <- bubble(data_stations_var_pred,"constant",main=paste0("Average Residuals by validation stations ",

  #p <- spplot(data_stations_var_pred,"constant",col.regions=NA, col="black",

  title_str <- paste0("Stations ",station_type_name," on ",date_processed, " for ",y_var_name," ", var_selected)

  plot(data_stations_var_pred,add=T,cex=0.5,col="blue")

  legend("topleft",legend=paste("n= ", nrow(data_stations_var_pred)),bty = "n")

  png_filename <- paste("Figure_ac_metrics_map_stations",station_type_name,"averaged","_fixed_intervals_",model_name,y_var_name,date_processed,

  class_intervals <- c(-20,-10,-5,-4,-3,-2,-1,0,1,2,3,4,5,10,20)

  p <- bubble(data_stations_var_pred,zcol="res_mod1",key.entries = class_intervals , 

              fill=F, #plot circle with filling

              #col= matlab.like(length(class_intervals)),

              main=title_str)

  p1 <- p + p_shp

  #### Add additional plot with quantiles and min max?...

  #library(classInt)

  #breaks.qt <- classIntervals(palo_alto$PrCpInc, n = 6, style = "quantile", intervalClosure = "right")

  #spplot(palo_alto, "PrCpInc", col = "transparent", col.regions = my.palette, 

  #  at = breaks.qt$brks)

  title_str <- paste0("Stations ",station_type_name," on ",date_processed, " for ",y_var_name," ", var_selected)

  h<- histogram(data_stations_var_pred$res_mod1,breaks=seq(-50,50,5),

                main=title_str)

  df_basic_stat <- fav_stats(data_stations_var_pred$res_mod1)

  df_basic_stat$date <- date_processed

  #df_basic_stat$reg <- reg

  #quantile(data_stations_var_pred$res_mod1,c(1,5,10,90,95,99))

  df_quantile_val <- quantile(data_stations_var_pred$res_mod1,c(0.01,0.05,0.10,0.45,0.50,0.55,0.90,0.95,0.99))

  #names(df_quantile_val)

  #as.list(df_quantile_val)

  #df_test <- data.frame(names(df_quantile_val))[numeric(0), ]

  #quantile(c(1,5,10,90,95,99),data_stations_var_pred$res_mod1,)

  #rmse(data_stations_var_pred$res_mod1)

  plot_obj <- list(data_stations_var_pred,df_basic_stat,df_quantile_val,freq_tile_id)

  names(plot_obj) <- c("data_stations_var_pred","df_basic_stat","df_quantile_val","freq_tile_id")

  return(plot_obj)

#data_stations_var_pred <- aggregate(id ~ date, data = data_stations, min)

#data_stations_var_pred <- aggregate(id ~ x + y + date + dailyTmax + mod1 + res_mod1 , data = data_stations, min)

##Add tile id here...and check if data stations was training or testing.

data_stations_var_pred <- aggregate(id ~ x + y + date + dailyTmax,data = data_stations, min ) #+ mod1 + res_mod1 , data = data_stations, min)

dim(data_stations_var_pred)

#> dim(data_stations_var_pred)

data_stations_var_pred$date_str <- data_stations_var_pred$date

data_stations_var_pred$date <- as.Date(strptime(data_stations_var_pred$date_str,"%Y%m%d"))

#data_stations_var_pred2 <- aggregate(id ~ training,data = data_stations, sum ) #+ mod1 + res_mod1 , data = data_stations, min)

#data_stations_var_pred2 <- aggregate(date ~ training,data = data_stations, sum ) #+ mod1 + res_mod1 , data = data_stations, min)

df2 <- data_stations_var_pred

#df_selected2 <- data_stations_var_pred

data_stations$date <- as.Date(strptime(data_stations_var_pred$date_str,"%Y%m%d"))