Revision f9023de7
Added by Benoit Parmentier almost 9 years ago
| climate/research/oregon/interpolation/global_run_scalingup_mosaicing.R | ||
|---|---|---|
| 5 | 5 |
#Analyses, figures, tables and data are also produced in the script. |
| 6 | 6 |
#AUTHOR: Benoit Parmentier |
| 7 | 7 |
#CREATED ON: 04/14/2015 |
| 8 |
#MODIFIED ON: 04/07/2016
|
|
| 8 |
#MODIFIED ON: 04/08/2016
|
|
| 9 | 9 |
#Version: 6 |
| 10 | 10 |
#PROJECT: Environmental Layers project |
| 11 | 11 |
#COMMENTS: analyses run for reg4 1991 for test of mosaicing using 1500x4500km and other tiles |
| ... | ... | |
| 79 | 79 |
#27) algorithm: python or R, if R use mosaic function for R, if python use modified gdal merge, PARAM 27 |
| 80 | 80 |
#28) match_extent : if "FALSE" try without matching geographic extent #PARAM 28 |
| 81 | 81 |
#29) list_models : if NULL use y~1 formula #PARAM 29 |
| 82 |
|
|
| 82 |
#30) layers_option: mosaic to create as a layer from var_pred (e.g. TMax), res_training, res_testing, ac_testing |
|
| 83 |
|
|
| 83 | 84 |
###OUTPUT |
| 84 | 85 |
# |
| 85 | 86 |
# |
| ... | ... | |
| 167 | 168 |
#list_models <- paste(var_pred,"~","1",sep=" ") #if null then this is the default... |
| 168 | 169 |
layers_option <- list_param_run_mosaicing_prediction$layers_option |
| 169 | 170 |
|
| 171 |
|
|
| 172 |
################################################################# |
|
| 170 | 173 |
####### PART 1: Read in data and process data ######## |
| 174 |
######################################################## |
|
| 171 | 175 |
|
| 172 | 176 |
#out_dir <- in_dir #PARAM 11 |
| 173 | 177 |
#in_dir_tiles <- file.path(in_dir,"tiles") #this is valid both for Atlas and NEX |
| ... | ... | |
| 230 | 234 |
in_dir_tiles_tmp <- file.path(in_dir, region_name) |
| 231 | 235 |
#fix this later and add the year.. |
| 232 | 236 |
#gam_CAI_dailyTmax_predicted_mod1 |
| 233 |
|
|
| 237 |
#this is very slow!!! it takes 8 minutes?! |
|
| 234 | 238 |
lf_mosaic <- lapply(1:length(day_to_mosaic),FUN=function(i){list.files(path=file.path(in_dir_tiles_tmp),
|
| 235 | 239 |
pattern=paste("gam_CAI_dailyTmax_predicted_",pred_mod_name,".*.",day_to_mosaic[i],".*.tif$",sep=""),full.names=T,recursive=T)})
|
| 236 | 240 |
|
| 237 |
#reg_lf_mosaic[[k]] <- list.files(path=file.path(in_dir_tiles_tmp),pattern=paste(".*.",day_to_mosaic[i],".*.tif$",sep=""),full.names=T,recursive=T)
|
|
| 241 |
#########################################################################
|
|
| 238 | 242 |
##################### PART 2: produce the mosaic ################## |
| 243 |
###################################################################### |
|
| 239 | 244 |
|
| 240 | 245 |
#This is is assuming a list of file for a region!! |
| 241 | 246 |
#this is where the main function for mosaicing region starts!! |
| ... | ... | |
| 292 | 297 |
### Now create accuracy surfaces from residuals... |
| 293 | 298 |
|
| 294 | 299 |
if(layers_option=="res_testing"){
|
| 295 |
|
|
| 300 |
#This part took 19 minutes and 45 seconds |
|
| 301 |
|
|
| 296 | 302 |
## Create accuracy surface by kriging |
| 297 | 303 |
num_cores_tmp <-num_cores |
| 298 | 304 |
lf_day_tiles <- lf_mosaic #list of raster files by dates |
| ... | ... | |
| 357 | 363 |
##Run for data_day_s |
| 358 | 364 |
## |
| 359 | 365 |
|
| 360 |
if(layers_option=="res_testing"){
|
|
| 366 |
if(layers_option=="res_training"){
|
|
| 367 |
#This part took 19 minutes and 40 seconds |
|
| 361 | 368 |
|
| 362 | 369 |
data_df <- data_day_s # data.frame table/spdf containing stations with residuals and variable |
| 363 | 370 |
|
| ... | ... | |
| 413 | 420 |
lf_accuracy_residuals_data_s_raster <- lapply(1:length(list_create_accuracy_residuals_raster_obj),FUN=function(i,x){as.character(unlist(extract_from_list_obj(x[[i]]$list_pred_res_obj,"raster_name")))},x=list_create_accuracy_residuals_raster_obj)
|
| 414 | 421 |
|
| 415 | 422 |
} |
| 416 |
##took 31 minutes to generate the residuals maps for each tiles (28) for region 4 |
|
| 417 |
##Revised on 04/07 for three dates, it took 40 minutes |
|
| 418 | 423 |
|
| 419 | 424 |
###################################################### |
| 420 |
#### PART 2: GENERATE MOSAIC FOR LIST OF FILES #####
|
|
| 425 |
#### PART 3: GENERATE MOSAIC FOR LIST OF FILES #####
|
|
| 421 | 426 |
################################# |
| 422 | 427 |
#### Mosaic tiles for the variable predicted and accuracy metric |
| 423 | 428 |
|
| ... | ... | |
| 427 | 432 |
list_mosaic_obj <- vector("list",length=length(day_to_mosaic))
|
| 428 | 433 |
for(i in 1:length(day_to_mosaic)){
|
| 429 | 434 |
# |
| 430 |
mosaic_method <- "use_edge_weights" #this is distance from edge |
|
| 431 |
out_suffix_tmp <- paste(interpolation_method,y_var_name,day_to_mosaic[i],out_suffix,sep="_") |
|
| 432 |
#debug(mosaicFiles) |
|
| 433 |
#can also loop through methods!!! |
|
| 434 |
#python_bin <- "/usr/bin/" #python gdal bin, on Atlas NCEAS |
|
| 435 |
#python_bin <- "/nobackupp6/aguzman4/climateLayers/sharedModules/bin" #on NEX |
|
| 436 |
#gdal_merge_sum_noDataTest.py |
|
| 435 |
|
|
| 436 |
if(layers_option=="var_pred"){
|
|
| 437 |
|
|
| 438 |
mosaic_method <- "use_edge_weights" #this is distance from edge |
|
| 439 |
out_suffix_tmp <- paste(interpolation_method,y_var_name,day_to_mosaic[i],out_suffix,sep="_") |
|
| 440 |
#debug(mosaicFiles) |
|
| 441 |
#can also loop through methods!!! |
|
| 442 |
#python_bin <- "/usr/bin/" #python gdal bin, on Atlas NCEAS |
|
| 443 |
#python_bin <- "/nobackupp6/aguzman4/climateLayers/sharedModules/bin" #on NEX |
|
| 444 |
#gdal_merge_sum_noDataTest.py |
|
| 437 | 445 |
|
| 438 |
mosaic_edge_obj_prediction <- mosaicFiles(lf_mosaic[[i]],
|
|
| 446 |
mosaic_obj <- mosaicFiles(lf_mosaic[[i]],
|
|
| 439 | 447 |
mosaic_method="use_edge_weights", |
| 440 | 448 |
num_cores=num_cores, |
| 441 | 449 |
r_mask_raster_name=infile_mask, |
| ... | ... | |
| 448 | 456 |
file_format=file_format, |
| 449 | 457 |
out_suffix=out_suffix_tmp, |
| 450 | 458 |
out_dir=out_dir) |
| 451 |
#runs in 15-16 minutes for 3 dates and mosaicing of 28 tiles... |
|
| 459 |
#runs in 15-16 minutes for 3 dates and mosaicing of 28 tiles... |
|
| 460 |
list_mosaic_obj[[i]] <- mosaic_obj |
|
| 461 |
} |
|
| 452 | 462 |
|
| 453 |
## Now accuracy based on center of centroids |
|
| 454 |
mosaic_method <- "use_edge_weights" #this is distance from edge |
|
| 455 |
#Adding metric name in the name... |
|
| 456 |
out_suffix_tmp <- paste(interpolation_method,metric_name,day_to_mosaic[i],out_suffix,sep="_") |
|
| 457 |
|
|
| 458 |
#debug(mosaicFiles) |
|
| 459 |
#can also loop through methods!!! |
|
| 460 |
mosaic_edge_obj_accuracy <- mosaicFiles(lf_accuracy_raster[[i]], |
|
| 463 |
if(layers_option=="ac_testing"){
|
|
| 464 |
## Now accuracy based on center of centroids |
|
| 465 |
mosaic_method <- "use_edge_weights" #this is distance from edge |
|
| 466 |
#Adding metric name in the name... |
|
| 467 |
out_suffix_tmp <- paste(interpolation_method,metric_name,day_to_mosaic[i],out_suffix,sep="_") |
|
| 468 |
|
|
| 469 |
#debug(mosaicFiles) |
|
| 470 |
#can also loop through methods!!! |
|
| 471 |
mosaic_obj <- mosaicFiles(lf_accuracy_raster[[i]], |
|
| 461 | 472 |
mosaic_method="use_edge_weights", |
| 462 | 473 |
num_cores=num_cores, |
| 463 | 474 |
r_mask_raster_name=infile_mask, |
| ... | ... | |
| 469 | 480 |
file_format=file_format, |
| 470 | 481 |
out_suffix=out_suffix_tmp, |
| 471 | 482 |
out_dir=out_dir) |
| 472 |
##Took 39 minutes for 28 tiles and one date...!!! |
|
| 473 |
list_mosaic_obj[[i]] <- list(prediction=mosaic_edge_obj_prediction,accuracy=mosaic_edge_obj_accuracy) |
|
| 474 |
|
|
| 483 |
##Took 29 minutes for 28 tiles and one date...!!! |
|
| 484 |
list_mosaic_obj[[i]] <- mosaic_obj |
|
| 485 |
} |
|
| 486 |
|
|
| 487 |
#list_mosaic_obj[[i]] <- list(prediction=mosaic_edge_obj_prediction,accuracy=mosaic_edge_obj_accuracy) |
|
| 488 |
|
|
| 475 | 489 |
### produce residuals mosaics |
| 476 |
#for now add data_day_s in the name!! |
|
| 477 |
mosaic_method <- "use_edge_weights" #this is distance from edge |
|
| 478 |
out_suffix_tmp <- paste(interpolation_method,"kriged_residuals","data_day_s",day_to_mosaic[i],out_suffix,sep="_") |
|
| 479 |
#lf_tmp<-list.files(pattern="*kriged_residuals.*.tif",full.names=T) |
|
| 480 |
lf_tmp <- lf_accuracy_residuals_raster[[i]] |
|
| 481 |
#lf_accuracy_residuals_raster[[i]] |
|
| 482 |
#debug(mosaicFiles) |
|
| 483 |
mosaic_edge_obj_residuals <- mosaicFiles(lf_tmp, |
|
| 490 |
if(layers_option=="res_testing"){
|
|
| 491 |
#for now add data_day_s in the name!! |
|
| 492 |
mosaic_method <- "use_edge_weights" #this is distance from edge |
|
| 493 |
out_suffix_tmp <- paste(interpolation_method,"kriged_residuals","data_day_v",day_to_mosaic[i],out_suffix,sep="_") |
|
| 494 |
#lf_tmp<-list.files(pattern="*kriged_residuals.*.tif",full.names=T) |
|
| 495 |
lf_tmp <- lf_accuracy_residuals_raster[[i]] |
|
| 496 |
#lf_accuracy_residuals_raster[[i]] |
|
| 497 |
#debug(mosaicFiles) |
|
| 498 |
mosaic_obj <- mosaicFiles(lf_tmp, |
|
| 484 | 499 |
mosaic_method="use_edge_weights", |
| 485 | 500 |
num_cores=num_cores, |
| 486 | 501 |
r_mask_raster_name=infile_mask, |
| ... | ... | |
| 493 | 508 |
file_format=file_format, |
| 494 | 509 |
out_suffix=out_suffix_tmp, |
| 495 | 510 |
out_dir=out_dir) |
| 496 |
#Took 11 to 12 minues for one day and 28 tiles in region 4 |
|
| 497 |
list_mosaic_obj[[i]] <- list(prediction=mosaic_edge_obj_prediction,accuracy=mosaic_edge_obj_accuracy,mosaic_edge_obj_residuals) |
|
| 498 |
#} |
|
| 499 |
|
|
| 500 |
##End of mosaicing function for region predictions |
|
| 501 |
} |
|
| 502 |
|
|
| 503 |
|
|
| 504 |
##################### |
|
| 505 |
###### PART 2: Analysis and figures for the outputs of mosaic function ##### |
|
| 506 |
|
|
| 507 |
#### compute and aspect and slope with figures |
|
| 508 |
#list_lf_mosaic_obj <- vector("list",length(day_to_mosaic))
|
|
| 509 |
#lf_mean_mosaic <- vector("list",length(mosaicing_method))#2methods only
|
|
| 510 |
#l_method_mosaic <- vector("list",length(mosaicing_method))
|
|
| 511 |
#list_out_suffix <- vector("list",length(mosaicing_method))
|
|
| 512 |
|
|
| 513 |
#for(i in 1:length(day_to_mosaic)){
|
|
| 514 |
# list_lf_mosaic_obj[[i]] <- list.files(path=out_dir,pattern=paste("*",day_to_mosaic[i],
|
|
| 515 |
# "_.*.RData",sep="")) |
|
| 516 |
# lf_mean_mosaic[[i]] <- unlist(lapply(list_lf_mosaic_obj[[i]],function(x){load_obj(x)[["mean_mosaic"]]}))
|
|
| 517 |
# l_method_mosaic[[i]] <- paste(unlist(lapply(list_lf_mosaic_obj[[i]],function(x){load_obj(x)[["method"]]})),day_to_mosaic[i],sep="_")
|
|
| 518 |
# list_out_suffix[[i]] <- unlist(paste(l_method_mosaic[[i]],day_to_mosaic[[i]],out_suffix,sep="_")) |
|
| 519 |
#} |
|
| 520 |
|
|
| 521 |
#if(plot_figures==TRUE){
|
|
| 522 |
|
|
| 523 |
#list_param_plot_mosaic <- list(lf_mosaic=unlist(lf_mean_mosaic), |
|
| 524 |
# method=unlist(l_method_mosaic), |
|
| 525 |
# out_suffix=unlist(list_out_suffix)) |
|
| 511 |
#Took 11 to 12 minues for one day and 28 tiles in region 4 |
|
| 512 |
list_mosaic_obj[[i]] <- mosaic_obj |
|
| 513 |
} |
|
| 526 | 514 |
|
| 527 |
#plot and produce png movie... |
|
| 528 |
#plot_mosaic(1,list_param=list_param_plot_mosaic) |
|
| 529 |
#num_cores <- 1 |
|
| 530 |
#l_png_files <- mclapply(1:length(unlist(lf_mean_mosaic)),FUN=plot_mosaic, |
|
| 531 |
# list_param= list_param_plot_mosaic, |
|
| 532 |
# mc.preschedule=FALSE,mc.cores = num_cores) |
|
| 533 |
|
|
| 534 |
#lf_plot<- list.files(pattern="r_m_use.*.mask.*.tif$") |
|
| 535 |
#lf_mean_mosaic <- lf_plot |
|
| 536 |
|
|
| 537 |
|
|
| 538 |
#list_param_plot_mosaic <- list(lf_raster_fname=unlist(lf_mean_mosaic[1:2]), |
|
| 539 |
# screenRast=TRUE, |
|
| 540 |
# l_dates=day_to_mosaic, |
|
| 541 |
# out_dir_str=out_dir, |
|
| 542 |
# out_prefix_str <- "dailyTmax_", |
|
| 543 |
# out_suffix_str=out_suffix) |
|
| 544 |
#debug(plot_screen_raster_val) |
|
| 545 |
#plot_screen_raster_val(1,list_param_plot_mosaic) |
|
| 546 |
|
|
| 547 |
|
|
| 548 |
#num_cores <- 2 |
|
| 549 |
#l_png_files <- mclapply(1:length(unlist(lf_mean_mosaic)[1:2]),FUN=plot_screen_raster_val, |
|
| 550 |
# list_param= list_param_plot_mosaic, |
|
| 551 |
# mc.preschedule=FALSE,mc.cores = num_cores) |
|
| 552 |
|
|
| 553 |
|
|
| 554 |
#list_param_plot_mosaic <- list(lf_raster_fname=unlist(lf_mean_mosaic[4:6]), |
|
| 555 |
# screenRast=FALSE, |
|
| 556 |
# l_dates=day_to_mosaic, |
|
| 557 |
# out_dir_str=out_dir, |
|
| 558 |
# out_prefix_str <- paste("rmse_",out_suffix,sep=""),
|
|
| 559 |
# out_suffix_str=paste("rmse_",out_suffix,sep=""))
|
|
| 560 |
#num_cores <- 3 |
|
| 561 |
#l_png_files_rmse <- mclapply(1:length(unlist(lf_mean_mosaic)[4:6]),FUN=plot_screen_raster_val, |
|
| 562 |
# list_param= list_param_plot_mosaic, |
|
| 563 |
# mc.preschedule=FALSE,mc.cores = num_cores) |
|
| 515 |
### produce residuals mosaics |
|
| 516 |
if(layers_option=="res_training"){
|
|
| 517 |
#for now add data_day_s in the name!! |
|
| 518 |
mosaic_method <- "use_edge_weights" #this is distance from edge |
|
| 519 |
out_suffix_tmp <- paste(interpolation_method,"kriged_residuals","data_day_s",day_to_mosaic[i],out_suffix,sep="_") |
|
| 520 |
#lf_tmp<-list.files(pattern="*kriged_residuals.*.tif",full.names=T) |
|
| 521 |
lf_tmp <- lf_accuracy_residuals_raster[[i]] |
|
| 522 |
#lf_accuracy_residuals_raster[[i]] |
|
| 523 |
#debug(mosaicFiles) |
|
| 524 |
mosaic_obj <- mosaicFiles(lf_tmp, |
|
| 525 |
mosaic_method="use_edge_weights", |
|
| 526 |
num_cores=num_cores, |
|
| 527 |
r_mask_raster_name=infile_mask, |
|
| 528 |
python_bin=python_bin, |
|
| 529 |
mosaic_python=mosaic_python, |
|
| 530 |
algorithm=algorithm, |
|
| 531 |
match_extent=match_extent, |
|
| 532 |
df_points=NULL, |
|
| 533 |
NA_flag=NA_flag_val, |
|
| 534 |
file_format=file_format, |
|
| 535 |
out_suffix=out_suffix_tmp, |
|
| 536 |
out_dir=out_dir) |
|
| 537 |
list_mosaic_obj[[i]] <- mosaic_obj |
|
| 538 |
#Took 11 to 12 minues for one day and 28 tiles in region 4 |
|
| 539 |
} |
|
| 564 | 540 |
|
| 565 |
#} |
|
| 541 |
##End of mosaicing function for region predictions |
|
| 542 |
}## end of day_to_mosaic loop |
|
| 566 | 543 |
|
| 567 | 544 |
##Create return object |
| 568 |
#list of mosaiced files |
|
| 569 |
|
|
| 545 |
#list of mosaiced files: get the list of files now to include in the output object!! |
|
| 546 |
mosaicing_prediction_obj <- list(list_mosaic_obj,layer_option) |
|
| 547 |
names(mosaicing_prediction_obj) <- c("list_mosaic_obj","layer_option")
|
|
| 570 | 548 |
return(run_mosaicing_prediction_obj) |
| 549 |
|
|
| 571 | 550 |
} |
| 572 | 551 |
|
| 573 | 552 |
############### |
Also available in: Unified diff
mosaicing script, clean up of code and debugging for options