Revision fc8d8661
Added by Benoit Parmentier about 10 years ago
| climate/research/oregon/interpolation/global_run_scalingup_assessment_part1.R | ||
|---|---|---|
| 5 | 5 |
#Part 1 create summary tables and inputs for figure in part 2 and part 3. |
| 6 | 6 |
#AUTHOR: Benoit Parmentier |
| 7 | 7 |
#CREATED ON: 03/23/2014 |
| 8 |
#MODIFIED ON: 08/14/2014
|
|
| 8 |
#MODIFIED ON: 08/25/2014
|
|
| 9 | 9 |
#Version: 3 |
| 10 | 10 |
#PROJECT: Environmental Layers project |
| 11 | 11 |
#TO DO: |
| ... | ... | |
| 88 | 88 |
extract_from_list_obj<-function(obj_list,list_name){
|
| 89 | 89 |
#extract object from list of list. This useful for raster_prediction_obj |
| 90 | 90 |
library(plyr) |
| 91 |
|
|
| 92 | 91 |
list_tmp<-vector("list",length(obj_list))
|
| 93 | 92 |
for (i in 1:length(obj_list)){
|
| 94 |
tmp<-obj_list[[i]][[list_name]] #double bracket to return data.frame |
|
| 95 |
list_tmp[[i]]<- as.data.frame(tmp) #if spdf |
|
| 93 |
tmp <- obj_list[[i]] |
|
| 94 |
if(inherits(tmp,"try-error")){
|
|
| 95 |
print(paste("no model generated or error in list",sep=" ")) #change message for any model type...
|
|
| 96 |
list_tmp[[i]] <- NULL #double bracket to return data.frame |
|
| 97 |
}else{
|
|
| 98 |
#tmp<-obj_list[[i]][[list_name]] #double bracket to return data.frame |
|
| 99 |
list_tmp[[i]] <- as.data.frame(tmp[[list_name]]) |
|
| 100 |
} |
|
| 101 |
# |
|
| 102 |
#tmp<-obj_list[[i]][[list_name]] #double bracket to return data.frame |
|
| 103 |
#list_tmp[[i]]<- as.data.frame(tmp) #if spdf |
|
| 96 | 104 |
} |
| 105 |
# |
|
| 106 |
#list_tmp <-list_tmp[!is.null(list_tmp)] |
|
| 107 |
list_tmp <- list_tmp[unlist(lapply(list_tmp,FUN=function(x){!is.null(x)}))]
|
|
| 108 |
|
|
| 97 | 109 |
tb_list_tmp<-do.call(rbind.fill,list_tmp) #long rownames |
| 98 |
#tb_list_tmp<-do.call(rbind,list_tmp) #long rownames |
|
| 99 |
|
|
| 110 |
#tb_list_tmp<-do.call(rbind,list_tmp) #long rownames |
|
| 100 | 111 |
return(tb_list_tmp) #this is a data.frame |
| 101 | 112 |
} |
| 102 | 113 |
|
| ... | ... | |
| 318 | 329 |
|
| 319 | 330 |
#/nobackupp4/aguzman4/climateLayers/output10Deg/reg1/finished.txt |
| 320 | 331 |
#in_dir_list <- list.dirs(path=in_dir1,recursive=FALSE) #get the list of directories with resutls by 10x10 degree tiles |
| 321 |
in_dir_list <- c("/nobackupp4/aguzman4/climateLayers/output20Deg/reg5/20.0_0.0/",
|
|
| 322 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg3/-20.0_-70.0/", |
|
| 323 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg5/20.0_30.0/", |
|
| 324 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg4/40.0_0.0/", |
|
| 325 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg5/20.0_-10.0/", |
|
| 326 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg4/50.0_0.0/", |
|
| 327 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg6/60.0_40.0/", |
|
| 328 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg6/30.0_40.0/") |
|
| 329 |
|
|
| 332 |
in_dir_list <- c( |
|
| 333 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg2//-10.0_-70.0/", |
|
| 334 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg4//40.0_0.0/", |
|
| 335 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg4//50.0_0.0/", |
|
| 336 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg6//60.0_40.0/", |
|
| 337 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg6//30.0_40.0/", |
|
| 338 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg8//40.0_130.0/") |
|
| 339 |
|
|
| 340 |
#Models used. |
|
| 341 |
#list_models<-c("y_var ~ s(lat,lon,k=4) + s(elev_s,k=3) + s(LST,k=3)",
|
|
| 342 |
# "y_var ~ s(lat,lon,k=5) + s(elev_s,k=3) + s(LST,k=3)", |
|
| 343 |
# "y_var ~ s(lat,lon,k=8) + s(elev_s,k=4) + s(LST,k=4)", |
|
| 344 |
|
|
| 330 | 345 |
#use subset for now: |
| 331 | 346 |
|
| 332 | 347 |
#in_dir_list <- c( |
| ... | ... | |
| 338 | 353 |
#in_dir_list <- in_dir_list[grep("bak",basename(basename(in_dir_list)),invert=TRUE)] #the first one is the in_dir1
|
| 339 | 354 |
#in_dir_shp <- in_dir_list[grep("shapefiles",basename(in_dir_list),invert=FALSE)] #select directory with shapefiles...
|
| 340 | 355 |
in_dir_shp <- c( |
| 341 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg3/subset/shapefiles/", |
|
| 342 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg5/subset/shapefiles/", |
|
| 356 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg2/subset/shapefiles/", |
|
| 343 | 357 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg4/subset/shapefiles/", |
| 344 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg5/subset/shapefiles/", |
|
| 345 | 358 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg4/subset/shapefiles/", |
| 346 | 359 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg6/subset/shapefiles/", |
| 347 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg6/subset/shapefiles/") |
|
| 360 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg6/subset/shapefiles/", |
|
| 361 |
"/nobackupp4/aguzman4/climateLayers/output20Deg/reg8/subset/shapefiles/") |
|
| 362 |
|
|
| 348 | 363 |
#in_dir_shp <- "/nobackupp4/aguzman4/climateLayers/output10Deg/reg1/subset/shapefiles/" |
| 349 | 364 |
#in_dir_shp <- "/nobackupp4/aguzman4/climateLayers/output20Deg/reg2/subset/shapefiles" |
| 350 | 365 |
in_dir_shp_list <- list.files(in_dir_shp,".shp",full.names=T) |
| ... | ... | |
| 354 | 369 |
# the last directory contains shapefiles |
| 355 | 370 |
y_var_name <- "dailyTmax" |
| 356 | 371 |
interpolation_method <- c("gam_CAI")
|
| 357 |
out_prefix<-"run4_global_analyses_08142014"
|
|
| 372 |
out_prefix<-"run5_global_analyses_08252014"
|
|
| 358 | 373 |
|
| 359 | 374 |
#out_dir<-"/data/project/layers/commons/NEX_data/" #On NCEAS Atlas |
| 360 | 375 |
out_dir <- "/nobackup/bparmen1/" #on NEX |
| ... | ... | |
| 398 | 413 |
lf_diagnostic_obj <- lapply(in_dir_list,FUN=function(x){list.files(path=x,pattern="diagnostics_.*.RData",full.names=T)})
|
| 399 | 414 |
lf_diagnostic_obj <- lf_diagnostic_obj[grep("lk_min",lf_diagnostic_obj,invert=T)] #remove object that have lk_min...
|
| 400 | 415 |
|
| 416 |
lf_validation_obj <- lapply(in_dir_list,FUN=function(x){list.files(path=x,pattern="gam_CAI_validation_mod_obj_dailyTmax.*.RData",full.names=T)})
|
|
| 417 |
#validation_mod_obj <-load_obj("/nobackupp4/aguzman4/climateLayers/output20Deg/reg6/60.0_40.0/gam_CAI_validation_mod_obj_dailyTmax60.0_40.0.RData")
|
|
| 418 |
debug(extract_from_list_obj) |
|
| 419 |
tb_diagnostic_v<-extract_from_list_obj(validation_mod_obj,"metrics_v") |
|
| 420 |
|
|
| 401 | 421 |
########################## START SCRIPT ############################## |
| 402 | 422 |
|
| 403 | 423 |
################################################################ |
| ... | ... | |
| 425 | 445 |
#Get the number of models predicted |
| 426 | 446 |
nb_mod <- length(unique(robj1$tb_diagnostic_v$pred_mod)) |
| 427 | 447 |
|
| 448 |
validation_mod_obj <-load_obj("/nobackupp4/aguzman4/climateLayers/output20Deg/reg6/60.0_40.0/gam_CAI_validation_mod_obj_dailyTmax60.0_40.0.RData")
|
|
| 449 |
tb_diagnostic_v<-extract_from_list_obj(validation_mod_obj,"metrics_v") |
|
| 450 |
|
|
| 451 |
|
|
| 452 |
#clim_method_mod_obj <- load_obj("/nobackupp4/aguzman4/climateLayers/output20Deg/reg6/60.0_40.0/gam_CAI_mod_dailyTmax60.0_40.0.RData")
|
|
| 453 |
#list_data_v <- extract_list_from_list_obj(clim_method_mod_obj,"data_month_v") #extract monthly testing/validation dataset |
|
| 454 |
#list_data_s <- extract_list_from_list_obj(clim_method_mod_obj,"data_month") #extract monthly training/fitting dataset |
|
| 455 |
#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 |
|
| 456 |
#list_sampling_dat <- extract_list_from_list_obj(clim_method_mod_obj,"sampling_month_dat") |
|
| 457 |
|
|
| 458 |
list_tb_diagnostic_v <- mclapply(lf_validation_obj,FUN=function(x){try( x<- load_obj(x)); try(extract_from_list_obj(x,"metrics_v"))},mc.preschedule=FALSE,mc.cores = 6)
|
|
| 459 |
names(list_tb_diagnostic_v) <- list_names_tile_id |
|
| 460 |
|
|
| 461 |
#undebug(extract_from_list_obj) |
|
| 462 |
#validation_mod_month_obj <-load_obj("/nobackupp4/aguzman4/climateLayers/output20Deg/reg6/60.0_40.0/gam_CAI_validation_mod_month_obj_dailyTmax60.0_40.0.RData")
|
|
| 463 |
#tb_diagnostic_v<-extract_from_list_obj(validation_mod_obj,"metrics_v") |
|
| 464 |
|
|
| 465 |
|
|
| 428 | 466 |
################ |
| 429 | 467 |
#### Table 1: Average accuracy metrics |
| 430 | 468 |
|
| ... | ... | |
| 508 | 546 |
write.table(gam_diagnostic_df, |
| 509 | 547 |
file=file.path(out_dir,paste("gam_diagnostic_df_",out_prefix,".txt",sep="")),sep=",")
|
| 510 | 548 |
|
| 549 |
|
|
| 550 |
#Now look at the 100 tiles of 10x10 |
|
| 551 |
#lf_test<-list.files("/nobackupp4/aguzman4/climateLayers/output10Deg/*/*/","diagnostics_obj_gam_fitting*")
|
|
| 552 |
lf_test <-list.files("/nobackupp4/aguzman4/climateLayers/output10Deg/","diagnostics_obj_gam_fitting.*.RData",recursive=T,full.names=T)
|
|
| 553 |
|
|
| 554 |
gam_diagnostic_10x10tb_list <- vector("list",length=length(lf_test))
|
|
| 555 |
lf_diagnostic_obj_tmp <- lf_test |
|
| 556 |
for(i in 1:length( lf_diagnostic_obj_tmp)){
|
|
| 557 |
l_diagnostic_obj_tmp <- lf_diagnostic_obj_tmp[[i]] |
|
| 558 |
tile_coord <- basename(dirname(lf_diagnostic_obj_tmp[i])) |
|
| 559 |
#l_diagnostic_obj_tmp <- l_diagnostic_obj_tmp[grep("lk_min",l_diagnostic_obj_tmp,invert=T)] #remove object that have lk_min...
|
|
| 560 |
l_diagnostic_obj_tmp_list <- lapply(l_diagnostic_obj_tmp,FUN=function(x){try(x<-load_obj(x));try(x[["df_diagnostics"]])})#,mc.preschedule=FALSE,mc.cores = 6)
|
|
| 561 |
gam_diagnostic_tb <- do.call(rbind.fill,l_diagnostic_obj_tmp_list)#create a df for NA tiles with all accuracy metrics |
|
| 562 |
gam_diagnostic_tb$tile_coord <- tile_coord |
|
| 563 |
gam_diagnostic_10x10tb_list[[i]] <- gam_diagnostic_tb |
|
| 564 |
} |
|
| 565 |
|
|
| 566 |
gam_diagnostic_10x10_df <- do.call(rbind.fill,gam_diagnostic_10x10tb_list) #create a df for NA tiles with all accuracy metrics |
|
| 567 |
|
|
| 568 |
list_tile_coord <- unique(gam_diagnostic_10x10_df$tile_coord) |
|
| 569 |
list_tile_id <- paste("tile_",1:length(list_tile_coord),sep="")
|
|
| 570 |
|
|
| 571 |
tile_id_df <- data.frame(tile_coord=list_tile_coord,tile_id=list_tile_id) |
|
| 572 |
gam_diagnostic_10x10_df <- merge(gam_diagnostic_10x10_df,tile_id_df,by="tile_coord") |
|
| 573 |
|
|
| 574 |
# write.table(gam_diagnostic_10x10_df, |
|
| 575 |
# file=file.path(out_dir,paste("gam_diagnostic_10x10_df_",out_prefix,".txt",sep="")),sep=",")
|
|
| 576 |
|
|
| 511 | 577 |
################# |
| 512 | 578 |
###Table 3: monthly station information with predictions for all tiles |
| 513 | 579 |
|
Also available in: Unified diff
run5 assessment NEX part1 using 3 different k sets of values, adding functions