library(rgeos)                                          # Polygon buffering and other vector operations

var <- "PRCP"                    #name of the variables to keep: TMIN, TMAX or PRCP

buffer=100                      #  size of buffer around shapefile to include in spatial query (in km), set to 0 to use no buffer

outpath=path                                                              # create different output path because we don't have write access to other's home dirs

#outpath="/home/wilson/data"

out_prefix<-"y2010_2010_OR_20110705"                                                 #User defined output prefix

#####  Connect to Station database

drv <- dbDriver("PostgreSQL")

db <- dbConnect(drv, dbname=db.name)#,options="statement_timeout = 1m")

#####  Buffer shapefile if desired

##     This is done to include stations from outside the region in the interpolation fitting process and reduce edge effects when stiching regions

if(buffer>0){  #only apply buffer if buffer >0

  interp_area=gUnionCascaded(interp_area)  #dissolve any subparts of roi (if there are islands, lakes, etc.)

  interp_areaC=gCentroid(interp_area)       # get centroid of region

  interp_areaB=spTransform(                # buffer roi (transform to azimuthal equidistant with centroid of region for most (?) accurate buffering, add buffer, then transform to WGS84)

    gBuffer(

      spTransform(interp_area,

                  CRS(paste("+proj=aeqd +lat_0=",interp_areaC@coords[2]," +lon_0=",interp_areaC@coords[1]," +ellps=WGS84 +datum=WGS84 +units=m +no_defs ",sep=""))),

      width=buffer*1000),                  # convert buffer (km) to meters

    CRS(CRS_interp))                       # reproject back to original projection

  interp_area=interp_areaB                 # replace original region with buffered region

}

## get bounding box of study area

bbox=bbox(interp_area)

### read in station location information from database

### use the bbox of the region to include only station in rectangular region to speed up overlay

dat_stat=dbGetQuery(db, paste("SELECT id,name,latitude,longitude

                  FROM stations

                  WHERE latitude>=",bbox[2,1]," AND latitude<=",bbox[2,2],"

                  AND longitude>=",bbox[1,1]," AND longitude<=",bbox[1,2],"

                  ;",sep=""))

coordinates(dat_stat)<-c("longitude","latitude")

##  Query to link station location information and observations

##  Concatenate date columns into single field for easy convert to date

##  Divide value by 10 to convert to degrees C and mm

##  Subset to years in year_start -> year_end

##  Drop missing values (-9999)

##  Drop observations that failed quality control (keep only qflag==NA)

system.time(

  data_table<<-dbGetQuery(db,

                          paste("SELECT station,year||'-'||month||'-'||day AS date,value / 10.0 as value,latitude,longitude,elevation

                                 FROM ghcn, stations

                                 WHERE station = id

                                 AND id IN ('",paste(stat_OR$id,collapse="','"),"')

                                 AND element='",var,"'

                                 AND year>=",year_start,"

                                 AND year<",year_end,"

                                 AND value<>-9999

                                 AND qflag IS NULL

                                 ;",sep=""))

  )  ### print used time in seconds  - only taking ~30 seconds....

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

write.table(as.data.frame(stat_OR), file=paste(outpath,"/","location_study_area_",out_prefix,".txt",sep=""),sep=",")

write.table(data_table, file= paste(outpath,"/","ghcn_data_",var,out_prefix,".txt",sep=""), sep=",")

outfile<-paste(outpath,"/ghcn_data_",var,out_prefix,sep="")         #Name of the file