#DATE: 10/04/2013                                                                                 

  #29)list_models2: models' formulas as string vector for daily devation

  #30)interp_method2: intepolation method for daily devation step

  interp_method2 <- list_param_raster_prediction$interp_method2

    list_param_run_prediction_daily_deviation <-list(i,clim_yearlist,daily_dev_sampling_dat,sampling_month_obj,sampling_obj,dst,list_models2,interp_method2,

    names(list_param_run_prediction_daily_deviation)<-c("list_index","clim_yearlist","daily_dev_sampling_dat","sampling_month_obj","sampling_obj","dst","list_models2","interp_method2",

    method_mod_obj<-mclapply(1:nrow(daily_dev_sampling_dat),list_param=list_param_run_prediction_daily_deviation,run_prediction_daily_deviation,mc.preschedule=FALSE,mc.cores = 11) #This is the end bracket from mclapply(...) statement

    #test_val2 <-calculate_accuracy_metrics(1,list_param_validation_month)

    validation_mod_month_obj <- mclapply(1:length(clim_method_mod_obj), list_param=list_param_validation_month, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 6) 

#DATE: 10/03/2013                                                                                 

#Function to glue all methods together...still need to separate fit and training for gwr and kriging, ok for now

interpolate_area_fun <- function(method_interp,list_models,s_raster,list_out_filename,data_df){

  ##Function to fit and predict an interpolation surface

  ##Author: Benoit Parmentier

  ##Function depends on other functions!!!

  #inpputs:

  #method_interp: interpolation method with value "gam","gwr","kriging"

  #list_models: models to fit and predict as string (i.e.vector char)

  #s_raster: stack with covariate variables, must match in name the data.frame input

  #data_df: spatial point data.frame with covariates, must be projected match names of covariates

  #list_out_filename: list of char containing output names for models

  #Conver to formula object

  list_formulas<-lapply(list_models,as.formula,env=.GlobalEnv) #mulitple arguments passed to lapply!!

  cname<-paste("mod",1:length(list_formulas),sep="") #change to more meaningful name?

  names(list_out_filename)<-cname  

  ##Now carry out prediction

  if(method_interp=="gam"){    

    #First fitting

    mod_list<-fit_models(list_formulas,data_df) #only gam at this stage

    names(mod_list)<-cname

    #if raster provided then predict surface

    if(!is.null(s_raster)){

      #Second predict values for raster image...by providing fitted model list, raster brick and list of output file names

      rast_pred_list<-predict_raster_model(mod_list,s_raster,list_out_filename)

      names(rast_pred_list)<-cname

    }

  }

  if(method_interp%in%c("gwr","kriging")){

    #Call funciton to fit and predict gwr and/or kriging

    #month_prediction_obj<-predict_auto_krige_raster_model(list_formulas,s_raster,data_month,list_out_filename)

    rast_prediction_obj<-predict_autokrige_gwr_raster_model(method_interp,list_formulas,s_raster,data_df,list_out_filename)

    mod_list <-rast_prediction_obj$list_fitted_models

    rast_pred_list <-rast_prediction_obj$list_rast_pred

    names(rast_pred_list)<-cname

  }

  #Now prepare to return object

  interp_area_obj <-list(mod_list,list_formulas,rast_pred_list)

  names(interp_area_obj) <- c("mod_list","list_formulas","rast_pred_list")

  return(interp_area_obj)

} 

  #9)list_models2 : interpolation model's formulas as string

  #10)interp_methods2: "gam","gwr","kriging"

  #11)s_raster: stack for covariates and toher variables

  list_models2 <-list_param$list_models2

  interp_method2 <- list_param$interp_method2

  s_raster <- list_param$s_raster

  #Adding layer LST to the raster stack  

  #names(s_raster)<-covar_names

  pos<-match("LST",names(s_raster)) #Find the position of the layer with name "LST", if not present pos=NA

  s_raster<-dropLayer(s_raster,pos)      # If it exists drop layer

  LST<-subset(s_raster,LST_month)

  names(LST)<-"LST"

  s_raster<-addLayer(s_raster,LST)            #Adding current month

    if(is.null(list_models2)){ #change here...

      list_daily_delta_rast <- vector("list",length=1) #only one delta surface in this case!!

      list_mod_krtmp2 <- vector("list",length=1) #only one delta model in this case!!

      model_name<-paste("mod_stat_kr",sep="_")

      daily_delta_xy<-as.matrix(cbind(dmoday$x,dmoday$y))

      fitdelta<-Krig(daily_delta_xy,daily_delta,theta=1e5) #use TPS or krige

      mod_krtmp2 <- fitdelta

      #names(mod_krtmp2)[k] <- model_name

      #data_s$daily_delta<-daily_delta

      #rast_clim_list<-rast_clim_yearlist[[index_m]]  #select relevant monthly climatology image ...

      rast_clim_list<-rast_clim_yearlist[[index_m]]  #select relevant monthly climatology image ...

      rast_clim_mod <- stack(rast_clim_list)

      names(rast_clim_mod) <- names(rast_clim_list)

      rast_clim_month <- subset(rast_clim_mod,1) #example layer to interpolate to

      daily_delta_rast<-interpolate(rast_clim_month,fitdelta) #Interpolation of the of the daily devation

      #there is only one daily devation (delta) sruface in this case

      #To many I/O out of swap memory on atlas

      #Saving kriged surface in raster images

      data_name<-paste("daily_delta_",y_var_name,"_",model_name,"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",

                       sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],sep="")

      raster_name_delta<-file.path(out_path,paste(interpolation_method,"_",var,"_",data_name,out_prefix,".tif", sep=""))

      writeRaster(daily_delta_rast, filename=raster_name_delta,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)

      list_daily_delta_rast[[1]] <- raster_name_delta

      list_mod_krtmp2[[1]] <- mod_krtmp2

    }

    if(!is.null(list_models2)){ #change here...

      list_daily_delta_rast <- vector("list",length=1) #several delta surfaces in this case but stored as one list!!

      list_mod_krtmp2 <- vector("list",length=1) #several delta model in this case but stored as one list!!

      dev_mod_name<-paste("dev_mod",1:length(list_models2),sep="") #change to more meaningful name?

      model_name<-paste("mod_stat_",sep="_")

      #Now generate file names for the predictions...

      list_out_filename<-vector("list",length(list_models2))

      names(list_out_filename)<- dev_mod_name  

      ##Change name...

      for (j in 1:length(list_out_filename)){

        #j indicate which month is predicted, var indicates TMIN or TMAX

        data_name<-paste("daily_delta_",y_var_name,"_",model_name,"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",

                         sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],

                         "_",interp_method2,"_",dev_mod_name[j],sep="")

        raster_name_delta<-file.path(out_path,paste(interpolation_method,"_",var,"_",data_name,out_prefix,".tif", sep=""))

        list_out_filename[[j]]<-raster_name_delta

      }

      #Now call function

      #for (j in 1:length(list_models2)){

      dmoday$y_var <- daily_delta

      #coordinates(data_s)<-cbind(data_s$x,data_s$y)

      #proj4string(data_s)<-proj_str

      #coordinates(data_v)<-cbind(data_v$x,data_v$y)

      #proj4string(data_v)<-proj_str

      interp_area_obj <-interpolate_area_fun(interp_method2,list_models2,s_raster,list_out_filename,dmoday)

      rast_pred_list <- interp_area_obj$rast_pred_list

      rast_pred_list <-rast_pred_list[!sapply(rast_pred_list,is.null)] #remove NULL elements in list

      list_daily_delta_rast[[1]] <-rast_pred_list

      #names(list_daily_delta_rast) <- names(daily_delta_df)

      list_mod_krtmp2[[1]] <-interp_area_obj$mod_list      

    }

    dmoday <-spCbind(dmoday,daily_delta_df) #contains the predicted clim at locations

    names(list_daily_delta_rast) <- names(daily_delta_df)

    names(list_mod_krtmp2) <- names(daily_delta_df)

      daily_delta <- daily_delta_df[[k]] #Current daily deviation being process: the reference monthly prediction varies...

      if(is.null(list_models2)){

        daily_delta_xy<-as.matrix(cbind(dmoday$x,dmoday$y))

        fitdelta<-Krig(daily_delta_xy,daily_delta,theta=1e5) #use TPS or krige

        list_mod_krtmp2[[k]] <-fitdelta

        names(list_mod_krtmp2)[k] <- model_name

        #data_s$daily_delta<-daily_delta

        #rast_clim_list<-rast_clim_yearlist[[index_m]]  #select relevant monthly climatology image ...

        rast_clim_month <- subset(rast_clim_mod,1) #example layer to interpolate to

        daily_delta_rast<-interpolate(rast_clim_month,fitdelta) #Interpolation of the bias surface...

        #list_daily_delta_rast[[k]] <- raster_name_delta

        data_name<-paste("daily_delta_",y_var_name,"_",model_name,"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",

                         sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],sep="")

        raster_name_delta<-file.path(out_path,paste(interpolation_method,"_",var,"_",data_name,out_prefix,".tif", sep=""))

        writeRaster(delta_rast_s, filename=raster_name_delta,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)

        #writeRaster(r_spat, NAflag=NA_flag_val,filename=raster_name,bylayer=TRUE,bandorder="BSQ",overwrite=TRUE)   

        #raster_name_delta <- list_daily_delta_rast

        #mod_krtmp2 <- list_mod_krtmp2

        list_daily_delta_rast[[k]] <- raster_name_delta

      }

      if (!is.null(list_models2)){

        #list_formulas<-lapply(list_models,as.formula,env=.GlobalEnv) #mulitple arguments passed to lapply!!

        dev_mod_name<-paste("dev_mod",1:length(list_models2),sep="") #change to more meaningful name?

        #Now generate file names for the predictions...

        list_out_filename<-vector("list",length(list_models2))

        names(list_out_filename)<- dev_mod_name  

        ##Change name...

        for (j in 1:length(list_out_filename)){

          #j indicate which month is predicted, var indicates TMIN or TMAX

          data_name<-paste("daily_delta_",y_var_name,"_",model_name,"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",

                           sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],

                           "_",interp_method2,"_",dev_mod_name[j],sep="")

          raster_name_delta<-file.path(out_path,paste(interpolation_method,"_",var,"_",data_name,out_prefix,".tif", sep=""))

          list_out_filename[[j]]<-raster_name_delta

        }

        #Now call function

        #for (j in 1:length(list_models2)){

        dmoday$y_var <- daily_delta

        #coordinates(data_s)<-cbind(data_s$x,data_s$y)

        #proj4string(data_s)<-proj_str

        #coordinates(data_v)<-cbind(data_v$x,data_v$y)

        #proj4string(data_v)<-proj_str

        interp_area_obj <-interpolate_area_fun(interp_method2,list_models2,s_raster,list_out_filename,dmoday)

        rast_pred_list <- interp_area_obj$rast_pred_list

        names(rast_pred_list) <- dev_mod_name

        rast_pred_list <-rast_pred_list[!sapply(rast_pred_list,is.null)] #remove NULL elements in list

        list_daily_delta_rast[[k]] <-rast_pred_list

        names(list_daily_delta_rast) <- names(daily_delta_df)

        mod_list <-interp_area_obj$mod_list

        names(mod_list) <- dev_mod_name

        list_mod_krtmp2[[k]] <-interp_area_obj$mod_list

      }

    #delta_rast_s <-stack(list_daily_delta_rast)

    #names(delta_rast_s) <- names(daily_delta_df)

    #data_name<-paste("daily_delta_",y_var_name,"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",

    #                 sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],sep="")

    #raster_name_delta<-file.path(out_path,paste(interpolation_method,"_",var,"_",data_name,out_prefix,".tif", sep=""))

    #writeRaster(delta_rast_s, filename=raster_name_delta,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)

    #mod_krtmp2 <- list_mod_krtmp2

      if (is.null(list_models2)){ 

        daily_delta_rast <- raster(list_daily_delta_rast[[k]]) #There is only one image of deviation per model if list_models2 is NULL

      }

      if (!is.null(list_models2)){ #then possible multiple daily dev predictions

        daily_delta_rast <- stack(unlist(list_daily_delta_rast[[k]]))

      }

      #daily_delta_rast <- subset(delta_rast_s,k)

    if(use_clim_image==FALSE){

      if (is.null(list_models2)){ 

        daily_delta_rast <- raster(list_daily_delta_rast[[1]]) #There is only one image of deviation per model if list_models2 is NULL

      }

      if (!is.null(list_models2)){ #then possible multiple daily dev predictions hence use stack

        daily_delta_rast <- stack(unlist(list_daily_delta_rast[[1]]))

      }

      #daily_delta_rast <- subset(delta_rast_s,k)

    }

    rast_clim_month <- subset(rast_clim_mod,k) #long term monthly prediction

    if (is.null(list_models2)){ 

      temp_predicted<-rast_clim_month + daily_delta_rast

      data_name<-paste(y_var_name,"_predicted_",names(rast_clim_mod)[k],"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",

                       sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],sep="")

      raster_name<-file.path(out_path,paste(interpolation_method,"_",data_name,out_prefix,".tif", sep=""))

      writeRaster(temp_predicted, filename=raster_name,overwrite=TRUE) 

      temp_list[[k]]<-raster_name

    }

    if (!is.null(list_models2)){ 

      dev_mod_name <- paste("dev_mod",1:length(list_models2),sep="") #change to more meaningful name?

      raster_name_list <- vector("list",length(list_models2))

      for(j in 1:length(list_models2)){

        temp_predicted <- rast_clim_month + subset(daily_delta_rast,j)

        data_name<-paste(y_var_name,"_predicted_",names(rast_clim_mod)[k],"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",

                         sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],"_",

                         interp_method2,"_",dev_mod_name[j],sep="")

        raster_name<-file.path(out_path,paste(interpolation_method,"_",data_name,out_prefix,".tif", sep=""))

        writeRaster(temp_predicted, filename=raster_name,overwrite=TRUE) 

        raster_name_list[[j]] <- raster_name

      }

      names(raster_name_list) <- dev_mod_name

      temp_list[[k]]<-raster_name_list #record names of the daily predictions 

    }

  data_s <- dmoday #put the 

  model_name<-paste("mod_","day",sep="_")

  delta_obj<-list(temp_list,rast_clim_list,list_daily_delta_rast, data_s, data_v,

                 modst,data_month_v,sampling_dat[index_d,],daily_dev_sampling_dat[i,],list_mod_krtmp2)

#DATE: 10/04/2013                                                                                 

source(file.path(script_path,"GAM_fusion_analysis_raster_prediction_multisampling_10042013.R"))

source(file.path(script_path,"GAM_fusion_function_multisampling_10042013.R")) #Includes Fusion and CAI methods

source(file.path(script_path,"GAM_fusion_function_multisampling_validation_metrics_10042013.R"))

out_prefix<-"_365d_gam_cai_lst_comb3_10042013"                #User defined output prefix

out_suffix<-"_OR_10042013"                                       #Regional suffix

#list_models<-c("y_var ~ s(lat,lon)",

#              "y_var ~ s(lat,lon) + s(elev_s)",

#                "y_var ~ s(lat,lon) + s(N_w)",

#                "y_var ~ s(lat,lon) + s(E_w)",

#               "y_var ~ s(lat,lon) + s(LST)",

#               "y_var ~ s(lat,lon) + s(DISTOC)",

#              "y_var ~ s(lat,lon) + s(LC1)",

#             "y_var ~ s(lat,lon) + s(CANHGHT)",

#            "y_var ~ s(lat,lon) + s(LST) + ti(LST,LC1)",

#           "y_var ~ s(lat,lon) + s(LST) + ti(LST,CANHGHT)")

#Combination 5: for paper baseline=s(lat,lon)+s(elev)

list_models<-c("y_var ~ s(lat,lon) + s(elev_s)",

                "y_var ~ s(lat,lon) + s(elev_s) + s(N_w,E_w)",

                "y_var ~ s(lat,lon) + s(elev_s) + s(LST)",

                "y_var ~ s(lat,lon) + s(elev_s) + s(DISTOC)",

                "y_var ~ s(lat,lon) + s(elev_s) + s(LC1)",

                "y_var ~ s(lat,lon) + s(elev_s) + s(CANHGHT)",

                "y_var ~ s(lat,lon) + s(elev_s) + s(LST) + ti(LST,LC1)")#,

                #"y_var ~ s(lat,lon) + s(elev_s) + s(LST) + ti(LST,CANHGHT)")

list_models2<-c("y_var ~ s(lat,lon) + s(elev_s)",

                "y_var ~ s(lat,lon) + s(LST)")

interp_method2 <- "gam" #other options are "gwr" and "kriging"

                                list_models,list_models2,interp_method2,lst_avg,out_path,script_path,use_clim_image,join_daily,

                                "list_models","list_models2","interp_method2","lst_avg","out_path","script_path","use_clim_image","join_daily",

#debug(debug_fun_test)

#debug_fun_test(list_param_raster_prediction)