Revision 54e5115e
Added by Benoit Parmentier over 9 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: 03/11/2015
|
|
| 8 |
#MODIFIED ON: 03/23/2015
|
|
| 9 | 9 |
#Version: 4 |
| 10 | 10 |
#PROJECT: Environmental Layers project |
| 11 | 11 |
#COMMENTS: analyses for run 10 global analyses,all regions 1500x4500km and other tiles |
| ... | ... | |
| 305 | 305 |
|
| 306 | 306 |
} |
| 307 | 307 |
|
| 308 |
plot_screen_raster_val<-function(i,list_param){
|
|
| 309 |
##USAGE### |
|
| 310 |
#Screen raster list and produced plot as png. |
|
| 311 |
fname <-list_param$lf_raster_fname[i] |
|
| 312 |
screenRast <- list_param$screenRast |
|
| 313 |
l_dates<- list_param$l_dates |
|
| 314 |
out_dir_str <- list_param$out_dir_str |
|
| 315 |
prefix_str <-list_param$prefix_str |
|
| 316 |
out_suffix_str <- list_param$out_suffix_str |
|
| 317 |
|
|
| 318 |
### START SCRIPT #### |
|
| 319 |
date_proc <- l_dates[i] |
|
| 320 |
|
|
| 321 |
if(screenRast==TRUE){
|
|
| 322 |
r_test <- raster(fname) |
|
| 323 |
|
|
| 324 |
m <- c(-Inf, -100, NA, |
|
| 325 |
-100, 100, 1, |
|
| 326 |
100, Inf, NA) #can change the thresholds later |
|
| 327 |
rclmat <- matrix(m, ncol=3, byrow=TRUE) |
|
| 328 |
rc <- reclassify(r_test, rclmat,filename=paste("rc_tmp_",i,".tif",sep=""),dataType="FLT4S",overwrite=T)
|
|
| 329 |
#file_name <- unlist(strsplit(basename(lf_m[i]),"_")) |
|
| 330 |
extension_str <- extension(filename(r_test)) |
|
| 331 |
raster_name_tmp <- gsub(extension_str,"",basename(filename(r_test))) |
|
| 332 |
raster_name <- file.path(out_dir_str,paste(raster_name_tmp,"_masked.tif",sep="")) |
|
| 333 |
|
|
| 334 |
r_pred <- mask(r_test,rc,filename=raster_name,overwrite=TRUE) |
|
| 335 |
}else{
|
|
| 336 |
r_pred <- raster(fname) |
|
| 337 |
} |
|
| 338 |
|
|
| 339 |
#date_proc <- file_name[7] #specific tot he current naming of files |
|
| 340 |
#date_proc<- "2010010101" |
|
| 341 |
|
|
| 342 |
#paste(raster_name[1:7],collapse="_") |
|
| 343 |
#add filename option later |
|
| 344 |
|
|
| 345 |
res_pix <- 960 |
|
| 346 |
#res_pix <- 480 |
|
| 347 |
|
|
| 348 |
col_mfrow <- 1 |
|
| 349 |
row_mfrow <- 1 |
|
| 350 |
|
|
| 351 |
# png(filename=paste("Figure10_clim_world_mosaics_day_","_",date_proc,"_",tile_size,"_",out_suffix,".png",sep=""),
|
|
| 352 |
# width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 353 |
png(filename=paste(prefix_str,"_",date_proc,"_",tile_size,"_",out_suffix_str,".png",sep=""), |
|
| 354 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 355 |
|
|
| 356 |
plot(r_pred,main=paste("Predicted on ",date_proc , " ", tile_size,sep=""),cex.main=1.5)
|
|
| 357 |
dev.off() |
|
| 358 |
|
|
| 359 |
} |
|
| 360 |
|
|
| 308 | 361 |
############################################ |
| 309 | 362 |
#### Parameters and constants |
| 310 | 363 |
|
| ... | ... | |
| 328 | 381 |
interpolation_method <- c("gam_CAI") #PARAM2
|
| 329 | 382 |
#out_suffix<-"run10_global_analyses_01282015" |
| 330 | 383 |
#out_suffix <- "output_run10_1000x3000_global_analyses_02102015" |
| 331 |
out_suffix <- "run10_1500x4500_global_analyses_03112015" #PARAM3
|
|
| 332 |
out_dir <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_03112015" #PARAM4
|
|
| 384 |
out_suffix <- "run10_1500x4500_global_analyses_03232015" #PARAM3
|
|
| 385 |
out_dir <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_03232015" #PARAM4
|
|
| 333 | 386 |
create_out_dir_param <- FALSE #PARAM 5 |
| 334 | 387 |
|
| 335 | 388 |
mosaic_plot <- FALSE #PARAM6 |
| ... | ... | |
| 355 | 408 |
mulitple_region <- TRUE #PARAM 12 |
| 356 | 409 |
|
| 357 | 410 |
region_name <- "world" #PARAM 13 |
| 411 |
plot_region <- FALSE |
|
| 358 | 412 |
|
| 359 | 413 |
########################## START SCRIPT ############################## |
| 360 | 414 |
|
| ... | ... | |
| 716 | 770 |
#plot(ac_mod1,cex=sqrt(ac_mod1$rmse),pch=1,add=T) |
| 717 | 771 |
#plot(ac_mod,cex=(ac_mod$rmse^2)/10,pch=1,col="red",add=T) |
| 718 | 772 |
|
| 719 |
coordinates(ac_mod) <- ac_mod[,c("lon","lat")]
|
|
| 720 |
#coordinates(ac_mod) <- ac_mod[,c("lon.x","lat.x")] #solve this later
|
|
| 773 |
#coordinates(ac_mod) <- ac_mod[,c("lon","lat")]
|
|
| 774 |
coordinates(ac_mod) <- ac_mod[,c("lon.x","lat.x")] #solve this later
|
|
| 721 | 775 |
p_shp <- layer(sp.polygons(reg_layer, lwd=1, col='black')) |
| 722 | 776 |
#title("(a) Mean for 1 January")
|
| 723 | 777 |
p <- bubble(ac_mod,"rmse",main=paste("Averrage RMSE per tile and by ",model_name[i]))
|
| ... | ... | |
| 739 | 793 |
sum(df_tile_processed$metrics_v)/length(df_tile_processed$metrics_v) #87.85% of tiles with info |
| 740 | 794 |
|
| 741 | 795 |
#coordinates |
| 742 |
coordinates(summary_metrics_v) <- c("lon","lat")
|
|
| 743 |
#coordinates(summary_metrics_v) <- c("lon.y","lat.y")
|
|
| 796 |
#coordinates(summary_metrics_v) <- c("lon","lat")
|
|
| 797 |
coordinates(summary_metrics_v) <- c("lon.y","lat.y")
|
|
| 744 | 798 |
|
| 745 | 799 |
threshold_missing_day <- c(367,365,300,200) |
| 746 | 800 |
|
| ... | ... | |
| 970 | 1024 |
|
| 971 | 1025 |
lf_world_pred <-list.files(path=file.path(out_dir,"mosaics"), |
| 972 | 1026 |
pattern=paste("^world_mosaics.*.tif$",sep=""),full.names=T)
|
| 973 |
) |
|
| 1027 |
|
|
| 974 | 1028 |
|
| 975 | 1029 |
#mosaic_list_mean <- test_list |
| 976 | 1030 |
#out_rastnames <- "world_test_mosaic_20100101" |
| ... | ... | |
| 985 | 1039 |
list_param_plot_screen_raster <- list(lf_raster_fname,screenRast,l_dates,out_dir,prefix_str,out_suffix) |
| 986 | 1040 |
names(list_param_plot_screen_raster) <- c("lf_raster_fname","screenRast","l_dates","out_dir_str","prefix_str","out_suffix_str")
|
| 987 | 1041 |
|
| 988 |
undebug(plot_screen_raster_val) |
|
| 1042 |
#undebug(plot_screen_raster_val)
|
|
| 989 | 1043 |
|
| 990 | 1044 |
#world_m_list1<- plot_screen_raster_val(1,list_param_plot_screen_raster) |
| 991 | 1045 |
#world_m_list <- mclapply(1:10, list_param=list_param_plot_screen_raster, plot_screen_raster_val,mc.preschedule=FALSE,mc.cores = 5) #This is the end bracket from mclapply(...) statement |
| 992 |
world_m_list <- mclapply(1:30, list_param=list_param_plot_screen_raster, plot_screen_raster_val,mc.preschedule=FALSE,mc.cores = 7) #This is the end bracket from mclapply(...) statement |
|
| 993 |
|
|
| 994 |
plot_screen_raster_val<-function(i,list_param){
|
|
| 995 |
##USAGE### |
|
| 996 |
#Screen raster list and produced plot as png. |
|
| 997 |
fname <-list_param$lf_raster_fname[i] |
|
| 998 |
screenRast <- list_param$screenRast |
|
| 999 |
l_dates<- list_param$l_dates |
|
| 1000 |
out_dir_str <- list_param$out_dir_str |
|
| 1001 |
prefix_str <-list_param$prefix_str |
|
| 1002 |
out_suffix_str <- list_param$out_suffix_str |
|
| 1003 |
|
|
| 1004 |
### START SCRIPT #### |
|
| 1005 |
date_proc <- l_dates[i] |
|
| 1006 |
|
|
| 1007 |
if(screenRast==TRUE){
|
|
| 1008 |
r_test <- raster(fname) |
|
| 1009 |
|
|
| 1010 |
m <- c(-Inf, -100, NA, |
|
| 1011 |
-100, 100, 1, |
|
| 1012 |
100, Inf, NA) #can change the thresholds later |
|
| 1013 |
rclmat <- matrix(m, ncol=3, byrow=TRUE) |
|
| 1014 |
rc <- reclassify(r_test, rclmat,filename=paste("rc_tmp_",i,".tif",sep=""),dataType="FLT4S",overwrite=T)
|
|
| 1015 |
#file_name <- unlist(strsplit(basename(lf_m[i]),"_")) |
|
| 1016 |
extension_str <- extension(filename(r_test)) |
|
| 1017 |
raster_name_tmp <- gsub(extension_str,"",basename(filename(r_test))) |
|
| 1018 |
raster_name <- file.path(out_dir_str,paste(raster_name_tmp,"_masked.tif",sep="")) |
|
| 1019 |
|
|
| 1020 |
r_pred <- mask(r_test,rc,filename=raster_name,overwrite=TRUE) |
|
| 1021 |
}else{
|
|
| 1022 |
r_pred <- raster(fname) |
|
| 1023 |
} |
|
| 1024 |
|
|
| 1025 |
#date_proc <- file_name[7] #specific tot he current naming of files |
|
| 1026 |
#date_proc<- "2010010101" |
|
| 1027 |
|
|
| 1028 |
#paste(raster_name[1:7],collapse="_") |
|
| 1029 |
#add filename option later |
|
| 1030 |
|
|
| 1031 |
res_pix <- 960 |
|
| 1032 |
#res_pix <- 480 |
|
| 1033 |
|
|
| 1034 |
col_mfrow <- 1 |
|
| 1035 |
row_mfrow <- 1 |
|
| 1036 |
|
|
| 1037 |
# png(filename=paste("Figure10_clim_world_mosaics_day_","_",date_proc,"_",tile_size,"_",out_suffix,".png",sep=""),
|
|
| 1038 |
# width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 1039 |
png(filename=paste(prefix_str,"_",date_proc,"_",tile_size,"_",out_suffix_str,".png",sep=""), |
|
| 1040 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 1041 |
|
|
| 1042 |
plot(r_pred,main=paste("Predicted on ",date_proc , " ", tile_size,sep=""),cex.main=1.5)
|
|
| 1043 |
dev.off() |
|
| 1044 |
|
|
| 1045 |
} |
|
| 1046 |
|
|
| 1046 |
world_m_list <- mclapply(1:length(l_dates), list_param=list_param_plot_screen_raster, plot_screen_raster_val,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement |
|
| 1047 | 1047 |
|
| 1048 | 1048 |
#lf_world_mask_reg <- vector("list",length=length(lf_world_pred))
|
| 1049 | 1049 |
#for(i in 1:length(lf_world_pred)){
|
Also available in: Unified diff
global assessment part2 1500x4500km production of figures, checking for missing tiles