Project

General

Profile

« Previous | Next » 

Revision c3bfbe4e

Added by Benoit Parmentier almost 9 years ago

stage6 testing assessment part2 figures and call from the shell

View differences:

climate/research/oregon/interpolation/master_script_stage_6.R
8 8 
#STAGE 3: Data preparation: meteorological station database query and extraction of covariates values from raster brick
9 9 
#STAGE 4: Raster prediction: run interpolation method (-- gam fusion, gam CAI, ...) and perform validation
10 10 
#STAGE 5: Output analyses: assessment of results for specific dates...
11 
#STAGE 6: Assessement of predictions by tiles and regions with mosaicing of predictions and accuracy
11 
#STAGE 6: Assessement of predictions by tiles and regions
12 
#STAGE 7: Mosaicing of predictions and accuracy layer productions
12 13 
#AUTHOR: Benoit Parmentier
13 14 
#CREATED ON: 12/29/2015
14 
#MODIFIED ON: 01/06/2016
15 
#MODIFIED ON: 01/22/2016
15 16 
#PROJECT: NCEAS-IPLANT-NASA: Environment Layers
16 17 

  		




  
  18
  
    
#First source these files:


  		




  
  19
  
    
#Resolved call issues from R.


  		




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


  		




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


  		




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


  		




  
  23
  
    

  		




  17
  24
  
    
## TODO:


  		




  18
  25
  
    
# Add  assessment part 2 (figures): still testing


  		




  19
  
  
    
# 


  		




  
  26
  
    
# Make this callable from the shell


  		




  
  27
  
    
# Adapt for python script


  		




  
  28
  
    
# Fix figure for part2


  		




  
  29
  
    
# Call also use library(optparse)


  		




  
  30
  
    

  		




  20
  31
  
    
##################################################################################################


  		




  
  32
  
    
#


  		




  
  33
  
    
### PARAMETERS DEFINED IN THE SCRIPT


  		




  
  34
  
    
#There are 21 parameters, 1 constant and 8 arguments (drawn from the parameters) for the Rscript call.


  		




  
  35
  
    
#The arguments are passed directly from Rscript:


  		




  
  36
  
    
#var <- args[1] # variable being interpolated #param 1, arg 1


  		




  
  37
  
    
#in_dir1 <- args[2] # This is the output directory containing global prediction e.g./nobackupp6/aguzman4/climateLayers/out/ param 5, arg 2


  		




  
  38
  
    
#region_name <- args[3] # region e.g. "reg4" param 6, arg 3


  		




  
  39
  
    
#out_prefix <- args[4] # this is used in creating an output directory,include region name? param 7, arg 4


  		




  
  40
  
    
#out_dir <- args[5] # output parent dir, can be home dir or other, param 8, arg 5)


  		




  
  41
  
    
#create_out_dir_param <- args[6] # if TRUE create a output from "output"+out_prefix param 9, arg 6


  		




  
  42
  
    
#list_year_predicted <- args[7] # enter as list but currently runs on the first element of the list, param 10, arg 7


  		




  
  43
  
    
#num_cores <- args[8] #number of cores used # param 13, arg 8


  		




  
  44
  
    
#max_mem <- args[9] # maximum memory, param 21


  		




  
  45
  
    

  		




  
  46
  
    
#var = "TMAX" # variable being interpolated #param 1, arg 1


  		




  
  47
  
    
#in_dir1 = "/nobackupp6/aguzman4/climateLayers/out/" #param 5, arg 2


  		




  
  48
  
    
#region_name = "reg4" #param 6, arg 3


  		




  
  49
  
    
#out_prefix = "run_global_analyses_pred_12282015" #param 7, arg 4


  		




  
  50
  
    
#out_dir = "/nobackupp8/bparmen1/" #param 8, arg 5


  		




  
  51
  
    
#create_out_dir_param = "TRUE" #param 9, arg 6


  		




  
  52
  
    
#list_year_predicted = c(2010) # param 10, arg 7


  		




  
  53
  
    
#num_cores = 6 #number of cores used # param 13, arg 8


  		




  
  54
  
    
#max_mem = 1e+07 #param 21, arg 9


  		




  
  55
  
    
  


  		




  
  56
  
    
#"Rscript /nobackupp8/bparmen1/env_layers_scripts/master_script_stage_6_01182016.R TMAX /nobackupp6/aguzman4/climateLayers/out/ reg4 run_global_analyses_pred_12282015 /nobackupp8/bparmen1/ TRUE 2010 6 1e+07


  		




  21
  57
  
    

  		




  22
  58
  
    
###Loading R library and packages  ou 


  		




  23
  59
  
    
library(RPostgreSQL)


  		




  ......




  41
  77
  
    
### Need to add documentation ###


  		




  42
  78
  
    

  		




  43
  79
  
    
#Adding command line arguments to use mpiexec


  		




  44
  
  
    
args<-commandArgs(TRUE)


  		




  
  80
  
    
args <- commandArgs(TRUE)


  		




  45
  81
  
    
#script_path<-"/nobackupp6/aguzman4/climateLayers/finalCode/environmental-layers/climate/research/oregon/interpolation"


  		




  46
  82
  
    
#dataHome<-"/nobackupp6/aguzman4/climateLayers/interp/testdata/"


  		




  47
  83
  
    
#script_path2<-"/nobackupp6/aguzman4/climateLayers/finalCode/environmental-layers/climate/research/world/interpolation"


  		




  
  84
  
    
#var <- args[1] # variable being interpolated #param 1, arg 1


  		




  
  85
  
    
#in_dir1 <- args[2] #param 5, arg 2


  		




  
  86
  
    
#region_name <- args[3] #param 6, arg 3


  		




  
  87
  
    
#out_prefix <- args[4] #param 7, arg 4


  		




  
  88
  
    
#out_dir <- args[5] #param 8, arg 5


  		




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


  		




  
  90
  
    
#create_out_dir_param <- args[6] #param 9, arg 6


  		




  
  91
  
    
#list_year_predicted <- args[7] # param 10, arg 7


  		




  
  92
  
    
#num_cores <- args[8] #number of cores used # param 13, arg 8


  		




  
  93
  
    
#max_mem <- args[9] #param 21


  		




  48
  94
  
    

  		




  49
  95
  
    
#CALLED FROM MASTER SCRIPT:


  		




  50
  96
  
    

  		




  ......




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


  		




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


  		




  60
  106
  
    

  		




  61
  
  
    
### Parameters and arguments ###


  		




  
  107
  
    
### Parameters, constants and arguments ###


  		




  62
  108
  
    
  


  		




  63
  
  
    
var<-"TMAX" # variable being interpolated


  		




  
  109
  
    
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #constant 1


  		




  
  110
  
    

  		




  
  111
  
    
#var<-"TMAX" # variable being interpolated #param 1, arg 1


  		




  
  112
  
    
var <- args[1] # variable being interpolated #param 1, arg 1


  		




  
  113
  
    

  		




  
  114
  
    
##Add for precip later...


  		




  64
  115
  
    
if (var == "TMAX") {


  		




  65
  116
  
    
  y_var_name <- "dailyTmax"


  		




  66
  117
  
    
  y_var_month <- "TMax"


  		




  ......




  71
  122
  
    
}


  		




  72
  123
  
    

  		




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


  		




  74
  
  
    
interpolation_method<-c("gam_CAI")


  		




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


  		




  
  125
  
    
interpolation_method<-c("gam_CAI") #param 2


  		




  
  126
  
    
CRS_interp <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0" #param 3


  		




  76
  127
  
    
#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";


  		




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


  		




  78
  128
  
    

  		




  79
  129
  
    
out_region_name<-""


  		




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


  		




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


  		




  81
  131
  
    

  		




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


  		




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


  		




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


  		




  85
  
  
    
#/nobackupp6/aguzman4/climateLayers/out_15x45/1982


  		




  
  134
  
    
#in_dir1 <- "/nobackupp6/aguzman4/climateLayers/out/" #param 5, arg 2


  		




  
  135
  
    
in_dir1 <- args[2] #param 5, arg 2


  		




  86
  136
  
    

  		




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


  		




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


  		




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


  		




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


  		




  91
  
  
    
out_prefix <- "run_global_analyses_pred_12282015" #PARAM5


  		




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


  		




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


  		




  94
  
  
    
create_out_dir_param <- TRUE #PARAM7


  		




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


  		




  
  139
  
    
#run assessment by region, this is a unique region only 


  		




  
  140
  
    
#region_name <- c("reg4") #param 6, arg 3


  		




  
  141
  
    
region_name <- args[3] #param 6, arg 3


  		




  95
  142
  
    

  		




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


  		




  97
  
  
    
#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";


  		




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


  		




  
  143
  
    
#out_prefix <- "run_global_analyses_pred_12282015" #param 7, arg 4


  		




  
  144
  
    
#out_dir <- "/nobackupp8/bparmen1/" #param 8, arg 5


  		




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


  		




  
  146
  
    
create_out_dir_param <- TRUE #param 9, arg 6


  		




  
  147
  
    
out_prefix <- args[4] #param 7, arg 4


  		




  
  148
  
    
out_dir <- args[5] #param 8, arg 5


  		




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


  		




  
  150
  
    
create_out_dir_param <- args[6] #param 9, arg 6


  		




  99
  151
  
    

  		




  100
  152
  
    
#list_year_predicted <- 1984:2004


  		




  101
  
  
    
list_year_predicted <- c("2014")


  		




  
  153
  
    
#list_year_predicted <- c("2014") # param 10, arg 7


  		




  102
  154
  
    
#year_predicted <- list_year_predicted[1]


  		




  
  155
  
    
list_year_predicted <- args[7] # param 10, arg 7


  		




  103
  156
  
    

  		




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


  		




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


  		




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


  		




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


  		




  
  157
  
    
file_format <- ".tif" #format for mosaiced files # param 11


  		




  
  158
  
    
NA_flag_val <- -9999  #No data value, # param 12


  		




  
  159
  
    
#num_cores <- 6 #number of cores used # param 13, arg 8


  		




  
  160
  
    
plotting_figures <- TRUE #running part2 of assessment to generate figures... # param 14


  		




  
  161
  
    
num_cores <- args[8] #number of cores used # param 13, arg 8


  		




  108
  162
  
    
  


  		




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


  		




  110
  
  
    
mosaic_plot <- FALSE #PARAM14


  		




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


  		




  112
  
  
    
multiple_region <- TRUE #PARAM16


  		




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


  		




  
  164
  
    
mosaic_plot <- FALSE #param 15


  		




  
  165
  
    
day_to_mosaic <- c("19920102","19920103","19920103") #param 16, not in use...


  		




  
  166
  
    
multiple_region <- TRUE #param 17


  		




  
  167
  
    
countries_shp <- "/nobackupp8/bparmen1/NEX_data/countries.shp" #param 18


  		




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


  		




  115
  
  
    
plot_region <- TRUE  #PARAM18


  		




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


  		




  
  169
  
    
plot_region <- TRUE  # param 19


  		




  
  170
  
    
threshold_missing_day <- c(367,365,300,200) # param 20


  		




  117
  171
  
    

  		




  118
  172
  
    
list_param_run_assessment_prediction <- list(in_dir1,region_name,y_var_name,interpolation_method,out_prefix,


  		




  119
  173
  
    
                                  out_dir,create_out_dir_param,CRS_locs_WGS84,list_year_predicted,


  		




  ......




  124
  178
  
    
                                  "file_format","NA_flag_val","num_cores","plotting_figures",


  		




  125
  179
  
    
                                  "mosaic_plot","day_to_mosaic","multiple_region","countries_shp","plot_region")


  		




  126
  180
  
    
                                      


  		




  127
  
  
    

  		




  128
  181
  
    
names(list_param_run_assessment_prediction)<-list_names


  		




  129
  182
  
    
  


  		




  130
  183
  
    
#max number of cells to read in memory


  		




  131
  
  
    
max_mem<-args[11]


  		




  
  184
  
    
max_mem <- args[9] #param 21


  		




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


  		




  133
  186
  
    

  		




  134
  
  
    

  		




  135
  187
  
    
#debug(run_assessment_prediction_fun)


  		




  136
  188
  
    
#debug(debug_fun_test)


  		




  137
  189
  
    
#debug_fun_test(list_param_raster_prediction)


  		




  ......




  163
  215
  
    
# LC11: Barren lands/sparse vegetation


  		




  164
  216
  
    
# LC12: Open water


  		




  165
  217
  
    

  		




  
  218
  
    
#"Rscript %s %s wgs84Grid %s %s %s %s/subset/mean_LST_%s_jan_6_wgs84.tif FALSE %s/%s/covar_obj_%s.RData %s/%s/%s/met_stations_outfiles_obj_gam_CAI_%s.RData 10 4800  %s %s > %s/outLogs/%s_stage4_%s.log 2>  %s/outLogs/%s_stage4_err_%s.log" % (scriptFn,ll,ll,outputFolder,b[0],outputFolder,ll,outputFolder,ll,ll,outputFolder,ll,year,ll,year,yearInt,outputFolder,ll,year,outputFolder,ll,year)


  		




  
  219
  
    
      #print outSt


  		




  
  220
  
    

  		












Also available in:  Unified diff