diff --git a/.github/workflows/ci_env.yml b/.github/workflows/ci_env.yml index 5fd9d04..9be0e90 100644 --- a/.github/workflows/ci_env.yml +++ b/.github/workflows/ci_env.yml @@ -24,4 +24,4 @@ jobs: - name: Test with pytest run: | conda install pytest - pytest + pytest \ No newline at end of file diff --git a/.gitignore b/.gitignore index 90821c7..63653d7 100644 --- a/.gitignore +++ b/.gitignore @@ -37,10 +37,8 @@ vignettes/*.pdf *.utf8.md *.knit.md -*.yml *.yaml *.json -!environment.yml dist/ *.egg-info/ **/build/ diff --git a/README.md b/README.md index dc2aa40..08f9f52 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. @@ -228,10 +230,10 @@ This tool cut netcdf files, using a mask, a bounding box or a list of stations a ### Usage: The tool accepts as input: -* a mask map (either PCRaster or netCDF format) or - - alternatively, matrix indices in the form xmini_xmaxi:ymini_ymaxi or - - alternatively, coordinates bounding box in the form xmin_xmax:ymin_ymax - - alternatively, list of stations with coordinates and a LDD map. +* a mask map (either PCRaster or netCDF format) using the -m argument or + - alternatively, using the -i argument, matrix indices in the form `imin imax jmin jmax` (imin, imax, jmin, jmax must be integer numbers) + - alternatively, using the -c argument, coordinates bounding box in the form `xmin xmax ymin ymax` (xmin, xmax, ymin, ymax can be integer or floating point numbers; x = longitude, y = latitude) + - alternatively, using the -N and -l arguments, list of stations with coordinates and a LDD map. * a path to a folder containing netCDF files to cut or a static dataset path like LISFLOOD static files. * a path to a folder where to write cut files. @@ -245,14 +247,14 @@ The mask can also be in PCRaster format. cutmaps -m /workarea/Madeira/maps/MaskMap/Bacia_madeira.nc -f /workarea/Madeira/lai/ -o ./ ``` -**Indices can also be passed as an argument (using -c argument instead of -m). Knowing your area of interest from your netCDF files, -you can determine indices of the array and you can pass in the form `imin imax jmin jmax`.** +**Indices can also be passed as an argument (using -i argument instead of -m). Knowing your area of interest from your netCDF files, +you can determine indices of the array and you can pass in the form `imin imax jmin jmax` (imin, imax, jmin, jmax must be integer numbers).** ```bash -cutmaps -c "150 350 80 180" -f /workarea/Madeira/lai/ -o ./ +cutmaps -i "150 350 80 180" -f /workarea/Madeira/lai/ -o ./ ``` -**Example with coordinates and path to EFAS/GloFAS static data (-S option), with -W to allow overwriting existing files in output directory:** +**Example with coordinates (using -c argument) `xmin xmax ymin ymax` (xmin, xmax, ymin, ymax can be integer or floating point numbers; x = longitude, y = latitude) and path to EFAS/GloFAS static data (-S option), with -W to allow overwriting existing files in output directory:** ```bash cutmaps -S /home/projects/lisflood-eu -c "4078546.12 4463723.85 811206.57 1587655.50" -o /Work/Tunisia/cutmaps -W @@ -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/environment.yml b/environment.yml index 0b1b89c..e076115 100644 --- a/environment.yml +++ b/environment.yml @@ -215,4 +215,4 @@ dependencies: - wcwidth==0.2.5 - webencodings==0.5.1 - xarray==0.20.2 - - zipp==3.8.0 + - zipp==3.8.0 \ No newline at end of file diff --git a/gfit/AUTHORS b/gfit/AUTHORS deleted file mode 100644 index 10f364e..0000000 --- a/gfit/AUTHORS +++ /dev/null @@ -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 -# 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 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=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/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/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/VERSION b/src/lisfloodutilities/VERSION index c018e45..330a1e9 100644 --- a/src/lisfloodutilities/VERSION +++ b/src/lisfloodutilities/VERSION @@ -1 +1 @@ -0.12.20 +0.12.21 diff --git a/src/lisfloodutilities/compare/main.py b/src/lisfloodutilities/compare/main.py index cb186bc..614c1fa 100644 --- a/src/lisfloodutilities/compare/main.py +++ b/src/lisfloodutilities/compare/main.py @@ -71,6 +71,9 @@ def main(cliargs): if errors!=None: for i, e in enumerate(errors): logger.error('%d - %s', i, e) + else: + logger.info('No differences in Datasets %s and %s', dataset_a, dataset_b) + def main_script(): diff --git a/src/lisfloodutilities/compare/nc.py b/src/lisfloodutilities/compare/nc.py index 08960a2..43d5b22 100644 --- a/src/lisfloodutilities/compare/nc.py +++ b/src/lisfloodutilities/compare/nc.py @@ -113,6 +113,7 @@ def compare_files(self, file_a, file_b, timestep=None): raise ValueError('timestep must be of type datetime.datetime or a range of dates, but type {} was found'.format(str(type(timestep)))) logger.info('Comparing %s and %s %s', file_a, file_b, '(from %s to %s)' % (min(timestep), max(timestep)) if timestep else '') + self.diff_timesteps = [] with Dataset(file_a) as nca, Dataset(file_b) as ncb: num_dims = 3 if 'time' in nca.variables else 2 @@ -197,7 +198,6 @@ def compare_arrays(self, vara_step, varb_step, varname=None, step=None, filepath filepath = os.path.basename(filepath) if filepath else '' varname = varname or '' message = '{}/{}@{} - {:3.2f}% of different values - max diff: {:3.6f}'.format(filepath, varname, step, perc_wrong, max_diff) - logger.error(message) if self.for_testing: assert False, message else: diff --git a/src/lisfloodutilities/cutmaps/cutlib.py b/src/lisfloodutilities/cutmaps/cutlib.py index 06c9f31..84305f0 100644 --- a/src/lisfloodutilities/cutmaps/cutlib.py +++ b/src/lisfloodutilities/cutmaps/cutlib.py @@ -33,17 +33,20 @@ encoding_netcdf_vars = {'zlib': False} -def cutmap(f, fileout, x_min, x_max, y_min, y_max): +def cutmap(f, fileout, x_min, x_max, y_min, y_max, use_coords = True): nc, num_dims = open_dataset(f) var = [v for v in nc.variables if len(nc.variables[v].dims) == num_dims][0] logger.info('Variable: %s', var) - - if isinstance(x_min, float): + + if use_coords: # bounding box input from user # FIXME weak isinstance test sliced_var = cut_from_coords(nc, var, x_min, x_max, y_min, y_max) else: - # user provides with indices directly (not coordinates) - sliced_var = cut_from_indices(nc, var, x_min, x_max, y_min, y_max) + if isinstance(x_min, float): + raise ValueError('box values must be integer when using cut_indices') + else: + # user provides with indices directly (not coordinates) + sliced_var = cut_from_indices(nc, var, x_min, x_max, y_min, y_max) if sliced_var is not None: if 'missing_value' in sliced_var.encoding: @@ -153,7 +156,7 @@ def get_filelist(input_folder=None, static_data_folder=None): return list_to_cut -def get_cuts(cuts=None, mask=None): +def get_cuts(cuts=None, cuts_indices=None, mask=None): if mask: if not os.path.isfile(mask): raise FileNotFoundError('Wrong input mask: %s not a file' % mask) @@ -174,13 +177,18 @@ def get_cuts(cuts=None, mask=None): else: logger.error('Mask map format not recognized. Must be either .map or .nc. Found %s', ext) sys.exit(1) + logger.info('MASK: \nmin x: %s \nmax x: %s \nmin y: %s \nmax y: %s', x_min, x_max, y_min, y_max) elif cuts: # user provided coordinates bounds x_min, x_max, y_min, y_max = cuts + logger.info('CUTS: \nmin x: %s \nmax x: %s \nmin y: %s \nmax y: %s', x_min, x_max, y_min, y_max) + elif cuts_indices: + # user provided indices bounds + x_min, x_max, y_min, y_max = cuts_indices + logger.info('CUTS_INDICES: \nmin x: %s \nmax x: %s \nmin y: %s \nmax y: %s', x_min, x_max, y_min, y_max) else: - logger.error('You must provide either cuts (in the format "lonmin lonmax latmin latmax") or a mask map') + logger.error('You must provide either cuts (in the format "lonmin lonmax latmin latmax") or cuts_indices (in the format "imin imax jmin jmax") or a mask map') sys.exit(1) - logger.info('CUTS: \nmin x: %s \nmax x: %s \nmin y: %s \nmax y: %s', x_min, x_max, y_min, y_max) return x_min, x_max, y_min, y_max diff --git a/src/lisfloodutilities/cutmaps/main.py b/src/lisfloodutilities/cutmaps/main.py index cda486c..034acc2 100644 --- a/src/lisfloodutilities/cutmaps/main.py +++ b/src/lisfloodutilities/cutmaps/main.py @@ -30,14 +30,14 @@ def parse_and_check_args(parser, cliargs): args = parser.parse_args(cliargs) - if args.mask and args.cuts: - parser.error('[--mask | --cuts] arguments are mutually exclusive') - if not (args.mask or args.cuts) and not (args.ldd and args.stations): - parser.error('(--mask | --cuts | [--ldd, --stations]) You need to pass mask path or cuts coordinates ' + if (args.mask!=None) + (args.cuts!=None) + (args.cuts_indices!=None) > 1: + parser.error('[--mask | --cuts | --cuts_indices] arguments are mutually exclusive') + if not (args.mask or args.cuts or args.cuts_indices) and not (args.ldd and args.stations): + parser.error('(--mask | --cuts | --cuts_indices | [--ldd, --stations]) You need to pass mask path or cuts coordinates ' 'or a list of stations along with LDD path') - if (args.mask or args.cuts) and (args.ldd or args.stations): - parser.error('(--mask | --cuts | [--ldd, --stations]) ' - '--mask, --cuts and --ldd and --stations arguments are mutually exclusive') + if (args.mask or args.cuts or args.cuts_indices) and (args.ldd or args.stations): + parser.error('(--mask | --cuts | --cuts_indices | [--ldd, --stations]) ' + '--mask, --cuts, --cuts_indices and --ldd and --stations arguments are mutually exclusive') return args def get_arg_coords(value): @@ -59,7 +59,8 @@ def add_args(self): 'create mask cookie-cutter on-fly from stations list and ldd map') group_filelist = self.add_mutually_exclusive_group(required=True) group_mask.add_argument("-m", "--mask", help='mask file cookie-cutter, .map if pcraster, .nc if netcdf') - group_mask.add_argument("-c", "--cuts", help='Cut coordinates in the form "lonmin lonmax latmin latmax"', type=get_arg_coords) + group_mask.add_argument("-c", "--cuts", help='Cut coordinates in the form "lonmin lonmax latmin latmax" using coordinates bounding box', type=get_arg_coords) + group_mask.add_argument("-i", "--cuts_indices", help='Cut coordinates in the form "imin imax jmin jmax" using matrix indices', type=get_arg_coords) group_mask.add_argument("-l", "--ldd", help='Path to LDD file') group_mask.add_argument("-N", "--stations", help='Path to stations.txt file.' @@ -84,6 +85,7 @@ def main(cliargs): args = parse_and_check_args(parser, cliargs) mask = args.mask cuts = args.cuts + cuts_indices = args.cuts_indices ldd = args.ldd stations = args.stations @@ -107,13 +109,13 @@ def main(cliargs): shutil.copy(mask, os.path.join(pathout, 'my_mask.map')) shutil.copy(mask_nc, os.path.join(pathout, 'my_mask.nc')) + x_min, x_max, y_min, y_max = get_cuts(cuts=cuts, cuts_indices=cuts_indices, mask=mask) logger.info('\n\nCutting using: %s\n Files to cut from: %s\n Output: %s\n Overwrite existing: %s\n\n', - mask or cuts, + mask or ([x_min, x_max, y_min, y_max if cuts or cuts_indices else None]), input_folder or static_data_folder, pathout, overwrite) list_to_cut = get_filelist(input_folder, static_data_folder) - x_min, x_max, y_min, y_max = get_cuts(cuts=cuts, mask=mask) # walk through list_to_cut for file_to_cut in list_to_cut: @@ -142,7 +144,7 @@ def main(cliargs): logger.warning('%s already existing. This file will not be overwritten', fileout) continue - cutmap(file_to_cut, fileout, x_min, x_max, y_min, y_max) + cutmap(file_to_cut, fileout, x_min, x_max, y_min, y_max, use_coords=(cuts_indices is None)) if ldd and stations: with Dataset(os.path.join(pathout, 'my_mask.nc'),'r',format='NETCDF4_CLASSIC') as mask_map: for k in mask_map.variables.keys(): diff --git a/src/lisfloodutilities/readers/nc.py b/src/lisfloodutilities/readers/nc.py index 2d88869..0d45a5e 100644 --- a/src/lisfloodutilities/readers/nc.py +++ b/src/lisfloodutilities/readers/nc.py @@ -1,4 +1,5 @@ import xarray as xr +from netCDF4 import default_fillvals x_coordinates_names = ('lon', 'x', 'longitude') y_coordinates_names = ('lat', 'y', 'latitude') @@ -19,7 +20,7 @@ def mv(self): for variable in self.ds.variables.values(): if len(variable.dims) < 2: continue - return variable.attrs['_FillValue'] if '_FillValue' in variable.attrs else -9999 + return variable.attrs['_FillValue'] if '_FillValue' in variable.attrs else default_fillvals[variable.dtype.str[1:]] @property def data(self): diff --git a/src/lisfloodutilities/readers/pcr.py b/src/lisfloodutilities/readers/pcr.py index fca8cf3..7556f47 100644 --- a/src/lisfloodutilities/readers/pcr.py +++ b/src/lisfloodutilities/readers/pcr.py @@ -197,12 +197,12 @@ def __init__(self, input_set): def get_metadata_from_set(self): """ Get metadata from first map of a set. - :return: A dictionary with keys ['rows', 'cols', 'lats', 'lons', 'dtype'] + :return: A dictionary with keys ['rows', 'cols', 'lats', 'lons', 'dtype', 'mv'] """ first_map, _ = next(self.fileset) return {'rows': first_map.rows, 'cols': first_map.cols, 'lats': first_map.lats, 'lons': first_map.lons, - 'dtype': first_map.data.dtype} + 'dtype': first_map.data.dtype, 'mv': first_map.mv} @property def fileset(self): 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() diff --git a/src/lisfloodutilities/writers/nc.py b/src/lisfloodutilities/writers/nc.py index a94bde3..3a768fb 100644 --- a/src/lisfloodutilities/writers/nc.py +++ b/src/lisfloodutilities/writers/nc.py @@ -25,7 +25,7 @@ from pathlib import Path import numpy as np -from netCDF4 import Dataset +from netCDF4 import Dataset, default_fillvals from ..readers import PCRasterMap from .. import logger @@ -84,11 +84,16 @@ def __init__(self, filename, is_mapstack=True, **metadata): self.hour = float(self.metadata.get('time', {}).get('hour') or 0) # you can pass the MV to set in netcdf files directly in yaml configuration, otherwuse np.nan is used - self.mv = self.metadata['variable'].get('mv') + self.mv = self.metadata['variable'].get('mv') + if self.mv is None: + self.mv = self.metadata.get('mv') if self.mv is not None: self.mv = int(self.mv) if np.issubdtype(self.metadata['dtype'], np.integer) else float(self.mv) else: - self.mv = np.nan + # take missing values from the default netCDF fillvals values or NaN + dtype = self.metadata['dtype'] + self.mv = default_fillvals[dtype.str[1:]] if np.issubdtype(self.metadata['dtype'], np.integer) else np.nan + self.current_idx1 = 0 self.current_idx2 = 0 self.time, self.variable = self._init_dataset() diff --git a/src/lisfloodutilities/writers/pcr.py b/src/lisfloodutilities/writers/pcr.py index 8e21df6..203dbca 100644 --- a/src/lisfloodutilities/writers/pcr.py +++ b/src/lisfloodutilities/writers/pcr.py @@ -54,12 +54,12 @@ def __init__(self, *args, **kwargs): rs = ma.masked_values(rs, no_data_value, copy=False) self._mask = ma.getmask(rs) rs = None - self.mv = no_data_value if not mv else float(mv) + self.mv = no_data_value if mv is None else float(mv) else: self.cols = self._coordinates['x'].size self.rows = self._coordinates['y'].size self.src_geo_transform = self._get_geo_transform_from_coords() - self.mv = np.nan if not mv else float(mv) + self.mv = np.nan if mv is None else float(mv) self.flipped_x = self._coordinates['x'][0] > self._coordinates['x'][1] self.flipped_y = self._coordinates['y'][0] < self._coordinates['y'][1] diff --git a/tests/test_cutmaps.py b/tests/test_cutmaps.py index 4d6c4cc..910320f 100644 --- a/tests/test_cutmaps.py +++ b/tests/test_cutmaps.py @@ -74,14 +74,14 @@ def test_get_cuts_withcoords(self): def test_get_cuts_indices(self): # "minxi maxxi minyi maxyi" - cuts = '3 7 1 2' - cuts = get_arg_coords(cuts) - ix_min, ix_max, iy_min, iy_max = get_cuts(cuts=cuts) + cuts_indices = '3 7 1 2' + cuts_indices = get_arg_coords(cuts_indices) + ix_min, ix_max, iy_min, iy_max = get_cuts(cuts_indices=cuts_indices) assert (ix_min, ix_max, iy_min, iy_max) == (3, 7, 1, 2) fin = 'tests/data/folder_a/ta.nc' fout = 'tests/data/folder_a/ta_cut.nc' self.cleanups.append((os.unlink, (fout,))) - cutmap(fin, fout, ix_min, ix_max, iy_min, iy_max) + cutmap(fin, fout, ix_min, ix_max, iy_min, iy_max, use_coords=False) with Dataset(fout) as nc: lons = nc.variables['lon'][:] lats = nc.variables['lat'][:] 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)