Revision e0928996
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/global_run_scalingup_assessment_part5.R | ||
|---|---|---|
| 363 | 363 |
df_pts <- as.data.frame(do.call(rbind,tmp_pts)) |
| 364 | 364 |
#df_pts <- cbind(df_pts,df_tiles_reg) |
| 365 | 365 |
df_tiles_reg <- cbind(df_pts,df_tiles_reg) #(shp_files=(list_shp_reg_files),tile_id=list_tile_id) |
| 366 |
df_tiles_reg$id <- unlist(lapply(strsplit(df_tiles_reg$tile_id,"_"),FUN=function(x){x[2]}))
|
|
| 366 |
df_tiles_reg$id <- as.numeric(unlist(lapply(strsplit(df_tiles_reg$tile_id,"_"),FUN=function(x){x[2]})))
|
|
| 367 | 367 |
|
| 368 |
#plot info: with labels |
|
| 368 |
coordinates(df_tiles_reg)<- cbind(df_tiles_reg$x,df_tiles_reg$y) |
|
| 369 |
|
|
| 370 |
#plot info: with labels |
|
| 369 | 371 |
res_pix <-1200 |
| 370 | 372 |
col_mfrow <- 1 |
| 371 | 373 |
row_mfrow <- 1 |
| ... | ... | |
| 397 | 399 |
|
| 398 | 400 |
dev.off() |
| 399 | 401 |
|
| 402 |
res_pix <-1200 |
|
| 403 |
col_mfrow <- 1 |
|
| 404 |
row_mfrow <- 1 |
|
| 400 | 405 |
|
| 406 |
png(filename=paste("Figure1b_tile_processed_centroids_region_",region_name,"_",out_suffix,".png",sep=""),
|
|
| 407 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 408 |
|
|
| 409 |
#plot(reg_layer) |
|
| 410 |
#Add polygon tiles... |
|
| 411 |
#title(paste("Tiles ", tile_size,region_name,sep=""))
|
|
| 412 |
#plot(df_tiles_reg,add=T,pch=2) |
|
| 413 |
#label_id <- df_tiles_reg$id |
|
| 414 |
#text(coordinates(df_tiles_reg)[1],coordinates(df_tiles_reg)[2],labels=i,cex=1.3,font=2,col=c("red"),family="HersheySerif")
|
|
| 415 |
|
|
| 416 |
p_shp <- spplot(reg_layer,"ISO3" ,col.regions=NA, col="black") #ok problem solved!! |
|
| 417 |
#title("(a) Mean for 1 January")
|
|
| 418 |
df_tiles_reg$lab <- 1 |
|
| 419 |
sl1 <- list('sp.points',df_tiles_reg, pch=19, cex=.8, col='red')
|
|
| 420 |
sl2 <- list('sp.pointLabel', df_tiles_reg, label=list_id,
|
|
| 421 |
cex=2.4, font=2,col='red',col.regions="red", |
|
| 422 |
fontfamily='Palatino') #Add labels at centroids |
|
| 423 |
|
|
| 424 |
p <- spplot(df_tiles_reg,"id",main=paste("Tile id processed",sep=""),sp.layout=list(sl1, sl2))
|
|
| 425 |
#spplot(meuse.grid["dist"], col.regions=myCols, sp.layout=list(sl1, sl2) |
|
| 426 |
#p <- spplot(df_tiles_reg,"lab",main=paste("Tile id processed",sep=""))
|
|
| 427 |
p1 <- p+p_shp |
|
| 428 |
try(print(p1)) #error raised if number of missing values below a threshold does not exist |
|
| 429 |
#ltext(coordinates(df_tiles_reg)[,1],coordinates(df_tiles_reg)[,2],labels=list_tile_id) |
|
| 430 |
|
|
| 431 |
#list_id <- df_tiles_reg$id |
|
| 432 |
dev.off() |
|
| 433 |
|
|
| 401 | 434 |
list_outfiles[[counter_fig+1]] <- paste("Figure1_tile_processed_region_",region_name,"_",out_suffix,".png",sep="")
|
| 402 | 435 |
counter_fig <- counter_fig+1 |
| 403 | 436 |
#this will be changed to be added to data.frame on the fly |
| ... | ... | |
| 420 | 453 |
|
| 421 | 454 |
tb_tmp <- try(subset(tb_subset,tb_subset[[metric_name]]>threshold_val[i])) |
| 422 | 455 |
hist(tb_tmp$rmse) |
| 423 |
|
|
| 456 |
fig_filename1 <-paste("Figure2a_barplot_extremes_val_centroids_tile_",model_name[j],"_",threshold_val[i],
|
|
| 457 |
"_",out_suffix,".png",sep="") |
|
| 458 |
fig_filename2 <- paste("Figure2b_ac_metrics_extremes_map_centroids_tile_",model_name[j],"_",threshold_val[i],
|
|
| 459 |
"_",out_suffix,".png",sep="") |
|
| 460 |
list_outfiles[[counter_fig+i]] <- fig_filename1 |
|
| 461 |
list_outfiles[[counter_fig+i]] <- fig_filename2 |
|
| 462 |
|
|
| 424 | 463 |
if(nrow(tb_subset)>0){
|
| 425 | 464 |
|
| 426 | 465 |
df_extremes <- as.data.frame(table(tb_tmp$tile_id)) |
| ... | ... | |
| 447 | 486 |
#barplot(table(tb_subset$tile_id),main=paste("Extremes threshold val for ",threshold_val[i],sep=""),
|
| 448 | 487 |
# ylab="frequency",xlab="tile_id",las=2) |
| 449 | 488 |
dev.off() |
| 450 |
list_outfiles[[counter_fig+i]] <- fig_filename |
|
| 451 |
|
|
| 452 | 489 |
|
| 453 | 490 |
#test<-(subset(tb,tb$tile_id==unique(tb_subset$tile_id))) |
| 454 | 491 |
#df_ac_mod <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("pred_mod","rmse","mae","tile_id")]
|
| ... | ... | |
| 480 | 517 |
|
| 481 | 518 |
} |
| 482 | 519 |
|
| 483 |
|
|
| 484 | 520 |
#list_outfiles[[counter_fig+i]] <- fig_filename |
| 485 | 521 |
} |
| 486 | 522 |
counter_fig <- counter_fig+length(threshold_missing_day) #currently 4 days... |
| ... | ... | |
| 490 | 526 |
r20 <-c("figure_7","Number of missing days threshold3 map at centroids","mod1",metric_name,region_name,year_predicted,list_outfiles[[20]])
|
| 491 | 527 |
r21 <-c("figure_7","Number of missing days threshold4 map at centroids","mod1",metric_name,region_name,year_predicted,list_outfiles[[21]])
|
| 492 | 528 |
|
| 529 |
##### Figure 3 ### |
|
| 530 |
|
|
| 531 |
test<- subset(tb_subset,tb_subset$tile_id=="tile_14") |
|
| 532 |
test2 <- aggregate(rmse~date,test,min) |
|
| 533 |
idx <- test$date #transform this format... |
|
| 534 |
d_z <- zoo(test,idx) #make a time series ... |
|
| 535 |
#add horizontal line... |
|
| 536 |
|
|
| 537 |
plot(test$rmse,type="b") |
|
| 538 |
unique(test$year_predicted) |
|
| 539 |
unique(tb$year_predicted) |
|
| 540 |
|
|
| 493 | 541 |
###################################################### |
| 494 | 542 |
##### Prepare objet to return #### |
| 495 | 543 |
|
Also available in: Unified diff
analyses of extremes, fixing figures for assessment part 5