#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)

  #3) daily_covar_ghcn_data: ghcn daily data with covariates 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...

  year_start <-list_param_prep$range_years[1] #"2010"               #starting year for the query (included)

  year_end <-list_param_prep$range_years[2] #"2011"                 #end year for the query (excluded)

  year_start_clim <-list_param_prep$range_years_clim[1] #right bound not included in the range!! starting year for monthly query to calculate clime

  infile1<- list_param_prep$infile1  #This is the shape file of outline of the study area                                                      #It is an input/output of the covariate script

  infile2<-list_param_prep$infile2      #"/home/layers/data/climate/ghcn/v2.92-upd-2012052822/ghcnd-stations.txt"                              #This is the textfile of station locations from GHCND

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

  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

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

  ### 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

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

                               FROM ghcn

                               WHERE element=",shQuote(var),

                               "AND year>=",year_start_clim,

                               "AND station IN ",list_s,";",sep=""))  #Selecting station using a SQL query

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

  time_minutes<-time_duration[3]/60

  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

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

  #d<-subset(data_m,year>=2000 & mflag=="0" ) #Selecting dataset 2000-2010 with good quality: 193 stations

  d<-subset(data_m,mflag=="0" | mflag=="S") #should be input arguments!!

  #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

  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

  ####

  #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 # 

}