Project

General

Profile

« Previous | Next » 

Revision 09b23c86

Added by Benoit Parmentier over 8 years ago

stage 7 mosaicing following assessment, script cleaned and integrated to workflow

View differences:

climate/research/oregon/interpolation/master_script_stage_7.R
10 10 
#STAGE 5: Output analyses: assessment of results for specific dates and tiles
11 11 
#STAGE 6: Assessement of predictions by tiles and regions
12 12 
#STAGE 7: Mosaicing of predicted surfaces and accuracy metrics (RMSE,MAE) by regions
13 
#STAGE 8: Comparison of predictions across regions and years with figures generation.
14 

  		




  13
  15
  
    
#AUTHOR: Benoit Parmentier                                                                        


  		




  14
  
  
    
#CREATED ON: 01/01/2015  


  		




  15
  
  
    
#MODIFIED ON: 01/05/2016  


  		




  
  16
  
    
#CREATED ON: 01/01/2016  


  		




  
  17
  
    
#MODIFIED ON: 04/05/2016  


  		




  16
  18
  
    
#PROJECT: NCEAS INPLANT: Environment and Organisms                                                                           


  		




  17
  19
  
    

  		




  
  20
  
    
#First source these files:


  		




  
  21
  
    
#Resolved call issues from R.


  		




  
  22
  
    
#source /nobackupp6/aguzman4/climateLayers/sharedModules2/etc/environ.sh  


  		




  
  23
  
    
#MODULEPATH=$MODULEPATH:/nex/modules/files


  		




  
  24
  
    
#module load pythonkits/gdal_1.10.0_python_2.7.3_nex


  		




  
  25
  
    

  		




  18
  26
  
    
## TODO:


  		




  19
  27
  
    
# 


  		




  20
  28
  
    
# Clean up temporary files


  		




  ......




  50
  58
  
    

  		




  51
  59
  
    
#CALLED FROM MASTER SCRIPT:


  		




  52
  60
  
    

  		




  
  61
  
    
#script_path <- "/home/parmentier/Data/IPLANT_project/env_layers_scripts"


  		




  53
  62
  
    
script_path <- "/nobackupp8/bparmen1/env_layers_scripts" #path to script


  		




  54
  63
  
    
function_mosaicing_functions <- "global_run_scalingup_mosaicing_function_12192015.R" #PARAM12


  		




  55
  
  
    
function_mosaicing <-"global_run_scalingup_mosaicing_01012016.Rglobal_run_scalingup_mosaicing_01052016.R"


  		




  56
  
  
    
source(file.path(function_mosaicing)) #source all functions used in this script 


  		




  57
  
  
    
source(file.path(function_mosaicing_functions)) #source all functions used in this script 


  		




  
  64
  
    
function_mosaicing <-"global_run_scalingup_mosaicing_04052016.R"


  		




  
  65
  
    
source(file.path(script_path,function_mosaicing)) #source all functions used in this script 


  		




  
  66
  
    
source(file.path(script_path,function_mosaicing_functions)) #source all functions used in this script 


  		




  
  67
  
    

  		




  
  68
  
    
function_assessment_part1_functions <- "global_run_scalingup_assessment_part1_functions_02112015.R" #PARAM12


  		




  
  69
  
    
function_assessment_part1a <-"global_run_scalingup_assessment_part1a_01042016.R"


  		




  
  70
  
    
function_assessment_part2 <- "global_run_scalingup_assessment_part2_02092016.R"


  		




  
  71
  
    
function_assessment_part2_functions <- "global_run_scalingup_assessment_part2_functions_01032016.R"


  		




  
  72
  
    
function_assessment_part3 <- "global_run_scalingup_assessment_part3_02102016.R"


  		




  
  73
  
    
source(file.path(script_path,function_assessment_part1_functions)) #source all functions used in this script 


  		




  
  74
  
    
source(file.path(script_path,function_assessment_part1a)) #source all functions used in this script 


  		




  
  75
  
    
source(file.path(script_path,function_assessment_part2)) #source all functions used in this script 


  		




  
  76
  
    
source(file.path(script_path,function_assessment_part2_functions)) #source all functions used in this script 


  		




  
  77
  
    
source(file.path(script_path,function_assessment_part3)) #source all functions used in this script 


  		




  58
  78
  
    

  		




  59
  79
  
    
### Parameters and arguments ###


  		




  60
  80
  
    
  


  		




  ......




  68
  88
  
    
  y_var_month <- "TMin"


  		




  69
  89
  
    
}


  		




  70
  90
  
    

  		




  71
  
  
    
#interpolation_method<-c("gam_fusion") #other otpions to be added later


  		




  72
  
  
    
interpolation_method<-c("gam_CAI")


  		




  73
  
  
    
CRS_interp<-"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs";


  		




  74
  
  
    
#CRS_interp <-"+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs";


  		




  75
  91
  
    
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84


  		




  76
  92
  
    

  		




  77
  
  
    
out_region_name<-""


  		




  78
  
  
    
list_models<-c("y_var ~ s(lat,lon,k=5) + s(elev_s,k=3) + s(LST,k=3)")


  		




  79
  
  
    

  		




  80
  
  
    
#reg1 (North Am), reg2(Europe),reg3(Asia), reg4 (South Am), reg5 (Africa), reg6 (Australia-Asia)


  		




  81
  
  
    
#master directory containing the definition of tile size and tiles predicted


  		




  82
  
  
    
in_dir1 <- "/nobackupp6/aguzman4/climateLayers/out/"


  		




  83
  
  
    
#/nobackupp6/aguzman4/climateLayers/out_15x45/1982


  		




  84
  
  
    

  		




  85
  
  
    
#region_names <- c("reg23","reg4") #selected region names, #PARAM2


  		




  86
  
  
    
region_name <- c("reg4") #run assessment by region, this is a unique region only


  		




  87
  
  
    
#region_names <- c("reg1","reg2","reg3","reg4","reg5","reg6") #selected region names, #PARAM2


  		




  88
  
  
    
interpolation_method <- c("gam_CAI") #PARAM4


  		




  89
  
  
    
out_prefix <- "run_global_analyses_pred_12282015" #PARAM5


  		




  90
  
  
    
out_dir <- "/nobackupp8/bparmen1/" #PARAM6


  		




  91
  
  
    
#out_dir <-paste(out_dir,"_",out_prefix,sep="")


  		




  92
  
  
    
create_out_dir_param <- TRUE #PARAM7


  		




  93
  
  
    

  		




  94
  
  
    
#CRS_interp<-"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs";


  		




  95
  
  
    
#CRS_interp <-"+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs";


  		




  96
  
  
    
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84


  		




  97
  
  
    

  		




  98
  
  
    
#list_year_predicted <- 1984:2004


  		




  99
  
  
    
list_year_predicted <- c("2014")


  		




  100
  
  
    
#year_predicted <- list_year_predicted[1]


  		




  101
  
  
    

  		




  102
  
  
    
file_format <- ".tif" #format for mosaiced files #PARAM10


  		




  103
  
  
    
NA_flag_val <- -9999  #No data value, #PARAM11


  		




  104
  
  
    
num_cores <- 6 #number of cores used #PARAM13


  		




  105
  
  
    
plotting_figures <- TRUE #running part2 of assessment to generate figures...


  		




  106
  
  
    
  


  		




  107
  
  
    
##Additional parameters used in part 2, some these may be removed as code is simplified


  		




  108
  
  
    
mosaic_plot <- FALSE #PARAM14


  		




  109
  
  
    
day_to_mosaic <- c("19920102","19920103","19920103") #PARAM15


  		




  110
  
  
    
multiple_region <- TRUE #PARAM16


  		




  111
  
  
    
countries_shp <- "/nobackupp8/bparmen1/NEX_data/countries.shp" #PARAM17


  		




  112
  
  
    
#countries_shp <-"/data/project/layers/commons/NEX_data/countries.shp" #Atlas


  		




  113
  
  
    
plot_region <- TRUE  #PARAM18


  		




  114
  
  
    
threshold_missing_day <- c(367,365,300,200)#PARAM19


  		




  115
  
  
    

  		




  116
  
  
    
list_param_run_assessment_prediction <- list(in_dir1,region_name,y_var_name,interpolation_method,out_prefix,


  		




  117
  
  
    
                                  out_dir,create_out_dir_param,CRS_locs_WGS84,list_year_predicted,


  		




  118
  
  
    
                                  file_format,NA_flag_val,num_cores,plotting_figures,


  		




  119
  
  
    
                                  mosaic_plot,day_to_mosaic,multiple_region,countries_shp,plot_region)


  		




  120
  
  
    
list_names <- c("in_dir1","region_name","y_var_name","interpolation_method","out_prefix",


  		




  121
  
  
    
                                  "out_dir","create_out_dir_param","CRS_locs_WGS84","list_year_predicted",


  		




  122
  
  
    
                                  "file_format","NA_flag_val","num_cores","plotting_figures",


  		




  123
  
  
    
                                  "mosaic_plot","day_to_mosaic","multiple_region","countries_shp","plot_region")


  		




  124
  
  
    
                                      


  		




  125
  
  
    

  		




  126
  
  
    
names(list_param_run_assessment_prediction)<-list_names


  		




  
  93
  
    
#####mosaicing parameters


  		




  
  94
  
    

  		




  
  95
  
    
#Data is on ATLAS or NASA NEX


  		




  
  96
  
    
#PARAM 1


  		




  
  97
  
    
#in_dir <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_pred_1992_12072015" #NCEAS


  		




  
  98
  
    
#in_dir <- "/nobackupp8/bparmen1/output_run10_1500x4500_global_analyses_pred_1992_12072015" #NEX


  		




  
  99
  
    
in_dir <- "/nobackupp6/aguzman4/climateLayers/out/"


  		




  
  100
  
    

  		




  
  101
  
    
y_var_name <- "dailyTmax" #PARAM2


  		




  
  102
  
    
interpolation_method <- c("gam_CAI") #PARAM3


  		




  
  103
  
    
region_name <- "reg4" #PARAM 4 #reg4 South America, Africa reg5,Europe reg2, North America reg1, Asia reg3


  		




  
  104
  
    
mosaicing_method <- c("unweighted","use_edge_weights") #PARAM5


  		




  
  105
  
    
#out_suffix <- paste(region_name,"_","run10_1500x4500_global_analyses_pred_1991_04052016",sep="") #PARAM 6


  		




  
  106
  
    
#out_suffix_str <- "run10_1500x4500_global_analyses_pred_1991_04052016" #PARAM 7


  		




  
  107
  
    

  		




  
  108
  
    
out_suffix <- region_name #PARAM 6


  		




  
  109
  
    
out_suffix_str <- region_name #PARAM 7


  		




  
  110
  
    

  		




  
  111
  
    
metric_name <- "rmse" #RMSE, MAE etc. #PARAM 8


  		




  
  112
  
    
pred_mod_name <- "mod1" #PARAM 9


  		




  
  113
  
    
var_pred <- "res_mod1" #used in residuals mapping #PARAM 10


  		




  
  114
  
    

  		




  
  115
  
    
#out_dir <- in_dir #PARAM 11


  		




  
  116
  
    
out_dir <- "/nobackupp8/bparmen1/climateLayers/out/reg4" #PARAM 11, use this location for now


  		




  
  117
  
    
create_out_dir_param <- TRUE #PARAM 12


  		




  
  118
  
    

  		




  
  119
  
    
#if daily mosaics NULL then mosaicas all days of the year #PARAM 13


  		




  
  120
  
    
day_to_mosaic <- c("19910101","19910102","19910103") #,"19920104","19920105") #PARAM9, two dates note in /tiles for now on NEX


  		




  
  121
  
    

  		




  
  122
  
    
#CRS_WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84 #CONSTANT1


  		




  
  123
  
    
#CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84


  		




  
  124
  
    
#proj_str<- CRS_WGS84 #PARAM 8 #check this parameter


  		




  
  125
  
    
 


  		




  
  126
  
    
file_format <- ".tif" #PARAM 14


  		




  
  127
  
    
NA_value <- -9999 #PARAM 15


  		




  
  128
  
    
NA_flag_val <- NA_value #PARAM 16


  		




  
  129
  
    
     


  		




  
  130
  
    
num_cores <- 6 #PARAM 17                 


  		




  
  131
  
    
region_names <- c("reg23","reg4") #selected region names, ##PARAM 18 


  		




  
  132
  
    
use_autokrige <- F #PARAM 19


  		




  
  133
  
    
proj_str <- CRS_locs_WGS84


  		




  
  134
  
    

  		




  
  135
  
    
###Separate folder for masks by regions, should be listed as just the dir!!... #PARAM 20


  		




  
  136
  
    
infile_mask <- "/nobackupp8/bparmen1/NEX_data/regions_input_files/r_mask_reg4.tif"


  		




  
  137
  
    
#infile_mask <- "/data/project/layers/commons/NEX_data/regions_input_files/r_mask_reg4.tif"


  		




  
  138
  
    
## All of this is interesting so use df_assessment!!


  		




  
  139
  
    

  		




  
  140
  
    
year_processed <- 1991 #PARAM 31


  		




  
  141
  
    
#path_assessment <- "/nobackupp6/aguzman4/climateLayers/out/reg4/assessment/output_reg4_1991"


  		




  
  142
  
    
#path_assessment <- file.path(in_dir,region_name,"assessment",paste("output_",region_name,year_processed,sep=""))


  		




  
  143
  
    
df_assessment_files_name <- "/nobackupp6/aguzman4/climateLayers/out/reg4/assessment/output_reg4_1991/df_assessment_files_reg4_1991_reg4_1991.txt" # data.frame with all files used in assessmnet, PARAM 21


  		




  
  144
  
    

  		




  
  145
  
    
#in_dir can be on NEX or Atlas


  		




  
  146
  
    

  		




  
  147
  
    
#python script and gdal on NEX NASA:


  		




  
  148
  
    
mosaic_python <- "/nobackupp6/aguzman4/climateLayers/sharedCode/"


  		




  
  149
  
    
python_bin <- "/nobackupp6/aguzman4/climateLayers/sharedModules2/bin"


  		




  
  150
  
    
#python script and gdal on Atlas NCEAS


  		




  
  151
  
    
#mosaic_python <- "/data/project/layers/commons/NEX_data/sharedCode" #PARAM 26


  		




  
  152
  
    
#python_bin <- "/usr/bin" #PARAM 27


  		




  
  153
  
    

  		




  
  154
  
    
algorithm <- "python" #PARAM 28 #if R use mosaic function for R, if python use modified gdalmerge script from Alberto Guzmann


  		




  
  155
  
    
#algorithm <- "R" #if R use mosaic function for R, if python use modified gdalmerge script from Alberto Guzmann


  		




  
  156
  
    
match_extent <- "FALSE" #PARAM 29 #try without matching!!!


  		




  
  157
  
    

  		




  
  158
  
    
#for residuals...


  		




  
  159
  
    
list_models <- NULL #PARAM 30


  		




  
  160
  
    
#list_models <- paste(var_pred,"~","1",sep=" ") #if null then this is the default...


  		




  127
  161
  
    
  


  		




  128
  162
  
    
#max number of cells to read in memory


  		




  129
  163
  
    
max_mem<-args[11]


  		




  130
  
  
    
#rasterOptions(maxmemory=1e+07,timer=TRUE)


  		




  
  164
  
    
#in_dir_tiles <- file.path(in_dir,"tiles") #this is valid both for Atlas and NEX


  		




  131
  165
  
    

  		




  132
  
  
    
#debug(run_assessment_prediction_fun)


  		




  
  166
  
    
#rasterOptions(maxmemory=1e+07,timer=TRUE)


  		




  
  167
  
    
list_param_run_mosaicing_prediction <- list(in_dir,y_var_name,interpolation_method,region_name,


  		




  
  168
  
    
                 mosaicing_method,out_suffix,out_suffix_str,metric_name,pred_mod_name,var_pred,


  		




  
  169
  
    
                 create_out_dir_param,day_to_mosaic,proj_str,file_format,NA_value,num_cores,


  		




  
  170
  
    
                 region_name,use_autokrige,infile_mask,df_assessment_files_name,mosaic_python,


  		




  
  171
  
    
                 python_bin,algorithm,match_extent,list_models)


  		




  
  172
  
    
param_names <- c("in_dir","y_var_name","interpolation_method","region_name",


  		




  
  173
  
    
                 "mosaicing_method","out_suffix","out_suffix_str","metric_name","pred_mod_name","var_pred",


  		




  
  174
  
    
                 "create_out_dir_param","day_to_mosaic","proj_str","file_format","NA_value","num_cores",


  		




  
  175
  
    
                 "region_name","use_autokrige","infile_mask","df_assessment_files_name","mosaic_python",


  		




  
  176
  
    
                 "python_bin","algorithm","match_extent","list_models")


  		




  
  177
  
    
names(list_param_run_mosaicing_prediction) <- param_names


  		




  
  178
  
    
#list_param_run_mosaicing_prediction


  		




  
  179
  
    
#debug(run_mosaicing_prediction_fun)


  		




  133
  180
  
    
#debug(debug_fun_test)


  		




  134
  181
  
    
#debug_fun_test(list_param_raster_prediction)


  		




  135
  182
  
    
i <- 1 #this select the first year of list_year_predicted


  		




  136
  
  
    
if (stages_to_run[6]==6){


  		




  137
  
  
    
  assessment_prediction_obj <- run_assessment_prediction_fun(i,list_param_run_assessment_prediction)


  		




  
  183
  
    
if (stages_to_run[7]==7){


  		




  
  184
  
    
  assessment_prediction_obj <- run_mosaicing_prediction_fun(i,list_param_run_mosaicing_prediction)


  		




  
  185
  
    
  #list_param_run_mosaicing_prediction


  		




  138
  186
  
    
}


  		




  139
  187
  
    

  		




  140
  
  
    
## Add stage 7 (mosaicing) here??


  		




  141
  
  
    
#i <- 1 #this select the first year of list_year_predicted


  		




  142
  
  
    
#if (stages_to_run[7]==7){


  		




  143
  
  
    
#  assessment_prediction_obj <- run_assessment_prediction_fun(i,list_param_run_assessment_prediction)


  		




  144
  
  
    
#}


  		




  145
  
  
    

  		




  146
  188
  
    
###############   END OF SCRIPT   ###################


  		




  147
  189
  
    
#####################################################


  		




  148
  190
  
    

  		












Also available in:  Unified diff