#DATE: 08/25/2013                                                                                 

  #inserted #

  sampling_month_obj<-list_param$sampling_month_obj

  ghcn.month.subsets<-sampling_month_obj$ghcn_data

  sampling_month_dat <- sampling_month_obj$sampling_dat

  sampling_month_index <- sampling_month_obj$sampling_index

  #prop_month<-0 #proportion retained for validation...

  #run_samp<-1 #sample number, can be introduced later...

  prop_month <- sampling_month_dat$prop[j] #proportion retained for validation...

  run_samp <- sampling_month_dat$run_samp[j] #sample number if multisampling...will need create mulitple prediction at daily!!! could be complicated

                                       #possibility is to average per proportion !!!

  date_month <-strptime(sampling_month_dat$date[j], "%Y%m%d")   # interpolation date being processed

  month_no <-strftime(date_month, "%m")          # current month of the date being processed

  LST_month<-paste("mm_",month_no,sep="") # name of LST month to be matched

  LST_name <-LST_month  

  #change here...use training data...

  ###Regression part 1: Creating a validation dataset by creating training and testing datasets

  #LST_name <-lst_avg[j] # name of LST month to be matched

  #data_month$LST<-data_month[[LST_name]]

  dataset_month <-ghcn.month.subsets[[j]]

  mod_LST <- ghcn.month.subsets[[j]][,match(LST_month, names(ghcn.month.subsets[[j]]))]  #Match interpolation date and monthly LST average

  dataset_month$LST <- as.data.frame(mod_LST)[,1] #Add the variable LST to the dataset

  #change here...

  dst$LST<-dst[[LST_month]] #Add the variable LST to the monthly dataset

  proj_str<-proj4string(dst) #get the local projection information from monthly data

    dataset_month$y_var<-dataset_month$TMax #Adding TMax as the variable modeled

    dataset_month$y_var<-dataset_month$TMin #Adding TMin as the variable modeled

  ind.training <- sampling_month_index[[j]]

  ind.testing  <- setdiff(1:nrow(dataset_month), ind.training)

  data_month_s <- dataset_month[ind.training, ]   #Training dataset currently used in the modeling

  data_month_v <- dataset_month[ind.testing, ]    #Testing/validation dataset using input sampling

  data_month <- data_month_s #training data for  monthhly predictions...

  #date_proc<-strptime(sampling_dat$date[i], "%Y%m%d")   # interpolation date being processed

  #mo<-as.integer(strftime(date_proc, "%m"))          # current month of the date being processed

  #day<-as.integer(strftime(date_proc, "%d"))

  #year<-as.integer(strftime(date_proc, "%Y"))

  ## end of pasted

  #end of insert...

    data_name<-paste(var,"_clim_month_",as.integer(month_no),"_",cname[k],"_",prop_month,

  data_name<-paste(var,"_clim_month_",as.integer(month_no),"_",model_name,"_",prop_month,

  clim_obj<-list(rast_clim_list,data_month,data_month_v,sampling_month_dat[j,],mod_list,list_formulas)

  names(clim_obj)<-c("clim","data_month","data_month_v","sampling_month_dat","mod","formulas")

  save(clim_obj,file= file.path(out_path,paste("clim_obj_CAI_month_",as.integer(month_no),"_",var,"_",prop_month,

                                               "_",run_samp,"_",out_prefix,".RData",sep="")))

  ghcn.subsets<-sampling_obj$ghcn_data

  sampling_month_obj <- list_param$sampling_month_obj

  daily_dev_sampling_dat <- list_param$daily_dev_sampling_dat

  index_d <- daily_dev_sampling_dat$index_d[i]

  index_m <- daily_dev_sampling_dat$index_m[i]

  use_clim_image <- list_param$use_clim_image # use predicted image as a base...rather than average Tmin at the station for delta

  join_daily <- list_param$join_daily # join monthly and daily station before calucating delta

  #use_clim_image

  #use index_d and index_m

  date<-strptime(daily_dev_sampling_dat$date[i], "%Y%m%d")   # interpolation date being processed

  data_day<-ghcn.subsets[[index_d]]

  mod_LST <- ghcn.subsets[[index_d]][,match(LST_month, names(ghcn.subsets[[index_d]]))]  #Match interpolation date and monthly LST average

  ind.training<-sampling[[index_d]]

  date_proc<-sampling_dat$date[index_d]

  date_proc<-strptime(sampling_dat$date[index_d], "%Y%m%d")   # interpolation date being processed

  #Now get monthly data...

  ghcn.month.subsets<-sampling_month_obj$ghcn_data

  sampling_month_dat <- sampling_month_obj$sampling_dat

  sampling_month_index <- sampling_month_obj$sampling_index

  dataset_month <-ghcn.month.subsets[[index_m]]

  mod_LST <- ghcn.month.subsets[[index_m]][,match(LST_month, names(ghcn.month.subsets[[index_m]]))]  #Match interpolation date and monthly LST average

  dataset_month$LST <- as.data.frame(mod_LST)[,1] #Add the variable LST to the dataset

  #change here...

  dst$LST<-dst[[LST_month]] #Add the variable LST to the monthly dataset

  proj_str<-proj4string(dst) #get the local projection information from monthly data

  ind.training_month <- sampling_month_index[[index_m]]

  ind.testing_month  <- setdiff(1:nrow(dataset_month), ind.training_month)

  data_month_s <- dataset_month[ind.training_month, ]   #Training dataset currently used in the modeling

  data_month_v <- dataset_month[ind.testing_month, ]    #Testing/validation dataset using input sampling

  modst <- data_month_s #training data for  monthhly predictions...

  # STEP 2 - CLEAN DATA AND JOIN DAILY TO MONTHLY STATION INFORMATION

  #if use join

  #modst<-dst[dst$month==mo,] #Subsetting dataset for the relevant month of the date being processed

  #if used images

  # extract from image 

  if(use_clim_image==FALSE){

    #must join daily and monthly data first...

    dmoday <-merge(modst,d,by="id",suffixes=c("",".y2"))  

    xmoday <-merge(modst,x,by="id",suffixes=c("",".y2"))  

    mod_pat<-glob2rx("*.y2")   #remove duplicate columns that have ".y2" in their names

    var_pat<-grep(mod_pat,names(dmoday),value=FALSE) # using grep with "value" extracts the matching names

    dmoday<-dmoday[,-var_pat] #dropping relevant columns

    mod_pat<-glob2rx("*.y2")   

    var_pat<-grep(mod_pat,names(xmoday),value=FALSE) # using grep with "value" extracts the matching names

    xmoday<-xmoday[,-var_pat] #Removing duplicate columns

    data_v<-xmoday

    #coords <-dmoday[,c("coords.x1","coords.x2")]

    coords <-dmoday[,c("x","y")]

    coordinates(dmoday)<-coords

    proj4string(dmoday)<-proj_str

    #dmoday contains the daily tmax values for training with TMax/TMin being the monthly station tmax/tmin mean

    #xmoday contains the daily tmax values for validation with TMax/TMin being the monthly station tmax/tmin mean

    if (var=="TMIN"){

      daily_delta <-dmoday$dailyTmin-dmoday$TMin #daily detl is the difference between monthly and daily temperatures

    }

    if (var=="TMAX"){

      daily_delta <- dmoday$dailyTmax-dmoday$TMax

    }  

    #daily_delta <- dmoday[[y_var_name]] - 

    #only one delta in this case!!!

    #list(mod)

    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 bias surface...

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

  if(use_clim_image==TRUE){

    # User can choose to join daily and monthly station before interpolation:

    #-this ensures that the delta difference is "more" exact since its starting point is basesd on average value but there is risk to loose some stations

    #may need to change this option later!!

    #if jion_Daily is true then daily station used as training will match monthly station used as training

    if(join_daily==TRUE){

      dmoday <-merge(modst,d,by="id",suffixes=c("",".y2"))  

      xmoday <-merge(modst,x,by="id",suffixes=c("",".y2"))  

      mod_pat<-glob2rx("*.y2")   #remove duplicate columns that have ".y2" in their names

      var_pat<-grep(mod_pat,names(dmoday),value=FALSE) # using grep with "value" extracts the matching names

      dmoday<-dmoday[,-var_pat] #dropping relevant columns

      mod_pat<-glob2rx("*.y2")   

      var_pat<-grep(mod_pat,names(xmoday),value=FALSE) # using grep with "value" extracts the matching names

      xmoday<-xmoday[,-var_pat] #Removing duplicate columns

      data_v<-xmoday

    }else{

      dmoday<-d

      data_v<-x

    }

    #dmoday contains the daily tmax values for training with TMax/TMin being the monthly station tmax/tmin mean

    #xmoday contains the daily tmax values for validation with TMax/TMin being the monthly station tmax/tmin mean

    #coords <-dmoday[,c("coords.x1","coords.x2")]

    coords <-dmoday[,c("x","y")]

    coordinates(dmoday)<-coords

    proj4string(dmoday)<-proj_str

    #Now compute daily delta deviation from climatology layer:

    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)

    extract_data_s <-extract(rast_clim_mod,dmoday,df=TRUE)

    #list_daily_delta

    daily_delta_df <- dmoday[[y_var_name]] - extract_data_s  

    daily_delta_df <- daily_delta_df[,-1]

    names(daily_delta_df) <- paste(names(daily_delta_df),"_del",sep="")

    names(extract_data_s) <- paste(names(extract_data_s),"_m",sep="") # "m" for monthly predictions...

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

    #Now krige  forevery model !! loop

    list_mod_krtmp2 <- vector("list",length=nlayers(rast_clim_mod)) 

    list_daily_delta_rast <- vector("list",length=nlayers(rast_clim_mod)) 

    for(k in 1:nlayers(rast_clim_mod)){

      daily_delta <- daily_delta_df[[k]]

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

      model_name<- names(daily_delta_df)[k]

      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...

      #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[[k]] <- raster_name_delta   

    }

    raster_name_delta <- list_daily_delta_rast

    mod_krtmp2 <- list_mod_krtmp2

  if(use_clim_image==FALSE){

    list_daily_delta_rast <- rep(raster_name_delta,length=nlayers(rast_clim_mod))

  }

  temp_list<-vector("list",nlayers(rast_clim_mod))  

  for (k in 1:nlayers(rast_clim_mod)){

    rast_clim_month<-raster(rast_clim_list[[k]])

    daily_delta_rast <- raster(list_daily_delta_rast[[k]])

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

    temp_list[[k]]<-raster_name

  data_s<-dmoday #put the 

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

  #proj4string(data_s)<-proj_str

  #mod_krtmp2<-list_mod_krtmp2    

  #mod_krtmp2<-fitdelta

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

  names(temp_list)<-names(rast_clim_list)

  #add data_month_s and data_month_v?

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

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

               "data_month_s","data_month_v","sampling_dat","daily_dev_sampling_dat",model_name)

  save(delta_obj,file= file.path(out_path,paste("delta_obj_",var,"_",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],

                                                out_prefix,".RData",sep="")))

combine_sampling_daily_monthly_for_daily_deviation_pred <- function(sampling_obj,sampling_month_obj){

  #This function combines sampling objects information at the daily and monthly time scale.

  #The data.frame created records propotion of hold out, sampling number as well as provides keys (index_d,index_m)

  #to access training and testing samples at the daily and monthly time scale (long term monthly!!!)

  #Inputs:

  #sampling_obj: contains daily sampling information and is output of "sampling_training_testing" function defined in sampling_script_functions.

  #sampling_month_obj: contains monthly sampling information and is output of "sampling_training_testing" function defined in sampling_script_functions

  sampling_month_dat <- sampling_month_obj$sampling_dat #extract monthly information about sampling (in data.frame format)

  sampling_dat <- sampling_obj$sampling_dat #extract daily information about samplint (in data.frame format)

  #Build new data.frame with combined information

  data_daily_dev <- sampling_dat 

  data_daily_dev$index_d <- 1:nrow(data_daily_dev) # index_d is a identifier key to the list of daily samples generated earlier in sampling_obj

  data_daily_dev$month_no <- as.integer(strftime(strptime(data_daily_dev$date, "%Y%m%d"),"%m")) #Extract month given specific format  "%Y%m%d"

  names(sampling_month_dat) <- c("dates_month","run_samp_month","prop_month") #rename before merging as table product

  sampling_month_dat$month_no <-as.integer(strftime(strptime(sampling_month_dat$dates_month, "%Y%m%d"),"%m")) #month of the record

  sampling_month_dat$index_m <-1:nrow(sampling_month_dat) # index_m is a identifier key to the list of monthly samples 

                                                          #generated earlier in sampling_month_obj

  data_daily_dev <-merge(data_daily_dev,sampling_month_dat,by = "month_no")

  return(data_daily_dev)

}