database_covariates_preparation<-function(list_param_prep){

  #This function performs queries on the Postgres ghcnd database for stations matching the             

  #interpolation area. It requires 11 inputs:                                           

  # 10) covar_names: names of covariates used for the interpolation --may be removed later? (should be stored in the brick)

  #

  #The output is a list of four shapefile names produced by the function:

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

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

  #AUTHOR: Benoit Parmentier                                                                       

  #PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#363--     

  #Comments and TODO

  #-Add buffer option...

  infile3<-list_param_prep$infile_covariates #"covariates__venezuela_region__VE_01292013.tif" #this is an output from covariate script

  CRS_locs_WGS84<-list_param_prep$CRS_locs_WGS84 #Station coords WGS84: same as earlier

  in_path <- list_param_prep$in_path #CRS_locs_WGS84"/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/input_data/"

  covar_names<-list_param_prep$covar_names # names should be written in the tif file!!!

  ## working directory is the same for input and output for this function  

  setwd(in_path) 

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

  time_minutes<-time_duration[3]/60

  dbDisconnect(db)

  data_table<-merge(data2,as.data.frame(stat_reg), by.x = "station", by.y = "STAT_ID")

  #Transform the subset data frame in a spatial data frame and reproject

  data_reg<-data_table                               #Make a copy of the data frame

  coords<- data_reg[c('lon','lat')]              #Define coordinates in a data frame: clean up here!!

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

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

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

  data_d <-data_reg  #data_d: daily data containing the query without screening

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

  ### STEP 3: Save results and outuput in textfile and a shape file

  #browser()

  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)

  #writeOGR(data_proj, paste(outfile, "shp", sep="."), outfile, driver ="ESRI Shapefile") #Note that the layer name is the file name without extension

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

  ### STEP 4: Extract values at stations from covariates stack of raster images

  }

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

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

  ######## STEP 5: Preparing monthly averages from the ProstGres database

  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 *

  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

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

  #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

  ####

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

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

  ### list of outputs return

  return(outfiles_obj)

  #END OF FUNCTION #