library(rasterVis)

#tinfile2<-"list_10_dates_04212012.txt"                     #List of 10 dates for the regression

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

path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_07192012_GAM"

#path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_GAM"

nmodels<-8   #number of models running

y_var_name<-"dailyTmax"

predval<-1

out_prefix<-"_07242012_365d_GAM_fusion5"                                                   #User defined output prefix

source("GAM_fusion_function_07192012d.R")

ghcn <- transform(ghcn,Northness = cos(ASPECT*pi/180)) #Adding a variable to the dataframe

ghcn <- transform(ghcn,Eastness = sin(ASPECT*pi/180))  #adding variable to the dataframe.

ghcn <- transform(ghcn,Northness_w = sin(slope*pi/180)*cos(ASPECT*pi/180)) #Adding a variable to the dataframe

ghcn <- transform(ghcn,Eastness_w = sin(slope*pi/180)*sin(ASPECT*pi/180))  #adding variable to the dataframe.

##Extracting the variables values from the raster files                                             

lines<-read.table(paste(path,"/",inlistf,sep=""), sep=" ")                  #Column 1 contains the names of raster files

inlistvar<-lines[,1]

inlistvar<-paste(path,"/",as.character(inlistvar),sep="")

covar_names<-as.character(lines[,2])                                         #Column two contains short names for covaraites

s_raster<- stack(inlistvar)                                                  #Creating a stack of raster images from the list of variables.

layerNames(s_raster)<-covar_names                                            #Assigning names to the raster layers

projection(s_raster)<-CRS

#stat_val<- extract(s_raster, ghcn3)                                          #Extracting values from the raster stack for every point location in coords data frame.

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

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

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

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

N<-cos(r1*pi/180)

E<-sin(r1*pi/180)

Nw<-sin(r2*pi/180)*cos(r1*pi/180)   #Adding a variable to the dataframe

Ew<-sin(r2*pi/180)*sin(r1*pi/180)   #Adding variable to the dataframe.

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

xy<-coordinates(r1)  #get x and y projected coordinates...

xy_latlon<-project(xy, CRS, inv=TRUE) # find lat long for projected coordinats (or pixels...)

lon<-raster(xy_latlon) #Transform a matrix into a raster object ncol=ncol(r1), nrow=nrow(r1))

ncol(lon)<-ncol(r1)

nrow(lon)<-nrow(r1)

extent(lon)<-extent(r1)

projection(lon)<-CRS  #At this stage this is still an empty raster with 536 nrow and 745 ncell 

lat<-lon

values(lon)<-xy_latlon[,1]

values(lat)<-xy_latlon[,2]

r<-stack(N,E,Nw,Ew,lon,lat,LC1,LC3)

rnames<-c("Northness","Eastness","Northness_w","Eastness_w", "lon","lat","LC1","LC3")

layerNames(r)<-rnames

s_raster<-addLayer(s_raster, r)

#s_sgdf<-as(s_raster,"SpatialGridDataFrame") #Conversion to spatial grid data frame

####### 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 tables for model assessment: specific diagnostic/metrics

#Model assessment: specific diagnostics/metrics

results_AIC<- matrix(1,1,nmodels+3)  

results_GCV<- matrix(1,1,nmodels+3)

results_DEV<- matrix(1,1,nmodels+3)

results_RMSE <- matrix(1,1,nmodels+4)

results_MAE <- matrix(1,1,nmodels+4)

results_ME <- matrix(1,1,nmodels+4)       #There are 8+1 models

results_R2 <- matrix(1,1,nmodels+4)       #Coef. of determination for the validation dataset

results_RMSE_f<- matrix(1,1,nmodels+4)    #RMSE fit, RMSE for the training dataset

results_MAE_f <- matrix(1,1,nmodels+4)

######## Preparing monthly averages from the ProstGres database

# 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

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

##Sampling: training and testing sites...

######## Prediction for the range of dates

#source("GAM_fusion_function_07192012d.R")

gam_fus_mod<-mclapply(1:length(dates), runGAMFusion,mc.preschedule=FALSE,mc.cores = 8) #This is the end bracket from mclapply(...) statement

tb<-gam_fus_mod[[1]][[3]][0,]  #empty data frame with metric table structure that can be used in rbinding...

tb_tmp<-gam_fus_mod #copy

save(gam_fus_mod,file= paste(path,"/","results2_fusion_Assessment_measure_all",out_prefix,".RData",sep=""))