Project

General

Profile

« Previous | Next » 

Revision 5ed5c7a3

Added by Benoit Parmentier over 11 years ago

raster prediction, modifications log file and output tracking for workflow run

View differences:

climate/research/oregon/interpolation/GAM_fusion_analysis_raster_prediction_multisampling.R
11 11 
#5)possibilty of running GAM+FUSION or GAM+CAI and other options added
12 12 
#The interpolation is done first at the monthly time scale then delta surfaces are added.
13 13 
#AUTHOR: Benoit Parmentier
14 
#DATE: 04/22/2013
14 
#DATE: 05/06/2013
15 15 
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#568--
16 16 
#
17 17 
# TO DO:
......
112 112 
  #source(file.path(script_path,"GAM_fusion_function_multisampling_03122013.R"))
113 113 
  #source(file.path(script_path,"GAM_fusion_function_multisampling_validation_metrics_03182013.R"))
114 114
115
116 115 
  ###################### START OF THE SCRIPT ########################
117 116
118 117
......
127 126
128 127 
  #create log file to keep track of details such as processing times and parameters.
129 128
130 
  log_fname<-paste("R_log_raster_prediction",out_prefix, ".log",sep="")
129 
  #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="")
131 
  #sink(log_fname) #create new log file
132 
  file.create(log_fname) #create new log file
131 133
132 
  if (file.exists(log_fname)){  #Stop the script???
133 
    file.remove(log_fname)
134 
    log_file<-file(log_fname,"w")
135 
  }
136 
  if (!file.exists(log_fname)){
137 
    log_file<-file(log_fname,"w")
138 
  }
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 
  #}
139 141
140 142 
  time1<-proc.time()    #Start stop watch
141 
  writeLines(paste("Starting script at this local Date and Time: ",as.character(Sys.time()),sep=""),
142 
             con=log_file,sep="\n")
143 
  writeLines("Starting script process time:",con=log_file,sep="\n")
144 
  writeLines(as.character(time1),con=log_file,sep="\n")
143
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 
  cat(paste("Starting script at this local Date and Time: ",as.character(Sys.time()),sep=""),
149 
             file=log_fname,sep="\n")
150 
  cat("Starting script process time:",file=log_fname,sep="\n",append=TRUE)
151 
  cat(as.character(time1),file=log_fname,sep="\n",append=TRUE)
145 152
146 153 
  ############### READING INPUTS: DAILY STATION DATA AND OTEHR DATASETS  #################
147 154
148 
  ghcn<-readOGR(dsn=in_path,layer=sub(".shp","",basename(infile_daily)))
155 
  ghcn<-readOGR(dsn=out_path,layer=sub(".shp","",basename(infile_daily)))
149 156 
  CRS_interp<-proj4string(ghcn)                       #Storing projection information (ellipsoid, datum,etc.)
150 
  stat_loc<-readOGR(dsn=in_path,layer=sub(".shp","",basename(infile_locs)))
157 
  stat_loc<-readOGR(dsn=out_path,layer=sub(".shp","",basename(infile_locs)))
151 158 
  #dates2 <-readLines(file.path(in_path,dates_selected)) #dates to be predicted, now read directly from the file
152 159 
  if (dates_selected==""){
153 160 
    dates<-as.character(sort(unique(ghcn$date))) #dates to be predicted
......
156 163 
    dates<-dates_selected #dates to be predicted
157 164 
  }
158 165
159 
  #Reading of covariate brick covariates can be changed...
166 
  #Reaccccding of covariate brick covariates can be changed...
160 167
161 168 
  s_raster<-brick(infile_covariates)                   #read in the data brck
162 169 
  names(s_raster)<-covar_names               #Assigning names to the raster layers: making sure it is included in the extraction
163 170
164 171 
  #Reading monthly data
165 
  dst<-readOGR(dsn=in_path,layer=sub(".shp","",basename(infile_monthly)))
172 
  dst<-readOGR(dsn=out_path,layer=sub(".shp","",basename(infile_monthly)))
166 173
167 174 
  ### TO DO -important ###
168 175 
  #SCREENING IN COVARIATE SCRIPT AND DATA PREP SCRIPT !!! Only perform predictions here
......
187 194 
  ########### PREDICT FOR MONTHLY SCALE  ##################
188 195
189 196 
  #First predict at the monthly time scale: climatology
190 
  writeLines("Predictions at monthly scale:",con=log_file,sep="\n")
197 
  cat("Predictions at monthly scale:",file=log_fname,sep="\n", append=TRUE)
198 
  cat(paste("Local Date and Time: ",as.character(Sys.time()),sep=""),
199 
      file=log_fname,sep="\n")
191 200 
  t1<-proc.time()
192 201 
  j=12
193 202 
  #browser()
......
202 211 
  }
203 212
204 213 
  if (interpolation_method=="gam_CAI"){
205 
    list_param_runClim_KGCAI<-list(j,s_raster,covar_names,lst_avg,list_models,dst,var,y_var_name, out_prefix)
206 
    names(list_param_runClim_KGCAI)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","var","y_var_name","out_prefix")
214 
    list_param_runClim_KGCAI<-list(j,s_raster,covar_names,lst_avg,list_models,dst,var,y_var_name, out_prefix,out_path)
215 
    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")
207 216 
    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
208 217 
    #test<-runClim_KGCAI(1,list_param=list_param_runClim_KGCAI)
209 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
210 219 
  }
211 220
212 221 
  #gamclim_fus_mod<-runClim_KGFusion(1,list_param=list_param_runClim_KGFusion) #This is the end bracket from mclapply(...) statement
213 
  save(clim_method_mod_obj,file= paste(interpolation_method,"_mod_",y_var_name,out_prefix,".RData",sep=""))
222 
  save(clim_method_mod_obj,file= file.path(out_path,paste(interpolation_method,"_mod_",y_var_name,out_prefix,".RData",sep="")))
214 223 
  t2<-proc.time()-t1
215 
  writeLines(as.character(t2),con=log_file,sep="\n")
224 
  cat(as.character(t2),file=log_fname,sep="\n", append=TRUE)
216 225
217 226 
  #now get list of raster clim layers
218 227
......
231 240 
  #Second predict at the daily time scale: delta
232 241
233 242 
  #method_mod_obj<-mclapply(1:1, runGAMFusion,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement
234 
  writeLines("Predictions at the daily scale:",con=log_file,sep="\n")
243 
  cat("Predictions at the daily scale:",file=log_fname,sep="\n", append=TRUE)
235 244 
  t1<-proc.time()
245 
  cat(paste("Local Date and Time: ",as.character(Sys.time()),sep=""),
246 
      file=log_fname,sep="\n")
236 247
237 248 
  if (interpolation_method=="gam_CAI" | interpolation_method=="gam_fusion"){
238 249 
    #input a list:note that ghcn.subsets is not sampling_obj$data_day_ghcn
239 
    list_param_run_prediction_daily_deviation <-list(i,clim_yearlist,sampling_obj,dst,var,y_var_name, interpolation_method,out_prefix)
240 
    names(list_param_run_prediction_daily_deviation)<-c("list_index","clim_yearlist","sampling_obj","dst","var","y_var_name","interpolation_method","out_prefix")
250 
    list_param_run_prediction_daily_deviation <-list(i,clim_yearlist,sampling_obj,dst,var,y_var_name, interpolation_method,out_prefix,out_path)
251 
    names(list_param_run_prediction_daily_deviation)<-c("list_index","clim_yearlist","sampling_obj","dst","var","y_var_name","interpolation_method","out_prefix","out_path")
241 252 
    #test<-mclapply(1:18, runGAMFusion,list_param=list_param_runGAMFusion,mc.preschedule=FALSE,mc.cores = 9)
242 253 
    #test<-runGAMFusion(1,list_param=list_param_runGAMFusion)
243 254
......
249 260
250 261 
  }
251 262
252 
  save(method_mod_obj,file= paste("method_mod_obj_",interpolation_method,"_",y_var_name,out_prefix,".RData",sep=""))
263 
  save(method_mod_obj,file= file.path(out_path,paste("method_mod_obj_",interpolation_method,"_",y_var_name,out_prefix,".RData",sep="")))
253 264 
  t2<-proc.time()-t1
254 
  writeLines(as.character(t2),con=log_file,sep="\n")
265 
  cat(as.character(t2),file=log_fname,sep="\n", append=TRUE)
255 266 
  #browser()
256 267
257 268 
  ############### NOW RUN VALIDATION #########################
......
264 275 
  }
265 276 
  rast_day_yearlist<-list_tmp #list of predicted images
266 277
267 
  writeLines("Validation step:",con=log_file,sep="\n")
278 
  cat("Validation step:",file=log_fname,sep="\n", append=TRUE)
268 279 
  t1<-proc.time()
269 
  #calculate_accuary_metrics<-function(i)
270 
  list_param_validation<-list(i,rast_day_yearlist,method_mod_obj,y_var_name, out_prefix)
271 
  names(list_param_validation)<-c("list_index","rast_day_year_list","method_mod_obj","y_var_name","out_prefix") #same names for any method
280 
  cat(paste("Local Date and Time: ",as.character(Sys.time()),sep=""),
281 
      file=log_fname,sep="\n")
272 282
273 
  #gam_fus_validation_mod<-mclapply(1:length(method_mod_obj), calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement
274 
  validation_mod_obj <-mclapply(1:length(method_mod_obj), list_param=list_param_validation, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement
283 
  list_param_validation<-list(i,rast_day_yearlist,method_mod_obj,y_var_name, out_prefix, out_path)
284 
  names(list_param_validation)<-c("list_index","rast_day_year_list","method_mod_obj","y_var_name","out_prefix", "out_path") #same names for any method
275 285
276 
  #gam_fus_validation_mod<-mclapply(1:1, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement
277 
  save(validation_mod_obj,file= paste(interpolation_method,"_validation_mod_obj_",y_var_name,out_prefix,".RData",sep=""))
286 
  validation_mod_obj <-mclapply(1:length(method_mod_obj), list_param=list_param_validation, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 9)
287 
  #test_val<-calculate_accuracy_metrics(1,list_param_validation)
288 
  save(validation_mod_obj,file= file.path(out_path,paste(interpolation_method,"_validation_mod_obj_",y_var_name,out_prefix,".RData",sep="")))
278 289 
  t2<-proc.time()-t1
279 
  writeLines(as.character(t2),con=log_file,sep="\n")
290 
  cat(as.character(t2),file=log_fname,sep="\n", append=TRUE)
280 291
281 292 
  #################### ASSESSMENT OF PREDICTIONS: PLOTS OF ACCURACY METRICS ###########
282 293
......
284 295 
  tb_diagnostic_v<-extract_from_list_obj(validation_mod_obj,"metrics_v")
285 296 
  rownames(tb_diagnostic_v)<-NULL #remove row names
286 297
287 
  #Call function to create plots of metrics for validation dataset
298 
  #Call functions to create plots of metrics for validation dataset
288 299 
  metric_names<-c("rmse","mae","me","r","m50")
289 
  summary_metrics_v<-boxplot_from_tb(tb_diagnostic_v,metric_names,out_prefix)
300 
  summary_metrics_v<- boxplot_from_tb(tb_diagnostic_v,metric_names,out_prefix,out_path)
290 301 
  names(summary_metrics_v)<-c("avg","median")
302 
  summary_month_metrics_v<- boxplot_month_from_tb(tb_diagnostic_v,metric_names,out_prefix,out_path)
291 303
292 304 
  #################### CLOSE LOG FILE  ####################
293 305
294 306 
  #close log_file connection and add meta data
295 
  writeLines("Finished script process time:",con=log_file,sep="\n")
307 
  cat("Finished script process time:",file=log_fname,sep="\n", append=TRUE)
296 308 
  time2<-proc.time()-time1
297 
  writeLines(as.character(time2),con=log_file,sep="\n")
298 
  #later on add all the paramters used in the script...
299 
  writeLines(paste("Finished script at this local Date and Time: ",as.character(Sys.time()),sep=""),
300 
             con=log_file,sep="\n")
301 
  writeLines("End of script",con=log_file,sep="\n")
302 
  close(log_file)
309 
  cat(as.character(time2),file=log_fname,sep="\n", append=TRUE)
310 
  #later on add all the parameters used in the script...
311 
  cat(paste("Finished script at this local Date and Time: ",as.character(Sys.time()),sep=""),
312 
             file=log_fname,sep="\n", append=TRUE)
313 
  cat("End of script",file=log_fname,sep="\n", append=TRUE)
314 
  #close(log_fname)
303 315
304 316 
  ################### PREPARE RETURN OBJECT ###############
305 317 
  #Will add more information to be returned
306 318
307 
  raster_prediction_obj<-list(clim_method_mod_obj,method_mod_obj,validation_mod_obj,tb_diagnostic_v,summary_metrics_v)
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)
308 321 
  names(raster_prediction_obj)<-c("clim_method_mod_obj","method_mod_obj","validation_mod_obj","tb_diagnostic_v",
309 
                                  "summary_metrics_v")
310 
  save(raster_prediction_obj,file= paste("raster_prediction_obj_",interpolation_method,"_", y_var_name,out_prefix,".RData",sep=""))
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="")))
311 324
312 325 
  return(raster_prediction_obj)
313 326 
}

Also available in: Unified diff