##################    MULTI SAMPLING GAM FUSION METHOD ASSESSMENT ####################################

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

out_path<-"/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/output_data"

script_path<-"/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/"

setwd(in_path)

y_var_name<-"dailyTmax"                                       

out_prefix<-"_365d_GAM_fus5_all_lstd_02202013"                #User defined output prefix

#projection used in the interpolation of the study area: should be read directly from the outline 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

data3<-readOGR(dsn=in_path,layer=sub(".shp","",infile_monthly))

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

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

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

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

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

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

##Extracting the variables values from the raster files                                             

#The names of covariates can be changed...

                    "nobs_09","nobs_10","nobs_11","nobs_12")

covar_names<-c(rnames,lc_names,lst_names)

s_raster<-stack(infile3)                   #read in the data stack

#Deal with no data value and zero      

#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

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

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

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

#ghcn$value<-as.numeric(ghcn$value)

#ghcn_all<-ghcn

#ghcn_test<-subset(ghcn,ghcn$value>-150 & ghcn$value<400)

#ghcn_test<-ghcn

#ghcn_test2<-subset(ghcn_test,ghcn_test$elev_1>0)

#ghcn<-ghcn_test2

#coords<- ghcn[,c('x_OR83M','y_OR83M')]

#Now clean and screen monthly values

#dst_all<-dst

#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

##### Create sampling: select training and testing sites ###