#Method is assedsed using constant sampling with variation  of validation sample with different  #

#hold out proportions.                                                                           #

#AUTHOR: Benoit Parmentier                                                                       #

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

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

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

y_var_name<-"dailyTmax"

climgam=1                                                     #if 1, then GAM is run on the climatology rather than the daily deviation surface...

predval<-1

seed_number<- 100                                             #Seed number for random sampling, if seed_number<0, no seed number is used..

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

bias_val<-0            #if value 1 then daily training data is used in the bias surface rather than the all monthly stations (added on 07/11/2012)

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

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

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

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

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

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

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

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

values(lon)<-xy_latlon[,1]

values(lat)<-xy_latlon[,2]

layerNames(r)<-rnames

s_raster<-addLayer(s_raster, r)

stations_val<-as.data.frame(stations_val)

dst_extract<-cbind(dst_month,stations_val)

dst<-dst_extract

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

#ghcn.subsets <-lapply(dates, function(d) subset(ghcn, date==d)) #this creates a list of 10 or 365 subsets dataset based on dates

ghcn.subsets <-lapply(as.character(sampling_dat$date), function(d) subset(ghcn, date==d)) #this creates a list of 10 or 365 subsets dataset based on dates

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

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

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

  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){

  }

}

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

# Save before plotting

#sampling_obj<-list(sampling_dat=sampling_dat,training=sampling)

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=",")

write.table(tb_diagnostic, file= paste(path,"/","results2_fusion_Assessment_measure",out_prefix,".txt",sep=""), sep=",")

write.table(tb, file= paste(path,"/","results2_fusion_Assessment_measure_all",out_prefix,".txt",sep=""), sep=",")

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

save(gam_CAI_mod,file= paste(path,"/","results2_CAI_Assessment_measure_all",out_prefix,".RData",sep=""))