Project

General

Profile

« Previous | Next » 

Revision b2563a44

Added by Benoit Parmentier over 11 years ago

raster prediction function, modifications related to gam daily method

View differences:

climate/research/oregon/interpolation/GAM_fusion_analysis_raster_prediction_multisampling.R
99 99 
  #6 additional parameters for monthly climatology and more
100 100 
  list_models<-list_param_raster_prediction$list_models
101 101 
  lst_avg<-list_param_raster_prediction$lst_avg
102 
  in_path<-list_param_raster_prediction$in_path
103 102 
  out_path<-list_param_raster_prediction$out_path
104 103 
  script_path<-list_param_raster_prediction$script_path
105 104 
  interpolation_method<-list_param_raster_prediction$interpolation_method
106 105
107 
  #setwd(in_path)
108 106 
  setwd(out_path)
109 107 
  #Sourcing in the master script to avoid confusion on the latest versions of scripts and functions!!!
110 108
......
113 111 
  #source(file.path(script_path,"GAM_fusion_function_multisampling_validation_metrics_03182013.R"))
114 112
115 113 
  ###################### START OF THE SCRIPT ########################
116
117
114
118 115 
  if (var=="TMAX"){
119 116 
    y_var_name<-"dailyTmax"
120 117 
  }
......
127 124 
  #create log file to keep track of details such as processing times and parameters.
128 125
129 126 
  #log_fname<-paste("R_log_raster_prediction",out_prefix, ".log",sep="")
130 
  log_fname<-paste(out_path,"R_log_raster_prediction",out_prefix, ".log",sep="")
127 
  log_fname<-paste("R_log_raster_prediction",out_prefix, ".log",sep="")
131 128 
  #sink(log_fname) #create new log file
132 
  file.create(log_fname) #create new log file
133
134 
  #if (file.exists(log_fname)){  #Stop the script???
135 
  #  file.remove(log_fname)
136 
  #  log_file<-file(log_fname,"w")
137 
  #}
138 
  #if (!file.exists(log_fname)){
139 
  #  log_file<-file(file.path(out_path,log_fname),"w")
140 
  #}
129 
  file.create(file.path(path,log_fname)) #create new log file
141 130
142 131 
  time1<-proc.time()    #Start stop watch
143 132
144 
  #writeLines(paste("Starting script at this local Date and Time: ",as.character(Sys.time()),sep=""),
145 
  #           con=log_file,sep="\n")
146 
  #writeLines("Starting script process time:",con=log_file,sep="\n")
147 
  #writeLines(as.character(time1),con=log_file,sep="\n")
148 133 
  cat(paste("Starting script at this local Date and Time: ",as.character(Sys.time()),sep=""),
149 134 
             file=log_fname,sep="\n")
150 135 
  cat("Starting script process time:",file=log_fname,sep="\n",append=TRUE)
......
207 192 
    clim_method_mod_obj<-mclapply(1:12, list_param=list_param_runClim_KGFusion, runClim_KGFusion,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement
208 193 
    #test<-runClim_KGFusion(1,list_param=list_param_runClim_KGFusion)
209 194 
    #gamclim_fus_mod<-mclapply(1:6, list_param=list_param_runClim_KGFusion, runClim_KGFusion,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement
210
195 
    save(clim_method_mod_obj,file= file.path(out_path,paste(interpolation_method,"_mod_",y_var_name,out_prefix,".RData",sep="")))
196 
    list_tmp<-vector("list",length(clim_method_mod_obj))
197 
    for (i in 1:length(clim_method_mod_obj)){
198 
      tmp<-clim_method_mod_obj[[i]]$clim
199 
      list_tmp[[i]]<-tmp
200 
    }
201 
    clim_yearlist<-list_tmp
211 202 
  }
212 203
213 204 
  if (interpolation_method=="gam_CAI"){
......
215 206 
    names(list_param_runClim_KGCAI)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","var","y_var_name","out_prefix","out_path")
216 207 
    clim_method_mod_obj<-mclapply(1:12, list_param=list_param_runClim_KGCAI, runClim_KGCAI,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement
217 208 
    #test<-runClim_KGCAI(1,list_param=list_param_runClim_KGCAI)
218 
    #gamclim_fus_mod<-mclapply(1:6, list_param=list_param_runClim_KGFusion, runClim_KGFusion,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement
209 
    #gamclim_fus_mod<-mclapply(1:6, list_param=list_param_runClim_KGFusion, runClim_KGFusion,mc.preschedule=FALSE,mc.cores = 6)
210 
    save(clim_method_mod_obj,file= file.path(out_path,paste(interpolation_method,"_mod_",y_var_name,out_prefix,".RData",sep="")))
211 
    list_tmp<-vector("list",length(clim_method_mod_obj))
212 
    for (i in 1:length(clim_method_mod_obj)){
213 
      tmp<-clim_method_mod_obj[[i]]$clim
214 
      list_tmp[[i]]<-tmp
215 
    }
216 
    clim_yearlist<-list_tmp
219 217 
  }
220 218
221 
  #gamclim_fus_mod<-runClim_KGFusion(1,list_param=list_param_runClim_KGFusion) #This is the end bracket from mclapply(...) statement
222 
  save(clim_method_mod_obj,file= file.path(out_path,paste(interpolation_method,"_mod_",y_var_name,out_prefix,".RData",sep="")))
223 
  t2<-proc.time()-t1
224 219 
  cat(as.character(t2),file=log_fname,sep="\n", append=TRUE)
225
226 
  #now get list of raster clim layers
227
228 
  list_tmp<-vector("list",length(clim_method_mod_obj))
229 
  for (i in 1:length(clim_method_mod_obj)){
230 
    tmp<-clim_method_mod_obj[[i]]$clim
231 
    list_tmp[[i]]<-tmp
232 
  }
220 
  t2<-proc.time()-t1
233 221
234 222 
  ################## PREDICT AT DAILY TIME SCALE #################
235 223 
  #Predict at daily time scale from single time scale or multiple time scale methods: 2 methods availabe now
236 224
237 225 
  #put together list of clim models per month...
238 226 
  #rast_clim_yearlist<-list_tmp
239 
  clim_yearlist<-list_tmp
227
240 228 
  #Second predict at the daily time scale: delta
241 229
242 230 
  #method_mod_obj<-mclapply(1:1, runGAMFusion,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement
......
253 241 
    #test<-runGAMFusion(1,list_param=list_param_runGAMFusion)
254 242
255 243 
    method_mod_obj<-mclapply(1:length(sampling_obj$ghcn_data_day),list_param=list_param_run_prediction_daily_deviation,run_prediction_daily_deviation,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement
256 
    #method_mod_obj<-mclapply(1:1,list_param=list_param_run_prediction_daily_deviation,run_prediction_daily_deviation,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement
244 
    save(method_mod_obj,file= file.path(out_path,paste("method_mod_obj_",interpolation_method,"_",y_var_name,out_prefix,".RData",sep="")))
245 
   }
246
247 
  if (interpolation_method=="gam_daily"){
248 
    #input a list:note that ghcn.subsets is not sampling_obj$data_day_ghcn
257 249
258 
    #method_mod_obj<-mclapply(1:length(sampling_obj$ghcn_data_day),runGAMFusion,list_param_runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement
259 
    #method_mod_obj<-mclapply(1:length(ghcn.subsets), runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement
250 
    list_param_run_prediction_gam_daily <-list(i,s_raster,covar_names,lst_avg,list_models,dst,var,y_var_name, sampling_obj,interpolation_method,out_prefix,out_path)
251 
    names(list_param_run_prediction_gam_daily)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","var","y_var_name","sampling_obj","interpolation_method","out_prefix","out_path")
252 
    #test <- runGAM_day_fun(1,list_param_run_prediction_gam_daily)
253 
    method_mod_obj<-mclapply(1:length(sampling_obj$ghcn_data_day),list_param=list_param_run_prediction_gam_daily,runGAM_day_fun,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement
254 
    #method_mod_obj<-mclapply(1:18,list_param=list_param_run_prediction_gam_daily,runGAM_day_fun,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement
260 255
261 
  }
256 
    save(method_mod_obj,file= file.path(out_path,paste("method_mod_obj_",interpolation_method,"_",y_var_name,out_prefix,".RData",sep="")))
262 257
263 
  save(method_mod_obj,file= file.path(out_path,paste("method_mod_obj_",interpolation_method,"_",y_var_name,out_prefix,".RData",sep="")))
258 
  }
259
264 260 
  t2<-proc.time()-t1
265 261 
  cat(as.character(t2),file=log_fname,sep="\n", append=TRUE)
266 262 
  #browser()
......
273 269 
    tmp<-method_mod_obj[[i]][[y_var_name]]  #y_var_name is the variable predicted (dailyTmax or dailyTmin)
274 270 
    list_tmp[[i]]<-tmp
275 271 
  }
276 
  rast_day_yearlist<-list_tmp #list of predicted images
272 
  rast_day_yearlist<-list_tmp #list of predicted images over full year...
277 273
278 274 
  cat("Validation step:",file=log_fname,sep="\n", append=TRUE)
279 275 
  t1<-proc.time()
......
316 312 
  ################### PREPARE RETURN OBJECT ###############
317 313 
  #Will add more information to be returned
318 314
319 
  raster_prediction_obj<-list(clim_method_mod_obj,method_mod_obj,validation_mod_obj,tb_diagnostic_v,
320 
                              summary_metrics_v,summary_month_metrics_v)
321 
  names(raster_prediction_obj)<-c("clim_method_mod_obj","method_mod_obj","validation_mod_obj","tb_diagnostic_v",
322 
                                  "summary_metrics_v","summary_month_metrics_v")
323 
  save(raster_prediction_obj,file= file.path(out_path,paste("raster_prediction_obj_",interpolation_method,"_", y_var_name,out_prefix,".RData",sep="")))
315 
  if (interpolation_method=="gam_CAI" | interpolation_method=="gam_fusion"){
316 
    raster_prediction_obj<-list(clim_method_mod_obj,method_mod_obj,validation_mod_obj,tb_diagnostic_v,
317 
                                summary_metrics_v,summary_month_metrics_v)
318 
    names(raster_prediction_obj)<-c("clim_method_mod_obj","method_mod_obj","validation_mod_obj","tb_diagnostic_v",
319 
                                    "summary_metrics_v","summary_month_metrics_v")
320 
    save(raster_prediction_obj,file= file.path(out_path,paste("raster_prediction_obj_",interpolation_method,"_", y_var_name,out_prefix,".RData",sep="")))
321
322 
  }
323
324 
  if (interpolation_method=="gam_daily"){
325 
    raster_prediction_obj<-list(method_mod_obj,validation_mod_obj,tb_diagnostic_v,
326 
                                summary_metrics_v,summary_month_metrics_v)
327 
    names(raster_prediction_obj)<-c("method_mod_obj","validation_mod_obj","tb_diagnostic_v",
328 
                                    "summary_metrics_v","summary_month_metrics_v")
329 
    save(raster_prediction_obj,file= file.path(out_path,paste("raster_prediction_obj_",interpolation_method,"_", y_var_name,out_prefix,".RData",sep="")))
330
331 
  }
324 332
325 333 
  return(raster_prediction_obj)
326 334 
}

Also available in: Unified diff