Revision da7bcc4a
Added by Benoit Parmentier about 11 years ago
| climate/research/oregon/interpolation/GAM_fusion_analysis_raster_prediction_multisampling.R | ||
|---|---|---|
| 54 | 54 |
#24)out_path |
| 55 | 55 |
#25)script_path: path to script |
| 56 | 56 |
#26)interpolation_method: c("gam_fusion","gam_CAI") #other otpions to be added later
|
| 57 |
#27) use_clim_image |
|
| 58 |
#28) join_daily |
|
| 57 | 59 |
|
| 58 | 60 |
###Loading R library and packages |
| 59 | 61 |
|
| ... | ... | |
| 116 | 118 |
interpolation_method<-list_param_raster_prediction$interpolation_method |
| 117 | 119 |
screen_data_training <-list_param_raster_prediction$screen_data_training |
| 118 | 120 |
|
| 121 |
use_clim_image <- list_param_raster_prediction$use_clim_image # use predicted image as a base...rather than average Tmin at the station for delta |
|
| 122 |
join_daily <- list_param_raster_prediction$join_daily # join monthly and daily station before calucating delta |
|
| 123 |
|
|
| 119 | 124 |
setwd(out_path) |
| 120 | 125 |
|
| 121 | 126 |
###################### START OF THE SCRIPT ######################## |
| ... | ... | |
| 189 | 194 |
dates_month<-as.character(sort(unique(dst$date))) |
| 190 | 195 |
list_param_sampling<-list(seed_number_month,nb_sample_month,step_month,constant_month,prop_minmax_month,dates_month,dst) |
| 191 | 196 |
#list_param_sampling<-list(seed_number,nb_sample,step,constant,prop_minmax,dates,ghcn_name) |
| 192 |
|
|
| 193 | 197 |
names(list_param_sampling)<-c("seed_number","nb_sample","step","constant","prop_minmax","dates","ghcn")
|
| 194 | 198 |
|
| 195 | 199 |
sampling_month_obj <- sampling_training_testing(list_param_sampling) |
| ... | ... | |
| 257 | 261 |
|
| 258 | 262 |
daily_dev_sampling_dat <- combine_sampling_daily_monthly_for_daily_deviation_pred(sampling_obj,sampling_month_obj) |
| 259 | 263 |
|
| 260 |
use_clim_image<- TRUE |
|
| 264 |
#use_clim_image<- TRUE
|
|
| 261 | 265 |
#use_clim_image<-FALSE |
| 262 |
join_daily <- FALSE |
|
| 266 |
#join_daily <- FALSE |
|
| 267 |
#join_daily <- TRUE |
|
| 263 | 268 |
|
| 264 | 269 |
if (interpolation_method %in% c("gam_CAI","kriging_CAI","gwr_CAI","gam_fusion","kriging_fusion","gwr_fusion")){
|
| 265 | 270 |
#input a list:note that ghcn.subsets is not sampling_obj$data_day_ghcn |
| ... | ... | |
| 268 | 273 |
use_clim_image,join_daily,var,y_var_name, interpolation_method,out_prefix,out_path) |
| 269 | 274 |
names(list_param_run_prediction_daily_deviation)<-c("list_index","clim_yearlist","daily_dev_sampling_dat","sampling_month_obj","sampling_obj","dst",
|
| 270 | 275 |
"use_clim_image","join_daily","var","y_var_name","interpolation_method","out_prefix","out_path") |
| 271 |
#test<-mclapply(1:18, runGAMFusion,list_param=list_param_runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) |
|
| 272 |
#test<-runGAMFusion(1,list_param=list_param_runGAMFusion) |
|
| 273 | 276 |
#method_mod_obj<-mclapply(1:length(sampling_obj$ghcn_data),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 |
| 277 |
#debug(run_prediction_daily_deviation) |
|
| 274 | 278 |
#test <- run_prediction_daily_deviation(1,list_param=list_param_run_prediction_daily_deviation) #This is the end bracket from mclapply(...) statement |
| 275 | 279 |
#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 = 9) #This is the end bracket from mclapply(...) statement |
| 276 | 280 |
|
| ... | ... | |
| 326 | 330 |
############### NOW RUN VALIDATION ######################### |
| 327 | 331 |
#SIMPLIFY THIS PART: one call |
| 328 | 332 |
|
| 329 |
list_tmp<-vector("list",length(method_mod_obj))
|
|
| 330 |
for (i in 1:length(method_mod_obj)){
|
|
| 331 |
tmp<-method_mod_obj[[i]][[y_var_name]] #y_var_name is the variable predicted (dailyTmax or dailyTmin) |
|
| 332 |
list_tmp[[i]]<-tmp |
|
| 333 |
} |
|
| 334 |
rast_day_yearlist<-list_tmp #list of predicted images over full year... |
|
| 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") |
|
| 335 | 344 |
|
| 336 | 345 |
cat("Validation step:",file=log_fname,sep="\n", append=TRUE)
|
| 337 | 346 |
t1<-proc.time() |
| 338 | 347 |
cat(paste("Local Date and Time: ",as.character(Sys.time()),sep=""),
|
| 339 | 348 |
file=log_fname,sep="\n") |
| 340 | 349 |
|
| 341 |
list_param_validation<-list(i,rast_day_yearlist,method_mod_obj,y_var_name, out_prefix, out_path) |
|
| 342 |
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
|
|
| 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) |
|
| 343 | 363 |
|
| 344 | 364 |
validation_mod_obj <-mclapply(1:length(method_mod_obj), list_param=list_param_validation, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 9) |
| 345 |
#test_val<-calculate_accuracy_metrics(1,list_param_validation) |
|
| 346 | 365 |
save(validation_mod_obj,file= file.path(out_path,paste(interpolation_method,"_validation_mod_obj_",y_var_name,out_prefix,".RData",sep=""))) |
| 347 | 366 |
t2<-proc.time()-t1 |
| 348 | 367 |
cat(as.character(t2),file=log_fname,sep="\n", append=TRUE) |
| 349 | 368 |
|
| 369 |
### Run monthly validation if multi-time scale methods and add information to daily... |
|
| 370 |
|
|
| 371 |
if (interpolation_method %in% c("gam_CAI","kriging_CAI","gwr_CAI","gam_fusion","kriging_fusion","gwr_fusion")){
|
|
| 372 |
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
|
|
| 375 |
|
|
| 376 |
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 |
#test_val<-calculate_accuracy_metrics(1,list_param_validation) |
|
| 378 |
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 |
|
|
| 380 |
|
|
| 381 |
} |
|
| 350 | 382 |
#################### ASSESSMENT OF PREDICTIONS: PLOTS OF ACCURACY METRICS ########### |
| 351 | 383 |
|
| 352 | 384 |
##Create data.frame with validation and fit metrics for a full year/full numbe of runs |
Also available in: Unified diff
raster script prediction modifications to deal with for monthly hold out