Revision 5182652d
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/global_run_scalingup_assessment_part5.R | ||
|---|---|---|
| 448 | 448 |
j<-1 #for model name 1,mod1 |
| 449 | 449 |
#model_name <- c("mod1","mod_kr")
|
| 450 | 450 |
#threshold_val <- c(5,10,20,50) |
| 451 |
|
|
| 451 |
list_tb_extremes <- vector("list",length=length(threshold_val))
|
|
| 452 | 452 |
for(i in 1:length(threshold_val)){
|
| 453 |
|
|
| 453 |
#formula_str <- paste0(metric_name,"~","date+tile_id+lat+lon") |
|
| 454 |
#test2 <- aggregate(formula_str,tb_subset,min) |
|
| 454 | 455 |
tb_tmp <- try(subset(tb_subset,tb_subset[[metric_name]]>threshold_val[i])) |
| 456 |
#test2 <- aggregate(rmse~date,test,min) |
|
| 457 |
|
|
| 455 | 458 |
hist(tb_tmp$rmse) |
| 456 | 459 |
fig_filename1 <-paste("Figure2a_barplot_extremes_val_centroids_tile_",model_name[j],"_",threshold_val[i],
|
| 457 | 460 |
"_",out_suffix,".png",sep="") |
| ... | ... | |
| 460 | 463 |
list_outfiles[[counter_fig+i]] <- fig_filename1 |
| 461 | 464 |
list_outfiles[[counter_fig+i]] <- fig_filename2 |
| 462 | 465 |
|
| 463 |
if(nrow(tb_subset)>0){
|
|
| 464 |
|
|
| 465 |
df_extremes <- as.data.frame(table(tb_tmp$tile_id)) |
|
| 466 |
if(nrow(tb_tmp)>0){
|
|
| 467 |
table(tb_sorted$date) |
|
| 468 |
df_extremes <- as.data.frame(table(tb_tmp$tile_id)/17)
|
|
| 466 | 469 |
names(df_extremes)<- c("tile_id","freq_extremes")
|
| 467 | 470 |
tb_sorted <- merge(tb_tmp,df_extremes,"tile_id") |
| 468 | 471 |
tb_sorted <- arrange(tb_sorted,desc(freq_extremes)) #[,c("pred_mod","rmse","mae","tile_id")]
|
| ... | ... | |
| 509 | 512 |
#p_shp <- layer(sp.polygons(reg_layer, lwd=1, col='black')) |
| 510 | 513 |
p_shp <- spplot(reg_layer,"ISO3" ,col.regions=NA, col="black") #ok problem solved!! |
| 511 | 514 |
#title("(a) Mean for 1 January")
|
| 515 |
#tb_sorted$freq_extremes2 <- tb_sorted$freq_extremes/17 |
|
| 512 | 516 |
p <- bubble(tb_sorted,"freq_extremes",main=paste("Extremes per tile and by ",model_name[j]," for ",
|
| 513 | 517 |
threshold_val[i])) |
| 514 | 518 |
p1 <- p+p_shp |
| 515 | 519 |
try(print(p1)) #error raised if number of missing values below a threshold does not exist |
| 516 | 520 |
dev.off() |
| 517 | 521 |
|
| 522 |
}else{
|
|
| 523 |
df_extremes <- NULL |
|
| 524 |
tb_sorted <- NULL |
|
| 518 | 525 |
} |
| 519 |
|
|
| 526 |
extremes_obj <- list(tb_tmp,df_extremes,tb_sorted) |
|
| 527 |
names(extremes_obj) <- c("tb_tmp","df_extremes","tb_sorted")
|
|
| 528 |
list_tb_extremes[[i]] <- extremes_obj |
|
| 529 |
|
|
| 520 | 530 |
#list_outfiles[[counter_fig+i]] <- fig_filename |
| 521 | 531 |
} |
| 522 | 532 |
counter_fig <- counter_fig+length(threshold_missing_day) #currently 4 days... |
| ... | ... | |
| 528 | 538 |
|
| 529 | 539 |
##### Figure 3 ### |
| 530 | 540 |
|
| 531 |
test<- subset(tb_subset,tb_subset$tile_id=="tile_14") |
|
| 541 |
|
|
| 542 |
for(i in 1:length(tile_selected_extremes)){
|
|
| 543 |
d |
|
| 544 |
d |
|
| 545 |
d |
|
| 546 |
test<- subset(tb_subset,tb_subset$tile_id=="tile_14") |
|
| 532 | 547 |
test2 <- aggregate(rmse~date,test,min) |
| 533 |
idx <- test$date #transform this format... |
|
| 534 |
d_z <- zoo(test,idx) #make a time series ... |
|
| 548 |
#idx <- test2$date #transform this format... |
|
| 549 |
|
|
| 550 |
idx <- as.Date(strptime(test2$date, "%Y%m%d")) # interpolation date being processed |
|
| 551 |
|
|
| 552 |
|
|
| 553 |
d_z <- zoo(test2,idx) #make a time series ... |
|
| 535 | 554 |
#add horizontal line... |
| 555 |
plot(d_z[[metric_name]]) |
|
| 556 |
plot(d_z$rmse,ylab=metric_name,xlab="dates") |
|
| 557 |
title(title_str) |
|
| 536 | 558 |
|
| 537 | 559 |
plot(test$rmse,type="b") |
| 538 |
unique(test$year_predicted) |
|
| 539 |
unique(tb$year_predicted) |
|
| 560 |
length(unique(test$year_predicted)) |
|
| 561 |
unique(tb$year_predicted) |
|
| 562 |
|
|
| 563 |
} |
|
| 540 | 564 |
|
| 541 | 565 |
###################################################### |
| 542 | 566 |
##### Prepare objet to return #### |
Also available in: Unified diff
analyses of extremes, debugging assessment part5