#DATE: 02/06/2013                                                                                 

infile1<- "ghcn_or_tmax_covariates_06262012_OR83M.shp"             #GHCN shapefile containing variables for modeling 2010                 

#infile3<-"LST_dates_var_names.txt"                        #LST dates name

#infile4<-"models_interpolation_05142012.txt"              #Interpolation model names

#infile5<-"mean_day244_rescaled.rst"                       #Raster or grid for the locations of predictions

#infile6<-"lst_climatology.txt"

#infile6<-"LST_files_monthly_climatology.txt"

#inlistf<-"list_files_05032012.txt"                        #Stack of images containing the Covariates

infile_monthly<-"monthly_covariates_ghcn_data_TMAXy2010_2010_VE_02062013.shp"

infile_daily<-"daily_covariates_ghcn_data_TMAXy2010_2010_VE_02062013.shp"

infile_locs<-"stations_venezuela_region_y2010_2010_VE_02062013.shp"

prop<-0.3             #Proportion of testing retained for validation   

#prop<-0.25

out_prefix<-"_10d_GAM_fus5_all_lstd_020632013"                #User defined output prefix

#out_prefix<-"_365d_GAM_12272012"                #User defined output prefix

CRS_interp<-"+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs";

source("GAM_fusion_function_multisampling_02062013.R")

#filename<-sub(".shp","",infile1)             #Removing the extension from file.

#ghcn<-readOGR(".", filename)                 #reading shapefile 

#mean_LST<- readGDAL(infile5)                 #Reading the whole raster in memory. This provides a grid for kriging

#proj4string(mean_LST)<-CRS_interp                   #Assigning coordinate information to prediction grid.

#Station location of the study area

#stat_loc<-read.table(paste(path,"/","location_study_area_OR_0602012.txt",sep=""),sep=",", header=TRUE)

#GHCN Database for 1980-2010 for study area (OR) 

#data3<-read.table(paste(path,"/","ghcn_data_TMAXy1980_2010_OR_0602012.txt",sep=""),sep=",", header=TRUE)

#data3<-file.path(in_path,infile_monthly)

#Remove NA for LC and CANHEIGHT: Need to check this part

dates <-readLines(paste(file.path(in_path,infile2))

#LST_dates <-readLines(file.path(in_path,infile3))

#models <-readLines(paste(path,"/",infile4, sep=""))

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

pos<-match("LC1",layerNames(s_raster)) #Find column with name "value"

LC1<-raster(s_raster,layer=pos)             #Select layer from stack

s_raster<-dropLayer(s_raster,pos)

LC1[is.na(LC1)]<-0

pos<-match("LC3",layerNames(s_raster)) #Find column with name "value"

LC3<-raster(s_raster,layer=pos)             #Select layer from stack

s_raster<-dropLayer(s_raster,pos)

LC3[is.na(LC3)]<-0

pos<-match("CANHEIGHT",layerNames(s_raster)) #Find column with name "value"

CANHEIGHT<-raster(s_raster,layer=pos)             #Select layer from stack

s_raster<-dropLayer(s_raster,pos)

CANHEIGHT[is.na(CANHEIGHT)]<-0

pos<-match("ELEV_SRTM",layerNames(s_raster)) #Find column with name "ELEV_SRTM"

ELEV_SRTM<-raster(s_raster,layer=pos)             #Select layer from stack on 10/30

s_raster<-dropLayer(s_raster,pos)

ELEV_SRTM[ELEV_SRTM <0]<-NA

####### Preparing LST stack of climatology...

#l=list.files(pattern="mean_month.*rescaled.rst")

l <-readLines(paste(path,"/",infile6, sep=""))

molst<-stack(l)  #Creating a raster stack...

#setwd(old)

molst<-molst-273.16  #K->C          #LST stack of monthly average...

idx <- seq(as.Date('2010-01-15'), as.Date('2010-12-15'), 'month')

molst <- setZ(molst, idx)

layerNames(molst) <- month.abb

######## Preparing monthly averages from the ProstGres database and extracting covarvariates from stack

# do this work outside of (before) this function

# to avoid making a copy of the data frame inside the function call

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

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

# data3$date<-date2

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

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

# id<-as.data.frame(unique(d1$station))     #Unique station in OR for year 2000-2010: 193 but 7 loss of monthly avg    

# 

# dst<-merge(d1, stat_loc, 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"

# names(dst)[pos]<-c("TMax")

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

# #dstjan=dst[dst$month==9,]  #dst contains the monthly averages for tmax for every station over 2000-2010

# 

# #Extracting covariates from stack

# coords<- dst[c('lon','lat')]              #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)  #extraction of the infomration at station location

# stations_val<-as.data.frame(stations_val)

# dst_extract<-cbind(dst_month,stations_val)

# dst<-dst_extract

#Now clean and screen monthly values

#dst_all<-dst

dst_all<-data3

dst<-data3

#dst<-subset(dst,dst$TMax>-15 & dst$TMax<45) #may choose different threshold??

#dst<-subset(dst,dst$ELEV_SRTM>0) #This will drop two stations...or 24 rows

######### Preparing daily values for training and testing

#Screening for bad values: value is tmax in this case