##################    Data preparation for interpolation  #######################################

############################ Extraction of station data##########################################

  # 14)

  #1) loc_stations: outfile1, locations of stations as shapefile in EPSG 4326 

  #2) loc_stations_ghcn: outfile2, ghcn daily data for the year range of interpolation (locally projected) 

  #3) daily_query_ghcn_data: outfile3, ghcn daily data from daily query before application of quality flag

  #4) daily_covar_ghcn_data: outfile4, ghcn daily data with covariates for the year range of interpolation (locally projected)

  #5) monthly_query_ghcn_data: outfile5, ghcn daily data from monthly query before application of quality flag

  #6) monthly_covar_ghcn_data: outfile6, ghcn monthly averaged data with covariates for the year range of interpolation (locally projected)

  #CREATED ON: 02/22/2013

  #MODIFIED ON: 12/21/2015

  #Comments:

  #Modifications to take into account the existence of clim query, to avoid querying many times the database at the climatology step.

  #TODO

  #-Add buffer option: solved by overlapping tiles

  # 

  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: create the output dir as defined 

    if(!file.exists(out_dir)){

      dir.create(out_dir)

    }

    return(out_dir)

  }

  ## New parameters added for  sub samplineg in areas with important density of meteo stations

  sub_sampling <- list_param_prep$sub_sampling  #if TRUE then monthly stations data are resampled

  sub_sample_rnd <- list_param_prep$sub_sample_rnd #if  TRUE use random sampling  in addition to spatial  sub-sampling

  target_range_nb <- list_param_prep$target_range_nb # number of stations desired as min and max, convergence to  min  for  now

  #needs to be added in master script!!	

  sub_sampling_day <- list_param_prep$sub_sampling_day

  target_range_daily_nb <- list_param_prep$target_range_daily_nb #desired number range of daily stations

  dist_range <- list_param_prep$dist_range #distance range  for pruning,  usually (0,5) in km or 0,0.009*5 for  degreee

  step_dist <- list_param_prep$step_dist #stepping distance used in pruning  spatially, use 1km or 0.009 for degree data

  ##### Create additional output dir for clim datasets

  out_path_clim <- file.path(out_path,paste("clim_",year_start_clim,"_",year_end_clim,sep=""))

  create_dir_fun(out_dir=out_path_clim,out_suffix=NULL) #create dir if it does not exists

  setwd(out_path) #thisvaries for clim datasets in the code

  out_path <- file.path(out_path,year_start) #if predictions for year 1992, then a 1992 dir is created

  create_dir_fun(out_dir=out_path,out_suffix=NULL) #create dir if it does not exists

  setwd(out_path) #thisvaries for clim datasets in the code

  ##TODO: Add buffer option? Not needed right now as tiles overlap

  ##This is the query for daily time step predictions

  #Also save in clim data

  outfile1_clim <-file.path(out_path_clim,paste("stations","_",out_prefix,".shp",sep=""))

  writeOGR(stat_reg,dsn= dirname(outfile1),layer= sub(".shp","",basename(outfile1)), driver="ESRI Shapefile",overwrite_layer=TRUE)

  #writeOGR(dst,dsn= ".",layer= sub(".shp","",outfile4), driver="ESRI Shapefile",overwrite_layer=TRUE)

  #Also save as rds

  outfile1_rds<-sub(".shp",".rds",outfile1)

  saveRDS(stat_reg,outfile1)

  ##This goes in the daily output dir (e.g. for year 1992)

  #/nobackupp6/aguzman4/climateLayers/out/reg4/20_-60/1992/

   ## Adding subsampling for daily stations...

    print("daily subsample")

    sub_sampling_obj <- sub_sampling_by_dist_nb_stat(target_range_nb=target_range_nb,dist_range=dist_range,step_dist=step_dist,data_in=data_RST_SDF,sampling=T,combined=F)

  ## Only query for monthly climatology data if data does not exist in the directory

  #This file contains the input data for climatology predictions

  outfile6 <- file.path(out_path_clim,paste("monthly_covariates_ghcn_data_",var,"_",year_start_clim,"_",year_end_clim,out_prefix,".shp",sep=""))  #Name of the file

  #Could also be replaced by run_clim_query parameter.

  if(!(file.exists(outfile6))){

    drv <- dbDriver("PostgreSQL")

    db <- dbConnect(drv, dbname=db.name)

    #year_start_clim: set at the start of the script

    time1<-proc.time()    #Start stop watch

    list_s<-format_s(stat_reg$STAT_ID)

    data_m <- dbGetQuery(db, paste("SELECT *

    time_duration<-proc.time()-time1             #Time for the query may be long given the size of the database

    time_minutes<-time_duration[3]/60

    print(time_minutes)

    dbDisconnect(db)

    #Clean out this section!!

    date1<-ISOdate(data_m$year,data_m$month,data_m$day) #Creating a date object from 3 separate column

    date2<-as.POSIXlt(as.Date(date1))

    data_m$date<-date2

    #Save the query data here...

    data_m<-merge(data_m, stat_reg, by.x="station", by.y="STAT_ID")   #Inner join all columns are retained

    #Extracting covariates from stack for the monthly dataset...

    coords<- data_m[c('longitude','latitude')]              #Define coordinates in a data frame

    coordinates(data_m)<-coords                      #Assign coordinates to the data frame

    proj4string(data_m)<-CRS_locs_WGS84                  #Assign coordinates reference system in PROJ4 format

    data_m<-spTransform(data_m,CRS(CRS_interp))     #Project from WGS84 to new coord. system

    outfile5<-file.path(out_path_clim,paste("monthly_query_ghcn_data_",var,"_",year_start_clim,"_",year_end_clim,out_prefix,".shp",sep=""))  #Name of the file

    #outfile5<-file.path(out_path,paste("monthly_query_ghcn_data_",var,"_",year_start_clim,"_",year_end_clim,out_prefix,".shp",sep=""))  #Name of the file

    writeOGR(data_m,dsn= dirname(outfile5),layer= sub(".shp","",basename(outfile5)), driver="ESRI Shapefile",overwrite_layer=TRUE)

    outfile5_rds<-sub(".shp",".rds",outfile5)

    saveRDS(data_m,outfile5_rds) 

    #In Venezuela and other regions where there are not many stations...mflag==S should be added..see Durenne etal.2010.

    #d<-subset(data_m,mflag==qc_flags_stations[1] | mflag==qc_flags_stations[2])

    d<-subset(data_m,mflag %in% qc_flags_stations) 

    #Add screening here ...May need some screeing??? i.e. range of temp and elevation...

    d1<-aggregate(value~station+month, data=d, mean)  #Calculate monthly mean for every station in OR

    #d2<-aggregate(value~station+month, data=d, length)  #Calculate monthly mean for every station in OR

    is_not_na_fun<-function(x) sum(!is.na(x)) #count the number of available observation

    d2<-aggregate(value~station+month, data=d, is_not_na_fun)  #Calculate monthly mean for every station in OR

    #names(d2)[names(d2)=="value"] <-"nobs_station"

    d1$nobs_station<-d2$value

    dst<-merge(d1, stat_reg, by.x="station", by.y="STAT_ID")   #Inner join all columns are retained

    #This allows to change only one name of the data.frame

    pos<-match("value",names(dst)) #Find column with name "value"

    if (var=="TMAX"){

      names(dst)[pos]<-c("TMax")           #renaming the column "value" extracted from the Postgres database

      dst$TMax<-dst$TMax/10                #TMax is the average max temp for monthy data

    }

    if (var=="TMIN"){

      names(dst)[pos]<-c("TMin")

      dst$TMin<-dst$TMin/10                #TMin is the average min temp for monthy data

    }

    #Extracting covariates from stack for the monthly dataset...

    #names(dst)[5:6] <-c('latitude','longitude')

    coords<- dst[c('longitude','latitude')]              #Define coordinates in a data frame

    coordinates(dst)<-coords                      #Assign coordinates to the data frame

    proj4string(dst)<-CRS_locs_WGS84                  #Assign coordinates reference system in PROJ4 format

    dst_month<-spTransform(dst,CRS(CRS_interp))     #Project from WGS84 to new coord. system

    stations_val<-extract(s_raster,dst_month,df=TRUE)  #extraction of the information at station location in a data frame

    #dst_extract<-spCbind(dst_month,stations_val) #this is in sinusoidal from the raster stack

    dst_extract<-cbind(dst_month,stations_val) #this is in sinusoidal from the raster stack

    dst<-dst_extract #problem!!! two column named elev!!! use elev_s??

    #browser()

    coords<- dst[c('x','y')]              #Define coordinates in a data frame, this is the local x,y

    index<-!is.na(coords$x)               #remove if NA, may need to revisit at a later stage

    dst<-dst[index,]

    coords<- dst[c('x','y')]              #Define coordinates in a data frame, this is the local x,y

    coordinates(dst)<-coords                    #Assign coordinates to the data frame

    proj4string(dst)<-CRS_interp        #Assign coordinates reference system in PROJ4 format

    ##Added 01-07-2015

    dst$id <- dst$station #the id field is needed for possible subsampling

    ### ADD SCREENING HERE BEFORE WRITING OUT DATA

    #Covariates ok since screening done in covariate script

    #screening on var i.e. value, TMIN, TMAX...

    #### Adding  subsampling for regions  that  have  too  many stations...

    #This must be set up in master script

    #target_max_nb <- 100,000 #this is not actually used yet in the current implementation,can be set to very high value...

    #target_min_nb <- 8,000 #this is the target number of stations we would like: to be set by Alberto...

    ##max_dist <- 1000 # the maximum distance used for pruning ie removes stations that are closer than 1000m, this in degree...? 

    #max_dist <- 0.009*5 #5km in degree

    #min_dist <- 0    #minimum distance to start with

    #step_dist <- 0.009 #iteration step to remove the stations

    #test5 <- sub_sampling_by_dist_nb_stat(target_range_nb=target_range_nb,dist_range=dist_range,step_dist=step_dist,data_in=data_month,sampling=T,combined=F)

    #note that  this is for monthly stations.

    if(sub_sampling==TRUE){

      print("monthly subsample")

      sub_sampling_obj <- sub_sampling_by_dist_nb_stat(target_range_nb=target_range_nb,dist_range=dist_range,step_dist=step_dist,data_in=dst,sampling=T,combined=F)

      dst <- sub_sampling_obj$data #get sub-sampled data...for monhtly stations

      #save the information for later use (validation at monthly step!!)

      #save(sub_sampling_obj,file= file.path(out_path,paste("sub_sampling_obj_",interpolation_method,"_", out_prefix,".RData",sep="")))

      save(sub_sampling_obj,file= file.path(out_path_clim,paste("sub_sampling_obj_",interpolation_method,"_", out_prefix,".RData",sep="")))

    }

    #Make sure this is still a shapefile...!! This might need to be uncommented...

    #coordinates(dst)<-cbind(dst$x,dst$y) #Transforming data_RST_SDF into a spatial point dataframe

    #CRS_reg<-proj4string(data_reg)

    #proj4string(dst)<-CRS_reg  #Need to assign coordinates...

    ####

    #write out a new shapefile (including .prj component)

    dst$OID<-1:nrow(dst) #need a unique ID?

    #outfile6<-file.path(out_path,paste("monthly_covariates_ghcn_data_",var,"_",year_start_clim,"_",year_end_clim,out_prefix,".shp",sep=""))  #Name of the file

    outfile6<-file.path(out_path_clim,paste("monthly_covariates_ghcn_data_",var,"_",year_start_clim,"_",year_end_clim,out_prefix,".shp",sep=""))  #Name of the file

    writeOGR(dst,dsn= dirname(outfile6),layer= sub(".shp","",basename(outfile6)), driver="ESRI Shapefile",overwrite_layer=TRUE)

    outfile6_rds<-sub(".shp",".rds",outfile6)

    saveRDS(dst,outfile6_rds)

  }

  ###If clim exists, just create name for outfile5 to object used later...

  #outfile6 <- file.path(out_path_clim,paste("monthly_covariates_ghcn_data_",var,"_",year_start_clim,"_",year_end_clim,out_prefix,".shp",sep=""))  #Name of the file

  #Could also be replaced by run_clim_query parameter.

  if(file.exists(outfile6)){

    outfile5<-file.path(out_path_clim,paste("monthly_query_ghcn_data_",var,"_",year_start_clim,"_",year_end_clim,out_prefix,".shp",sep=""))  #Name of the file

  }