Revision d31ebea7
Added by Benoit Parmentier about 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: 08/28/2014
|
|
| 8 |
#MODIFIED ON: 09/16/2014
|
|
| 9 | 9 |
#Version: 3 |
| 10 | 10 |
#PROJECT: Environmental Layers project |
| 11 | 11 |
#COMMENTS: analyses for run 5 global using 6 specific tiles |
| ... | ... | |
| 69 | 69 |
#on ATLAS |
| 70 | 70 |
#in_dir1 <- "/data/project/layers/commons/NEX_data/test_run1_03232014/output" #On Atlas |
| 71 | 71 |
#parent output dir : contains subset of the data produced on NEX |
| 72 |
in_dir1 <- "/data/project/layers/commons/NEX_data/output_run5_global_analyses_08252014/output20Deg"
|
|
| 72 |
in_dir1 <- "/data/project/layers/commons/NEX_data/output_run6_global_analyses_09162014/output20Deg2"
|
|
| 73 | 73 |
# parent output dir for the curent script analyes |
| 74 | 74 |
#out_dir <-"/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/" #On NCEAS Atlas |
| 75 |
out_dir <-"/data/project/layers/commons/NEX_data/output_run5_global_analyses_08252014/"
|
|
| 75 |
out_dir <-"/data/project/layers/commons/NEX_data/output_run6_global_analyses_09162014/"
|
|
| 76 | 76 |
# input dir containing shapefiles defining tiles |
| 77 | 77 |
#in_dir_shp <- "/data/project/layers/commons/NEX_data/output_run5_global_analyses_08252014/output/subset/shapefiles" |
| 78 | 78 |
|
| ... | ... | |
| 83 | 83 |
#out_dir <- "/nobackup/bparmen1/" #on NEX |
| 84 | 84 |
#in_dir_shp <- "/nobackupp4/aguzman4/climateLayers/output4/subset/shapefiles/" |
| 85 | 85 |
|
| 86 |
y_var_nay_var_name <- "dailyTmax"
|
|
| 86 |
y_var_name <- "dailyTmax" |
|
| 87 | 87 |
interpolation_method <- c("gam_CAI")
|
| 88 |
out_prefix<-"run5_global_analyses_08252014"
|
|
| 88 |
out_prefix<-"run6_global_analyses_09162014"
|
|
| 89 | 89 |
|
| 90 | 90 |
|
| 91 | 91 |
#out_dir <-paste(out_dir,"_",out_prefix,sep="") |
| ... | ... | |
| 113 | 113 |
tb_s <- read.table(file=file.path(out_dir,paste("tb_diagnostic_s_NA","_",out_prefix,".txt",sep="")),sep=",")
|
| 114 | 114 |
|
| 115 | 115 |
tb_month_s_<- read.table(file=file.path(out_dir,paste("tb_month_diagnostic_s_NA","_",out_prefix,".txt",sep="")),sep=",")
|
| 116 |
tb_month_s <- read.table("tb_month_diagnostic_s_NA_run5_global_analyses_08252014.txt",sep=",")
|
|
| 116 |
#tb_month_s <- read.table("tb_month_diagnostic_s_NA_run5_global_analyses_08252014.txt",sep=",")
|
|
| 117 | 117 |
pred_data_month_info <- read.table(file=file.path(out_dir,paste("pred_data_month_info_",out_prefix,".txt",sep="")),sep=",")
|
| 118 | 118 |
pred_data_day_info <- read.table(file=file.path(out_dir,paste("pred_data_day_info_",out_prefix,".txt",sep="")),sep=",")
|
| 119 | 119 |
df_tile_processed <- read.table(file=file.path(out_dir,paste("df_tile_processed_",out_prefix,".txt",sep="")),sep=",")
|
| 120 | 120 |
#in_dir_list <- file.path(in_dir1,read.table(file.path(in_dir1,"processed.txt"))$V1) |
| 121 |
gam_diagnostic_df <- read.table(file=file.path(out_dir,"gam_diagnostic_df_run4_global_analyses_08142014.txt"),sep=",") |
|
| 121 |
#gam_diagnostic_df <- read.table(file=file.path(out_dir,"gam_diagnostic_df_run4_global_analyses_08142014.txt"),sep=",")
|
|
| 122 | 122 |
|
| 123 | 123 |
########################## START SCRIPT ############################## |
| 124 | 124 |
|
| ... | ... | |
| 128 | 128 |
#df_tiled_processed <- na.omit(df_tile_processed) #remove other list of folders irrelevant |
| 129 | 129 |
#list_shp_reg_files <- df_tiled_processed$shp_files |
| 130 | 130 |
list_shp_reg_files<- as.character(df_tile_processed$shp_files) |
| 131 |
list_shp_reg_files <- file.path("/data/project/layers/commons/NEX_data/output_run4_global_analyses_08142014/output20Deg",
|
|
| 132 |
as.character(df_tile_processed$tile_coord),"shapefiles",basename(list_shp_reg_files)) |
|
| 131 |
#list_shp_reg_files <- file.path("/data/project/layers/commons/NEX_data/",out_dir,
|
|
| 132 |
# as.character(df_tile_processed$tile_coord),"shapefiles",basename(list_shp_reg_files)) |
|
| 133 |
list_shp_reg_files <- file.path("/data/project/layers/commons/NEX_data/",out_dir,
|
|
| 134 |
"shapefiles",basename(list_shp_reg_files)) |
|
| 133 | 135 |
|
| 134 | 136 |
### First get background map to display where study area is located |
| 135 | 137 |
#can make this more general later on.. |
| ... | ... | |
| 154 | 156 |
shps_tiles <- vector("list",length(list_shp_reg_files))
|
| 155 | 157 |
#collect info: read in all shapfiles |
| 156 | 158 |
#This is slow...make a function and use mclapply?? |
| 157 |
|
|
| 159 |
/data/project/layers/commons/NEX_data/output_run6_global_analyses_09162014/shapefiles |
|
| 158 | 160 |
for(i in 1:length(list_shp_reg_files)){
|
| 159 |
path_to_shp <- dirname(list_shp_reg_files[[i]]) |
|
| 160 |
path_to_shp <- in_dir1
|
|
| 161 |
#path_to_shp <- dirname(list_shp_reg_files[[i]])
|
|
| 162 |
path_to_shp <- "/data/project/layers/commons/NEX_data/output_run6_global_analyses_09162014/shapefiles"
|
|
| 161 | 163 |
layer_name <- sub(".shp","",basename(list_shp_reg_files[[i]]))
|
| 162 | 164 |
shp1 <- readOGR(path_to_shp, layer_name) |
| 163 | 165 |
#shp1<-readOGR(dirname(list_shp_reg_files[[i]]),sub(".shp","",basename(list_shp_reg_files[[i]])))
|
| ... | ... | |
| 168 | 170 |
|
| 169 | 171 |
|
| 170 | 172 |
#plot info: with labels |
| 171 |
res_pix <- 480
|
|
| 173 |
res_pix <- 1200
|
|
| 172 | 174 |
col_mfrow <- 1 |
| 173 | 175 |
row_mfrow <- 1 |
| 174 | 176 |
|
| ... | ... | |
| 183 | 185 |
plot(shp1,add=T,border="blue") |
| 184 | 186 |
#plot(pt,add=T,cex=2,pch=5) |
| 185 | 187 |
label_id <- df_tile_processed$tile_id[i] |
| 186 |
text(coordinates(pt)[1],coordinates(pt)[2],labels=i,cex=1,col=c("red"))
|
|
| 188 |
text(coordinates(pt)[1],coordinates(pt)[2],labels=i,cex=1.3,font=2,col=c("red"))
|
|
| 187 | 189 |
} |
| 188 | 190 |
title(paste("Tiles location 20x20 degrees for ", region_name,sep=""))
|
| 189 | 191 |
|
| ... | ... | |
| 292 | 294 |
|
| 293 | 295 |
##### Diagnostic gam |
| 294 | 296 |
#### |
| 295 |
gam_diag_10x10 <- read.table("gam_diagnostic_10x10_df_run4_global_analyses_08142014.txt",sep=",")
|
|
| 296 |
gam_diag_10x10$month <- as.factor(gam_diag_10x10$month) |
|
| 297 |
#gam_diag_10x10 <- read.table("gam_diagnostic_10x10_df_run4_global_analyses_08142014.txt",sep=",")
|
|
| 298 |
#gam_diag_10x10$month <- as.factor(gam_diag_10x10$month)
|
|
| 297 | 299 |
|
| 298 |
xyplot(rmse~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 299 |
pred_mod!="mod_kr"),type="b") |
|
| 300 |
#xyplot(rmse~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v),
|
|
| 301 |
# pred_mod!="mod_kr"),type="b")
|
|
| 300 | 302 |
|
| 301 |
xyplot(n~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 302 |
pred_mod!="mod_kr"),type="h") |
|
| 303 |
#xyplot(n~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v),
|
|
| 304 |
# pred_mod!="mod_kr"),type="h")
|
|
| 303 | 305 |
|
| 304 |
xyplot(n~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 305 |
pred_mod!="mod_kr"),type="h") |
|
| 306 |
#xyplot(n~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v),
|
|
| 307 |
# pred_mod!="mod_kr"),type="h")
|
|
| 306 | 308 |
|
| 307 |
xyplot(n~month | tile_id + pred_mod,data=subset(as.data.frame(tb_month_s), |
|
| 308 |
pred_mod!="mod_kr"),type="h") |
|
| 309 |
#xyplot(n~month | tile_id + pred_mod,data=subset(as.data.frame(tb_month_s),
|
|
| 310 |
# pred_mod!="mod_kr"),type="h")
|
|
| 309 | 311 |
|
| 310 | 312 |
# |
| 311 | 313 |
## Figure 3b |
| ... | ... | |
| 349 | 351 |
####### Figure 5... |
| 350 | 352 |
### Adding tiles do a plot of mod1 with tiles |
| 351 | 353 |
|
| 352 |
r_pred <- raster(lf_pred_list[i]) |
|
| 354 |
#r_pred <- raster(lf_pred_list[i])
|
|
| 353 | 355 |
|
| 354 | 356 |
res_pix <- 480 |
| 355 | 357 |
col_mfrow <- 1 |
| 356 | 358 |
row_mfrow <- 1 |
| 357 |
|
|
| 359 |
model_name <- as.character(model_name) |
|
| 360 |
|
|
| 361 |
i<-2 |
|
| 358 | 362 |
png(filename=paste("Figure5_tiles_with_models_predicted_surfaces_",model_name[i],"_",name_method_var,out_prefix,".png",sep=""),
|
| 359 | 363 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
| 360 | 364 |
|
| 361 |
plot(r_pred) |
|
| 362 |
title(paste("Mosaiced",model_name[i],name_method_var,date_selected,"with tiles",sep=" "))
|
|
| 365 |
#plot(r_pred) |
|
| 366 |
plot(reg_layer) |
|
| 367 |
|
|
| 368 |
#title(paste("Mosaiced",model_name[i],name_method_var,date_selected,"with tiles",sep=" "))
|
|
| 363 | 369 |
|
| 364 | 370 |
#Add polygon tiles... |
| 365 | 371 |
for(i in 1:length(list_shp_reg_files)){
|
| ... | ... | |
| 415 | 421 |
|
| 416 | 422 |
coordinates(summary_metrics_v) <- cbind(summary_metrics_v$lon,summary_metrics_v$lat) |
| 417 | 423 |
proj4string(summary_metrics_v) <- CRS_WGS84 |
| 418 |
lf_list <- lf_pred_list |
|
| 424 |
#lf_list <- lf_pred_list
|
|
| 419 | 425 |
list_df_ac_mod <- vector("list",length=length(lf_pred_list))
|
| 420 | 426 |
for (i in 1:length(lf_list)){
|
| 421 | 427 |
|
| ... | ... | |
| 430 | 436 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
| 431 | 437 |
|
| 432 | 438 |
plot(r_pred) |
| 433 |
|
|
| 439 |
|
|
| 434 | 440 |
#plot(ac_mod1,cex=sqrt(ac_mod1$rmse),pch=1,add=T) |
| 435 | 441 |
plot(ac_mod,cex=(ac_mod$rmse^2)/10,pch=1,add=T) |
| 436 | 442 |
#plot(ac_mod1,cex=(ac_mod1$rmse1)*2,pch=1,add=T) |
| ... | ... | |
| 446 | 452 |
#quick kriging... |
| 447 | 453 |
#autokrige(rmse~1,r2,) |
| 448 | 454 |
|
| 455 |
|
|
| 456 |
### Without |
|
| 457 |
|
|
| 458 |
list_df_ac_mod <- vector("list",length=length(lf_pred_list))
|
|
| 459 |
list_df_ac_mod <- vector("list",length=3)
|
|
| 460 |
|
|
| 461 |
for (i in 1:length(model_name)){
|
|
| 462 |
|
|
| 463 |
ac_mod <- summary_metrics_v[summary_metrics_v$pred_mod==model_name[i],] |
|
| 464 |
#r_pred <- raster(lf_list[i]) |
|
| 465 |
|
|
| 466 |
res_pix <- 1200 |
|
| 467 |
#res_pix <- 480 |
|
| 468 |
|
|
| 469 |
col_mfrow <- 1 |
|
| 470 |
row_mfrow <- 1 |
|
| 471 |
|
|
| 472 |
png(filename=paste("Figure6_ac_metrics_map_centroids_tile_",model_name[i],"_",out_prefix,".png",sep=""),
|
|
| 473 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 474 |
|
|
| 475 |
#plot(r_pred) |
|
| 476 |
#plot(reg_layer) |
|
| 477 |
#plot(ac_mod1,cex=sqrt(ac_mod1$rmse),pch=1,add=T) |
|
| 478 |
#plot(ac_mod,cex=(ac_mod$rmse^2)/10,pch=1,col="red",add=T) |
|
| 479 |
|
|
| 480 |
p_shp <- layer(sp.polygons(reg_layer, lwd=1, col='black')) |
|
| 481 |
#title("(a) Mean for 1 January")
|
|
| 482 |
p <- bubble(ac_mod,"rmse",main=paste("Averrage RMSE per tile and by ",model_name[i]))
|
|
| 483 |
p1 <- p+p_shp |
|
| 484 |
print(p1) |
|
| 485 |
#plot(ac_mod1,cex=(ac_mod1$rmse1)*2,pch=1,add=T) |
|
| 486 |
#title(paste("Averrage RMSE per tile and by ",model_name[i]))
|
|
| 487 |
|
|
| 488 |
dev.off() |
|
| 489 |
|
|
| 490 |
### Ranking by tile... |
|
| 491 |
#df_ac_mod <- |
|
| 492 |
list_df_ac_mod[[i]] <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("rmse","mae","tile_id")]
|
|
| 493 |
} |
|
| 494 |
|
|
| 449 | 495 |
###################### |
| 450 | 496 |
### Figure 7: Delta and clim... |
| 451 | 497 |
|
Also available in: Unified diff
run6 assessment NEX part2: add new plots of assessment of first global run