Revision 298461c0
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/global_run_scalingup_assessment_part5.R | ||
|---|---|---|
| 283 | 283 |
stat_val <- c(std_dev_val,mean_val,median_val,max_val,min_val,n_val) |
| 284 | 284 |
df_stat <- data.frame(stat_name=stat_name,stat_val=stat_val) |
| 285 | 285 |
|
| 286 |
threshold_val <- c(5,6,10) |
|
| 287 |
n_threshold_val <- sum((tb_subset[[metric_name]]) > threshold_val[1]) |
|
| 288 |
100*n_threshold_val/n_val |
|
| 289 |
|
|
| 290 |
n_threshold_val <- sum((tb_subset[[metric_name]]) > threshold_val[2]) |
|
| 291 |
100*n_threshold_val/n_val |
|
| 292 |
|
|
| 293 |
n_threshold_val <- sum((tb_subset[[metric_name]]) > threshold_val[3]) |
|
| 294 |
100*n_threshold_val/n_val |
|
| 295 |
|
|
| 296 |
#Find out where these values are located...by mapping extremes! |
|
| 297 |
|
|
| 298 |
|
|
| 286 |
model_name <- c("mod1","mod_kr")
|
|
| 287 |
threshold_val <- c(5,10,20,50) |
|
| 299 | 288 |
|
| 300 | 289 |
########################### |
| 301 | 290 |
### Figure 1: plot location of the study area with tiles processed |
| ... | ... | |
| 304 | 293 |
#list_shp_reg_files <- df_tiled_processed$shp_files |
| 305 | 294 |
|
| 306 | 295 |
#list_shp_reg_files<- as.character(df_tile_processed$shp_files) #this could be the solution!! |
| 296 |
#Use melt!! quick solution |
|
| 307 | 297 |
list_shp_reg_files <- as.character(basename(unique(df_tile_processed$shp_files))) #this could be the solution!! |
| 308 |
#the shapefiles must be copied in the proper folder!!! |
|
| 309 |
#list_shp_reg_files<- file.path(in_dir,in_dir_list[1],"shapefile",list_shp_reg_files) |
|
| 310 |
#list_shp_reg_files <- file.path("/data/project/layers/commons/NEX_data/",out_dir,
|
|
| 311 |
# as.character(df_tile_processed$tile_coord),"shapefiles",basename(list_shp_reg_files)) |
|
| 312 |
#list_shp_reg_files <- file.path("/data/project/layers/commons/NEX_data/",out_dir,
|
|
| 313 |
#"shapefiles",basename(list_shp_reg_files)) |
|
| 314 |
|
|
| 298 |
list_tile_id <- as.character((unique(df_tile_processed$tile_id))) #this is in the order of appearance |
|
| 299 |
#list_tile_lat <- as.character((unique(df_tile_processed$lat))) #this is in the order of appearance |
|
| 300 |
#list_tile_lon <- as.character((unique(df_tile_processed$lon))) #this is in the order of appearance |
|
| 301 |
|
|
| 302 |
#df_tile_processed$shp_files2 <- basename(df_tile_processed$shp_files) |
|
| 303 |
df_tiles_reg <- data.frame(shp_files=(list_shp_reg_files),tile_id=list_tile_id) |
|
| 304 |
df_tiles_reg$tile_id <- as.character(df_tiles_reg$tile_id) |
|
| 305 |
#,lat=list_tile_lat,lon=list_tile_lon) |
|
| 306 |
#cast(df_tile_processed, shp_files ~ tile_id+lat+lon, mean, value = 'income') |
|
| 307 |
|
|
| 315 | 308 |
### Potential function starts here: |
| 316 | 309 |
#function(in_dir,out_dir,list_shp_reg_files,title_str,region_name,num_cores,out_suffix,out_suffix) |
| 317 | 310 |
|
| ... | ... | |
| 367 | 360 |
centroids_pts <- remove_errors_list(centroids_pts) #[[!inherits(shps_tiles,"try-error")]] |
| 368 | 361 |
#centroids_pts <- tmp_pts |
| 369 | 362 |
|
| 370 |
#plot info: with labels |
|
| 363 |
df_pts <- as.data.frame(do.call(rbind,tmp_pts)) |
|
| 364 |
#df_pts <- cbind(df_pts,df_tiles_reg) |
|
| 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]}))
|
|
| 367 |
|
|
| 368 |
#plot info: with labels |
|
| 371 | 369 |
res_pix <-1200 |
| 372 | 370 |
col_mfrow <- 1 |
| 373 | 371 |
row_mfrow <- 1 |
| 374 | 372 |
|
| 375 |
png(filename=paste("Figure1_tile_processed_region_",region_name,"_",out_suffix,".png",sep=""),
|
|
| 373 |
png(filename=paste("Figure1a_tile_processed_region_",region_name,"_",out_suffix,".png",sep=""),
|
|
| 376 | 374 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
| 377 | 375 |
|
| 378 | 376 |
plot(reg_layer) |
| ... | ... | |
| 399 | 397 |
|
| 400 | 398 |
dev.off() |
| 401 | 399 |
|
| 400 |
|
|
| 402 | 401 |
list_outfiles[[counter_fig+1]] <- paste("Figure1_tile_processed_region_",region_name,"_",out_suffix,".png",sep="")
|
| 403 | 402 |
counter_fig <- counter_fig+1 |
| 404 | 403 |
#this will be changed to be added to data.frame on the fly |
| ... | ... | |
| 411 | 410 |
## Figure 2a |
| 412 | 411 |
|
| 413 | 412 |
#Plot location of extremes and select them for further analyses? |
| 414 |
|
|
| 415 | 413 |
|
| 416 |
## Number of tiles with information: |
|
| 417 |
sum(df_tile_processed$metrics_v) #26,number of tiles with raster object |
|
| 418 |
length(df_tile_processed$metrics_v) #26,number of tiles in the region |
|
| 419 |
sum(df_tile_processed$metrics_v)/length(df_tile_processed$metrics_v) #80 of tiles with info |
|
| 420 |
|
|
| 421 |
#coordinates |
|
| 422 |
try(coordinates(summary_metrics_v) <- c("lon","lat"))
|
|
| 423 |
#try(coordinates(summary_metrics_v) <- c("lon.y","lat.y"))
|
|
| 424 |
|
|
| 425 |
#threshold_missing_day <- c(367,365,300,200) |
|
| 426 |
|
|
| 427 |
nb<-nrow(subset(summary_metrics_v,model_name=="mod1")) |
|
| 428 |
sum(subset(summary_metrics_v,model_name=="mod1")$n_missing)/nb #33/35 |
|
| 429 |
|
|
| 430 | 414 |
|
| 431 | 415 |
j<-1 #for model name 1,mod1 |
| 432 |
model_name <- c("mod1","mod_kr")
|
|
| 416 |
#model_name <- c("mod1","mod_kr")
|
|
| 417 |
#threshold_val <- c(5,10,20,50) |
|
| 418 |
|
|
| 433 | 419 |
for(i in 1:length(threshold_val)){
|
| 434 | 420 |
|
| 421 |
tb_tmp <- try(subset(tb_subset,tb_subset[[metric_name]]>threshold_val[i])) |
|
| 422 |
hist(tb_tmp$rmse) |
|
| 435 | 423 |
|
| 436 |
tb_subset <- subset(tb_subset,tb_subset[[metric_name]]>threshold_val[i]) |
|
| 437 |
|
|
| 438 |
df_extremes <- as.data.frame(table(tb_subset$tile_id)) |
|
| 439 |
names(df_extremes)<- c("tile_id","freq_extremes")
|
|
| 440 |
tb_sorted <- merge(tb_subset,df_extremes,"tile_id") |
|
| 441 |
tb_sorted <- arrange(tb_sorted,desc(freq_extremes)) #[,c("pred_mod","rmse","mae","tile_id")]
|
|
| 442 |
coordinates(tb_sorted) <- c("lon","lat")
|
|
| 424 |
if(nrow(tb_subset)>0){
|
|
| 425 |
|
|
| 426 |
df_extremes <- as.data.frame(table(tb_tmp$tile_id)) |
|
| 427 |
names(df_extremes)<- c("tile_id","freq_extremes")
|
|
| 428 |
tb_sorted <- merge(tb_tmp,df_extremes,"tile_id") |
|
| 429 |
tb_sorted <- arrange(tb_sorted,desc(freq_extremes)) #[,c("pred_mod","rmse","mae","tile_id")]
|
|
| 430 |
coordinates(tb_sorted) <- c("lon","lat")
|
|
| 431 |
df_extremes <- arrange(df_extremes,desc(freq_extremes)) |
|
| 443 | 432 |
|
| 444 |
fig_filename <- paste("Figure2a_barplot_extremes_val_centroids_tile_",model_name[j],"_",threshold_val[i],
|
|
| 433 |
fig_filename <- paste("Figure2a_barplot_extremes_val_centroids_tile_",model_name[j],"_",threshold_val[i],
|
|
| 445 | 434 |
"_",out_suffix,".png",sep="") |
| 435 |
|
|
| 436 |
#res_pix <- 1200 |
|
| 437 |
res_pix <- 960 |
|
| 438 |
col_mfrow <- 1 |
|
| 439 |
row_mfrow <- 1 |
|
| 440 |
#only mod1 right now |
|
| 441 |
png(filename=fig_filename, |
|
| 442 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 443 |
|
|
| 444 |
#barplot(df_ac_mod$rmse, names.arg=x) |
|
| 445 |
barplot(df_extremes$freq_extremes,main=paste("Extremes threshold val for ",threshold_val[i],"deg C",sep=""),
|
|
| 446 |
ylab="frequency",xlab="tile_id",names.arg=df_extremes$tile_id,las=2) |
|
| 447 |
#barplot(table(tb_subset$tile_id),main=paste("Extremes threshold val for ",threshold_val[i],sep=""),
|
|
| 448 |
# ylab="frequency",xlab="tile_id",las=2) |
|
| 449 |
dev.off() |
|
| 450 |
list_outfiles[[counter_fig+i]] <- fig_filename |
|
| 446 | 451 |
|
| 447 |
barplot(table(tb_subset$tile_id),main=paste("Extremes threshold val for ",threshold_val[i],sep=""),
|
|
| 448 |
ylab="frequency",xlab="tile_id",las=2) |
|
| 449 |
dev.off() |
|
| 450 | 452 |
|
| 451 |
test<-(subset(tb,tb$tile_id==unique(tb_subset$tile_id))) |
|
| 452 |
#df_ac_mod <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("pred_mod","rmse","mae","tile_id")]
|
|
| 453 |
#test<-(subset(tb,tb$tile_id==unique(tb_subset$tile_id)))
|
|
| 454 |
#df_ac_mod <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("pred_mod","rmse","mae","tile_id")]
|
|
| 453 | 455 |
|
| 454 |
#plot top three, then all,and histogram...make this a function... |
|
| 455 |
list_df_ac_mod[[i]] <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("rmse","mae","tile_id")]
|
|
| 456 |
#plot top three, then all,and histogram...make this a function...
|
|
| 457 |
#list_df_ac_mod[[i]] <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("rmse","mae","tile_id")]
|
|
| 456 | 458 |
|
| 457 |
#summary_metrics_v$n_missing <- summary_metrics_v$n == 365 |
|
| 458 |
#summary_metrics_v$n_missing <- summary_metrics_v$n < 365 |
|
| 459 |
#summary_metrics_v$n_missing <- as.numeric(summary_metrics_v$n < threshold_missing_day[i]) |
|
| 460 |
#summary_metrics_v_subset <- subset(summary_metrics_v,model_name=="mod1") |
|
| 461 |
|
|
| 462 |
fig_filename <- paste("Figure7a_ac_metrics_extremes_map_centroids_tile_",model_name[j],"_",threshold_val[i],
|
|
| 459 |
fig_filename <- paste("Figure2b_ac_metrics_extremes_map_centroids_tile_",model_name[j],"_",threshold_val[i],
|
|
| 463 | 460 |
"_",out_suffix,".png",sep="") |
| 464 |
list_outfiles[[counter_fig+i]] <- fig_filename |
|
| 465 | 461 |
|
| 466 |
if(sum(summary_metrics_v_subset$n_missing) > 0){#then there are centroids to plot!!!
|
|
| 467 |
|
|
| 468 | 462 |
#res_pix <- 1200 |
| 469 | 463 |
res_pix <- 960 |
| 470 | 464 |
col_mfrow <- 1 |
| 471 | 465 |
row_mfrow <- 1 |
| 472 | 466 |
#only mod1 right now |
| 473 |
png(filename=paste("Figure7a_ac_metrics_map_centroids_tile_",model_name[j],"_","missing_day_",threshold_missing_day[i],
|
|
| 474 |
"_",out_suffix,".png",sep=""), |
|
| 467 |
png(filename=fig_filename, |
|
| 475 | 468 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
| 476 | 469 |
|
| 477 | 470 |
model_name[j] |
| ... | ... | |
| 485 | 478 |
try(print(p1)) #error raised if number of missing values below a threshold does not exist |
| 486 | 479 |
dev.off() |
| 487 | 480 |
|
| 488 |
} |
|
| 489 |
|
|
| 481 |
} |
|
| 482 |
|
|
| 483 |
|
|
| 490 | 484 |
#list_outfiles[[counter_fig+i]] <- fig_filename |
| 491 | 485 |
} |
| 492 | 486 |
counter_fig <- counter_fig+length(threshold_missing_day) #currently 4 days... |
Also available in: Unified diff
analyses of extreme values part 5 assessment, examining frequency of extremes by tiles