Revision ba485f60
Added by Benoit Parmentier about 11 years ago
| climate/research/oregon/interpolation/GAM_fusion_analysis_raster_prediction_multisampling.R | ||
|---|---|---|
| 125 | 125 |
|
| 126 | 126 |
###################### START OF THE SCRIPT ######################## |
| 127 | 127 |
|
| 128 |
#This should not be set here...? master script, modify for precip |
|
| 128 | 129 |
if (var=="TMAX"){
|
| 129 |
y_var_name<-"dailyTmax" |
|
| 130 |
y_var_name<-"dailyTmax" |
|
| 131 |
y_var_month<-"TMax" |
|
| 130 | 132 |
} |
| 131 |
|
|
| 132 | 133 |
if (var=="TMIN"){
|
| 133 |
y_var_name<-"dailyTmin" |
|
| 134 |
y_var_name<-"dailyTmin" |
|
| 135 |
y_var_month <-"TMin" |
|
| 134 | 136 |
} |
| 135 | 137 |
|
| 136 | 138 |
################# CREATE LOG FILE ##################### |
| ... | ... | |
| 329 | 331 |
|
| 330 | 332 |
############### NOW RUN VALIDATION ######################### |
| 331 | 333 |
#SIMPLIFY THIS PART: one call |
| 332 |
|
|
| 333 |
#list_tmp<-vector("list",length(method_mod_obj))
|
|
| 334 |
#for (i in 1:length(method_mod_obj)){
|
|
| 335 |
# tmp<-method_mod_obj[[i]][[y_var_name]] #y_var_name is the variable predicted (dailyTmax or dailyTmin) |
|
| 336 |
# list_tmp[[i]]<-tmp |
|
| 337 |
#} |
|
| 338 |
#rast_day_yearlist<-list_tmp #list of predicted images over full year... |
|
| 339 |
|
|
| 340 |
list_data_v <- extract_list_from_list_obj(method_mod_obj,"data_v") |
|
| 341 |
list_data_s <- extract_list_from_list_obj(method_mod_obj,"data_s") |
|
| 342 |
rast_day_yearlist <- extract_list_from_list_obj(method_mod_obj,y_var_name) #list_tmp #list of predicted images over full year... |
|
| 343 |
list_sampling_dat <- extract_list_from_list_obj(method_mod_obj,"sampling_dat") |
|
| 344 |
|
|
| 334 |
|
|
| 345 | 335 |
cat("Validation step:",file=log_fname,sep="\n", append=TRUE)
|
| 346 | 336 |
t1<-proc.time() |
| 347 | 337 |
cat(paste("Local Date and Time: ",as.character(Sys.time()),sep=""),
|
| 348 | 338 |
file=log_fname,sep="\n") |
| 349 | 339 |
|
| 350 |
multi_time_scale <- FALSE |
|
| 351 |
#list_param_validation<-list(i,rast_day_yearlist,method_mod_obj,y_var_name, multi_time_scale,out_prefix, out_path) |
|
| 352 |
#names(list_param_validation)<-c("list_index","rast_day_year_list","method_mod_obj","y_var_name","multi_time_scale","out_prefix", "out_path") #same names for any method
|
|
| 353 |
#debug(calculate_accuracy_metrics) |
|
| 354 |
#test_val<-calculate_accuracy_metrics(1,list_param_validation) |
|
| 355 |
#list_param_validation<-list(i,rast_day_yearlist,method_mod_obj,y_var_name, multi_time_scale,out_prefix, out_path) |
|
| 356 |
#names(list_param_validation)<-c("list_index","rast_day_year_list","method_mod_obj","y_var_name","multi_time_scale","out_prefix", "out_path") #same names for any method
|
|
| 357 |
|
|
| 358 |
list_param_validation<-list(i,rast_day_yearlist,list_data_v,list_data_s,list_sampling_dat,y_var_name, multi_time_scale,out_prefix, out_path) |
|
| 359 |
names(list_param_validation)<-c("list_index","rast_day_year_list",
|
|
| 360 |
"list_data_v","list_data_s","list_sampling_dat","y_var_name","multi_time_scale","out_prefix", "out_path") #same names for any method |
|
| 361 |
debug(calculate_accuracy_metrics) |
|
| 362 |
test_val2 <-calculate_accuracy_metrics(1,list_param_validation) |
|
| 363 |
|
|
| 364 |
validation_mod_obj <-mclapply(1:length(method_mod_obj), list_param=list_param_validation, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 9) |
|
| 365 |
save(validation_mod_obj,file= file.path(out_path,paste(interpolation_method,"_validation_mod_obj_",y_var_name,out_prefix,".RData",sep=""))) |
|
| 366 |
t2<-proc.time()-t1 |
|
| 367 |
cat(as.character(t2),file=log_fname,sep="\n", append=TRUE) |
|
| 368 |
|
|
| 340 |
if (interpolation_method=="gam_daily" | interpolation_method=="kriging_daily" | interpolation_method=="gwr_daily"){
|
|
| 341 |
multi_time_scale <- FALSE |
|
| 342 |
|
|
| 343 |
list_data_v <- extract_list_from_list_obj(method_mod_obj,"data_v") |
|
| 344 |
list_data_s <- extract_list_from_list_obj(method_mod_obj,"data_s") |
|
| 345 |
rast_day_yearlist <- extract_list_from_list_obj(method_mod_obj,y_var_name) #list_tmp #list of predicted images over full year... |
|
| 346 |
list_sampling_dat <- extract_list_from_list_obj(method_mod_obj,"sampling_dat") |
|
| 347 |
|
|
| 348 |
list_param_validation<-list(i,rast_day_yearlist,list_data_v,list_data_s,list_sampling_dat,y_var_name, multi_time_scale,out_prefix, out_path) |
|
| 349 |
names(list_param_validation)<-c("list_index","rast_day_year_list",
|
|
| 350 |
"list_data_v","list_data_s","list_sampling_dat","y_ref","multi_time_scale","out_prefix", "out_path") #same names for any method |
|
| 351 |
#debug(calculate_accuracy_metrics) |
|
| 352 |
#test_val2 <-calculate_accuracy_metrics(1,list_param_validation) |
|
| 353 |
|
|
| 354 |
validation_mod_obj <-mclapply(1:length(method_mod_obj), list_param=list_param_validation, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 9) |
|
| 355 |
save(validation_mod_obj,file= file.path(out_path,paste(interpolation_method,"_validation_mod_obj_",y_var_name,out_prefix,".RData",sep=""))) |
|
| 356 |
t2<-proc.time()-t1 |
|
| 357 |
cat(as.character(t2),file=log_fname,sep="\n", append=TRUE) |
|
| 358 |
} |
|
| 359 |
|
|
| 369 | 360 |
### Run monthly validation if multi-time scale methods and add information to daily... |
| 370 | 361 |
|
| 371 | 362 |
if (interpolation_method %in% c("gam_CAI","kriging_CAI","gwr_CAI","gam_fusion","kriging_fusion","gwr_fusion")){
|
| 372 | 363 |
multi_time_scale <- TRUE |
| 373 |
list_param_validation_month <-list(i,clim_yearlist,clim_method_mod_obj,y_var_name, multi_time_scale ,out_prefix, out_path) |
|
| 374 |
names(list_param_validation_month)<-c("list_index","rast_day_year_list","method_mod_obj","y_var_name","multi_time_scale","out_prefix", "out_path") #same names for any method
|
|
| 364 |
i<-1 |
|
| 365 |
|
|
| 366 |
## daily time scale |
|
| 367 |
list_data_v <- extract_list_from_list_obj(method_mod_obj,"data_v") |
|
| 368 |
list_data_s <- extract_list_from_list_obj(method_mod_obj,"data_s") |
|
| 369 |
rast_day_yearlist <- extract_list_from_list_obj(method_mod_obj,y_var_name) #list_tmp #list of predicted images over full year... |
|
| 370 |
list_sampling_dat <- extract_list_from_list_obj(method_mod_obj,"daily_dev_sampling_dat") |
|
| 371 |
|
|
| 372 |
list_param_validation<-list(i,rast_day_yearlist,list_data_v,list_data_s,list_sampling_dat,y_var_name, multi_time_scale,out_prefix, out_path) |
|
| 373 |
names(list_param_validation)<-c("list_index","rast_day_year_list",
|
|
| 374 |
"list_data_v","list_data_s","list_sampling_dat","y_ref","multi_time_scale","out_prefix", "out_path") #same names for any method |
|
| 375 |
#debug(calculate_accuracy_metrics) |
|
| 376 |
#test_val2 <-calculate_accuracy_metrics(1,list_param_validation) |
|
| 377 |
|
|
| 378 |
validation_mod_obj <-mclapply(1:length(method_mod_obj), list_param=list_param_validation, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 9) |
|
| 379 |
save(validation_mod_obj,file= file.path(out_path,paste(interpolation_method,"_validation_mod_obj_",y_var_name,out_prefix,".RData",sep=""))) |
|
| 380 |
|
|
| 381 |
### monthly time scale |
|
| 382 |
list_data_v <- extract_list_from_list_obj(clim_method_mod_obj,"data_month_v") #extract monthly testing/validation dataset |
|
| 383 |
list_data_s <- extract_list_from_list_obj(clim_method_mod_obj,"data_month") #extract monthly training/fitting dataset |
|
| 384 |
rast_day_yearlist <- extract_list_from_list_obj(clim_method_mod_obj,"clim") #list_tmp #list of predicted images over full year at monthly time scale |
|
| 385 |
list_sampling_dat <- extract_list_from_list_obj(clim_method_mod_obj,"sampling_month_dat") |
|
| 386 |
|
|
| 387 |
#list_param_validation_month <-list(i,clim_yearlist,clim_method_mod_obj,y_var_name, multi_time_scale ,out_prefix, out_path) |
|
| 388 |
#names(list_param_validation_month)<-c("list_index","rast_day_year_list","method_mod_obj","y_var_name","multi_time_scale","out_prefix", "out_path") #same names for any method
|
|
| 389 |
|
|
| 390 |
list_param_validation_month <-list(i,rast_day_yearlist,list_data_v,list_data_s,list_sampling_dat,y_var_month, multi_time_scale,out_prefix, out_path) |
|
| 391 |
names(list_param_validation_month)<-c("list_index","rast_day_year_list",
|
|
| 392 |
"list_data_v","list_data_s","list_sampling_dat","y_ref","multi_time_scale","out_prefix", "out_path") #same names for any method |
|
| 393 |
#debug(calculate_accuracy_metrics) |
|
| 394 |
#test_val2 <-calculate_accuracy_metrics(2,list_param_validation) |
|
| 375 | 395 |
|
| 376 | 396 |
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 = 9) |
| 377 | 397 |
#test_val<-calculate_accuracy_metrics(1,list_param_validation) |
| 378 | 398 |
save(validation_mod_month_obj,file= file.path(out_path,paste(interpolation_method,"_validation_mod_month_obj_",y_var_name,out_prefix,".RData",sep=""))) |
| 379 | 399 |
|
| 400 |
##Create data.frame with validation and fit metrics for a full year/full numbe of runs |
|
| 401 |
tb_month_diagnostic_v<-extract_from_list_obj(validation_mod_month_obj,"metrics_v") |
|
| 402 |
#tb_diagnostic_v contains accuracy metrics for models sample and proportion for every run...if full year then 365 rows maximum |
|
| 403 |
rownames(tb_month_diagnostic_v)<-NULL #remove row names |
|
| 404 |
tb_month_diagnostic_v$method_interp <- interpolation_method |
|
| 405 |
tb_month_diagnostic_s<-extract_from_list_obj(validation_mod_month_obj,"metrics_s") |
|
| 406 |
rownames(tb_month_diagnostic_s)<-NULL #remove row names |
|
| 407 |
tb_month_diagnostic_s$method_interp <- interpolation_method #add type of interpolation...out_prefix too?? |
|
| 380 | 408 |
|
| 381 | 409 |
} |
| 382 | 410 |
#################### ASSESSMENT OF PREDICTIONS: PLOTS OF ACCURACY METRICS ########### |
| ... | ... | |
| 412 | 440 |
#Will add more information to be returned |
| 413 | 441 |
|
| 414 | 442 |
if (interpolation_method %in% c("gam_CAI","kriging_CAI","gwr_CAI","gam_fusion","kriging_fusion","gwr_fusion")){
|
| 415 |
raster_prediction_obj<-list(clim_method_mod_obj,method_mod_obj,validation_mod_obj,tb_diagnostic_v, |
|
| 416 |
tb_diagnostic_s,summary_metrics_v,summary_month_metrics_v) |
|
| 417 |
names(raster_prediction_obj)<-c("clim_method_mod_obj","method_mod_obj","validation_mod_obj","tb_diagnostic_v",
|
|
| 418 |
"tb_diagnostic_s","summary_metrics_v","summary_month_metrics_v") |
|
| 443 |
raster_prediction_obj<-list(clim_method_mod_obj,method_mod_obj,validation_mod_obj,validation_mod_month_obj, tb_diagnostic_v,
|
|
| 444 |
tb_diagnostic_s,tb_month_diagnostic_v,tb_month_diagnostic_s,summary_metrics_v,summary_month_metrics_v)
|
|
| 445 |
names(raster_prediction_obj)<-c("clim_method_mod_obj","method_mod_obj","validation_mod_obj","validation_mod_month_obj","tb_diagnostic_v",
|
|
| 446 |
"tb_diagnostic_s","tb_month_diagnostic_v","tb_month_diagnostic_s","summary_metrics_v","summary_month_metrics_v")
|
|
| 419 | 447 |
save(raster_prediction_obj,file= file.path(out_path,paste("raster_prediction_obj_",interpolation_method,"_", y_var_name,out_prefix,".RData",sep="")))
|
| 420 | 448 |
|
| 421 | 449 |
} |
Also available in: Unified diff
raster prediction script, validation call for monthly and daily hold out and outputs objects