Merge pull request #5 from ec-jrc/ec_variables

Added variable varNumber to handle different variables position in ne…

gfit/Gfit2.sh

@@ -17,6 +17,28 @@ source $settFile
 eval echo "source $settFile"             #Import settings file
 set +o allexport
 
+#######################################################################
+# check default values for input variables
+if [ -z "$warmup_yrs" ]; then
+    echo "warmup_yrs is unset"
+    warmup_yrs=0
+fi
+
+if [ -z "$MaxGTavgDis" ]; then
+    echo "MaxGTavgDis is unset"
+    MaxGTavgDis=1
+fi
+
+if [ -z "$Nyears" ]; then
+    echo "Nyears is unset"
+    Nyears="NaN"
+fi
+if [ -z "$varNumber" ]; then
+    echo "varNumber is unset"
+    varNumber=1
+fi
+#######################################################################
+
 
 mkdir -p $outDir
 
@@ -46,6 +68,8 @@ cdo yearmax $outDir/${varName}2.nc $outDir/${varName}AMax.nc
 echo "Extract the average map"
 cdo timmean $outDir/${varName}2.nc $outDir/${varName}Avg.nc
 
+# Copy pcraster clone map
+echo "Copy pcraster clone map"
 cp $cloneMap $outDir
 ########################################################################
 # Fitting of Gumbel extreme value distributions with L-moments

gfit/gfit2.r

@@ -22,14 +22,14 @@
 ##    Jan 2018  - fixed bug in the threshold estimation in presence of NA and replaced for loops with parallel computing
 ##    Sep 2018  - fixed bug in the estimation of the gumbel parameters
 
-# USAGE: /usr/bin/Rscript gfit2.r <input.file> <output.file.path> <Return.Periods> <warmup_years> <N of years> <MaxGTavgDis>
+# USAGE: /usr/bin/Rscript gfit2.r <input.file> <output.file.path> <Return.Periods> <warmup_years> <N of years> <MaxGTavgDis> <VarNumber>
 # e.g. : /usr/bin/Rscript gfit2.r "/myInputPath/disAMax.nc" "/myOutputPath" 2-5-10-20-50-100-200-500 1 
 
 
 ######################################################
 ############   Input arguments   #####################
-
 args <- commandArgs(TRUE)
+
 input.file<-as.character(args[1])      #File with annual maxima (from cdo yearmax command)
 out.file.path<-as.character(args[2])   #Output folder
 rp1<-try(as.character(args[3]))        #Vector of output return periods
@@ -38,10 +38,10 @@ rp1<-try(as.character(args[3]))        #Vector of output return periods
 warmup_yrs<-try(as.numeric(args[4]))   #default=0         Warm-up years are excluded from the calculations of return levels. Normally this value is zero as Warmup years are removed already in the Gfit2.sh
 MaxGTavgDis<-try(as.numeric(args[5]))  #default=1 (Yes)   Option (1 or 0) to use only annual maxima larger than the long term average discharge. The latter must be put in out.file.path and named <whatever>Avg.nc (e.g., disAvg.nc)
 Nyears<-try(as.numeric(args[6]))       #default=all       Calculations of return levels is done on this number of years (normally 20-30 years)
+varNumber<-try(as.numeric(args[7]))    #default=1         Position of variable to extract in the netcdf file
 
 ptm <- proc.time() # used for estimating processing time
 dir.create(out.file.path, showWarnings = FALSE, recursive = TRUE)
-
 ####################################################################################################
 
 rp<-as.numeric(unlist(strsplit(rp1, split='-', fixed=TRUE)))
@@ -79,20 +79,25 @@ lmom.ub.LA <- function (x)
     return(z)
 }
 
+######################################
+#extract data from source nc file
+# varNumber<-2  #Change this if the variable to extract is in a different position in the netcdf file
+if(!exists("varNumber")){varNumber<-1} #default=1  Change this if the variable to extract is in a different position in the netcdf file
+######################################
 
 #extract data from source nc file
 nc = nc_open(filename = input.file,write=FALSE,suppress_dimvals=FALSE)	#open file with annual maxima from cdo
-xDim = nc$var[[1]]$dim[[1]]
-yDim = nc$var[[1]]$dim[[2]]
+xDim = nc$var[[varNumber]]$dim[[1]]
+yDim = nc$var[[varNumber]]$dim[[2]]
 data = ncvar_get(nc)
 
-Units<-nc$var[[1]]$units
-nCols<-nc$var[[1]]$size[1]   
-nRows<-nc$var[[1]]$size[2]
+Units<-nc$var[[varNumber]]$units
+nCols<-nc$var[[varNumber]]$size[1]   
+nRows<-nc$var[[varNumber]]$size[2]
 xllCenter<-min(nc$dim[[1]]$vals)
 yllCenter<-min(nc$dim[[2]]$vals)
 cellSize<-abs(nc$dim[[1]]$vals[1]-nc$dim[[1]]$vals[2])
-MV<-nc$var[[1]]$missval
+MV<-nc$var[[varNumber]]$missval
 
 Header<-list(paste0("ncols      ",nCols),paste0("nrows      ",nRows),paste0("xllcenter  ",xllCenter),paste0("yllcenter  ",yllCenter),paste0("cellsize   ",cellSize),paste0("NODATA_value  ",MV))
 

gfit/settingFile_test.sh

@@ -1,12 +1,13 @@
 #!/bin/bash
 # Settings file for the Gumbel fit
 
-inpDir="/climateRun4/HELIX/global/dis/LF05/out"   			#Input file with daily maps to analyze and from which to extract the annual maxima
-varName="dis"                                     			#Output filename 
+inpDir="/climateRun4/HELIX/global/dis/LF05/out"   			#Input path with daily maps to analyze and from which to extract the annual maxima
+varName="dis"                                     			#Input filename 
 outDir="/H07_Global/deleteMe"    							#Output directory
 cloneMap="/climateRun4/HELIX/global/maps/area.map"  		#Clone map (Only in case you want the output maps also in PCRaster format, otherwise omit)
 returnPeriods="2-5-10-20-50-100-200-500"              		#Return periods [years]. Use the dash (-) as separator
 rootdir="/nahaUsers/alfielo/CA/GloFAS/scripts/GumbelFit_v2/GFIT_tool.v2/" 	#Directory where gfit2.r is stored, with "/" at the end (omit if current directory) 
-warmup_yrs="0"   											#Number of years to discard in the beginning 
+warmup_yrs="0"   											#Number of years to discard in the beginning (omit if 0) 
 # Nyears="30"      										#Number of years to use in the fitting (omit if all)
 MaxGTavgDis=0												#if 1 then annual maxima which are below the long term average are ignored in the fitting (it may result in missing values in the return levels)
+varNumber=1										#Change this if the variable to extract is in a different position in the netcdf file (omit if 1)