#This script interpolates tmax values using MODIS LST and GHCND station data                      

#interpolation area. It requires the text file of stations and a shape file of the study area.           

#Note that the projection for both GHCND and study area is lonlat WGS84.       

#Options to run this program are:

#1) Multisampling: vary the porportions of hold out and use random samples for each run

#2)Constant sampling: use the same sample over the runs

#3)over dates: run over for example 365 dates without mulitsampling

#4)use seed number: use seed if random samples must be repeatable

#5)GAM fusion: possibilty of running GAM+FUSION or GAM separately 

#AUTHOR: Benoit Parmentier                                                                        

#DATE: 10/27/2012                                                                                 

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

library(reshape)

library(plotrix)

infile2<-"list_365_dates_04212012.txt"

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

prop<-0.3             #Proportion of testing retained for validation   

seed_number<- 100  #if seed zero then no seed?                                                                 #Seed number for random sampling

out_prefix<-"_365d_GAM_fusion_const_all_lstd_10282012"                #User defined output prefix

bias_prediction<-1     #if value 1 then use GAM for the BIAS prediction otherwise GAM direct repdiction for y_var (daily tmax)

nb_sample<-1           #number of time random sampling must be repeated for every hold out proportion

prop_min<-0.3          #if prop_min=prop_max and step=0 then predicitons are done for the number of dates...

prop_max<-0.3

step<-0         

constant<-1             #if value 1 then use the same samples as date one for the all set of dates

#projection used in the interpolation of the study area

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";

CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84

source("GAM_fusion_function_multisampling_10272012.R")

###################### START OF THE SCRIPT ########################

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

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

s_raster<-dropLayer(s_raster,pos)

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

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

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

s_raster<-dropLayer(s_raster,pos)

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

LC_s<-stack(LC1,LC3,LC4,LC6)

layerNames(LC_s)<-c("LC1_forest","LC3_grass","LC4_crop","LC6_urban")

plot(LC_s)

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

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

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

#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

if (seed_number>0) {

  set.seed(seed_number)                        #Using a seed number allow results based on random number to be compared...

}

prop_range<-(seq(from=prop_min,to=prop_max,by=step))*100     #range of proportion to run

sn<-length(dates)*nb_sample*length(prop_range)               #Number of samples to run

## adding choice of constant sample 

if (seed_number>0) {

  set.seed(seed_number)                        #Using a seed number allow results based on random number to be compared...

}

sampling<-vector("list",length(ghcn.subsets))

sampling_station_id<-vector("list",length(ghcn.subsets))

  #Find the corresponding 

  data_sampled<-ghcn.subsets[[i]][ind.training,] #selected the randomly sampled stations

  station_id.training<-data_sampled$station     #selected id for the randomly sampled stations (115)

  #Save the information

  sampling_station_id[[i]]<- station_id.training

}

## Use same samples across the year...

if (constant==1){

  sampled<-sampling[[1]]

  data_sampled<-ghcn.subsets[[1]][sampled,] #selected the randomly sampled stations

  station_sampled<-data_sampled$station     #selected id for the randomly sampled stations (115)

  list_const_sampling<-vector("list",sn)

  list_const_sampling_station_id<-vector("list",sn)

  for(i in 1:sn){

    station_id.training<-intersect(station_sampled,ghcn.subsets[[i]]$station)

    ind.training<-match(station_id.training,ghcn.subsets[[i]]$station)

    list_const_sampling[[i]]<-ind.training

    list_const_sampling_station_id[[i]]<-station_id.training

  }

  sampling<-list_const_sampling 

  sampling_station_id<-list_const_sampling_station_id

# sampling<-vector("list",length(ghcn.subsets))

# 

# for(i in 1:length(ghcn.subsets)){

#   n<-nrow(ghcn.subsets[[i]])

#   prop<-(sampling_dat$prop[i])/100

#   ns<-n-round(n*prop)   #Create a sample from the data frame with 70% of the rows

#   nv<-n-ns              #create a sample for validation with prop of the rows

#   ind.training <- sample(nrow(ghcn.subsets[[i]]), size=ns, replace=FALSE) #This selects the index position for 70% of the rows taken randomly

#   ind.testing <- setdiff(1:nrow(ghcn.subsets[[i]]), ind.training)

#   sampling[[i]]<-ind.training

# }

# 

# if (constant==1){

#   sampled<-sampling[[1]]

#   list_const_sampling<-vector("list",sn)

#   for(i in 1:sn){

#     list_const_sampling[[i]]<-sampled

#   }

#   sampling<-list_const_sampling  

# }

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

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

#gam_fus_mod2<-mclapply(4:4, runGAMFusion,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement

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

metrics<-as.character(unique(tb$metric))            #Name of accuracy metrics (RMSE,MAE etc.)

for(i in 1:length(metrics)){            # Reorganizing information in terms of metrics 

#tb_diagnostic[["prop"]]<-as.factor(tb_diagnostic[["prop"]])

sampling_obj<-list(sampling_dat=sampling_dat,training=sampling, training_id=sampling_station_id, tb=tb_diagnostic)

#write.table(avg_tb, file= paste(path,"/","results2_fusion_Assessment_measure_avg_",out_prefix,".txt",sep=""), sep=",")

#write.table(median_tb, file= paste(path,"/","results2_fusion_Assessment_measure_median_",out_prefix,".txt",sep=""), sep=",")

save(sampling_obj, file= paste(path,"/","results2_fusion_sampling_obj",out_prefix,".RData",sep=""))

## Summary of number of files and results

#l_f<-list.files(pattern="GAM_bias_tmax_predicted_mod1_20100103_30_1_365d_GAM_fusion_lstd_10062012.rst$")

#r6<-raster(l_f)

#plot(r6)

#results_list_metrics_objects_*_20101231_30_1_365d_GAM_fusion_lstd_10062012.RData

#l_f<-list.files(pattern=paste(".*","tmax_predicted_mod1","*outprefix,".rst",sep="""))