Revision 91e41f72
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: 12/23/2014
|
|
| 9 |
#Version: 3
|
|
| 8 |
#MODIFIED ON: 01/02/2015
|
|
| 9 |
#Version: 4
|
|
| 10 | 10 |
#PROJECT: Environmental Layers project |
| 11 |
#COMMENTS: analyses for run 10 global analyses, Europe 1000x300km |
|
| 11 |
#COMMENTS: analyses for run 10 global analyses, Europe, Australia, 1000x300km
|
|
| 12 | 12 |
################################################################################################# |
| 13 | 13 |
|
| 14 | 14 |
### Loading R library and packages |
| ... | ... | |
| 206 | 206 |
return(combined_spdf) |
| 207 | 207 |
} |
| 208 | 208 |
|
| 209 |
############################## |
|
| 209 |
plot_daily_mosaics <- function(i,list_param_plot_daily_mosaics){
|
|
| 210 |
#Purpose: |
|
| 211 |
#This functions mask mosaics files for a default range (-100,100 deg). |
|
| 212 |
#It produces a masked tif in a given dataType format (FLT4S) |
|
| 213 |
#It creates a figure of mosaiced region being interpolated. |
|
| 214 |
#Author: Benoit Parmentier |
|
| 215 |
#Parameters: |
|
| 216 |
#lf_m: list of files |
|
| 217 |
#reg_name:region name with tile size included |
|
| 218 |
#To do: |
|
| 219 |
#Add filenames |
|
| 220 |
#Add range |
|
| 221 |
#Add output dir |
|
| 222 |
#Add dataType_val option |
|
| 223 |
|
|
| 224 |
##### BEGIN ######## |
|
| 225 |
|
|
| 226 |
#Parse the list of parameters |
|
| 227 |
lf_m <- list_param_plot_daily_mosaics$lf_m |
|
| 228 |
reg_name <- list_param_plot_daily_mosaics$reg_name |
|
| 229 |
out_dir_str <- list_param_plot_daily_mosaics$out_dir_str |
|
| 230 |
out_suffix <- list_param_plot_daily_mosaics$out_suffix |
|
| 231 |
|
|
| 232 |
|
|
| 233 |
#list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str) |
|
| 234 |
|
|
| 235 |
|
|
| 236 |
#rast_list <- vector("list",length=length(lf_m))
|
|
| 237 |
r_test<- raster(lf_m[i]) |
|
| 238 |
|
|
| 239 |
m <- c(-Inf, -100, NA, |
|
| 240 |
-100, 100, 1, |
|
| 241 |
100, Inf, NA) #can change the thresholds later |
|
| 242 |
rclmat <- matrix(m, ncol=3, byrow=TRUE) |
|
| 243 |
rc <- reclassify(r_test, rclmat,filename=paste("rc_tmp_",i,".tif",sep=""),dataType="FLT4S",overwrite=T)
|
|
| 244 |
file_name <- unlist(strsplit(basename(lf_m[i]),"_")) |
|
| 245 |
date_proc <- file_name[7] #specific tot he current naming of files |
|
| 246 |
#paste(raster_name[1:7],collapse="_") |
|
| 247 |
#add filename option later |
|
| 248 |
extension_str <- extension(filename(r_test)) |
|
| 249 |
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="")) |
|
| 251 |
r_pred <- mask(r_test,rc,filename=raster_name,overwrite=TRUE) |
|
| 252 |
|
|
| 253 |
res_pix <- 1200 |
|
| 254 |
#res_pix <- 480 |
|
| 255 |
|
|
| 256 |
col_mfrow <- 1 |
|
| 257 |
row_mfrow <- 1 |
|
| 258 |
|
|
| 259 |
png(filename=paste("Figure9_clim_mosaics_day_test","_",date_proc,"_",reg_name,"_",out_suffix,".png",sep=""),
|
|
| 260 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 261 |
|
|
| 262 |
plot(r_pred,main=paste("climatology month ",date_proc , " ", reg_name,sep=""),cex.main=1.5)
|
|
| 263 |
dev.off() |
|
| 264 |
|
|
| 265 |
return(raster_name) |
|
| 266 |
|
|
| 267 |
} |
|
| 268 |
|
|
| 269 |
############################################ |
|
| 210 | 270 |
#### Parameters and constants |
| 211 | 271 |
|
| 212 | 272 |
#on ATLAS |
| ... | ... | |
| 215 | 275 |
#in_dir1 <- "/data/project/layers/commons/NEX_data/output_run6_global_analyses_09162014/output20Deg2" |
| 216 | 276 |
# parent output dir for the curent script analyes |
| 217 | 277 |
#out_dir <-"/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/" #On NCEAS Atlas |
| 218 |
out_dir <-"/data/project/layers/commons/NEX_data/output_run10_global_analyses_12152014/"
|
|
| 278 |
out_dir <-"/data/project/layers/commons/NEX_data/output_run10_global_analyses_12232014/"
|
|
| 219 | 279 |
# input dir containing shapefiles defining tiles |
| 220 | 280 |
#in_dir_shp <- "/data/project/layers/commons/NEX_data/output_run5_global_analyses_08252014/output/subset/shapefiles" |
| 221 | 281 |
|
| ... | ... | |
| 228 | 288 |
|
| 229 | 289 |
y_var_name <- "dailyTmax" |
| 230 | 290 |
interpolation_method <- c("gam_CAI")
|
| 231 |
out_prefix<-"run10_global_analyses_12152014"
|
|
| 291 |
out_prefix<-"run10_global_analyses_12232014"
|
|
| 232 | 292 |
mosaic_plot <- FALSE |
| 233 | 293 |
|
| 234 | 294 |
#CRS_locs_WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| ... | ... | |
| 277 | 337 |
tb$pred_mod <- as.character(tb$pred_mod) |
| 278 | 338 |
summary_metrics_v$pred_mod <- as.character(summary_metrics_v$pred_mod) |
| 279 | 339 |
|
| 280 |
mulitple_region <- FALSE
|
|
| 340 |
mulitple_region <- TRUE
|
|
| 281 | 341 |
|
| 282 | 342 |
#multiple regions? |
| 283 | 343 |
if(mulitple_region==TRUE){
|
| 284 | 344 |
df_tile_processed$reg <- basename(dirname(as.character(df_tile_processed$path_NEX))) |
| 285 | 345 |
} |
| 286 | 346 |
|
| 287 |
#This part is specifically related to this run : dropping model with elev |
|
| 288 |
#tb <- merge(tb,df_tile_processed,by="tile_id") |
|
| 289 |
|
|
| 290 |
#tb_tmp_reg5 <- subset(tb,reg=="reg5") |
|
| 291 |
#tb_tmp_reg4 <- subset(tb,reg=="reg4") |
|
| 292 |
#tb_tmp_reg4 <- subset(tb_tmp_reg4,pred_mod!="mod1") #remove mod1 because it is not in reg5 |
|
| 293 |
#tb_tmp_reg4$pred_mod <- replace(tb_tmp_reg4$pred_mod, tb_tmp_reg4$pred_mod=="mod2", "mod1") |
|
| 294 |
#tb <- rbind(tb_tmp_reg4,tb_tmp_reg5) |
|
| 295 |
|
|
| 296 |
#ac_mod <- summary_metrics_v[summary_metrics_v$pred_mod==model_name[i],] |
|
| 297 |
#summary_metrics_v <- merge(summary_metrics_v,df_tile_processed,by="tile_id") |
|
| 298 |
#table(summary_metrics_v$reg) |
|
| 299 |
|
|
| 300 |
#summary_metrics_v_reg5 <- subset(summary_metrics_v,reg=="reg5") |
|
| 301 |
#summary_metrics_v_reg4 <- subset(summary_metrics_v,reg=="reg4") |
|
| 302 |
#summary_metrics_v_reg4 <- subset(summary_metrics_v_reg4,pred_mod!="mod1") #remove mod1 because it is not in reg5 |
|
| 303 |
#summary_metrics_v_reg4$pred_mod <- replace(summary_metrics_v_reg4$pred_mod, |
|
| 304 |
# summary_metrics_v_reg4$pred_mod=="mod2", "mod1") |
|
| 305 |
#summary_metrics_v <- rbind(summary_metrics_v_reg4,summary_metrics_v_reg5) |
|
| 306 | 347 |
|
| 307 | 348 |
############### |
| 308 | 349 |
### Figure 1: plot location of the study area with tiles processed |
| ... | ... | |
| 666 | 707 |
coordinates(summary_metrics_v) <- c("lon","lat")
|
| 667 | 708 |
#coordinates(summary_metrics_v) <- c("lon.y","lat.y")
|
| 668 | 709 |
|
| 669 |
summary_metrics_v$n_missing <- summary_metrics_v$n == 365 |
|
| 670 |
summary_metrics_v$n_missing <- summary_metrics_v$n < 365 |
|
| 671 |
summary_metrics_v$n_missing <- summary_metrics_v$n < 300 |
|
| 710 |
threshold_missing_day <- c(367,365,300,200) |
|
| 672 | 711 |
|
| 673 | 712 |
nb<-nrow(subset(summary_metrics_v,model_name=="mod1")) |
| 674 | 713 |
sum(subset(summary_metrics_v,model_name=="mod1")$n_missing)/nb #33/35 |
| ... | ... | |
| 676 | 715 |
## Make this a figure... |
| 677 | 716 |
|
| 678 | 717 |
#plot(summary_metrics_v) |
| 679 |
i <- 1 |
|
| 680 |
model_name[1] |
|
| 681 |
p_shp <- layer(sp.polygons(reg_layer, lwd=1, col='black')) |
|
| 682 |
#title("(a) Mean for 1 January")
|
|
| 683 |
p <- bubble(summary_metrics_v,"n_missing",main=paste("Missing per tile and by ",model_name[i]))
|
|
| 684 |
p1 <- p+p_shp |
|
| 685 |
print(p1) |
|
| 718 |
#Make this a function later so that we can explore many thresholds... |
|
| 719 |
|
|
| 720 |
j<-1 #for model name 1 |
|
| 721 |
for(i in 1:length(threshold_missing_day)){
|
|
| 722 |
|
|
| 723 |
#summary_metrics_v$n_missing <- summary_metrics_v$n == 365 |
|
| 724 |
#summary_metrics_v$n_missing <- summary_metrics_v$n < 365 |
|
| 725 |
summary_metrics_v$n_missing <- summary_metrics_v$n < threshold_missing_day[i] |
|
| 726 |
summary_metrics_v_subset <- subset(summary_metrics_v,model_name=="mod1") |
|
| 727 |
|
|
| 728 |
#res_pix <- 1200 |
|
| 729 |
res_pix <- 960 |
|
| 730 |
|
|
| 731 |
col_mfrow <- 1 |
|
| 732 |
row_mfrow <- 1 |
|
| 733 |
|
|
| 734 |
png(filename=paste("Figure7a_ac_metrics_map_centroids_tile_",model_name[j],"_","missing_day_",threshold_missing_day[i],
|
|
| 735 |
"_",out_prefix,".png",sep=""), |
|
| 736 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 737 |
|
|
| 738 |
model_name[j] |
|
| 739 |
|
|
| 740 |
p_shp <- layer(sp.polygons(reg_layer, lwd=1, col='black')) |
|
| 741 |
#title("(a) Mean for 1 January")
|
|
| 742 |
p <- bubble(summary_metrics_v_subset,"n_missing",main=paste("Missing per tile and by ",model_name[j]," for ",
|
|
| 743 |
threshold_missing_day[i])) |
|
| 744 |
p1 <- p+p_shp |
|
| 745 |
print(p1) |
|
| 746 |
#plot(ac_mod1,cex=(ac_mod1$rmse1)*2,pch=1,add=T) |
|
| 747 |
#title(paste("Averrage RMSE per tile and by ",model_name[i]))
|
|
| 748 |
|
|
| 749 |
dev.off() |
|
| 750 |
} |
|
| 686 | 751 |
|
| 687 | 752 |
###################### |
| 688 | 753 |
### Figure 7: Number of predictions: daily and monthly |
| ... | ... | |
| 738 | 803 |
#dev.off() |
| 739 | 804 |
# |
| 740 | 805 |
|
| 741 |
## plot mosaics... |
|
| 742 | 806 |
|
| 743 |
lf_mosaics_reg5 <- mixedsort(list.files(path="/data/project/layers/commons/NEX_data/output_run10_global_analyses_11302014/mosaics/reg5", |
|
| 744 |
pattern="CAI_TMAX_clim_month_.*_mod1_all.tif", full.names=T)) |
|
| 807 |
########################################################## |
|
| 808 |
##### Figure 8: Breaking down accuaracy by regions!! ##### |
|
| 809 |
|
|
| 810 |
|
|
| 811 |
#This part is specifically related to this run : dropping model with elev |
|
| 812 |
#tb <- merge(tb,df_tile_processed,by="tile_id") |
|
| 813 |
|
|
| 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) |
|
| 819 |
|
|
| 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) |
|
| 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) |
|
| 830 |
|
|
| 831 |
|
|
| 832 |
##################################################### |
|
| 833 |
#### Figure 9: plotting mosaics of regions ########### |
|
| 834 |
## plot mosaics... |
|
| 835 |
l_reg_name <- unique(df_tile_processed$reg) |
|
| 836 |
#lf_mosaics_reg5 <- mixedsort(list.files(path="/data/project/layers/commons/NEX_data/output_run10_global_analyses_11302014/mosaics/reg5", |
|
| 837 |
# pattern="CAI_TMAX_clim_month_.*_mod1_all.tif", full.names=T)) |
|
| 838 |
lf_mosaics_reg <- vector("list",length=length(l_reg_name))
|
|
| 839 |
for(i in 1:length(l_reg_name)){
|
|
| 840 |
lf_mosaics_reg[[i]] <- mixedsort(list.files( |
|
| 841 |
path=file.path(out_dir,"mosaics"), |
|
| 842 |
#pattern="reg6_.*._CAI_TMAX_clim_month_.*._mod1_all_mean.tif", |
|
| 843 |
pattern=paste(l_reg_name[i],".*._CAI_TMAX_clim_month_.*._mod1_all_mean.tif",sep=""), |
|
| 844 |
full.names=T)) |
|
| 845 |
} |
|
| 745 | 846 |
|
| 746 | 847 |
#r_reg5 <- stack(lf_mosaics_reg5) |
| 747 | 848 |
lf_m <- lf_mosaics_reg5[1:12] |
| ... | ... | |
| 845 | 946 |
dev.off() |
| 846 | 947 |
} |
| 847 | 948 |
|
| 848 |
#################### IMAGE DIFFERENCING BETWEEN 1000x3000 and 1500x4500 predictions for reg4 and reg5 |
|
| 849 |
|
|
| 850 |
#### Use previous results for differencing |
|
| 851 |
### for reg4_1500x4500: use "mod2 in this case... |
|
| 852 |
|
|
| 853 |
out_prefix_str <- "reg4_1500x4500" |
|
| 854 |
lf_mosaics_reg4_1500x4500 <- mixedsort( |
|
| 855 |
list.files(path= |
|
| 856 |
"/data/project/layers/commons/NEX_data/output_run10_global_analyses_12152014/mosaics/reg4_1500x4500", |
|
| 857 |
pattern=".*._mod2_all_mean.tif$",full.names=T) |
|
| 858 |
) |
|
| 859 |
lf_m <- lf_mosaics_reg4_1500x4500 |
|
| 860 |
out_dir_str <- "/data/project/layers/commons/NEX_data/output_run10_global_analyses_12152014/mosaics/" |
|
| 861 |
|
|
| 862 |
reg_name <- "reg4_1500x4500" |
|
| 863 |
#lapply() |
|
| 864 |
#list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir=out_dir_str,out_suffix=out_prefix) |
|
| 865 |
list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str,out_suffix=out_prefix) |
|
| 866 |
|
|
| 867 |
#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) |
|
| 868 |
#debug(plot_daily_mosaics) |
|
| 869 |
#test<- plot_daily_mosaics(1,list_param=list_param_plot_daily_mosaics)#,mc.preschedule)#=FALSE,mc.cores = 6) |
|
| 870 |
|
|
| 871 |
lf_m_mask_reg4_1500x4500 <- mclapply(1:length(lf_m),FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 6) |
|
| 872 |
|
|
| 873 |
### for reg5_1500x4500: use "mod1 in this case |
|
| 874 |
|
|
| 875 |
out_prefix_str <- "reg5_1500x4500" |
|
| 876 |
lf_mosaics_reg5 <- mixedsort( |
|
| 877 |
list.files(path= |
|
| 878 |
"/data/project/layers/commons/NEX_data/output_run10_global_analyses_12152014/mosaics/reg5_1500x4500", |
|
| 879 |
pattern=".*._mod1_all_mean.tif$",full.names=T) |
|
| 880 |
) |
|
| 881 |
lf_m <- lf_mosaics_reg5 |
|
| 882 |
out_dir_str <- "/data/project/layers/commons/NEX_data/output_run10_global_analyses_12152014/mosaics/" |
|
| 883 |
reg_name <- "reg5_1500x4500" |
|
| 884 |
#lapply() |
|
| 885 |
list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str,out_suffix=out_prefix) |
|
| 886 |
#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) |
|
| 887 |
|
|
| 888 |
lf_m_mask_reg5_1500x4500 <- mclapply(1:length(lf_m),FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 6) |
|
| 889 |
|
|
| 890 |
### for reg5_1500x4500: use "mod1 in this case, this is Africa |
|
| 891 |
|
|
| 892 |
out_prefix_str <- "reg5_1000x3000" |
|
| 893 |
lf_mosaics_reg5 <- mixedsort( |
|
| 894 |
list.files(path= |
|
| 895 |
"/data/project/layers/commons/NEX_data/output_run10_global_analyses_12152014/mosaics/reg5_1000x4500", |
|
| 896 |
pattern=".*._mod1_all_mean.tif$",full.names=T) |
|
| 897 |
) |
|
| 898 |
lf_m <- lf_mosaics_reg5 |
|
| 899 |
out_dir_str <- "/data/project/layers/commons/NEX_data/output_run10_global_analyses_12152014/mosaics/" |
|
| 900 |
reg_name <- "reg5_1500x4500" |
|
| 901 |
#lapply() |
|
| 902 |
list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str,out_suffix=out_prefix) |
|
| 903 |
#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) |
|
| 904 |
|
|
| 905 |
lf_m_mask_reg5_1500x4500 <- mclapply(1:length(lf_m),FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 6) |
|
| 906 |
|
|
| 907 |
|
|
| 908 |
|
|
| 909 |
|
|
| 910 |
### for reg4_1000x3000: use "mod1 in this case |
|
| 911 |
|
|
| 912 |
out_prefix_str <- "reg4_1000x3000" |
|
| 913 |
lf_mosaics_reg5 <- mixedsort( |
|
| 914 |
list.files(path= |
|
| 915 |
"/data/project/layers/commons/NEX_data/output_run10_global_analyses_12152014/mosaics/reg5_1000x3000", |
|
| 916 |
pattern=".*._mod1_all_mean.tif$",full.names=T) |
|
| 917 |
) |
|
| 918 |
lf_m <- lf_mosaics_reg5 |
|
| 919 |
out_dir_str <- "/data/project/layers/commons/NEX_data/output_run10_global_analyses_12152014/mosaics/" |
|
| 920 |
reg_name <- "reg5_1000x4500" |
|
| 921 |
#lapply() |
|
| 922 |
list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str,out_suffix=out_prefix) |
|
| 923 |
#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) |
|
| 924 |
|
|
| 925 |
lf_m_mask_reg5_1000x4500 <- mclapply(1:length(lf_m),FUN=plot_daily_mosaics,list_param=list_param_plot_daily_mosaics,mc.preschedule=FALSE,mc.cores = 6) |
|
| 926 |
|
|
| 927 |
|
|
| 928 |
plot_daily_mosaics <- function(i,list_param_plot_daily_mosaics){
|
|
| 929 |
#Purpose: |
|
| 930 |
#This functions mask mosaics files for a default range (-100,100 deg). |
|
| 931 |
#It produces a masked tif in a given dataType format (FLT4S) |
|
| 932 |
#It creates a figure of mosaiced region being interpolated. |
|
| 933 |
#Author: Benoit Parmentier |
|
| 934 |
#Parameters: |
|
| 935 |
#lf_m: list of files |
|
| 936 |
#reg_name:region name with tile size included |
|
| 937 |
#To do: |
|
| 938 |
#Add filenames |
|
| 939 |
#Add range |
|
| 940 |
#Add output dir |
|
| 941 |
#Add dataType_val option |
|
| 942 |
|
|
| 943 |
##### BEGIN ######## |
|
| 944 |
|
|
| 945 |
#Parse the list of parameters |
|
| 946 |
lf_m <- list_param_plot_daily_mosaics$lf_m |
|
| 947 |
reg_name <- list_param_plot_daily_mosaics$reg_name |
|
| 948 |
out_dir_str <- list_param_plot_daily_mosaics$out_dir_str |
|
| 949 |
out_suffix <- list_param_plot_daily_mosaics$out_suffix |
|
| 950 |
|
|
| 951 |
|
|
| 952 |
#list_param_plot_daily_mosaics <- list(lf_m=lf_m,reg_name=reg_name,out_dir_str=out_dir_str) |
|
| 953 |
|
|
| 954 |
|
|
| 955 |
#rast_list <- vector("list",length=length(lf_m))
|
|
| 956 |
r_test<- raster(lf_m[i]) |
|
| 957 |
|
|
| 958 |
m <- c(-Inf, -100, NA, |
|
| 959 |
-100, 100, 1, |
|
| 960 |
100, Inf, NA) #can change the thresholds later |
|
| 961 |
rclmat <- matrix(m, ncol=3, byrow=TRUE) |
|
| 962 |
rc <- reclassify(r_test, rclmat,filename=paste("rc_tmp_",i,".tif",sep=""),dataType="FLT4S",overwrite=T)
|
|
| 963 |
file_name <- unlist(strsplit(basename(lf_m[i]),"_")) |
|
| 964 |
date_proc <- file_name[7] #specific tot he current naming of files |
|
| 965 |
#paste(raster_name[1:7],collapse="_") |
|
| 966 |
#add filename option later |
|
| 967 |
extension_str <- extension(filename(r_test)) |
|
| 968 |
raster_name_tmp <- gsub(extension_str,"",basename(filename(r_test))) |
|
| 969 |
raster_name <- file.path(out_dir_str,paste(raster_name_tmp,"_",date_proc,"_masked.tif",sep="")) |
|
| 970 |
r_pred <- mask(r_test,rc,filename=raster_name,overwrite=TRUE) |
|
| 971 |
|
|
| 972 |
res_pix <- 1200 |
|
| 973 |
#res_pix <- 480 |
|
| 974 |
|
|
| 975 |
col_mfrow <- 1 |
|
| 976 |
row_mfrow <- 1 |
|
| 977 |
|
|
| 978 |
png(filename=paste("Figure9_clim_mosaics_day_test","_",date_proc,"_",reg_name,"_",out_suffix,".png",sep=""),
|
|
| 979 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 980 |
|
|
| 981 |
plot(r_pred,main=paste("climatology month ",date_proc , " ", reg_name,sep=""),cex.main=1.5)
|
|
| 982 |
dev.off() |
|
| 983 |
|
|
| 984 |
return(raster_name) |
|
| 985 |
|
|
| 986 |
} |
|
| 987 |
|
|
| 988 |
|
|
| 989 |
###################### |
|
| 990 |
### Figure 10: Plot the difference in mosaics for processing tiles at 1500x4500 and 1000x3500 |
|
| 991 |
|
|
| 992 |
|
|
| 993 |
reg_name <- "reg4" |
|
| 994 |
for(i in 1:length(lf_m)){
|
|
| 995 |
r_test<- raster(lf_m[i]) |
|
| 996 |
|
|
| 997 |
m <- c(-Inf, -100, NA, |
|
| 998 |
-100, 100, 1, |
|
| 999 |
100, Inf, NA) |
|
| 1000 |
rclmat <- matrix(m, ncol=3, byrow=TRUE) |
|
| 1001 |
rc <- reclassify(r_test, rclmat,filename="rc.tif",dataType="FLT4S",overwrite=T) |
|
| 1002 |
raster_name <- unlist(strsplit(basename(lf_m[i]),"_")) |
|
| 1003 |
date_proc <- raster_name[5] |
|
| 1004 |
#paste(raster_name[1:7],collapse="_") |
|
| 1005 |
r_pred <- mask(r_test,rc,filename=paste("CAI_TMAX_clim_month_mod1_all_",reg_name,"_",date_proc,"_masked.tif",sep=""),overwrite=TRUE)
|
|
| 1006 |
|
|
| 1007 |
res_pix <- 1200 |
|
| 1008 |
#res_pix <- 480 |
|
| 1009 |
|
|
| 1010 |
col_mfrow <- 1 |
|
| 1011 |
row_mfrow <- 1 |
|
| 1012 |
|
|
| 1013 |
png(filename=paste("Figure9_clim_mosaics_day_test","_",date_proc,"_",reg_name,"_",out_prefix,".png",sep=""),
|
|
| 1014 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 1015 |
|
|
| 1016 |
plot(r_pred,main=paste("climatology month ",date_proc , " ", reg_name,sep=""),cex.main=1.5)
|
|
| 1017 |
dev.off() |
|
| 1018 |
} |
|
| 1019 |
|
|
| 1020 | 949 |
##################### END OF SCRIPT ###################### |
Also available in: Unified diff
NEX assessment part2 adding function to plot daily mosaics in parrallel