Revision 8e7d2304
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: 09/16/2014
|
|
| 8 |
#MODIFIED ON: 09/30/2014
|
|
| 9 | 9 |
#Version: 3 |
| 10 | 10 |
#PROJECT: Environmental Layers project |
| 11 | 11 |
#COMMENTS: analyses for run 5 global using 6 specific tiles |
| ... | ... | |
| 106 | 106 |
return(unlist(rast_list)) |
| 107 | 107 |
} |
| 108 | 108 |
|
| 109 |
plot_raster_tb_diagnostic <- function(reg_layer,tb_dat,df_tile_processed,date_selected,mod_selected,var_selected,out_suffix){
|
|
| 110 |
|
|
| 111 |
test <- subset(tb_dat,pred_mod==mod_selected & date==date_selected,select=c("tile_id",var_selected))
|
|
| 112 |
|
|
| 113 |
test_data_tb <- merge(df_tile_processed,test,by="tile_id",all=T) #keep all |
|
| 114 |
test_r <- subset(test_data_tb,select=c("lat","lon","tile_id",var_selected))
|
|
| 115 |
out_suffix_str <- paste(var_selected,mod_selected,date_selected,out_suffix,sep="_") |
|
| 116 |
coord_names <- c("lon","lat")
|
|
| 117 |
l_rast <- rasterize_df_fun(test_r,coord_names,proj_str,out_suffix_str,out_dir=".",file_format=".tif",NA_flag_val=-9999,tolerance_val=0.000120005) |
|
| 118 |
#mod_kr_stack <- stack(mod_kr_l_rast) |
|
| 119 |
d_tb_rast <- stack(l_rast) |
|
| 120 |
names(d_tb_rast) <- names(test_r) |
|
| 121 |
#plot(d_tb_rast) |
|
| 122 |
r <- subset(d_tb_rast,"rmse") |
|
| 123 |
names(r) <- paste(mod_selected,var_selected,date_selected,sep="_") |
|
| 124 |
#plot info: with labels |
|
| 125 |
|
|
| 126 |
res_pix <- 1200 |
|
| 127 |
col_mfrow <- 1 |
|
| 128 |
row_mfrow <- 1 |
|
| 129 |
|
|
| 130 |
png(filename=paste("Figure9_",names(r),"_map_processed_region_",region_name,"_",out_prefix,".png",sep=""),
|
|
| 131 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 132 |
|
|
| 133 |
#plot(reg_layer) |
|
| 134 |
#p1 <- spplot(reg_layer,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color? |
|
| 135 |
title_str <- paste(names(r),"for ", region_name,sep="") |
|
| 136 |
|
|
| 137 |
p0 <- levelplot(r,col.regions=matlab.like(25),margin=F,main=title_str) |
|
| 138 |
p_shp <- layer(sp.polygons(reg_layer, lwd=1, col='black')) |
|
| 139 |
|
|
| 140 |
p <- p0 + p_shp |
|
| 141 |
print(p) |
|
| 142 |
|
|
| 143 |
dev.off() |
|
| 144 |
|
|
| 145 |
} |
|
| 146 |
|
|
| 147 |
create_raster_from_tb_diagnostic <- function(i,list_param){
|
|
| 148 |
#create a raster image using tile centroids and given fields from tb diagnostic data |
|
| 149 |
tb_dat <- list_param$tb_dat |
|
| 150 |
df_tile_processed <- list_param$df_tile_processed |
|
| 151 |
date_selected <- list_param$date_selected[i] |
|
| 152 |
mod_selected <- list_param$mod_selected |
|
| 153 |
var_selected <- list_param$var_selected |
|
| 154 |
out_suffix <- list_param$out_suffix |
|
| 155 |
|
|
| 156 |
test <- subset(tb_dat,pred_mod==mod_selected & date==date_selected,select=c("tile_id",var_selected))
|
|
| 157 |
|
|
| 158 |
test_data_tb <- merge(df_tile_processed,test,by="tile_id",all=T) #keep all |
|
| 159 |
test_r <- subset(test_data_tb,select=c("lat","lon","tile_id",var_selected))
|
|
| 160 |
out_suffix_str <- paste(var_selected,mod_selected,date_selected,out_suffix,sep="_") |
|
| 161 |
coord_names <- c("lon","lat")
|
|
| 162 |
l_rast <- rasterize_df_fun(test_r,coord_names,proj_str,out_suffix_str,out_dir=".",file_format,NA_flag_val,tolerance_val=0.000120005) |
|
| 163 |
#mod_kr_stack <- stack(mod_kr_l_rast) |
|
| 164 |
#d_tb_rast <- stack(l_rast) |
|
| 165 |
#r <- subset(d_tb_rast,var_selected) |
|
| 166 |
#names(d_tb_rast) <- names(test_r) |
|
| 167 |
return(l_rast[4]) |
|
| 168 |
} |
|
| 109 | 169 |
|
| 110 | 170 |
############################## |
| 111 | 171 |
#### Parameters and constants |
| ... | ... | |
| 130 | 190 |
y_var_name <- "dailyTmax" |
| 131 | 191 |
interpolation_method <- c("gam_CAI")
|
| 132 | 192 |
out_prefix<-"run6_global_analyses_09162014" |
| 193 |
mosaic_plot <- FALSE |
|
| 133 | 194 |
|
| 134 | 195 |
proj_str<- CRS_WGS84 |
| 135 | 196 |
file_format <- ".rst" |
| ... | ... | |
| 161 | 222 |
tb <- read.table(file=file.path(out_dir,paste("tb_diagnostic_v_NA","_",out_prefix,".txt",sep="")),sep=",")
|
| 162 | 223 |
tb_s <- read.table(file=file.path(out_dir,paste("tb_diagnostic_s_NA","_",out_prefix,".txt",sep="")),sep=",")
|
| 163 | 224 |
|
| 164 |
tb_month_s_<- read.table(file=file.path(out_dir,paste("tb_month_diagnostic_s_NA","_",out_prefix,".txt",sep="")),sep=",")
|
|
| 225 |
tb_month_s <- read.table(file=file.path(out_dir,paste("tb_month_diagnostic_s_NA","_",out_prefix,".txt",sep="")),sep=",")
|
|
| 165 | 226 |
#tb_month_s <- read.table("tb_month_diagnostic_s_NA_run5_global_analyses_08252014.txt",sep=",")
|
| 166 |
pred_data_month_info <- read.table(file=file.path(out_dir,paste("pred_data_month_info_",out_prefix,".txt",sep="")),sep=",")
|
|
| 167 |
pred_data_day_info <- read.table(file=file.path(out_dir,paste("pred_data_day_info_",out_prefix,".txt",sep="")),sep=",")
|
|
| 227 |
#pred_data_month_info <- read.table(file=file.path(out_dir,paste("pred_data_month_info_",out_prefix,".txt",sep="")),sep=",")
|
|
| 228 |
#pred_data_day_info <- read.table(file=file.path(out_dir,paste("pred_data_day_info_",out_prefix,".txt",sep="")),sep=",")
|
|
| 168 | 229 |
df_tile_processed <- read.table(file=file.path(out_dir,paste("df_tile_processed_",out_prefix,".txt",sep="")),sep=",")
|
| 169 | 230 |
#in_dir_list <- file.path(in_dir1,read.table(file.path(in_dir1,"processed.txt"))$V1) |
| 170 | 231 |
#gam_diagnostic_df <- read.table(file=file.path(out_dir,"gam_diagnostic_df_run4_global_analyses_08142014.txt"),sep=",") |
| ... | ... | |
| 216 | 277 |
centroids_pts[[i]] <-pt |
| 217 | 278 |
shps_tiles[[i]] <- shp1 |
| 218 | 279 |
} |
| 219 |
coord_names <- c("lon","lat")
|
|
| 220 |
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) |
|
| 280 |
#coord_names <- c("lon","lat")
|
|
| 281 |
#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)
|
|
| 221 | 282 |
|
| 222 | 283 |
|
| 223 | 284 |
#plot info: with labels |
| ... | ... | |
| 313 | 374 |
dev.off() |
| 314 | 375 |
|
| 315 | 376 |
|
| 316 |
###################### |
|
| 317 |
### Figure 5: plot accuracy ranked |
|
| 318 |
|
|
| 319 |
#Turn summary table to a point shp |
|
| 320 |
|
|
| 321 |
list_df_ac_mod <- vector("list",length=length(model_name))
|
|
| 322 |
for (i in 1:length(model_name)){
|
|
| 323 |
|
|
| 324 |
ac_mod <- summary_metrics_v[summary_metrics_v$pred_mod==model_name[i],] |
|
| 325 |
### Ranking by tile... |
|
| 326 |
df_ac_mod <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("pred_mod","rmse","mae","tile_id")]
|
|
| 327 |
list_df_ac_mod[[i]] <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("rmse","mae","tile_id")]
|
|
| 328 |
|
|
| 329 |
res_pix <- 480 |
|
| 330 |
col_mfrow <- 1 |
|
| 331 |
row_mfrow <- 1 |
|
| 332 |
|
|
| 333 |
png(filename=paste("Figure5_ac_metrics_ranked_",model_name[i],"_",out_prefix,".png",sep=""),
|
|
| 334 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 335 |
x<- as.character(df_ac_mod$tile_id) |
|
| 336 |
barplot(df_ac_mod$rmse, names.arg=x) |
|
| 337 |
#plot(ac_mod1,cex=sqrt(ac_mod1$rmse),pch=1,add=T) |
|
| 338 |
#plot(ac_mod1,cex=(ac_mod1$rmse1)*2,pch=1,add=T) |
|
| 339 |
title(paste("RMSE ranked by tile for ",model_name[i],sep=""))
|
|
| 340 |
|
|
| 341 |
dev.off() |
|
| 342 |
|
|
| 343 |
} |
|
| 344 |
|
|
| 345 |
|
|
| 346 |
##### Diagnostic gam |
|
| 347 |
#### |
|
| 348 |
#gam_diag_10x10 <- read.table("gam_diagnostic_10x10_df_run4_global_analyses_08142014.txt",sep=",")
|
|
| 349 |
#gam_diag_10x10$month <- as.factor(gam_diag_10x10$month) |
|
| 350 |
|
|
| 351 |
#xyplot(rmse~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 352 |
# pred_mod!="mod_kr"),type="b") |
|
| 353 |
|
|
| 354 |
#xyplot(n~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 355 |
# pred_mod!="mod_kr"),type="h") |
|
| 356 |
|
|
| 357 |
xyplot(n~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 358 |
pred_mod!="mod_kr"),type="h") |
|
| 359 |
|
|
| 360 |
#xyplot(n~month | tile_id + pred_mod,data=subset(as.data.frame(tb_month_s_), |
|
| 361 |
# pred_mod!="mod_kr"),type="h") |
|
| 362 |
|
|
| 363 |
# |
|
| 364 |
## Figure 3b |
|
| 365 |
png(filename=paste("Figure3b_boxplot_overall_region_scaling_",model_name[i],"_",out_prefix,".png",sep=""),
|
|
| 366 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 367 |
|
|
| 368 |
boxplot(rmse~pred_mod,data=tb,ylim=c(0,5),outline=FALSE)#,names=tb$pred_mod) |
|
| 369 |
title("RMSE per model over all tiles")
|
|
| 370 |
|
|
| 371 |
dev.off() |
|
| 372 |
|
|
| 373 | 377 |
################ |
| 374 | 378 |
### Figure 4: plot predicted tiff for specific date per model |
| 375 | 379 |
|
| 376 | 380 |
#y_var_name <-"dailyTmax" |
| 377 | 381 |
#index <-244 #index corresponding to Sept 1 |
| 378 |
index <- 1 #index corresponding to Jan 1 |
|
| 379 |
date_selected <- "20100901" |
|
| 380 |
name_method_var <- paste(interpolation_method,"_",y_var_name,"_",sep="") |
|
| 381 | 382 |
|
| 382 |
pattern_str <- paste("mosaiced","_",name_method_var,"predicted",".*.",date_selected,".*.tif",sep="")
|
|
| 383 |
lf_pred_list <- list.files(pattern=pattern_str) |
|
| 383 |
if (mosaic_plot==TRUE){
|
|
| 384 |
index <- 1 #index corresponding to Jan 1 |
|
| 385 |
date_selected <- "20100901" |
|
| 386 |
name_method_var <- paste(interpolation_method,"_",y_var_name,"_",sep="") |
|
| 387 |
|
|
| 388 |
pattern_str <- paste("mosaiced","_",name_method_var,"predicted",".*.",date_selected,".*.tif",sep="")
|
|
| 389 |
lf_pred_list <- list.files(pattern=pattern_str) |
|
| 384 | 390 |
|
| 385 |
for(i in 1:length(lf_pred_list)){
|
|
| 391 |
for(i in 1:length(lf_pred_list)){
|
|
| 392 |
|
|
| 386 | 393 |
|
| 387 |
r_pred <- raster(lf_pred_list[i]) |
|
| 394 |
r_pred <- raster(lf_pred_list[i])
|
|
| 388 | 395 |
|
| 389 |
res_pix <- 480 |
|
| 390 |
col_mfrow <- 1 |
|
| 391 |
row_mfrow <- 1 |
|
| 392 |
|
|
| 393 |
png(filename=paste("Figure4_models_predicted_surfaces_",model_name[i],"_",name_method_var,"_",data_selected,"_",out_prefix,".png",sep=""),
|
|
| 394 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 396 |
res_pix <- 480 |
|
| 397 |
col_mfrow <- 1 |
|
| 398 |
row_mfrow <- 1 |
|
| 395 | 399 |
|
| 396 |
plot(r_pred) |
|
| 397 |
title(paste("Mosaiced",model_name[i],name_method_var,date_selected,sep=" "))
|
|
| 398 |
dev.off() |
|
| 400 |
png(filename=paste("Figure4_models_predicted_surfaces_",model_name[i],"_",name_method_var,"_",data_selected,"_",out_prefix,".png",sep=""),
|
|
| 401 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 399 | 402 |
|
| 403 |
plot(r_pred) |
|
| 404 |
title(paste("Mosaiced",model_name[i],name_method_var,date_selected,sep=" "))
|
|
| 405 |
dev.off() |
|
| 406 |
} |
|
| 407 |
#Plot Delta and clim... |
|
| 408 |
|
|
| 409 |
## plotting of delta and clim for later scripts... |
|
| 410 |
|
|
| 400 | 411 |
} |
| 401 | 412 |
|
| 402 | 413 |
####### Figure 5... |
| ... | ... | |
| 404 | 415 |
|
| 405 | 416 |
#r_pred <- raster(lf_pred_list[i]) |
| 406 | 417 |
|
| 407 |
res_pix <- 480 |
|
| 408 |
col_mfrow <- 1 |
|
| 409 |
row_mfrow <- 1 |
|
| 410 |
model_name <- as.character(model_name) |
|
| 411 |
|
|
| 412 |
i<-2 |
|
| 413 |
png(filename=paste("Figure5_tiles_with_models_predicted_surfaces_",model_name[i],"_",name_method_var,out_prefix,".png",sep=""),
|
|
| 414 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 415 |
|
|
| 418 |
#res_pix <- 480
|
|
| 419 |
#col_mfrow <- 1
|
|
| 420 |
#row_mfrow <- 1
|
|
| 421 |
#model_name <- as.character(model_name)
|
|
| 422 |
|
|
| 423 |
#i<-2
|
|
| 424 |
#png(filename=paste("Figure5_tiles_with_models_predicted_surfaces_",model_name[i],"_",name_method_var,out_prefix,".png",sep=""),
|
|
| 425 |
# width=col_mfrow*res_pix,height=row_mfrow*res_pix)
|
|
| 426 |
# |
|
| 416 | 427 |
#plot(r_pred) |
| 417 |
plot(reg_layer) |
|
| 428 |
#plot(reg_layer)
|
|
| 418 | 429 |
|
| 419 | 430 |
#title(paste("Mosaiced",model_name[i],name_method_var,date_selected,"with tiles",sep=" "))
|
| 420 | 431 |
|
| 421 | 432 |
#Add polygon tiles... |
| 422 |
for(i in 1:length(list_shp_reg_files)){
|
|
| 423 |
shp1 <- shps_tiles[[i]] |
|
| 424 |
pt <- centroids_pts[[i]] |
|
| 425 |
plot(shp1,add=T,border="blue") |
|
| 433 |
#for(i in 1:length(list_shp_reg_files)){
|
|
| 434 |
# shp1 <- shps_tiles[[i]]
|
|
| 435 |
# pt <- centroids_pts[[i]]
|
|
| 436 |
# plot(shp1,add=T,border="blue")
|
|
| 426 | 437 |
#plot(pt,add=T,cex=2,pch=5) |
| 427 |
label_id <- df_tile_processed$tile_id[i] |
|
| 428 |
text(coordinates(pt)[1],coordinates(pt)[2],labels=i,cex=1,col=c("red"))
|
|
| 429 |
} |
|
| 438 |
# label_id <- df_tile_processed$tile_id[i]
|
|
| 439 |
# text(coordinates(pt)[1],coordinates(pt)[2],labels=i,cex=1,col=c("red"))
|
|
| 440 |
#}
|
|
| 430 | 441 |
#title(paste("Prediction with tiles location 10x10 degrees for ", region_name,sep=""))
|
| 431 |
dev.off() |
|
| 442 |
#dev.off()
|
|
| 432 | 443 |
|
| 433 | 444 |
### |
| 434 | 445 |
|
| ... | ... | |
| 465 | 476 |
#print(p) |
| 466 | 477 |
#dev.off() |
| 467 | 478 |
|
| 479 |
|
|
| 468 | 480 |
###################### |
| 469 |
### Figure 6: plot map of average RMSE per tile at centroids
|
|
| 481 |
### Figure 5: plot accuracy ranked
|
|
| 470 | 482 |
|
| 471 | 483 |
#Turn summary table to a point shp |
| 472 | 484 |
|
| 473 |
coordinates(summary_metrics_v) <- cbind(summary_metrics_v$lon,summary_metrics_v$lat) |
|
| 474 |
proj4string(summary_metrics_v) <- CRS_WGS84 |
|
| 475 |
#lf_list <- lf_pred_list |
|
| 476 |
list_df_ac_mod <- vector("list",length=length(lf_pred_list))
|
|
| 477 |
for (i in 1:length(lf_list)){
|
|
| 485 |
list_df_ac_mod <- vector("list",length=length(model_name))
|
|
| 486 |
for (i in 1:length(model_name)){
|
|
| 478 | 487 |
|
| 479 | 488 |
ac_mod <- summary_metrics_v[summary_metrics_v$pred_mod==model_name[i],] |
| 480 |
r_pred <- raster(lf_list[i]) |
|
| 489 |
### Ranking by tile... |
|
| 490 |
df_ac_mod <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("pred_mod","rmse","mae","tile_id")]
|
|
| 491 |
list_df_ac_mod[[i]] <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("rmse","mae","tile_id")]
|
|
| 481 | 492 |
|
| 482 | 493 |
res_pix <- 480 |
| 483 | 494 |
col_mfrow <- 1 |
| 484 | 495 |
row_mfrow <- 1 |
| 485 | 496 |
|
| 486 |
png(filename=paste("Figure6_ac_metrics_map_centroids_tile_",model_name[i],"_",out_prefix,".png",sep=""),
|
|
| 497 |
png(filename=paste("Figure5_ac_metrics_ranked_",model_name[i],"_",out_prefix,".png",sep=""),
|
|
| 487 | 498 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
| 488 |
|
|
| 489 |
plot(r_pred) |
|
| 490 |
|
|
| 499 |
x<- as.character(df_ac_mod$tile_id) |
|
| 500 |
barplot(df_ac_mod$rmse, names.arg=x) |
|
| 491 | 501 |
#plot(ac_mod1,cex=sqrt(ac_mod1$rmse),pch=1,add=T) |
| 492 |
plot(ac_mod,cex=(ac_mod$rmse^2)/10,pch=1,add=T) |
|
| 493 | 502 |
#plot(ac_mod1,cex=(ac_mod1$rmse1)*2,pch=1,add=T) |
| 494 |
title(paste("Averrage RMSE per tile and by ",model_name[i]))
|
|
| 503 |
title(paste("RMSE ranked by tile for ",model_name[i],sep=""))
|
|
| 495 | 504 |
|
| 496 | 505 |
dev.off() |
| 497 | 506 |
|
| 498 |
### Ranking by tile... |
|
| 499 |
#df_ac_mod <- |
|
| 500 |
list_df_ac_mod[[i]] <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("rmse","mae","tile_id")]
|
|
| 501 | 507 |
} |
| 508 |
|
|
| 509 |
###################### |
|
| 510 |
### Figure 6: plot map of average RMSE per tile at centroids |
|
| 511 |
|
|
| 512 |
#Turn summary table to a point shp |
|
| 513 |
|
|
| 514 |
# coordinates(summary_metrics_v) <- cbind(summary_metrics_v$lon,summary_metrics_v$lat) |
|
| 515 |
# proj4string(summary_metrics_v) <- CRS_WGS84 |
|
| 516 |
# #lf_list <- lf_pred_list |
|
| 517 |
# list_df_ac_mod <- vector("list",length=length(lf_pred_list))
|
|
| 518 |
# for (i in 1:length(lf_list)){
|
|
| 519 |
# |
|
| 520 |
# ac_mod <- summary_metrics_v[summary_metrics_v$pred_mod==model_name[i],] |
|
| 521 |
# r_pred <- raster(lf_list[i]) |
|
| 522 |
# |
|
| 523 |
# res_pix <- 480 |
|
| 524 |
# col_mfrow <- 1 |
|
| 525 |
# row_mfrow <- 1 |
|
| 526 |
# |
|
| 527 |
# png(filename=paste("Figure6_ac_metrics_map_centroids_tile_",model_name[i],"_",out_prefix,".png",sep=""),
|
|
| 528 |
# width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 529 |
# |
|
| 530 |
# plot(r_pred) |
|
| 531 |
# |
|
| 532 |
# #plot(ac_mod1,cex=sqrt(ac_mod1$rmse),pch=1,add=T) |
|
| 533 |
# plot(ac_mod,cex=(ac_mod$rmse^2)/10,pch=1,add=T) |
|
| 534 |
# #plot(ac_mod1,cex=(ac_mod1$rmse1)*2,pch=1,add=T) |
|
| 535 |
# title(paste("Averrage RMSE per tile and by ",model_name[i]))
|
|
| 536 |
# |
|
| 537 |
# dev.off() |
|
| 538 |
# |
|
| 539 |
# ### Ranking by tile... |
|
| 540 |
# #df_ac_mod <- |
|
| 541 |
# list_df_ac_mod[[i]] <- arrange(as.data.frame(ac_mod),desc(rmse))[,c("rmse","mae","tile_id")]
|
|
| 542 |
# } |
|
| 502 | 543 |
|
| 503 | 544 |
#quick kriging... |
| 504 | 545 |
#autokrige(rmse~1,r2,) |
| ... | ... | |
| 544 | 585 |
} |
| 545 | 586 |
|
| 546 | 587 |
###################### |
| 547 |
### Figure 7: Delta and clim... |
|
| 588 |
### Figure 7: Number of predictions: daily and monthly |
|
| 589 |
|
|
| 590 |
#xyplot(rmse~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 591 |
# pred_mod!="mod_kr"),type="b") |
|
| 592 |
|
|
| 593 |
#xyplot(n~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 594 |
# pred_mod!="mod_kr"),type="h") |
|
| 595 |
|
|
| 596 |
|
|
| 548 | 597 |
|
| 549 |
## plotting of delta and clim for later scripts... |
|
| 598 |
# |
|
| 599 |
## Figure 7a |
|
| 600 |
png(filename=paste("Figure7a_number_daily_predictions_per_models","_",out_prefix,".png",sep=""),
|
|
| 601 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 602 |
|
|
| 603 |
xyplot(n~pred_mod | tile_id,data=subset(as.data.frame(summary_metrics_v), |
|
| 604 |
pred_mod!="mod_kr"),type="h") |
|
| 605 |
dev.off() |
|
| 606 |
|
|
| 607 |
## Figure 7b |
|
| 608 |
#png(filename=paste("Figure7b_number_daily_predictions_per_models","_",out_prefix,".png",sep=""),
|
|
| 609 |
# width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 610 |
|
|
| 611 |
#xyplot(n~month | tile_id + pred_mod,data=subset(as.data.frame(tb_month_s), |
|
| 612 |
# pred_mod!="mod_kr"),type="h") |
|
| 613 |
#xyplot(n~month | tile_id,data=subset(as.data.frame(tb_month_s), |
|
| 614 |
# pred_mod="mod_1"),type="h") |
|
| 615 |
#test=subset(as.data.frame(tb_month_s),pred_mod="mod_1") |
|
| 616 |
#table(tb_month_s$month) |
|
| 617 |
#dev.off() |
|
| 618 |
# |
|
| 550 | 619 |
|
| 551 | 620 |
###################### |
| 552 | 621 |
### Figure 9: Plot the number of stations in a processing tile |
| ... | ... | |
| 564 | 633 |
coord_names <- c("lon.x","lat.x")
|
| 565 | 634 |
#l_r_ast <- rasterize_df_fun(df_tile_processed,coord_names,proj_str,out_suffix,out_dir=".",file_format,NA_flag_val,tolerance_val=0.000120005) |
| 566 | 635 |
out_suffix_str <- paste("mod1_",out_suffix)
|
| 567 |
mod1_l_rast <- rasterize_df_fun(mod1_data_tb,coord_names,proj_str,out_suffix_str,out_dir=".",file_format,NA_flag_val,tolerance_val=0.000120005) |
|
| 636 |
mod1_l_rast <- rasterize_df_fun(mod1_data_tb,coord_names,proj_str=CRS_WGS84,out_suffix_str,out_dir=".",file_format,NA_flag_val,tolerance_val=0.000120005)
|
|
| 568 | 637 |
mod1_stack <- stack(mod1_l_rast) |
| 569 | 638 |
names(mod1_stack) <- names(mod1_data_tb) |
| 570 | 639 |
|
| 571 | 640 |
out_suffix_str <- paste("mod2_",out_suffix)
|
| 572 |
mod2_l_rast <- rasterize_df_fun(mod2_data_tb,coord_names,proj_str,out_suffix_str,out_dir=".",file_format,NA_flag_val,tolerance_val=0.000120005) |
|
| 641 |
mod2_l_rast <- rasterize_df_fun(mod2_data_tb,coord_names,proj_str=CRS_WGS84,out_suffix_str,out_dir=".",file_format,NA_flag_val,tolerance_val=0.000120005)
|
|
| 573 | 642 |
mod2_stack <- stack(mod2_l_rast) |
| 574 | 643 |
names(mod2_stack) <- names(mod2_data_tb) |
| 575 | 644 |
|
| 576 | 645 |
out_suffix_str <- paste("mod_kr_",out_suffix)
|
| 577 |
mod_kr_l_rast <- rasterize_df_fun(mod_kr_data_tb,coord_names,proj_str,out_suffix_str,out_dir=".",file_format,NA_flag_val,tolerance_val=0.000120005) |
|
| 646 |
mod_kr_l_rast <- rasterize_df_fun(mod_kr_data_tb,coord_names,proj_str=CRS_WGS84,out_suffix_str,out_dir=".",file_format,NA_flag_val,tolerance_val=0.000120005)
|
|
| 578 | 647 |
mod_kr_stack <- stack(mod_kr_l_rast) |
| 579 | 648 |
names(mod_kr_stack) <- names(mod_kr_data_tb) |
| 580 | 649 |
|
| 581 | 650 |
#Number of daily predictions |
| 582 |
p0 <- levelplot(mod1_stack,layer=21,margin=F) |
|
| 651 |
p0 <- levelplot(mod1_stack,layer=21,margin=F, |
|
| 652 |
main=paste("number_daily_predictions_for_","mod1",sep=""))
|
|
| 583 | 653 |
p<- p0+p_shp |
| 654 |
|
|
| 655 |
########### |
|
| 656 |
## Figure 9a: number of daily predictions |
|
| 657 |
|
|
| 658 |
res_pix <- 1200 |
|
| 659 |
#res_pix <- 480 |
|
| 660 |
col_mfrow <- 1 |
|
| 661 |
row_mfrow <- 1 |
|
| 662 |
|
|
| 663 |
png(filename=paste("Figure9a_raster_map_number_daily_predictions_forr_","mod1","_",out_prefix,".png",sep=""),
|
|
| 664 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 665 |
|
|
| 584 | 666 |
print(p) |
| 585 |
p0<- levelplot(mod2_stack,layer=21,margin=F) |
|
| 586 |
p<- p0+p_shp |
|
| 667 |
|
|
| 668 |
dev.off() |
|
| 669 |
|
|
| 670 |
## Figure 9a |
|
| 671 |
p0 <- levelplot(mod2_stack,layer=21,margin=F) |
|
| 672 |
p <- p0+p_shp |
|
| 673 |
|
|
| 674 |
res_pix <- 1200 |
|
| 675 |
#res_pix <- 480 |
|
| 676 |
col_mfrow <- 1 |
|
| 677 |
row_mfrow <- 1 |
|
| 678 |
png(filename=paste("Figure9a_raster_map_number_daily_predictions_for_","mod2","_",out_prefix,".png",sep=""),
|
|
| 679 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 680 |
|
|
| 587 | 681 |
print(p) |
| 588 |
plot(mod_kr_stack,y=21) |
|
| 589 | 682 |
|
| 590 |
#RMSE |
|
| 591 |
p0 <- levelplot(mod1_stack,layer=10,margin=F,col.regions=matlab.like(25),main="Average RMSE for mod1") |
|
| 683 |
dev.off() |
|
| 684 |
|
|
| 685 |
#plot(mod_kr_stack,y=21) |
|
| 686 |
|
|
| 687 |
######## |
|
| 688 |
###Figure 9b: RMSE plots |
|
| 689 |
## |
|
| 690 |
p0 <- levelplot(mod1_stack,layer=10,margin=F,col.regions=matlab.like(25), |
|
| 691 |
main="Average RMSE for mod1") |
|
| 592 | 692 |
p<- p0+p_shp |
| 693 |
|
|
| 694 |
res_pix <- 1200 |
|
| 695 |
#res_pix <- 480 |
|
| 696 |
col_mfrow <- 1 |
|
| 697 |
row_mfrow <- 1 |
|
| 698 |
png(filename=paste("Figure9b_raster_map_rmse_predictions_for_","mod1","_",out_prefix,".png",sep=""),
|
|
| 699 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 700 |
|
|
| 593 | 701 |
print(p) |
| 594 | 702 |
|
| 595 |
p0 <- levelplot(mod1_stack,layer=10,margin=F,col.regions=matlab.like(25),main="Average RMSE for mod2") |
|
| 703 |
dev.off() |
|
| 704 |
|
|
| 705 |
#print(p) |
|
| 706 |
p0 <- levelplot(mod2_stack,layer=10,margin=F,col.regions=matlab.like(25), |
|
| 707 |
main="Average RMSE for mod2") |
|
| 596 | 708 |
p<- p0+p_shp |
| 709 |
#print(p) |
|
| 710 |
|
|
| 711 |
res_pix <- 1200 |
|
| 712 |
#res_pix <- 480 |
|
| 713 |
col_mfrow <- 1 |
|
| 714 |
row_mfrow <- 1 |
|
| 715 |
png(filename=paste("Figure9b_raster_map_rmse_predictions_for_","mod2","_",out_prefix,".png",sep=""),
|
|
| 716 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 717 |
|
|
| 597 | 718 |
print(p) |
| 598 | 719 |
|
| 720 |
dev.off() |
|
| 721 |
|
|
| 599 | 722 |
#plot(mod1_stack,y=10) |
| 600 | 723 |
#plot(mod2_stack,y=10) |
| 601 | 724 |
#plot(mod_kr_stack,y=10) |
| 602 | 725 |
|
| 726 |
########## Figure 10: map of daily predictions accuracy |
|
| 727 |
|
|
| 603 | 728 |
## Can make maps of daily and accuracy metrics!!! |
| 604 | 729 |
|
| 605 | 730 |
#mod1_tb <- subset(tb,pred_mod=="mod1") |
| ... | ... | |
| 622 | 747 |
#idx <- seq(as.Date('20100101'), as.Date('20101231'), 'day')
|
| 623 | 748 |
#date_l <- strptime(idx[1], "%Y%m%d") # interpolation date being processed |
| 624 | 749 |
l_date_selected <- format(idx, "%Y%m%d") # interpolation date being processed |
| 625 |
|
|
| 750 |
tb_dat <- tb |
|
| 751 |
proj_str <- CRS_WGS84 |
|
| 626 | 752 |
list_param_raster_from_tb <- list(tb_dat,df_tile_processed,l_date_selected, |
| 627 | 753 |
mod_selected,var_selected,out_suffix, |
| 628 | 754 |
proj_str,file_format,NA_value,NA_flag_val) |
| ... | ... | |
| 633 | 759 |
#undebug(create_raster_from_tb_diagnostic) |
| 634 | 760 |
rast <- create_raster_from_tb_diagnostic (i,list_param=list_param_raster_from_tb) |
| 635 | 761 |
l_rast <- lapply(1:length(l_date_selected),FUN=create_raster_from_tb_diagnostic,list_param=list_param_raster_from_tb) |
| 636 |
r_rmse <- stack(l_rast) |
|
| 762 |
r_mod1_rmse <- stack(l_rast) |
|
| 763 |
list_param_raster_from_tb$var_selected <- "n" |
|
| 764 |
l_rast_n<- lapply(1:length(l_date_selected),FUN=create_raster_from_tb_diagnostic,list_param=list_param_raster_from_tb) |
|
| 765 |
r_mod1_n <- stack(l_rast_n) |
|
| 766 |
|
|
| 767 |
##now stack for mod2 |
|
| 768 |
|
|
| 769 |
list_param_raster_from_tb$mod_selected <- "mod2" |
|
| 770 |
list_param_raster_from_tb$var_selected <- "rmse" |
|
| 771 |
l_rast <- lapply(1:length(l_date_selected),FUN=create_raster_from_tb_diagnostic,list_param=list_param_raster_from_tb) |
|
| 772 |
r_mod1_rmse <- stack(l_rast) |
|
| 773 |
list_param_raster_from_tb$var_selected <- "n" |
|
| 774 |
l_rast_n<- lapply(1:length(l_date_selected),FUN=create_raster_from_tb_diagnostic,list_param=list_param_raster_from_tb) |
|
| 775 |
r_mod1_n <- stack(l_rast_n) |
|
| 776 |
|
|
| 777 |
|
|
| 778 |
|
|
| 637 | 779 |
#Make this a function... |
| 638 | 780 |
#test <- subset(tb,pred_mod=="mod1" & date=="20100101",select=c("tile_id","n","mae","rmse","me","r","m50"))
|
| 639 | 781 |
|
| ... | ... | |
| 664 | 806 |
plot_raster_tb_diagnostic(reg_layer,tb_dat,df_tile_processed,date_selected,mod_selected,var_selected,out_suffix) |
| 665 | 807 |
} |
| 666 | 808 |
|
| 667 |
plot_raster_tb_diagnostic <- function(reg_layer,tb_dat,df_tile_processed,date_selected,mod_selected,var_selected,out_suffix){
|
|
| 668 |
|
|
| 669 |
test <- subset(tb_dat,pred_mod==mod_selected & date==date_selected,select=c("tile_id",var_selected))
|
|
| 670 |
|
|
| 671 |
test_data_tb <- merge(df_tile_processed,test,by="tile_id",all=T) #keep all |
|
| 672 |
test_r <- subset(test_data_tb,select=c("lat","lon","tile_id",var_selected))
|
|
| 673 |
out_suffix_str <- paste(var_selected,mod_selected,date_selected,out_suffix,sep="_") |
|
| 674 |
coord_names <- c("lon","lat")
|
|
| 675 |
l_rast <- rasterize_df_fun(test_r,coord_names,proj_str,out_suffix_str,out_dir=".",file_format=".tif",NA_flag_val=-9999,tolerance_val=0.000120005) |
|
| 676 |
#mod_kr_stack <- stack(mod_kr_l_rast) |
|
| 677 |
d_tb_rast <- stack(l_rast) |
|
| 678 |
names(d_tb_rast) <- names(test_r) |
|
| 679 |
#plot(d_tb_rast) |
|
| 680 |
r <- subset(d_tb_rast,"rmse") |
|
| 681 |
names(r) <- paste(mod_selected,var_selected,date_selected,sep="_") |
|
| 682 |
#plot info: with labels |
|
| 683 |
|
|
| 684 |
res_pix <- 1200 |
|
| 685 |
col_mfrow <- 1 |
|
| 686 |
row_mfrow <- 1 |
|
| 687 |
|
|
| 688 |
png(filename=paste("Figure9_",names(r),"_map_processed_region_",region_name,"_",out_prefix,".png",sep=""),
|
|
| 689 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 690 |
|
|
| 691 |
#plot(reg_layer) |
|
| 692 |
#p1 <- spplot(reg_layer,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color? |
|
| 693 |
title_str <- paste(names(r),"for ", region_name,sep="") |
|
| 694 |
|
|
| 695 |
p0 <- levelplot(r,col.regions=matlab.like(25),margin=F,main=title_str) |
|
| 696 |
p_shp <- layer(sp.polygons(reg_layer, lwd=1, col='black')) |
|
| 697 |
|
|
| 698 |
p <- p0 + p_shp |
|
| 699 |
print(p) |
|
| 700 |
|
|
| 701 |
dev.off() |
|
| 702 |
|
|
| 703 |
} |
|
| 704 |
|
|
| 705 |
create_raster_from_tb_diagnostic <- function(i,list_param){
|
|
| 706 |
#create a raster image using tile centroids and given fields from tb diagnostic data |
|
| 707 |
tb_dat <- list_param$tb_dat |
|
| 708 |
df_tile_processed <- list_param$df_tile_processed |
|
| 709 |
date_selected <- list_param$date_selected[i] |
|
| 710 |
mod_selected <- list_param$mod_selected |
|
| 711 |
var_selected <- list_param$var_selected |
|
| 712 |
out_suffix <- list_param$out_suffix |
|
| 713 |
|
|
| 714 |
test <- subset(tb_dat,pred_mod==mod_selected & date==date_selected,select=c("tile_id",var_selected))
|
|
| 715 |
|
|
| 716 |
test_data_tb <- merge(df_tile_processed,test,by="tile_id",all=T) #keep all |
|
| 717 |
test_r <- subset(test_data_tb,select=c("lat","lon","tile_id",var_selected))
|
|
| 718 |
out_suffix_str <- paste(var_selected,mod_selected,date_selected,out_suffix,sep="_") |
|
| 719 |
coord_names <- c("lon","lat")
|
|
| 720 |
l_rast <- rasterize_df_fun(test_r,coord_names,proj_str,out_suffix_str,out_dir=".",file_format,NA_flag_val,tolerance_val=0.000120005) |
|
| 721 |
#mod_kr_stack <- stack(mod_kr_l_rast) |
|
| 722 |
#d_tb_rast <- stack(l_rast) |
|
| 723 |
#r <- subset(d_tb_rast,var_selected) |
|
| 724 |
#names(d_tb_rast) <- names(test_r) |
|
| 725 |
return(l_rast[4]) |
|
| 726 |
} |
|
| 727 | 809 |
|
| 728 | 810 |
dd <- merge(df_tile_processed,pred_data_month_info,by="tile_id") |
| 729 | 811 |
coordinates(dd) <- c(dd$x,dd$y) |
Also available in: Unified diff
run6 assessment NEX part2: generating raster map of accuracy and number of stations