Revision 599839d9
Added by Benoit Parmentier almost 10 years ago
| climate/research/oregon/interpolation/global_run_scalingup_assessment_part2.R | ||
|---|---|---|
| 5 | 5 |
#Analyses, figures, tables and data are also produced in the script. |
| 6 | 6 |
#AUTHOR: Benoit Parmentier |
| 7 | 7 |
#CREATED ON: 03/23/2014 |
| 8 |
#MODIFIED ON: 01/02/2015
|
|
| 8 |
#MODIFIED ON: 01/20/2015
|
|
| 9 | 9 |
#Version: 4 |
| 10 | 10 |
#PROJECT: Environmental Layers project |
| 11 | 11 |
#COMMENTS: analyses for run 10 global analyses, Europe, Australia, 1000x300km |
| ... | ... | |
| 64 | 64 |
return(out_dir) |
| 65 | 65 |
} |
| 66 | 66 |
|
| 67 |
#Remove models that were not fitted from the list |
|
| 68 |
#All modesl that are "try-error" are removed |
|
| 69 |
remove_errors_list<-function(list_items){
|
|
| 70 |
|
|
| 71 |
#This function removes "error" items in a list |
|
| 72 |
list_tmp<-list_items |
|
| 73 |
if(is.null(names(list_tmp))){
|
|
| 74 |
names(list_tmp) <- paste("l",1:length(list_tmp),sep="_")
|
|
| 75 |
names(list_items) <- paste("l",1:length(list_tmp),sep="_")
|
|
| 76 |
} |
|
| 77 |
|
|
| 78 |
for(i in 1:length(list_items)){
|
|
| 79 |
if(inherits(list_items[[i]],"try-error")){
|
|
| 80 |
list_tmp[[i]]<-0 |
|
| 81 |
}else{
|
|
| 82 |
list_tmp[[i]]<-1 |
|
| 83 |
} |
|
| 84 |
} |
|
| 85 |
cnames<-names(list_tmp[list_tmp>0]) |
|
| 86 |
x <- list_items[match(cnames,names(list_items))] |
|
| 87 |
return(x) |
|
| 88 |
} |
|
| 89 |
|
|
| 90 |
#turn term from list into data.frame |
|
| 91 |
#name_col<-function(i,x){
|
|
| 92 |
#x_mat<-x[[i]] |
|
| 93 |
#x_df<-as.data.frame(x_mat) |
|
| 94 |
#x_df$mod_name<-rep(names(x)[i],nrow(x_df)) |
|
| 95 |
#x_df$term_name <-row.names(x_df) |
|
| 96 |
#return(x_df) |
|
| 97 |
#} |
|
| 67 | 98 |
#Function to rasterize a table with coordinates and variables...,maybe add option for ref image?? |
| 68 | 99 |
rasterize_df_fun <- function(data_tb,coord_names,proj_str,out_suffix,out_dir=".",file_format=".rst",NA_flag_val=-9999,tolerance_val= 0.000120005){
|
| 69 | 100 |
data_spdf <- data_tb |
| ... | ... | |
| 228 | 259 |
reg_name <- list_param_plot_daily_mosaics$reg_name |
| 229 | 260 |
out_dir_str <- list_param_plot_daily_mosaics$out_dir_str |
| 230 | 261 |
out_suffix <- list_param_plot_daily_mosaics$out_suffix |
| 262 |
l_dates <- list_param_plot_daily_mosaics$l_dates |
|
| 231 | 263 |
|
| 232 | 264 |
|
| 233 | 265 |
#list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str) |
| ... | ... | |
| 242 | 274 |
rclmat <- matrix(m, ncol=3, byrow=TRUE) |
| 243 | 275 |
rc <- reclassify(r_test, rclmat,filename=paste("rc_tmp_",i,".tif",sep=""),dataType="FLT4S",overwrite=T)
|
| 244 | 276 |
file_name <- unlist(strsplit(basename(lf_m[i]),"_")) |
| 245 |
date_proc <- file_name[7] #specific tot he current naming of files |
|
| 277 |
|
|
| 278 |
#date_proc <- file_name[7] #specific tot he current naming of files |
|
| 279 |
date_proc <- l_dates[i] |
|
| 246 | 280 |
#paste(raster_name[1:7],collapse="_") |
| 247 | 281 |
#add filename option later |
| 248 | 282 |
extension_str <- extension(filename(r_test)) |
| 249 | 283 |
raster_name_tmp <- gsub(extension_str,"",basename(filename(r_test))) |
| 250 |
raster_name <- file.path(out_dir_str,paste(raster_name_tmp,"_",date_proc,"_masked.tif",sep=""))
|
|
| 284 |
raster_name <- file.path(out_dir_str,paste(raster_name_tmp,"_masked.tif",sep="")) |
|
| 251 | 285 |
r_pred <- mask(r_test,rc,filename=raster_name,overwrite=TRUE) |
| 252 | 286 |
|
| 253 | 287 |
res_pix <- 1200 |
| ... | ... | |
| 259 | 293 |
png(filename=paste("Figure9_clim_mosaics_day_test","_",date_proc,"_",reg_name,"_",out_suffix,".png",sep=""),
|
| 260 | 294 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
| 261 | 295 |
|
| 262 |
plot(r_pred,main=paste("climatology month ",date_proc , " ", reg_name,sep=""),cex.main=1.5)
|
|
| 296 |
plot(r_pred,main=paste("Predicted on ",date_proc , " ", reg_name,sep=""),cex.main=1.5)
|
|
| 263 | 297 |
dev.off() |
| 264 | 298 |
|
| 265 | 299 |
return(raster_name) |
| ... | ... | |
| 275 | 309 |
#in_dir1 <- "/data/project/layers/commons/NEX_data/output_run6_global_analyses_09162014/output20Deg2" |
| 276 | 310 |
# parent output dir for the curent script analyes |
| 277 | 311 |
#out_dir <-"/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/" #On NCEAS Atlas |
| 278 |
out_dir <-"/data/project/layers/commons/NEX_data/output_run10_global_analyses_12232014/" |
|
| 279 | 312 |
# input dir containing shapefiles defining tiles |
| 280 | 313 |
#in_dir_shp <- "/data/project/layers/commons/NEX_data/output_run5_global_analyses_08252014/output/subset/shapefiles" |
| 281 | 314 |
|
| ... | ... | |
| 288 | 321 |
|
| 289 | 322 |
y_var_name <- "dailyTmax" |
| 290 | 323 |
interpolation_method <- c("gam_CAI")
|
| 291 |
out_prefix<-"run10_global_analyses_12232014"
|
|
| 324 |
out_prefix<-"run10_global_analyses_01202015"
|
|
| 292 | 325 |
mosaic_plot <- FALSE |
| 293 | 326 |
|
| 327 |
day_to_mosaic <- c("20100101","20100102","20100103","20100104","20100105",
|
|
| 328 |
"20100301","20100302","20100303","20100304","20100305", |
|
| 329 |
"20100501","20100502","20100503","20100504","20100505", |
|
| 330 |
"20100701","20100702","20100703","20100704","20100705", |
|
| 331 |
"20100901","20100902","20100903","20100904","20100905", |
|
| 332 |
"20101101","20101102","20101103","20101104","20101105") |
|
| 333 |
|
|
| 294 | 334 |
#CRS_locs_WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| 295 | 335 |
CRS_WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| 296 | 336 |
|
| ... | ... | |
| 302 | 342 |
|
| 303 | 343 |
#out_dir <-paste(out_dir,"_",out_prefix,sep="") |
| 304 | 344 |
create_out_dir_param <- FALSE |
| 345 |
out_dir <-"/data/project/layers/commons/NEX_data/output_run10_global_analyses_01202015/" |
|
| 305 | 346 |
|
| 306 | 347 |
if(create_out_dir_param==TRUE){
|
| 307 | 348 |
out_dir <- create_dir_fun(out_dir,out_prefix) |
| ... | ... | |
| 316 | 357 |
CRS_WGS84<-c("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| 317 | 358 |
|
| 318 | 359 |
region_name <- "world" |
| 319 |
# |
|
| 360 |
|
|
| 320 | 361 |
###Table 1: Average accuracy metrics |
| 321 | 362 |
###Table 2: daily accuracy metrics for all tiles |
| 322 |
#lf_tables <- list.files(out_dir,) |
|
| 363 |
|
|
| 323 | 364 |
summary_metrics_v <- read.table(file=file.path(out_dir,paste("summary_metrics_v2_NA_",out_prefix,".txt",sep="")),sep=",")
|
| 324 | 365 |
tb <- read.table(file=file.path(out_dir,paste("tb_diagnostic_v_NA","_",out_prefix,".txt",sep="")),sep=",")
|
| 325 | 366 |
#tb_diagnostic_s_NA_run10_global_analyses_11302014.txt |
| 326 | 367 |
tb_s <- read.table(file=file.path(out_dir,paste("tb_diagnostic_s_NA","_",out_prefix,".txt",sep="")),sep=",")
|
| 327 | 368 |
|
| 328 | 369 |
tb_month_s <- read.table(file=file.path(out_dir,paste("tb_month_diagnostic_s_NA","_",out_prefix,".txt",sep="")),sep=",")
|
| 329 |
#tb_month_s <- read.table("tb_month_diagnostic_s_NA_run5_global_analyses_08252014.txt",sep=",")
|
|
| 330 | 370 |
pred_data_month_info <- read.table(file=file.path(out_dir,paste("pred_data_month_info_",out_prefix,".txt",sep="")),sep=",")
|
| 331 | 371 |
pred_data_day_info <- read.table(file=file.path(out_dir,paste("pred_data_day_info_",out_prefix,".txt",sep="")),sep=",")
|
| 332 | 372 |
df_tile_processed <- read.table(file=file.path(out_dir,paste("df_tile_processed_",out_prefix,".txt",sep="")),sep=",")
|
| 333 |
#in_dir_list <- file.path(in_dir1,read.table(file.path(in_dir1,"processed.txt"))$V1) |
|
| 334 |
#gam_diagnostic_df <- read.table(file=file.path(out_dir,"gam_diagnostic_df_run4_global_analyses_08142014.txt"),sep=",") |
|
| 335 | 373 |
|
| 336 | 374 |
########################## START SCRIPT ############################## |
| 375 |
|
|
| 337 | 376 |
tb$pred_mod <- as.character(tb$pred_mod) |
| 338 | 377 |
summary_metrics_v$pred_mod <- as.character(summary_metrics_v$pred_mod) |
| 378 |
summary_metrics_v$tile_id <- as.character(summary_metrics_v$tile_id) |
|
| 379 |
df_tile_processed$tile_id <- as.character(df_tile_processed$tile_id) |
|
| 339 | 380 |
|
| 340 | 381 |
mulitple_region <- TRUE |
| 341 | 382 |
|
| 342 | 383 |
#multiple regions? |
| 343 | 384 |
if(mulitple_region==TRUE){
|
| 344 | 385 |
df_tile_processed$reg <- basename(dirname(as.character(df_tile_processed$path_NEX))) |
| 386 |
tb <- merge(tb,df_tile_processed,by="tile_id") |
|
| 387 |
|
|
| 345 | 388 |
} |
| 346 | 389 |
|
| 347 | 390 |
|
| ... | ... | |
| 384 | 427 |
#path_to_shp <- dirname(list_shp_reg_files[[i]]) |
| 385 | 428 |
path_to_shp <- file.path(out_dir,"/shapefiles") |
| 386 | 429 |
layer_name <- sub(".shp","",basename(list_shp_reg_files[[i]]))
|
| 387 |
shp1 <- readOGR(path_to_shp, layer_name) |
|
| 430 |
shp1 <- try(readOGR(path_to_shp, layer_name)) #use try to resolve error below |
|
| 431 |
#shp_61.0_-160.0 |
|
| 432 |
#Geographical CRS given to non-conformant data: -186.331747678 |
|
| 433 |
|
|
| 388 | 434 |
#shp1<-readOGR(dirname(list_shp_reg_files[[i]]),sub(".shp","",basename(list_shp_reg_files[[i]])))
|
| 389 |
pt <- gCentroid(shp1) |
|
| 390 |
centroids_pts[[i]] <-pt |
|
| 435 |
if (!inherits(shp1,"try-error")) {
|
|
| 436 |
pt <- gCentroid(shp1) |
|
| 437 |
centroids_pts[[i]] <- pt |
|
| 438 |
}else{
|
|
| 439 |
centroids <- shp1 |
|
| 440 |
} |
|
| 391 | 441 |
shps_tiles[[i]] <- shp1 |
| 392 | 442 |
} |
| 443 |
|
|
| 393 | 444 |
#coord_names <- c("lon","lat")
|
| 394 | 445 |
#l_rast <- rasterize_df_fun(test,coord_names,proj_str,out_suffix=out_prefix,out_dir=".",file_format,NA_flag_val,tolerance_val=0.000120005) |
| 395 | 446 |
|
| 447 |
#remove try-error polygons...we loose three tiles because they extend beyond -180 deg |
|
| 448 |
tmp <- shps_tiles |
|
| 449 |
shps_tiles <- remove_errors_list(shps_tiles) #[[!inherits(shps_tiles,"try-error")]] |
|
| 450 |
#shps_tiles <- tmp |
|
| 451 |
|
|
| 452 |
tmp_pts <- centroids_pts |
|
| 453 |
centroids_pts <- remove_errors_list(centroids_pts) #[[!inherits(shps_tiles,"try-error")]] |
|
| 454 |
#centroids_pts <- tmp_pts |
|
| 455 |
|
|
| 396 | 456 |
#plot info: with labels |
| 397 | 457 |
res_pix <- 1200 |
| 398 |
col_mfrow <- 1 |
|
| 458 |
col_mfrow <- 1
|
|
| 399 | 459 |
row_mfrow <- 1 |
| 400 | 460 |
|
| 401 | 461 |
png(filename=paste("Figure1_tile_processed_region_",region_name,"_",out_prefix,".png",sep=""),
|
| ... | ... | |
| 403 | 463 |
|
| 404 | 464 |
plot(reg_layer) |
| 405 | 465 |
#Add polygon tiles... |
| 406 |
for(i in 1:length(list_shp_reg_files)){
|
|
| 466 |
for(i in 1:length(shps_tiles)){
|
|
| 407 | 467 |
shp1 <- shps_tiles[[i]] |
| 408 | 468 |
pt <- centroids_pts[[i]] |
| 409 |
plot(shp1,add=T,border="blue") |
|
| 410 |
#plot(pt,add=T,cex=2,pch=5) |
|
| 411 |
label_id <- df_tile_processed$tile_id[i] |
|
| 412 |
text(coordinates(pt)[1],coordinates(pt)[2],labels=i,cex=1.3,font=2,col=c("red"))
|
|
| 469 |
if(!inherits(shp1,"try-error")){
|
|
| 470 |
plot(shp1,add=T,border="blue") |
|
| 471 |
#plot(pt,add=T,cex=2,pch=5) |
|
| 472 |
label_id <- df_tile_processed$tile_id[i] |
|
| 473 |
text(coordinates(pt)[1],coordinates(pt)[2],labels=i,cex=1.3,font=2,col=c("red"))
|
|
| 474 |
} |
|
| 413 | 475 |
} |
| 414 |
title(paste("Tiles location 20x20 degrees for ", region_name,sep=""))
|
|
| 476 |
title(paste("Tiles 1000x3000 ", region_name,sep=""))
|
|
| 415 | 477 |
|
| 416 | 478 |
dev.off() |
| 417 | 479 |
|
| ... | ... | |
| 494 | 556 |
#y_var_name <-"dailyTmax" |
| 495 | 557 |
#index <-244 #index corresponding to Sept 1 |
| 496 | 558 |
|
| 497 |
if (mosaic_plot==TRUE){
|
|
| 498 |
index <- 1 #index corresponding to Jan 1 |
|
| 499 |
date_selected <- "20100901" |
|
| 500 |
name_method_var <- paste(interpolation_method,"_",y_var_name,"_",sep="") |
|
| 501 |
|
|
| 502 |
pattern_str <- paste("mosaiced","_",name_method_var,"predicted",".*.",date_selected,".*.tif",sep="")
|
|
| 503 |
lf_pred_list <- list.files(pattern=pattern_str) |
|
| 504 |
|
|
| 505 |
for(i in 1:length(lf_pred_list)){
|
|
| 506 |
|
|
| 507 |
|
|
| 508 |
r_pred <- raster(lf_pred_list[i]) |
|
| 509 |
|
|
| 510 |
res_pix <- 480 |
|
| 511 |
col_mfrow <- 1 |
|
| 512 |
row_mfrow <- 1 |
|
| 513 |
|
|
| 514 |
png(filename=paste("Figure4_models_predicted_surfaces_",model_name[i],"_",name_method_var,"_",data_selected,"_",out_prefix,".png",sep=""),
|
|
| 515 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 516 |
|
|
| 517 |
plot(r_pred) |
|
| 518 |
title(paste("Mosaiced",model_name[i],name_method_var,date_selected,sep=" "))
|
|
| 519 |
dev.off() |
|
| 520 |
} |
|
| 521 |
#Plot Delta and clim... |
|
| 522 |
|
|
| 523 |
## plotting of delta and clim for later scripts... |
|
| 524 |
|
|
| 525 |
} |
|
| 526 |
|
|
| 527 |
####### Figure 5... |
|
| 528 |
### Adding tiles do a plot of mod1 with tiles |
|
| 529 |
|
|
| 530 |
#r_pred <- raster(lf_pred_list[i]) |
|
| 531 |
|
|
| 532 |
#res_pix <- 480 |
|
| 533 |
#col_mfrow <- 1 |
|
| 534 |
#row_mfrow <- 1 |
|
| 535 |
#model_name <- as.character(model_name) |
|
| 536 |
|
|
| 537 |
#i<-2 |
|
| 538 |
#png(filename=paste("Figure5_tiles_with_models_predicted_surfaces_",model_name[i],"_",name_method_var,out_prefix,".png",sep=""),
|
|
| 539 |
# width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 540 |
# |
|
| 541 |
#plot(r_pred) |
|
| 542 |
#plot(reg_layer) |
|
| 543 |
|
|
| 544 |
#title(paste("Mosaiced",model_name[i],name_method_var,date_selected,"with tiles",sep=" "))
|
|
| 545 |
|
|
| 546 |
#Add polygon tiles... |
|
| 547 |
#for(i in 1:length(list_shp_reg_files)){
|
|
| 548 |
# shp1 <- shps_tiles[[i]] |
|
| 549 |
# pt <- centroids_pts[[i]] |
|
| 550 |
# plot(shp1,add=T,border="blue") |
|
| 551 |
#plot(pt,add=T,cex=2,pch=5) |
|
| 552 |
# label_id <- df_tile_processed$tile_id[i] |
|
| 553 |
# text(coordinates(pt)[1],coordinates(pt)[2],labels=i,cex=1,col=c("red"))
|
|
| 554 |
#} |
|
| 555 |
#title(paste("Prediction with tiles location 10x10 degrees for ", region_name,sep=""))
|
|
| 556 |
#dev.off() |
|
| 559 |
# if (mosaic_plot==TRUE){
|
|
| 560 |
# index <- 1 #index corresponding to Jan 1 |
|
| 561 |
# date_selected <- "20100901" |
|
| 562 |
# name_method_var <- paste(interpolation_method,"_",y_var_name,"_",sep="") |
|
| 563 |
# |
|
| 564 |
# pattern_str <- paste("mosaiced","_",name_method_var,"predicted",".*.",date_selected,".*.tif",sep="")
|
|
| 565 |
# lf_pred_list <- list.files(pattern=pattern_str) |
|
| 566 |
# |
|
| 567 |
# for(i in 1:length(lf_pred_list)){
|
|
| 568 |
# |
|
| 569 |
# |
|
| 570 |
# r_pred <- raster(lf_pred_list[i]) |
|
| 571 |
# |
|
| 572 |
# res_pix <- 480 |
|
| 573 |
# col_mfrow <- 1 |
|
| 574 |
# row_mfrow <- 1 |
|
| 575 |
# |
|
| 576 |
# png(filename=paste("Figure4_models_predicted_surfaces_",model_name[i],"_",name_method_var,"_",data_selected,"_",out_prefix,".png",sep=""),
|
|
| 577 |
# width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 578 |
# |
|
| 579 |
# plot(r_pred) |
|
| 580 |
# title(paste("Mosaiced",model_name[i],name_method_var,date_selected,sep=" "))
|
|
| 581 |
# dev.off() |
|
| 582 |
# } |
|
| 583 |
# #Plot Delta and clim... |
|
| 584 |
# |
|
| 585 |
# ## plotting of delta and clim for later scripts... |
|
| 586 |
# |
|
| 587 |
# } |
|
| 557 | 588 |
|
| 558 |
### |
|
| 559 |
|
|
| 560 |
#### Now combined plot... |
|
| 561 |
#Use the function to match extent... |
|
| 562 |
|
|
| 563 |
#pred_s <- stack(lf_list) #problem different extent!! |
|
| 564 |
#methods_names <-c("gam","kriging","gwr")
|
|
| 565 |
#methods_names <- interpolation_method |
|
| 566 |
|
|
| 567 |
#names_layers<-methods_names[1] |
|
| 568 |
#names_layers <-c("mod1 = lat*long + elev","mod2 = lat*long + elev + LST",
|
|
| 569 |
# "mod3 = lat*long + elev + LST*FOREST")#, "mod_kr") |
|
| 570 |
#nb_fig<- c("4a","4b")
|
|
| 571 |
#list_plots_spt <- vector("list",length=length(lf))
|
|
| 572 |
#png(filename=paste("Figure4_models_predicted_surfaces_",date_selected,"_",out_prefix,".png",sep=""),
|
|
| 573 |
# height=480*layout_m[1],width=480*layout_m[2]) |
|
| 574 |
# max_val <- 40 |
|
| 575 |
# min_val <- -40 |
|
| 576 |
# layout_m <- c(1,3) #one row two columns |
|
| 577 |
# no_brks <- length(seq(min_val,max_val,by=0.1))-1 |
|
| 578 |
# #temp.colors <- colorRampPalette(c('blue', 'white', 'red'))
|
|
| 579 |
# #temp.colors <- colorRampPalette(c('blue', 'lightgoldenrodyellow', 'red'))
|
|
| 580 |
# #temp.colors <- matlab.like(no_brks) |
|
| 581 |
# temp.colors <- colorRampPalette(c('blue', 'khaki', 'red'))
|
|
| 582 |
# |
|
| 583 |
# p <- levelplot(pred_s,main=methods_names[i], |
|
| 584 |
# ylab=NULL,xlab=NULL, |
|
| 585 |
# par.settings = list(axis.text = list(font = 2, cex = 2),layout=layout_m, |
|
| 586 |
# par.main.text=list(font=2,cex=2.5),strip.background=list(col="white")),par.strip.text=list(font=2,cex=2), |
|
| 587 |
# names.attr=names_layers, |
|
| 588 |
# col.regions=temp.colors(no_brks),at=seq(min_val,max_val,by=0.1)) |
|
| 589 |
##col.regions=temp.colors(25)) |
|
| 590 |
#print(p) |
|
| 591 |
#dev.off() |
|
| 592 | 589 |
|
| 593 | 590 |
###################### |
| 594 | 591 |
### Figure 5: plot accuracy ranked |
| ... | ... | |
| 700 | 697 |
} |
| 701 | 698 |
|
| 702 | 699 |
## Number of tiles with information: |
| 703 |
|
|
| 704 |
sum(df_tile_processed$metrics_v)/length(df_tile_processed$metrics_v) #81.40% |
|
| 700 |
sum(df_tile_processed$metrics_v) |
|
| 701 |
length(df_tile_processed$metrics_v) |
|
| 702 |
sum(df_tile_processed$metrics_v)/length(df_tile_processed$metrics_v) #70.69% |
|
| 705 | 703 |
|
| 706 | 704 |
#coordinates |
| 707 | 705 |
coordinates(summary_metrics_v) <- c("lon","lat")
|
| ... | ... | |
| 807 | 805 |
########################################################## |
| 808 | 806 |
##### Figure 8: Breaking down accuaracy by regions!! ##### |
| 809 | 807 |
|
| 808 |
summary_metrics_v <- merge(summary_metrics_v,df_tile_processed,by="tile_id") |
|
| 809 |
table(summary_metrics_v$reg) |
|
| 810 | 810 |
|
| 811 |
#This part is specifically related to this run : dropping model with elev |
|
| 812 |
#tb <- merge(tb,df_tile_processed,by="tile_id") |
|
| 811 |
## Figure 8a |
|
| 812 |
res_pix <- 480 |
|
| 813 |
col_mfrow <- 1 |
|
| 814 |
row_mfrow <- 1 |
|
| 813 | 815 |
|
| 814 |
#tb_tmp_reg5 <- subset(tb,reg=="reg5") |
|
| 815 |
#tb_tmp_reg4 <- subset(tb,reg=="reg4") |
|
| 816 |
#tb_tmp_reg4 <- subset(tb_tmp_reg4,pred_mod!="mod1") #remove mod1 because it is not in reg5 |
|
| 817 |
#tb_tmp_reg4$pred_mod <- replace(tb_tmp_reg4$pred_mod, tb_tmp_reg4$pred_mod=="mod2", "mod1") |
|
| 818 |
#tb <- rbind(tb_tmp_reg4,tb_tmp_reg5) |
|
| 816 |
png(filename=paste("Figure8a_boxplot_overall_separated_by_region_with_oultiers_",model_name[i],"_",out_prefix,".png",sep=""),
|
|
| 817 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 819 | 818 |
|
| 820 |
#ac_mod <- summary_metrics_v[summary_metrics_v$pred_mod==model_name[i],] |
|
| 821 |
#summary_metrics_v <- merge(summary_metrics_v,df_tile_processed,by="tile_id") |
|
| 822 |
#table(summary_metrics_v$reg) |
|
| 819 |
p<- bwplot(rmse~pred_mod | reg, data=tb, |
|
| 820 |
main="RMSE per model and region over all tiles") |
|
| 821 |
print(p) |
|
| 822 |
dev.off() |
|
| 823 | 823 |
|
| 824 |
#summary_metrics_v_reg5 <- subset(summary_metrics_v,reg=="reg5") |
|
| 825 |
#summary_metrics_v_reg4 <- subset(summary_metrics_v,reg=="reg4") |
|
| 826 |
#summary_metrics_v_reg4 <- subset(summary_metrics_v_reg4,pred_mod!="mod1") #remove mod1 because it is not in reg5 |
|
| 827 |
#summary_metrics_v_reg4$pred_mod <- replace(summary_metrics_v_reg4$pred_mod, |
|
| 828 |
# summary_metrics_v_reg4$pred_mod=="mod2", "mod1") |
|
| 829 |
#summary_metrics_v <- rbind(summary_metrics_v_reg4,summary_metrics_v_reg5) |
|
| 824 |
## Figure 8b |
|
| 825 |
png(filename=paste("Figure8b_boxplot_overall_separated_by_region_scaling_",model_name[i],"_",out_prefix,".png",sep=""),
|
|
| 826 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 830 | 827 |
|
| 828 |
boxplot(rmse~pred_mod,data=tb,ylim=c(0,5),outline=FALSE)#,names=tb$pred_mod) |
|
| 829 |
title("RMSE per model over all tiles")
|
|
| 830 |
p<- bwplot(rmse~pred_mod | reg, data=tb,ylim=c(0,5), |
|
| 831 |
main="RMSE per model and region over all tiles") |
|
| 832 |
print(p) |
|
| 833 |
dev.off() |
|
| 831 | 834 |
|
| 832 | 835 |
##################################################### |
| 833 | 836 |
#### Figure 9: plotting mosaics of regions ########### |
| ... | ... | |
| 844 | 847 |
full.names=T)) |
| 845 | 848 |
} |
| 846 | 849 |
|
| 847 |
#r_reg5 <- stack(lf_mosaics_reg5) |
|
| 848 |
lf_m <- lf_mosaics_reg5[1:12] |
|
| 849 |
reg_name <- "reg5" |
|
| 850 |
for(i in 1:length(lf_m)){
|
|
| 851 |
r_test<- raster(lf_m[i]) |
|
| 852 |
#r_test <- subset(r_reg5,1) |
|
| 853 |
|
|
| 854 |
r_test_mask_high <- r_test < 100 |
|
| 855 |
NAvalue(r_test_mask_high) <- 0 |
|
| 856 |
r_test_mask_low <- r_test > -100 |
|
| 857 |
NAvalue(r_test_mask_low) <- 0 |
|
| 858 |
r3 <- overlay(r_test, r_test_mask_high, r_test_mask_low, fun=function(x,y,z){return(x*y*z)})
|
|
| 859 |
#writeRaster(r3,file=paste("CAI_TMAX_clim_month_",i,"_mod1_all_",reg_name,"_masked.tif",sep=""),overwrite=TRUE)
|
|
| 860 |
|
|
| 861 |
res_pix <- 1200 |
|
| 862 |
#res_pix <- 480 |
|
| 863 |
|
|
| 864 |
col_mfrow <- 1 |
|
| 865 |
row_mfrow <- 1 |
|
| 866 |
|
|
| 867 |
png(filename=paste("Figure9_clim_mosaics_month","_",i,"_",reg_name,"_",out_prefix,".png",sep=""),
|
|
| 868 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 869 |
|
|
| 870 |
plot(r3,main=paste("climatology month ", i, " ", reg_name,sep=""),cex.main=1.5)
|
|
| 871 |
dev.off() |
|
| 872 |
} |
|
| 873 |
|
|
| 874 |
#### |
|
| 875 |
|
|
| 876 |
#Monthly |
|
| 877 |
lf_m <- lf_mosaics_reg2[13:15] |
|
| 878 |
lf_m <- lf_mosaics_reg2[1:12] |
|
| 879 |
|
|
| 880 |
reg_name <- "reg2" |
|
| 881 |
for(i in 1:length(lf_m)){
|
|
| 882 |
r_test<- raster(lf_m[i]) |
|
| 883 |
#r_test <- subset(r_reg5,1) |
|
| 884 |
|
|
| 885 |
#r_test_mask_high <- r_test < 100 |
|
| 886 |
#r_test_mask_high <- r_test[r_test < 100] |
|
| 887 |
|
|
| 888 |
#r_test_mask_high <- r_test < 100 |
|
| 889 |
|
|
| 890 |
m <- c(-Inf, -100, NA, |
|
| 891 |
-100, 100, 1, |
|
| 892 |
100, Inf, NA) |
|
| 893 |
rclmat <- matrix(m, ncol=3, byrow=TRUE) |
|
| 894 |
rc <- reclassify(r_test, rclmat,filename="rc.tif",dataType="FLT4S",overwrite=T) |
|
| 895 |
r_pred <- mask(r_test,rc,filename=paste("CAI_TMAX_clim_month_",i,"_mod1_all_",reg_name,"_masked.tif",sep=""),overwrite=TRUE)
|
|
| 896 |
#r <- raster(r_test_mask_high,dataType="INT2U") |
|
| 897 |
#r3 <- clamp(r_test,-100,100) |
|
| 898 |
#NAvalue(r_test_mask_high) <- 0 |
|
| 899 |
|
|
| 900 |
#r_test_mask_low <- r_test > -100 |
|
| 901 |
#NAvalue(r_test_mask_low) <- 0 |
|
| 902 |
#r3 <- overlay(r_test, r_test_mask_high, r_test_mask_low, fun=function(x,y,z){return(x*y*z)})
|
|
| 903 |
#writeRaster(r3,file=paste("CAI_TMAX_clim_month_",i,"_mod1_all_",reg_name,"_masked.tif",sep=""),overwrite=TRUE)
|
|
| 904 |
|
|
| 905 |
res_pix <- 1200 |
|
| 906 |
#res_pix <- 480 |
|
| 907 |
|
|
| 908 |
col_mfrow <- 1 |
|
| 909 |
row_mfrow <- 1 |
|
| 910 |
|
|
| 911 |
png(filename=paste("Figure9_clim_mosaics_day_test","_",i,"_",reg_name,"_",out_prefix,".png",sep=""),
|
|
| 912 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 913 |
|
|
| 914 |
plot(r_pred,main=paste("climatology month ", i, " ", reg_name,sep=""),cex.main=1.5)
|
|
| 915 |
dev.off() |
|
| 916 |
} |
|
| 917 |
|
|
| 918 |
#daily |
|
| 919 |
lf_m <- lf_mosaics_reg2[13:length(lf_mosaics_reg2)] |
|
| 920 |
#lf_m <- lf_mosaics_reg2[1:12] |
|
| 921 |
|
|
| 922 |
reg_name <- "reg2" |
|
| 923 |
for(i in 1:length(lf_m)){
|
|
| 924 |
r_test<- raster(lf_m[i]) |
|
| 925 |
|
|
| 926 |
m <- c(-Inf, -100, NA, |
|
| 927 |
-100, 100, 1, |
|
| 928 |
100, Inf, NA) |
|
| 929 |
rclmat <- matrix(m, ncol=3, byrow=TRUE) |
|
| 930 |
rc <- reclassify(r_test, rclmat,filename="rc.tif",dataType="FLT4S",overwrite=T) |
|
| 931 |
raster_name <- unlist(strsplit(basename(lf_m[i]),"_")) |
|
| 932 |
date_proc <- raster_name[5] |
|
| 933 |
#paste(raster_name[1:7],collapse="_") |
|
| 934 |
r_pred <- mask(r_test,rc,filename=paste("CAI_TMAX_clim_month_mod1_all_",reg_name,"_",date_proc,"_masked.tif",sep=""),overwrite=TRUE)
|
|
| 935 |
|
|
| 936 |
res_pix <- 1200 |
|
| 937 |
#res_pix <- 480 |
|
| 938 |
|
|
| 939 |
col_mfrow <- 1 |
|
| 940 |
row_mfrow <- 1 |
|
| 941 |
|
|
| 942 |
png(filename=paste("Figure9_clim_mosaics_day_test","_",date_proc,"_",reg_name,"_",out_prefix,".png",sep=""),
|
|
| 943 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 944 |
|
|
| 945 |
plot(r_pred,main=paste("climatology month ",date_proc , " ", reg_name,sep=""),cex.main=1.5)
|
|
| 946 |
dev.off() |
|
| 947 |
} |
|
| 850 |
#plot Australia |
|
| 851 |
lf_m <- lf_mosaics_reg[[2]] |
|
| 852 |
out_dir_str <- out_dir |
|
| 853 |
reg_name <- "reg6_1000x3000" |
|
| 854 |
#lapply() |
|
| 855 |
#list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str,out_suffix=out_prefix) |
|
| 856 |
list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str,out_suffix=out_prefix,l_dates=day_to_mosaic) |
|
| 857 |
|
|
| 858 |
#lf_m_mask_reg4_1500x4500 <- mclapply(1:2,FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 6) |
|
| 859 |
#debug(plot_daily_mosaics) |
|
| 860 |
#lf_m_mask_reg6_1000x3000 <- plot_daily_mosaics(1,list_param=list_param_plot_daily_mosaics) |
|
| 861 |
|
|
| 862 |
lf_m_mask_reg6_1000x3000 <- mclapply(1:length(lf_m),FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 10) |
|
| 863 |
|
|
| 864 |
#### North America |
|
| 865 |
lf_m <- lf_mosaics_reg[[1]] |
|
| 866 |
out_dir_str <- out_dir |
|
| 867 |
reg_name <- "reg1_1000x3000" |
|
| 868 |
#lapply() |
|
| 869 |
list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str,out_suffix=out_prefix,l_dates=day_to_mosaic) |
|
| 870 |
#lf_m_mask_reg4_1500x4500 <- mclapply(1:2,FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 6) |
|
| 871 |
|
|
| 872 |
lf_m_mask_reg1_1000x3000 <- mclapply(1:length(lf_m),FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 10) |
|
| 948 | 873 |
|
| 949 | 874 |
##################### END OF SCRIPT ###################### |
Also available in: Unified diff
NEX assessment script part 2, major clean up