# 5)the location of raser covariate stack.                                                                                             

#The outputs are text files and a shape file of a time subset of the database                    

#AUTHOR: Benoit Parmentier                                                                       

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

#Comments and TODO

#-Add buffer option...

#-Add calculation of monthly mean...

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

###Loading R library and packages   

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

year_start<-"2010"               #starting year for the query (included)

year_end<-"2011"                 #end year for the query (excluded)

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

new_proj<-"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"

##Paths to inputs and output

in_path <- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/input_data/"

setwd(in_path) 

### Functions used in the script

interp_area <- readOGR(".",filename)

CRS_interp<-proj4string(interp_area)         #Storing the coordinate information: geographic coordinates longlat WGS84

colnames(dat_stat)<-c("STAT_ID","lat","lon","elev","state","name","GSNF","HCNF","WMOID")

coords<- dat_stat[,c('lon','lat')]

coordinates(dat_stat)<-coords

#proj4string(dat_stat)<-CRS_interp

dat_stat2<-spTransform(dat_stat,CRS(new_proj))         # Project from WGS84 to new coord. system

#only 357 station for Venezuela??

####

####

#### STEP 2: Connecting to the database and query for relevant data 

                                                   #Wrong label...it is in fact projected...

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

#proj4string(data3)<-locs_coord                  #Assign coordinates reference system in PROJ4 format

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

plot(interp_area, axes =TRUE)

plot(stat_reg, pch=1, col="red", cex= 0.7, add=TRUE)

plot(data_reg,pch=2,col="blue",cex=2,add=TRUE)

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

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

#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

#Eventually this step may be skipped if the covariates information is stored in the database...

rnames<-c("x","y","lon","lat","N","E","N_w","E_w","elev","slope","aspect","CANHEIGHT","DISTOC")

lc_names<-c("LC1","LC2","LC3","LC4","LC5","LC6","LC7","LC8","LC9","LC10","LC11","LC12")

lst_names<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12",

}

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

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

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

drv <- dbDriver("PostgreSQL")

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

year_end<-2011

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

list_s<-format_s(stat_reg$STAT_ID)

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

                            FROM ghcn

                            WHERE element=",shQuote(var),

                            "AND year<",year_end,

                            "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

####

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

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

##### END OF SCRIPT ##########