Revision 6a0591ae
Added by Benoit Parmentier over 8 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/11/2016
|
|
| 8 |
#MODIFIED ON: 04/20/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 |
| ... | ... | |
| 26 | 26 |
#setfacl -Rmd user:aguzman4:rwx /nobackupp8/bparmen1/output_run10_1500x4500_global_analyses_pred_1992_10052015 |
| 27 | 27 |
# |
| 28 | 28 |
#reg1 : North America |
| 29 |
#reg23" : Europe + Asia
|
|
| 30 |
#reg4" : South America
|
|
| 29 |
#reg23 : Europe + Asia
|
|
| 30 |
#reg4 : South America
|
|
| 31 | 31 |
#reg5 : Africa |
| 32 | 32 |
#reg6 : Oceania+ South East Asia |
| 33 | 33 |
# |
| ... | ... | |
| 46 | 46 |
|
| 47 | 47 |
run_mosaicing_prediction_fun <-function(i,list_param_run_mosaicing_prediction){
|
| 48 | 48 |
|
| 49 |
##Function to predict temperature interpolation with 12 input parameters
|
|
| 50 |
#12 parameters used in the data preparation stage and input in the current script
|
|
| 49 |
##This is a general function to mosaic predicted tiles and accuracy layers with 34 input parameters.
|
|
| 50 |
# |
|
| 51 | 51 |
# |
| 52 | 52 |
#Input Parameters: |
| 53 |
#1) in_dir <- list_param_run_mosaicing_prediction$in_dir #PARAM1
|
|
| 54 |
#2) y_var_name <- list_param_run_mosaicing_prediction$y_var_name # #PARAM2
|
|
| 55 |
#3) interpolation_method <- list_param_run_mosaicing_prediction$interpolation_method ##PARAM3
|
|
| 53 |
#1) in_dir: input directory, parent directory containing predictions for all regions #PARAM1
|
|
| 54 |
#2) y_var_name: climate variable predicted (dailyTmax, dailyTmin, dailyPrecip) # #PARAM2
|
|
| 55 |
#3) interpolation_method: interpolation method being used, gam_CAI currently ##PARAM3
|
|
| 56 | 56 |
#4) region_name: region_name e.g. reg4 South America #PARAM4 |
| 57 | 57 |
#5) mosaicing_method: options include "unweighted","use_edge_weights" #PARAM5 |
| 58 | 58 |
#6) out_suffix : output suffix #PARAM 6 |
| ... | ... | |
| 82 | 82 |
#30) list_models : if NULL use y~1 formula #PARAM 29 |
| 83 | 83 |
#31) layers_option: mosaic to create as a layer from var_pred (e.g. TMax), res_training, res_testing, ac_testing |
| 84 | 84 |
#32) tmp_files: if TRUE keep temporary files generated during mosaicing |
| 85 |
#33) use_int: if TRUE, use int32 in the final output |
|
| 86 |
#34) scaling: scaling factor to multiply the original variable before conversation to int |
|
| 85 | 87 |
|
| 86 | 88 |
###OUTPUT |
| 87 | 89 |
# |
| ... | ... | |
| 132 | 134 |
create_out_dir_param <- list_param_run_mosaicing_prediction$create_out_dir_param # FALSE #PARAM 12 |
| 133 | 135 |
|
| 134 | 136 |
#if daily mosaics NULL then mosaicas all days of the year #PARAM 13 |
| 135 |
day_to_mosaic_range <- list_param_run_mosaicing_prediction$day_to_mosaic_range # c("19920101","19920102","19920103") #,"19920104","19920105") #PARAM9, two dates note in /tiles for now on NEX
|
|
| 136 |
year_processed <- list_param_run_mosaicing_prediction$year_predicted |
|
| 137 |
day_to_mosaic_range <- list_param_run_mosaicing_prediction$day_to_mosaic_range # c("19920101","19920102","19920103") #,"19920104","19920105") #PARAM14, two dates note in /tiles for now on NEX
|
|
| 138 |
year_processed <- list_param_run_mosaicing_prediction$year_predicted #PARAM 15
|
|
| 137 | 139 |
#CRS_WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84 #CONSTANT1
|
| 138 | 140 |
#CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| 139 |
proj_str <- list_param_run_mosaicing_prediction$proj_str# CRS_WGS84 #PARAM 8 #check this parameter
|
|
| 141 |
proj_str <- list_param_run_mosaicing_prediction$proj_str# CRS_WGS84 #PARAM 16 #check this parameter
|
|
| 140 | 142 |
|
| 141 |
file_format <- list_param_run_mosaicing_prediction$file_format # ".tif" #PARAM 14
|
|
| 142 |
NA_value <- list_param_run_mosaicing_prediction$NA_value # -9999 #PARAM 15
|
|
| 143 |
file_format <- list_param_run_mosaicing_prediction$file_format # ".tif" #PARAM 17
|
|
| 144 |
NA_value <- list_param_run_mosaicing_prediction$NA_value # -9999 #PARAM 18
|
|
| 143 | 145 |
|
| 144 |
num_cores <- list_param_run_mosaicing_prediction$num_cores #6 #PARAM 17
|
|
| 146 |
num_cores <- list_param_run_mosaicing_prediction$num_cores #6 #PARAM 19
|
|
| 145 | 147 |
#region_names <- list_param_run_mosaicing_prediction$region_names # c("reg23","reg4") #selected region names, ##PARAM 18
|
| 146 |
use_autokrige <- list_param_run_mosaicing_prediction$use_autokrige # F #PARAM 19
|
|
| 148 |
use_autokrige <- list_param_run_mosaicing_prediction$use_autokrige # F #PARAM 20
|
|
| 147 | 149 |
|
| 148 |
###Separate folder for masks by regions, should be listed as just the dir!!... #PARAM 20
|
|
| 150 |
###Separate folder for masks by regions, should be listed as just the dir!!... #PARAM 21
|
|
| 149 | 151 |
#infile_mask <- "/nobackupp8/bparmen1/regions_input_files/r_mask_reg4.tif" |
| 150 | 152 |
infile_mask <- list_param_run_mosaicing_prediction$infile_mask # input mask used in defining the region |
| 151 | 153 |
|
| 152 | 154 |
#in_dir can be on NEX or Atlas |
| 153 | 155 |
|
| 154 | 156 |
##skip this for now |
| 155 |
#df_assessment_files_name <- list_param_run_mosaicing_prediction$df_assessment_files_name # data.frame with all files used in assessmnet, PARAM 21
|
|
| 157 |
df_assessment_files_name <- list_param_run_mosaicing_prediction$df_assessment_files_name # data.frame with all files used in assessmnet, PARAM 21 |
|
| 156 | 158 |
|
| 157 | 159 |
#python script and gdal on NEX NASA: |
| 158 | 160 |
#mosaic_python <- "/nobackupp6/aguzman4/climateLayers/sharedCode/" |
| ... | ... | |
| 169 | 171 |
#for residuals... |
| 170 | 172 |
list_models <- list_param_run_mosaicing_prediction$list_models # NULL #PARAM 26 |
| 171 | 173 |
#list_models <- paste(var_pred,"~","1",sep=" ") #if null then this is the default... |
| 172 |
layers_option <- list_param_run_mosaicing_prediction$layers_option |
|
| 173 |
tmp_files <- list_param_run_mosaicing_prediction$tmp_files |
|
| 174 |
layers_option <- list_param_run_mosaicing_prediction$layers_option #PARAM 27 |
|
| 175 |
tmp_files <- list_param_run_mosaicing_prediction$tmp_files #PARAM 28 |
|
| 176 |
use_int <- list_param_run_mosaicing_prediction$use_int #PARAM 29 |
|
| 177 |
scaling <- list_param_run_mosaicing_prediction$scaling |
|
| 174 | 178 |
|
| 175 | 179 |
################################################################# |
| 176 | 180 |
####### PART 1: Read in data and process data ######## |
| ... | ... | |
| 238 | 242 |
#fix this later and add the year.. |
| 239 | 243 |
#gam_CAI_dailyTmax_predicted_mod1 |
| 240 | 244 |
#this is very slow!!! it takes 8 minutes?! |
| 241 |
lf_mosaic <- lapply(1:length(day_to_mosaic),FUN=function(i){list.files(path=file.path(in_dir_tiles_tmp),
|
|
| 242 |
pattern=paste("gam_CAI_dailyTmax_predicted_",pred_mod_name,".*.",day_to_mosaic[i],".*.tif$",sep=""),full.names=T,recursive=T)})
|
|
| 245 |
#lf_mosaic <- lapply(1:length(day_to_mosaic),FUN=function(i){
|
|
| 246 |
# list.files(path=file.path(in_dir_tiles_tmp), |
|
| 247 |
# pattern=paste("gam_CAI_dailyTmax_predicted_",pred_mod_name,".*.",day_to_mosaic[i],".*.tif$",sep=""),
|
|
| 248 |
# full.names=T,recursive=T)}) |
|
| 243 | 249 |
|
| 250 |
#lf_mosaic <- lapply(1:length(day_to_mosaic),FUN=function(i){list.files(path=file.path(in_dir_tiles_tmp),
|
|
| 251 |
# pattern=paste("gam_CAI_dailyTmax_predicted_",pred_mod_name,".*.",day_to_mosaic[i],".*.tif$",sep=""),full.names=T,recursive=F)})
|
|
| 252 |
#Using changes from Alberto |
|
| 253 |
lf_mosaic <- lapply(1:length(day_to_mosaic),FUN=function(i){
|
|
| 254 |
searchStr = paste(in_dir_tiles_tmp,"/*/",year_processed,"/gam_CAI_dailyTmax_predicted_",pred_mod_name,"*",day_to_mosaic[i],"*.tif",sep="") |
|
| 255 |
#print(searchStr) |
|
| 256 |
Sys.glob(searchStr)}) |
|
| 257 |
|
|
| 258 |
browser() |
|
| 259 |
|
|
| 244 | 260 |
######################################################################### |
| 245 | 261 |
##################### PART 2: produce the mosaic ################## |
| 246 | 262 |
###################################################################### |
| ... | ... | |
| 255 | 271 |
|
| 256 | 272 |
####################################### |
| 257 | 273 |
################### PART I: Accuracy layers by tiles ############# |
| 258 |
#first generate accuracy layers using tiles definitions and output from the accuracy assessment |
|
| 274 |
###Depending on value of layers_option, create accuracy surfaces based on day testing metric (e.g. rmse)... |
|
| 275 |
#Generate accuracy layers using tiles definitions and output from the accuracy assessment |
|
| 259 | 276 |
|
| 260 | 277 |
if(layers_option=="ac_testing"){
|
| 261 | 278 |
#this takes about 1 minute and 35 seconds for 3 days and 28 tiles... |
| ... | ... | |
| 297 | 314 |
} |
| 298 | 315 |
|
| 299 | 316 |
#################################### |
| 300 |
### Now create accuracy surfaces from residuals...
|
|
| 317 |
###Depending on value of layers_option, create accuracy surfaces based on testing residuals...
|
|
| 301 | 318 |
|
| 302 | 319 |
if(layers_option=="res_testing"){
|
| 303 | 320 |
#This part took 19 minutes and 45 seconds |
| ... | ... | |
| 363 | 380 |
} |
| 364 | 381 |
|
| 365 | 382 |
######################################### |
| 366 |
##Run for data_day_s
|
|
| 383 |
###Depending on value of layers_option, create accuracy surfaces based on training residuals...
|
|
| 367 | 384 |
## |
| 368 | 385 |
|
| 369 | 386 |
if(layers_option=="res_training"){
|
| ... | ... | |
| 429 | 446 |
###################################################### |
| 430 | 447 |
#### PART 3: GENERATE MOSAIC FOR LIST OF FILES ##### |
| 431 | 448 |
################################# |
| 432 |
#### Mosaic tiles for the variable predicted and accuracy metric |
|
| 449 |
#### Mosaic tiles for the variable predicted and accuracy metrics, residuals surfaces or other options
|
|
| 433 | 450 |
|
| 434 | 451 |
#methods availbable:use_sine_weights,use_edge,use_linear_weights |
| 435 | 452 |
#only use edge method for now |
| 436 | 453 |
#loop to dates..., make this a function... |
| 454 |
#This is a loop but uses multicores when calling the mosaic function |
|
| 437 | 455 |
list_mosaic_obj <- vector("list",length=length(day_to_mosaic))
|
| 438 | 456 |
for(i in 1:length(day_to_mosaic)){
|
| 439 | 457 |
# |
| ... | ... | |
| 461 | 479 |
file_format=file_format, |
| 462 | 480 |
out_suffix=out_suffix_tmp, |
| 463 | 481 |
out_dir=out_dir, |
| 464 |
tmp_files=tmp_files) |
|
| 482 |
tmp_files=tmp_files, |
|
| 483 |
use_int=use_int, |
|
| 484 |
scaling=scaling) |
|
| 465 | 485 |
#runs in 15-16 minutes for 3 dates and mosaicing of 28 tiles... |
| 466 | 486 |
list_mosaic_obj[[i]] <- mosaic_obj |
| 467 | 487 |
} |
| ... | ... | |
| 486 | 506 |
file_format=file_format, |
| 487 | 507 |
out_suffix=out_suffix_tmp, |
| 488 | 508 |
out_dir=out_dir, |
| 489 |
tmp_files=tmp_files) |
|
| 509 |
tmp_files=tmp_files, |
|
| 510 |
use_int=use_int, |
|
| 511 |
scaling=scaling) |
|
| 490 | 512 |
##Took 29 minutes for 28 tiles and one date...!!! |
| 491 | 513 |
list_mosaic_obj[[i]] <- mosaic_obj |
| 492 | 514 |
} |
| ... | ... | |
| 503 | 525 |
#lf_accuracy_residuals_raster[[i]] |
| 504 | 526 |
#debug(mosaicFiles) |
| 505 | 527 |
mosaic_obj <- mosaicFiles(lf_tmp, |
| 506 |
mosaic_method="use_edge_weights", |
|
| 507 |
num_cores=num_cores, |
|
| 508 |
r_mask_raster_name=infile_mask, |
|
| 509 |
python_bin=python_bin, |
|
| 510 |
mosaic_python=mosaic_python, |
|
| 511 |
algorithm=algorithm, |
|
| 512 |
match_extent=match_extent, |
|
| 513 |
df_points=NULL, |
|
| 514 |
NA_flag=NA_flag_val, |
|
| 515 |
file_format=file_format, |
|
| 516 |
out_suffix=out_suffix_tmp, |
|
| 517 |
out_dir=out_dir, |
|
| 518 |
tmp_files=tmp_files) |
|
| 528 |
mosaic_method="use_edge_weights", |
|
| 529 |
num_cores=num_cores, |
|
| 530 |
r_mask_raster_name=infile_mask, |
|
| 531 |
python_bin=python_bin, |
|
| 532 |
mosaic_python=mosaic_python, |
|
| 533 |
algorithm=algorithm, |
|
| 534 |
match_extent=match_extent, |
|
| 535 |
df_points=NULL, |
|
| 536 |
NA_flag=NA_flag_val, |
|
| 537 |
file_format=file_format, |
|
| 538 |
out_suffix=out_suffix_tmp, |
|
| 539 |
out_dir=out_dir, |
|
| 540 |
tmp_files=tmp_files, |
|
| 541 |
use_int, |
|
| 542 |
scaling) |
|
| 519 | 543 |
#Took 11 to 12 minues for one day and 28 tiles in region 4 |
| 520 | 544 |
list_mosaic_obj[[i]] <- mosaic_obj |
| 521 | 545 |
} |
| ... | ... | |
| 530 | 554 |
#lf_accuracy_residuals_raster[[i]] |
| 531 | 555 |
#debug(mosaicFiles) |
| 532 | 556 |
mosaic_obj <- mosaicFiles(lf_tmp, |
| 533 |
mosaic_method="use_edge_weights", |
|
| 534 |
num_cores=num_cores, |
|
| 535 |
r_mask_raster_name=infile_mask, |
|
| 536 |
python_bin=python_bin, |
|
| 537 |
mosaic_python=mosaic_python, |
|
| 538 |
algorithm=algorithm, |
|
| 539 |
match_extent=match_extent, |
|
| 540 |
df_points=NULL, |
|
| 541 |
NA_flag=NA_flag_val, |
|
| 542 |
file_format=file_format, |
|
| 543 |
out_suffix=out_suffix_tmp, |
|
| 544 |
out_dir=out_dir, |
|
| 545 |
tmp_files=tmp_files) |
|
| 557 |
mosaic_method="use_edge_weights", |
|
| 558 |
num_cores=num_cores, |
|
| 559 |
r_mask_raster_name=infile_mask, |
|
| 560 |
python_bin=python_bin, |
|
| 561 |
mosaic_python=mosaic_python, |
|
| 562 |
algorithm=algorithm, |
|
| 563 |
match_extent=match_extent, |
|
| 564 |
df_points=NULL, |
|
| 565 |
NA_flag=NA_flag_val, |
|
| 566 |
file_format=file_format, |
|
| 567 |
out_suffix=out_suffix_tmp, |
|
| 568 |
out_dir=out_dir, |
|
| 569 |
tmp_files=tmp_files, |
|
| 570 |
use_int=use_int, |
|
| 571 |
scaling=scaling) |
|
| 572 |
|
|
| 546 | 573 |
list_mosaic_obj[[i]] <- mosaic_obj |
| 547 | 574 |
#Took 11 to 12 minues for one day and 28 tiles in region 4 |
| 548 | 575 |
} |
Also available in: Unified diff
adding use int and scaling option to function mosicing script