Repository: ec-jrc/lisflood-utilities (language: python)

Issue:
Fix compare tool, lat/lon variable, updated readme.



Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index dc2aa40..ab2049b 100644
--- a/README.md
+++ b/README.md
@@ -37,6 +37,8 @@ Here's a list of utilities you can find in lisflood-utilities package.
   - an existing boolean area mask
   - a list of stations and a LDD (in netCDF or PCRaster format) **Note: PCRaster must be installed in the conda env**
  
+* __thresholds__ is a tool to compute the discharge return period thresholds from netCDF4 file containing a discharge time series.
+
 * __compare__ is a package containing a set of simple Python classes that helps to compare 
 netCDF, PCRaster and TSS files.
 
@@ -323,45 +325,30 @@ optional arguments:
                         threshold for large diffs percentage
 ```
 
-## gfit
-
-### Introduction
-
-The tool Gfit2 was created to extract flood warning thresholds from a map stack of daily discharge for several years, given in NetCDF format.
-It is designed to work with any other variable as well as with different sampling (sub- or super- daily)
-
-The tool is made of three files:
-
-1. Gfit2.sh: the main file, used to run the tool (e.g., ./Gfit2.sh ./settingFile_test.sh)
-2. settingFile_test.sh: a setting file, where paths and values of the input variables are defined. It can be renamed, as long as the correct name is called in the script (e.g.: set1.sh --> ./Gfit2.sh ./set1.sh)
-3. gfit2.r: the r script that  performs the extreme value analysis
-
-### Requirements
+## thresholds
 
-You need to have installed the following:
+The thresholds tool computes the discharge return period thresholds using the method of L-moments.
+It is used to post-process the discharge from the LISFLOOD long term run.
+The resulting thresholds can be used in a flood forecasting system to define the flood warning levels.
 
-- CDO
-- R
-
-### How to run
-The tool operates in a Linux environment (also as a job with qsub). To run the tool you'll need:
-- R (the script was tested with R version 3.5.0)
-- the following R packages: ncdf4, lmomco, ismev
-- CDO (the script was tested with CDO version 1.6.5.1)
-
-### What does the tool do
+### Usage:
+The tool takes as input a Netcdf file containing the annual maxima of the discharge signal. LISFLOOD computes time series of discharge values (average value over the selected computational time step). The users are therefore required to compute the annual maxima. As an example, this step can be achieved by using CDO (cdo yearmax), for all the details please refer to [https://code.mpimet.mpg.de/projects/cdo/embedded/index.html#x1-190001.2.5](https://code.mpimet.mpg.de/projects/cdo/embedded/index.html#x1-190001.2.5)
 
-The Gfit2 tool performs the following steps:
-- It takes the input file in NetCDF and extract the series of annual maxima. an optional number of warm-up years ($warmup_yrs) are excluded at the beginning of the data, to remove potential spin-up effect. Also the map of average value is computed
-- An extreme value fitting is performed on each pixel of the map, using L-moments and a 2-parameter Gumbel distribution. As option, the user can limit the analysis to a specific number of years. Also, one can use the option to remove from the fitting all values smaller than the long-term average. The fitting is performed only where at least 5 annual maxima are available, otherwise a NA is returned on the output map
-- Output return level maps (corresponding to user-selected years of recurrence interval) are saved in a netcdf file (return_levels.nc) and as ascii files. If a clone map in PCRaster format is provided, maps are also saved in PCRaster (return level maps and parameters of the Gumbel distribution)
-- After the EV fitting, the tool estimates some further statistics from the long-term input file, including minimum, maximum, and different percentile maps. This part can be commented out if not of interest
+The output NetCDF file contains the following return period thresholds [1.5, 2, 5, 10, 20, 50, 100, 200, 500], together with the Gumbel parameters (sigma and mu).
 
-### Reference
-[L-Moments: Analysis and Estimation of Distributions Using Linear Combinations of Order Statistics](https://www.jstor.org/stable/2345653)
+```text
+usage: thresholds [-h] [-d DISCHARGE] [-o OUTPUT]
 
-Hosking, J.R.M., 1990. L-Moments: Analysis and Estimation of Distributions Using Linear Combinations of Order Statistics. J. R. Stat. Soc. Ser. B Methodol. 52, 105124.
+Utility to compute the discharge return period thresholds using the method of L-moments.
+Thresholds computed: [1.5, 2, 5, 10, 20, 50, 100, 200, 500]
 
+options:
+  -h, --help            show this help message and exit
+  -d DISCHARGE, --discharge DISCHARGE
+                        Input discharge files (annual maxima)
+  -o OUTPUT, --output OUTPUT
+                        Output thresholds file
+```
 
 ## waterregions
 
diff --git a/bin/thresholds b/bin/thresholds
new file mode 100755
index 0000000..3db6d45
--- /dev/null
+++ b/bin/thresholds
@@ -0,0 +1,14 @@
+#!python
+
+import os
+import sys
+
+current_dir = os.path.dirname(os.path.abspath(__file__))
+src_dir = os.path.normpath(os.path.join(current_dir, '../src/'))
+if os.path.exists(src_dir):
+    sys.path.append(src_dir)
+
+from lisfloodutilities.thresholds import main_script
+
+if __name__ == '__main__':
+    main_script()
diff --git a/gfit/AUTHORS b/gfit/AUTHORS
deleted file mode 100644
index 10f364e..0000000
--- a/gfit/AUTHORS
+++ /dev/null
@@ -1,1 +0,0 @@
-Alfieri Lorenzo lorenzo.alfieri@ec.europa.eu
diff --git a/gfit/Gfit2.sh b/gfit/Gfit2.sh
deleted file mode 100644
index b0c021c..0000000
--- a/gfit/Gfit2.sh
+++ /dev/null
@@ -1,119 +0,0 @@
-#!/bin/bash
-#$ -S /bin/sh
-
-
-# L.A., 1/3/2019
-# Take discharge files (in netCDF format), extract annual maxima and then fit an extreme value distribution. Produce output maps in PCRaster and netcdf format
-
-# usage: ./Gfit2.sh ./settingFile_test.sh
-# usage (e.g.): qsub -V -j y -N GFit -q all.q -o /nahaUsers/alfielo/CA/GloFAS/scripts/GumbelFit_v2/GFIT_tool.v2/GumFit_nc.log /nahaUsers/alfielo/CA/GloFAS/scripts/GumbelFit_v2/GFIT_tool.v2/Gfit2.sh /nahaUsers/alfielo/CA/GloFAS/scripts/GumbelFit_v2/GFIT_tool.v2/settingFile_test.sh
-
-settFile=${1}
-echo $settFile
-source $settFile
-
-set -o allexport
-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
-
-########################################################################
-# cdo operations
-# Delete first ( number of ) warmup_yrs from the cdo file
-
-echo "Check warm-up years"
-if [ -v warmup_yrs ] ; then
-  if [ $warmup_yrs -ge 1 ] ; then
-    years=`cdo showyear $inpDir/$varName.nc`
-    arr=($years)
-    echo cdo delete,year=$(seq -s, ${arr[0]} ${arr[warmup_yrs-1]}) $inpDir/$varName.nc $outDir/${varName}2.nc
-    cdo delete,year=$(seq -s, ${arr[0]} ${arr[warmup_yrs-1]}) $inpDir/$varName.nc $outDir/${varName}2.nc
-  else
-    ln -sf $inpDir/$varName.nc $outDir/${varName}2.nc
-  fi
-else
-  ln -sf $inpDir/$varName.nc $outDir/${varName}2.nc
-fi
-
-# Extract the series of annual maxima from a netCDF file of the time series
-echo "Extract the series of annual maxima"
-cdo yearmax $outDir/${varName}2.nc $outDir/${varName}AMax.nc
-
-# Extract the average map over the entire time series
-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
-echo "Fit of Gumbel extreme value distributions with L-moments"
-
-/usr/bin/Rscript ${rootdir}gfit2.r $outDir/${varName}AMax.nc $outDir $returnPeriods 0 $MaxGTavgDis $Nyears 
-
-########################################################################
-# Transform ArcView readable ASCII files into PCRaster maps
-if [ -v cloneMap ]; then 
-  echo "create PCRaster maps"
-  cd $outDir
-  Files=$(ls *.txt | sed -e 's/\..*$//')
-  arr=($Files)
-  let nFiles=${#arr[@]}-1 #number of models minus 1 (because the first is #0)
-  
-  for iFile in $(seq 0 $nFiles)  
-  do
-    M=${arr[iFile]}
-    echo ${arr[iFile]}
-    asc2map --clone $cloneMap -S -a -m 1e31 ${arr[iFile]}.txt ${arr[iFile]}.map
-  done
-fi
-
-########################################################################
-# Some more cdo operations to extract statistics
-
-# Extract min and maximum
-echo "extract min, max, and percentile maps"
-cdo timmin $outDir/${varName}2.nc $outDir/${varName}Min.nc
-cdo timmax $outDir/${varName}2.nc $outDir/${varName}Max.nc
-
-# Extract percentile maps
-cdo timpctl,1 $outDir/${varName}2.nc $outDir/${varName}Min.nc $outDir/${varName}Max.nc $outDir/${varName}01p.nc   #1 %
-cdo timpctl,2 $outDir/${varName}2.nc $outDir/${varName}Min.nc $outDir/${varName}Max.nc $outDir/${varName}02p.nc   #2 %
-cdo timpctl,5 $outDir/${varName}2.nc $outDir/${varName}Min.nc $outDir/${varName}Max.nc $outDir/${varName}05p.nc   #5 %
-cdo timpctl,10 $outDir/${varName}2.nc $outDir/${varName}Min.nc $outDir/${varName}Max.nc $outDir/${varName}10p.nc  #10%
-cdo timpctl,25 $outDir/${varName}2.nc $outDir/${varName}Min.nc $outDir/${varName}Max.nc $outDir/${varName}25p.nc  #25%
-cdo timpctl,50 $outDir/${varName}2.nc $outDir/${varName}Min.nc $outDir/${varName}Max.nc $outDir/${varName}50p.nc  #50%
-cdo timpctl,75 $outDir/${varName}2.nc $outDir/${varName}Min.nc $outDir/${varName}Max.nc $outDir/${varName}75p.nc  #75%
-cdo timpctl,90 $outDir/${varName}2.nc $outDir/${varName}Min.nc $outDir/${varName}Max.nc $outDir/${varName}90p.nc  #90%
-cdo timpctl,95 $outDir/${varName}2.nc $outDir/${varName}Min.nc $outDir/${varName}Max.nc $outDir/${varName}95p.nc  #95%
-cdo timpctl,98 $outDir/${varName}2.nc $outDir/${varName}Min.nc $outDir/${varName}Max.nc $outDir/${varName}98p.nc  #98%
-cdo timpctl,99 $outDir/${varName}2.nc $outDir/${varName}Min.nc $outDir/${varName}Max.nc $outDir/${varName}99p.nc  #99%
-
-rm $outDir/${varName}2.nc
-########################################################################
diff --git a/gfit/gfit2.r b/gfit/gfit2.r
deleted file mode 100644
index e77ab43..0000000
--- a/gfit/gfit2.r
+++ /dev/null
@@ -1,223 +0,0 @@
-## Copyright 2019 European Union
-##
-## Licensed under the EUPL, Version 1.2 or as soon they will be approved by the European Commission  subsequent versions of the EUPL (the "Licence");
-##
-## You may not use this work except in compliance with the Licence.
-## You may obtain a copy of the Licence at:
-##
-## https://joinup.ec.europa.eu/sites/default/files/inline-files/EUPL%20v1_2%20EN(1).txt
-##
-## Unless required by applicable law or agreed to in writing, software distributed under the Licence is distributed on an "AS IS" basis,
-## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-## See the Licence for the specific language governing permissions and limitations under the Licence.
-
-
-##purpose: estimating return levels from a Gumbel extreme value distribution functions fit to annual maxima 
-##Created: Zuzanna, August 2015; 
-##Modified: Lorenzo:
-##	  Nov 2015  - user input is now read inline as argument to be run as batch script;
-##	  Nov 2015  - bug fix for warm_up year and addition of Nyears;
-##	  Nov 2015  - header for Ascii file is added on top of the output files. Output files can be directly converted to PCRaster with: asc2map --clone clone.map -S -a ASCIImap.txt ASCIImap.map
-##    Jan 2018  - add option to remove values smaller than the long term average discharge
-##    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> <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
-
-# Following options are not mandatory. If not specified, default values are used
-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)))
-print(rp)
-
-library(ncdf4) ###reading and saving netcdf files 
-library(lmomco) ### L-moments for fitting distributions
-library(ismev) ### for fitting Gumbel distribution, when L-moments is not valid, in such case GLM method is used within - gum.fit
-
-####define function for lmom (to optimize calculation of lmoments, lmom.ub calculates up L5 and other stuff that we don't need, we need only L1 and L2)
-lmom.ub.LA <- function (x) 
-{
-    x2<-x[!is.na(x)]
-	n <- length(x2)
-    if (n == 1) 
-        stop("use mean() for data with one value")
-    if (n < 5) 
-        stop("a minimum of 5 data values are required\n                   because 5 lmoments are to be computed")
-    #if (length(unique(x2)) == 1) 
-    #    stop("all values are equal--lmoments can not be computed")
-    x2 <- sort(x2)
-    L1 = 0
-    L2 = 0
-    for (i in seq(1, n)) {
-        CL1 <- i - 1
-        CR1 <- n - i
-        L1 <- L1 + x2[i]
-        L2 <- L2 + x2[i] * (CL1 - CR1)
-    }
-    C1 <- n
-    C2 <- C1 * (n - 1)/2
-    L1 <- L1/C1
-    L2 <- L2/C2/2
-     z <- list(L1 = L1, L2 = L2, source = "lmom.ub")
-    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[[varNumber]]$dim[[1]]
-yDim = nc$var[[varNumber]]$dim[[2]]
-data = ncvar_get(nc)
-
-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[[varNumber]]$missval
-
-Header<-list(paste0("ncols      ",nCols),paste0("nrows      ",nRows),paste0("xllcenter  ",xllCenter),paste0("yllcenter  ",yllCenter),paste0("cellsize   ",cellSize),paste0("NODATA_value  ",MV))
-
-nc_close(nc)
-
-
-# Remove warm up years from the database
-if (!exists("warmup_yrs")){warmup_yrs<-0}
-if(warmup_yrs>=1){
-  data<-data[,,-c(1:warmup_yrs)]
-}
-
-# Limit the analysis in the time slice from year 1 to Nyears 
-if (!exists("Nyears")){Nyears<-NA}
-if(!is.na(Nyears)){
-  data<-data[,,c(1:Nyears)]
-}
-
-# Remove values smaller than the climatological average discharge
-data2use<-data
-if(!exists("MaxGTavgDis")){MaxGTavgDis<-1}    #default=1 (Yes)   Option (1 or 0) to use only annual maxima larger than the long term average discharge.
-
-if(MaxGTavgDis==0){
-
-}else{
-
-avgFile<-dir(path=out.file.path, pattern="Avg.nc")[1]
-avgDisnc = nc_open(filename = paste0(out.file.path,"/",avgFile),write=FALSE,suppress_dimvals=FALSE)	#open file with annual maxima from cdo
-avgData = ncvar_get(avgDisnc)
-nc_close(nc)
-Ny<-dim(data)[3]
-for (k in 1:Ny){
-  dx<-data2use[,,k]
-  dx[dx<avgData]<-NA
-  data2use[,,k]<-dx
-}
-}
-
-
-# import grid dimensions from the source nc file
-num_grid_x = xDim$len
-num_grid_y = yDim$len
-num__grid = num_grid_x * num_grid_y
-
-out_return_level = array(MV,dim=c(num_grid_x,num_grid_y,length(rp)))	# define array used for storing calculated return levels
-
-
-FunGumFit <- function(yrMax){
-  if(sum(!is.na(yrMax))>=5)	#process only for cells within a mask (unlike observations that may have some missing data, simulation results have only !NA values)
-    {
-      try(xmom<-pargum(lmom.ub.LA(yrMax)))	# calculate L-moments,with modified function; only first 2 - this is all we need.
-	}
-}   
-
-gumParams<-apply(data2use,c(1,2),FunGumFit)
-mu<-sapply(sapply(gumParams, "[[", 2),"[[",1)
-sig<-sapply(sapply(gumParams, "[[", 2),"[[",2)
-
-mu[sapply(mu, is.null)] <- NA
-mu2<-matrix(unlist(mu),nrow=num_grid_x, ncol=num_grid_y)
-
-
-sig[sapply(sig, is.null)] <- NA
-sig2<-matrix(unlist(sig),nrow=num_grid_x, ncol=num_grid_y)
-
-
-q = 1/rp
-for (i in 1:length(rp)){
-
-z = mu2 - sig2 * log(-log(1-q[i])) # Gumbel return level
-z[is.na(z)]<-MV
-out_return_level[,,i] = z
-		
-}
-
-mu2[is.na(mu2)]<-MV 
-sig2[is.na(sig2)]<-MV
-
-### saving results as .nc file 
-# specify netCDF variables
-rpVars = vector("list", length(rp))
-for (k in 1:length(rp))
-	{
-	rpVars[[k]] = ncvar_def(paste("rl", rp[k], sep=""), Units, list(xDim, yDim), MV, longname=paste("Return Level", rp[k]), prec="double")
-	}
-	
-# create a netCDF file with all variables
-rpNcdf <- nc_create(paste0(out.file.path,"/return_levels.nc"), rpVars)
-
-# write values to each variable on disk.
-for (k in 1:length(rp))
-	{
-	ncvar_put( rpNcdf, rpVars[[k]], out_return_level[, , k], verbose=FALSE)
-	}
-
-# add source info metadata to file
-ncatt_put( rpNcdf, 0, "source", "return levels from cdo annual maxima")
-
-nc_close(rpNcdf)
-
-############################################################
-### option of saving results as .txt files (Comment out if not of interest)
- for (k in 1:length(rp))
-	 {
-	 if(rp[k]%%1!=0){
-	 write.table(sapply(Header, "[[", 1),file=paste0(out.file.path,"/rl_", sprintf("%03.1f",rp[k]), ".txt"), row.names = FALSE, col.names = FALSE, quote = FALSE)
-	 write.table(t(out_return_level[, , k]), file=paste0(out.file.path,"/rl_", sprintf("%03.1f",rp[k]), ".txt"), row.names = FALSE, col.names = FALSE, append=TRUE)
-	 }else{
-	 write.table(sapply(Header, "[[", 1),file=paste0(out.file.path,"/rl_", sprintf("%03.0f",rp[k]), ".txt"), row.names = FALSE, col.names = FALSE, quote = FALSE)
-	 write.table(t(out_return_level[, , k]), file=paste0(out.file.path,"/rl_", sprintf("%03.0f",rp[k]), ".txt"), row.names = FALSE, col.names = FALSE, append=TRUE)
-	 }
-	 }
-
-### save gumbel parameters as ascii files 
-write.table(sapply(Header, "[[", 1),file=paste0(out.file.path,"/gumbel_xi.txt"), row.names = FALSE, col.names = FALSE, quote = FALSE)
-write.table(t(mu2), file=paste0(out.file.path,"/gumbel_xi.txt"), row.names = FALSE, col.names = FALSE, append=TRUE)
-
-write.table(sapply(Header, "[[", 1),file=paste0(out.file.path,"/gumbel_alpha.txt"), row.names = FALSE, col.names = FALSE, quote = FALSE)
-write.table(t(sig2), file=paste0(out.file.path,"/gumbel_alpha.txt"), row.names = FALSE, col.names = FALSE, append=TRUE)
-############################################################
-
-print ("time [sec]") # display processing time
-print (proc.time() - ptm)
diff --git a/setup.py b/setup.py
index be55c89..e3fe02c 100644
--- a/setup.py
+++ b/setup.py
@@ -144,7 +144,7 @@ def run(self):
     keywords=['netCDF4', 'PCRaster', 'mapstack', 'lisflood', 'efas', 'glofas', 'ecmwf', 'copernicus'],
     license='EUPL 1.2',
     url='https://github.com/ec-jrc/lisflood-utilities',
-    scripts=['bin/pcr2nc', 'bin/cutmaps', 'bin/compare', 'bin/nc2pcr', ],
+    scripts=['bin/pcr2nc', 'bin/cutmaps', 'bin/compare', 'bin/nc2pcr', 'bin/thresholds', ],
     zip_safe=True,
     classifiers=[
           # complete classifier list: http://pypi.python.org/pypi?%3Aaction=list_classifiers
diff --git a/src/lisfloodutilities/thresholds/__init__.py b/src/lisfloodutilities/thresholds/__init__.py
new file mode 100644
index 0000000..1aca94f
--- /dev/null
+++ b/src/lisfloodutilities/thresholds/__init__.py
@@ -0,0 +1,18 @@
+"""
+
+Copyright 2019-2023 European Union
+
+Licensed under the EUPL, Version 1.2 or as soon they will be approved by the European Commission  subsequent versions of the EUPL (the "Licence");
+
+You may not use this work except in compliance with the Licence.
+You may obtain a copy of the Licence at:
+
+https://joinup.ec.europa.eu/sites/default/files/inline-files/EUPL%20v1_2%20EN(1).txt
+
+Unless required by applicable law or agreed to in writing, software distributed under the Licence is distributed on an "AS IS" basis,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the Licence for the specific language governing permissions and limitations under the Licence.
+
+"""
+
+from .thresholds import *
diff --git a/src/lisfloodutilities/thresholds/thresholds.py b/src/lisfloodutilities/thresholds/thresholds.py
new file mode 100755
index 0000000..37b160e
--- /dev/null
+++ b/src/lisfloodutilities/thresholds/thresholds.py
@@ -0,0 +1,194 @@
+"""
+Copyright 2019-2023 European Union
+Licensed under the EUPL, Version 1.2 or as soon they will be approved by the European Commission  subsequent versions of the EUPL (the "Licence");
+You may not use this work except in compliance with the Licence.
+You may obtain a copy of the Licence at:
+https://joinup.ec.europa.eu/sites/default/files/inline-files/EUPL%20v1_2%20EN(1).txt
+Unless required by applicable law or agreed to in writing, software distributed under the Licence is distributed on an "AS IS" basis,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the Licence for the specific language governing permissions and limitations under the Licence.
+
+"""
+
+import argparse
+import os
+import sys
+
+import numpy as np
+import xarray as xr
+
+
+def lmoments(values: np.ndarray) -> np.ndarray:
+    """
+    Compute first 2 L-moments of dataset (on first axis)
+    :param values: N-D array of values
+    :return: an estimate of the first two sample L-moments
+    """
+
+    nmoments = 3
+
+    # we need to have at least four values in order
+    # to make a sample L-moments estimation
+    nvalues = values.shape[0]
+    if nvalues < 4:
+        raise ValueError(
+            "Insufficient number of values to perform sample L-moments estimation"
+        )
+
+    # sort the values into ascending order
+    values = np.sort(values, axis=0)
+
+    sums = np.zeros((nmoments, *(values.shape[1:])))
+
+    for i in range(1, nvalues + 1):
+        z = i
+        term = values[i - 1]
+        sums[0] = sums[0] + term
+        for j in range(1, nmoments):
+            z -= 1
+            term = term * z
+            sums[j] = sums[j] + term
+
+    y = float(nvalues)
+    z = float(nvalues)
+    sums[0] = sums[0] / z
+    for j in range(1, nmoments):
+        y = y - 1.0
+        z = z * y
+        sums[j] = sums[j] / z
+
+    k = nmoments
+    p0 = -1.0
+    for _ in range(2):
+        ak = float(k)
+        p0 = -p0
+        p = p0
+        temp = p * sums[0]
+        for i in range(1, k):
+            ai = i
+            p = -p * (ak + ai - 1.0) * (ak - ai) / (ai * ai)
+            temp = temp + (p * sums[i])
+        sums[k - 1] = temp
+        k = k - 1
+
+    lmoments = np.zeros((2, *(values.shape[1:])))
+    lmoments[0] = sums[0]
+    lmoments[1] = sums[1]
+
+    return lmoments
+
+
+def gumbel_parameters_moments(dis):
+    sigma = np.sqrt(6) * np.nanstd(dis, ddof=1, axis=0) / np.pi
+    mu = np.nanmean(dis, axis=0) - 0.5772 * sigma
+    return sigma, mu
+
+
+def gumbel_parameters_lmoments(dis):
+    lambda_coef = lmoments(dis)
+    sigma = lambda_coef[1] / np.log(2)
+    mu = lambda_coef[0] - sigma * 0.5772
+    return sigma, mu
+
+
+def gumbel_function(y, sigma, mu):
+    return mu - sigma * np.log(np.log(y / (y - 1)))
+
+
+def find_main_var(ds, path):
+    variable_names = [k for k in ds.variables if len(ds.variables[k].dims) == 3]
+    if len(variable_names) > 1:
+        raise Exception("More than one variable in dataset {}".format(path))
+    elif len(variable_names) == 0:
+        raise Exception("Could not find a valid variable in dataset {}".format(path))
+    else:
+        var_name = variable_names[0]
+    return var_name
+
+
+def read_discharge(in_files):
+    ds = xr.open_dataset(in_files)
+    var = find_main_var(ds, in_files)
+    da = ds[var]
+    return da
+
+
+def unmask_array(mask, template, data):
+    data_unmask = np.empty_like(template)
+    data_unmask[...] = np.NaN
+    data_unmask[mask] = data
+    return data_unmask
+
+
+def create_dataset(dis_max, return_periods, mask, thresholds, sigma, mu):
+    print("Creating dataset")
+    ds_rp = xr.Dataset(
+        coords={"lat": dis_max.coords["lat"], "lon": dis_max.coords["lon"]}
+    )
+    for i, rp in enumerate(return_periods):
+        thres = unmask_array(mask, dis_max.isel(time=0).values, thresholds[i])
+        ds_rp[f"rl_{rp}"] = (["lat", "lon"], thres)
+
+    s = unmask_array(mask, dis_max.isel(time=0).values, sigma)
+    print(s.shape)
+    ds_rp[f"sigma"] = (["lat", "lon"], s)
+    m = unmask_array(mask, dis_max.isel(time=0).values, mu)
+    print(m.shape)
+    ds_rp[f"mu"] = (["lat", "lon"], m)
+
+    print(ds_rp)
+
+    return ds_rp
+
+
+def compute_thresholds_gumbel(dis_max, return_periods):
+    mask = np.isfinite(dis_max.isel(time=0).values)
+    dis_max_masked = dis_max.values[:, mask]
+
+    print("Computing Gumbel coefficients")
+    sigma, mu = gumbel_parameters_lmoments(dis_max_masked)
+
+    print("Computing return periods")
+    thresholds = []
+    for y in return_periods:
+        dis = gumbel_function(y, sigma, mu)
+        thresholds.append(dis)
+
+    ds_rp = create_dataset(dis_max, return_periods, mask, thresholds, sigma, mu)
+
+    return ds_rp
+
+
+def main(argv=sys.argv):
+    prog = os.path.basename(argv[0])
+    parser = argparse.ArgumentParser(
+        description="""
+        Utility to compute the discharge return period thresholds
+        using the method of L-moments.
+        Thresholds computed: [1.5, 2, 5, 10, 20, 50, 100, 200, 500]
+        """,
+        prog=prog,
+    )
+    parser.add_argument(
+        "-d", "--discharge", help="Input discharge files (annual maxima)"
+    )
+    parser.add_argument("-o", "--output", help="Output thresholds file")
+
+    args = parser.parse_args()
+
+    dis = read_discharge(args.discharge)
+    print(dis)
+
+    return_periods = np.array([1.5, 2, 5, 10, 20, 50, 100, 200, 500])
+
+    thresholds = compute_thresholds_gumbel(dis, return_periods)
+
+    thresholds.to_netcdf(args.output)
+
+
+def main_script():
+    sys.exit(main())
+
+
+if __name__ == "__main__":
+    main_script()


Test patch:
diff --git a/gfit/settingFile_test.sh b/gfit/settingFile_test.sh
deleted file mode 100644
index bb53b8d..0000000
--- a/gfit/settingFile_test.sh
+++ /dev/null
@@ -1,13 +0,0 @@
-#!/bin/bash
-# Settings file for the Gumbel fit
-
-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 (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) 
diff --git a/tests/test_thresholds.py b/tests/test_thresholds.py
new file mode 100644
index 0000000..a0bf3f5
--- /dev/null
+++ b/tests/test_thresholds.py
@@ -0,0 +1,72 @@
+import numpy as np
+import unittest
+
+from lisfloodutilities.thresholds import lmoments, gumbel_parameters_moments, \
+    gumbel_parameters_lmoments, gumbel_function
+
+
+class TestLMoments(unittest.TestCase):
+    def test_lmoments_shapes(self):
+        # Test with input array of shape (10, 5)
+        random = np.random.randn(100, 5)
+        lambda_coef = lmoments(random)
+        self.assertEqual(lambda_coef.shape, (2, 5))
+
+        # Test with input array of shape (2, 2)
+        random = np.random.randn(200, 2)
+        lambda_coef = lmoments(random)
+        self.assertEqual(lambda_coef.shape, (2, 2))
+
+        # Test with input array of shape (3,)
+        random = np.random.randn(3)
+        with self.assertRaises(ValueError):
+            lambda_coef = lmoments(random)
+
+    def test_lmoments(self):
+        values = np.arange(1, 100)
+        lambda_coef = lmoments(values)
+        print(lambda_coef)
+        self.assertEqual(lambda_coef[0], 50.)
+        self.assertEqual(lambda_coef[1], 16.66666666666667)
+
+
+class TestGumbelParameters(unittest.TestCase):
+    def setUp(self):
+        self.random = np.random.randn(100, 5)
+        self.values = np.arange(1, 100)
+
+    def test_moments_random(self):
+        # Test output shapes
+        sigma, mu = gumbel_parameters_moments(self.random)
+        self.assertEqual(sigma.shape, (5,))
+        self.assertEqual(mu.shape, (5,))
+
+    def test_lmoments_random(self):
+        # Test output shapes
+        sigma, mu = gumbel_parameters_lmoments(self.random)
+        self.assertEqual(sigma.shape, (5,))
+        self.assertEqual(mu.shape, (5,))
+
+    def test_moments(self):
+        sigma, mu = gumbel_parameters_moments(self.values)
+        # Test output values
+        self.assertEqual(sigma, 22.395085599960808)
+        self.assertEqual(mu, 37.07355659170262)
+
+    def test_lmoments(self):
+        sigma, mu = gumbel_parameters_lmoments(self.values)
+        # Test output values
+        self.assertEqual(sigma, 24.044917348149397)
+        self.assertEqual(mu, 36.12127370664817)
+
+
+class TestGumbelFunction(unittest.TestCase):
+    def test_gumbel_function(self):
+        y = np.array([2, 3, 4, 5])
+        sigma = 1.0
+        mu = 0.0
+        result = gumbel_function(y, sigma, mu)
+        self.assertEqual(result.shape, (4,))
+        print(result)
+        truth = np.array([0.36651292, 0.90272046, 1.24589932, 1.49993999])
+        np.testing.assert_allclose(result, truth)