Add ensemble functionality

makefile

@@ -19,8 +19,7 @@ coverage: test
 dist:
 	echo $(VERSION)
 	rm -rf dist
-	python3 setup.py sdist
-	python3 setup.py bdist_wheel
+	python3 -m build --sdist --wheel
 	@ echo "Next, run 'twine upload dist/*'"
 
 clean:

other/verif.sh

@@ -9,9 +9,9 @@ if [ "$cur" = "" ] || [[ "$cur" =~ -* ]]; then
    COMPREPLY=( $( compgen -f -W '-m --config --list-times --list-dates --list-locations --list-quantiles --list-thresholds --help -d -elevrange -l -lx -latrange -lonrange -o -obsrange -r -q -t -tod -x -acc -agg -b -c -C -fcst -hist -obs -sort -T -Tagg -Tx -a -af -afs -aspect -bottom -clabel -clim -cmap -dpi -gc -gs -gw -f -fs -labfs -lc -left -leg -legfs -legloc -ls -lw -maptype -ma -ms -nogrid -nomargin -obsleg -right -simple -sp -tickfs -title -titlefs -top -type -xlabel -xlim -xlog -xrot -xticks -xticklabels -ylabel -ylim -ylog -yrot -yticks -yticklabels' -- $cur ) )
 fi
 if [ "$prev" = "-m" ]; then
-   COMPREPLY=( $( compgen -W 'a against alphaindex autocorr autocov b baserate bias biasfreq bs bsdecomp bsrel bsres bss bssrel bssres bsunc c change cmae cond corr d derror diff discrimination dmb droc droc0 dscore economicvalue edi eds ef error ets fa far fcst fcstrate freq fss hit hss ign0 igncontrib invreliability kendallcorr kge kss leps lor mae marginal marginalratio mbias meteo miss murphy n nnsec nsec obs obsfcst or pc performance pit pithist pithistdev pithistshape pithistslope qq quantile quantilecoverage quantilescore rankcorr ratio reliability rmse rmsf roc scatter sedi seds spherical spread spreadskill spreadskillratio stderror obsstddev fcststddev taylor threat threshold timeseries within yulesq ' -- $cur ) )
+   COMPREPLY=( $( compgen -W 'a against alphaindex autocorr autocov b baserate bias biasfreq bs bsdecomp bsrel bsres bss bssrel bssres bsunc c change cmae cond corr d derror diff discrimination dmb droc droc0 dscore economicvalue edi eds ef error ets fa far fcrps fcst fcstrate fcststddev freq fss hit hss ign0 igncontrib invreliability kendallcorr kge kss leps lor mae marginal marginalratio mbias meteo miss murphy n nnsec nsec obs obsfcst obsstddev or pc performance pit pithist pithistdev pithistshape pithistslope qq quantile quantilecoverage quantilescore rankcorr rankhist ratio reliability rmse rmsf roc scatter sedi seds spherical spread spreadskill spreadskillratio stderror taylor threat threshold timeseries within yulesq ' -- $cur ) )
 elif [ "$prev" = "-cmap" ]; then
-   COMPREPLY=( $( compgen -W 'Accent Accent_r Blues Blues_r BrBG BrBG_r BuGn BuGn_r BuPu BuPu_r CMRmap CMRmap_r Dark2 Dark2_r GnBu GnBu_r Grays Grays_r Greens Greens_r Greys Greys_r OrRd OrRd_r Oranges Oranges_r PRGn PRGn_r Paired Paired_r Pastel1 Pastel1_r Pastel2 Pastel2_r PiYG PiYG_r PuBu PuBuGn PuBuGn_r PuBu_r PuOr PuOr_r PuRd PuRd_r Purples Purples_r RdBu RdBu_r RdGy RdGy_r RdPu RdPu_r RdYlBu RdYlBu_r RdYlGn RdYlGn_r Reds Reds_r Set1 Set1_r Set2 Set2_r Set3 Set3_r Spectral Spectral_r Wistia Wistia_r YlGn YlGnBu YlGnBu_r YlGn_r YlOrBr YlOrBr_r YlOrRd YlOrRd_r afmhot afmhot_r autumn autumn_r berlin berlin_r binary binary_r bone bone_r brg brg_r bwr bwr_r cividis cividis_r cool cool_r coolwarm coolwarm_r copper copper_r cubehelix cubehelix_r flag flag_r gist_earth gist_earth_r gist_gray gist_gray_r gist_grey gist_grey_r gist_heat gist_heat_r gist_ncar gist_ncar_r gist_rainbow gist_rainbow_r gist_stern gist_stern_r gist_yarg gist_yarg_r gist_yerg gist_yerg_r gnuplot gnuplot2 gnuplot2_r gnuplot_r gray gray_r grey grey_r hot hot_r hsv hsv_r inferno inferno_r jet jet_r magma magma_r managua managua_r nipy_spectral nipy_spectral_r ocean ocean_r pink pink_r plasma plasma_r prism prism_r rainbow rainbow_r seismic seismic_r spring spring_r summer summer_r tab10 tab10_r tab20 tab20_r tab20b tab20b_r tab20c tab20c_r terrain terrain_r turbo turbo_r twilight twilight_r twilight_shifted twilight_shifted_r vanimo vanimo_r viridis viridis_r winter winter_r' -- $cur ) )
+   COMPREPLY=( $( compgen -W 'Accent Accent_r Blues Blues_r BrBG BrBG_r BuGn BuGn_r BuPu BuPu_r CMRmap CMRmap_r Dark2 Dark2_r GnBu GnBu_r Grays Greens Greens_r Greys Greys_r OrRd OrRd_r Oranges Oranges_r PRGn PRGn_r Paired Paired_r Pastel1 Pastel1_r Pastel2 Pastel2_r PiYG PiYG_r PuBu PuBuGn PuBuGn_r PuBu_r PuOr PuOr_r PuRd PuRd_r Purples Purples_r RdBu RdBu_r RdGy RdGy_r RdPu RdPu_r RdYlBu RdYlBu_r RdYlGn RdYlGn_r Reds Reds_r Set1 Set1_r Set2 Set2_r Set3 Set3_r Spectral Spectral_r Wistia Wistia_r YlGn YlGnBu YlGnBu_r YlGn_r YlOrBr YlOrBr_r YlOrRd YlOrRd_r afmhot afmhot_r autumn autumn_r binary binary_r bone bone_r brg brg_r bwr bwr_r cividis cividis_r cool cool_r coolwarm coolwarm_r copper copper_r cubehelix cubehelix_r flag flag_r gist_earth gist_earth_r gist_gray gist_gray_r gist_grey gist_heat gist_heat_r gist_ncar gist_ncar_r gist_rainbow gist_rainbow_r gist_stern gist_stern_r gist_yarg gist_yarg_r gist_yerg gnuplot gnuplot2 gnuplot2_r gnuplot_r gray gray_r grey hot hot_r hsv hsv_r inferno inferno_r jet jet_r magma magma_r nipy_spectral nipy_spectral_r ocean ocean_r pink pink_r plasma plasma_r prism prism_r rainbow rainbow_r seismic seismic_r spring spring_r summer summer_r tab10 tab10_r tab20 tab20_r tab20b tab20b_r tab20c tab20c_r terrain terrain_r turbo turbo_r twilight twilight_r twilight_shifted twilight_shifted_r viridis viridis_r winter winter_r' -- $cur ) )
 elif [ "$prev" = "-agg" ] || [ "$prev" = "-Tagg" ]; then
    COMPREPLY=( $( compgen -W 'abschange absmean change count iqr max mean meanabs median min quantile range std sum variance' -- $cur ) )
 elif [ "$prev" = "-type" ]; then

pyproject.toml

@@ -10,12 +10,12 @@ authors = [
 description='A verification program for meteorological forecasts and observations'
 keywords=['meteorology', 'verification', 'weather prediction']
 
-license = {file = 'BSD-3'}
+license = "BSD-3-Clause"
+license-files = ["LICENSE"]
 readme = "README.md"
 
 requires-python = ">=3.9"
 classifiers = [
-      "License :: OSI Approved :: BSD License",
       'Intended Audience :: Science/Research',
       'Topic :: Scientific/Engineering :: Atmospheric Science',
       'Topic :: Scientific/Engineering :: Information Analysis',

scripts/time_interpolation.py

@@ -0,0 +1,179 @@
+import os
+import sys
+import numpy as np
+import matplotlib.pylab as mpl
+import argparse
+import verif.util
+import netCDF4
+import shutil
+import scipy.interpolate
+
+
+def main():
+    """
+    It takes an existing Verif file, extracts 6h forecasts, and fills in the intermediate hours with
+    either linearly or splin interpolated values.
+
+    Ir interpolates fcst, ensemble, ensemble_mean, ensmble_variable fields and creates
+    ensemble_crps.
+    """
+    parser = argparse.ArgumentParser(description='This program downsamples a verif file and interpolates the data back using an interpolation method')
+    parser.add_argument('ifile', help='')
+    parser.add_argument('ofile', help='')
+    parser.add_argument('-i', choices=['linear', 'spline'], dest='interpolator', required=True)
+    parser.add_argument('-f', type=int, dest="frequency", required=True)
+
+    if len(sys.argv) == 1:
+        parser.print_help()
+        sys.exit(1)
+
+    args = parser.parse_args()
+    frequency = args.frequency
+
+    if frequency < 2:
+        raise ValueError("frequency must be > 1")
+
+    shutil.copy(args.ifile, args.ofile)
+
+    with netCDF4.Dataset(args.ofile, 'a') as file:
+        leadtimes = file.variables["leadtime"][:]
+        num_leadtimes = len(leadtimes)
+
+        has_ens = "ensemble" in file.variables
+        if has_ens:
+            ensemble = file.variables["ensemble"][:]
+            print(np.sum(np.isnan(ensemble)))
+
+        if frequency >= num_leadtimes:
+            raise ValueError(f"frequency must be < {num_leadtimes}")
+
+        fcst = file.variables["fcst"][:]
+        print("Interpolating... ")
+
+        if args.interpolator == "linear":
+            # Linear interpolation works, even if there are NaNs in the dataset, the result will
+            # just be NaNs, as expected.
+            ends = range(frequency, num_leadtimes+1, frequency)
+            for end in ends:
+                start = end - frequency
+
+                fcst_start = fcst[:, start, :]
+                fcst_end = fcst[:, end, :]
+                fcst[:, (start+1):end, :] = np.nan  # Erase everything in the middle
+                for i in range(1, frequency):
+                    w1 = i / frequency
+                    w0 = 1 - w1
+                    fcst[:, start + i, :] = w0 * fcst_start + w1 * fcst_end
+
+                if has_ens:
+                    ensemble[:, (start+1):end, :, :] = np.nan
+                    fcst_start = ensemble[:, start, :, :]
+                    fcst_end = ensemble[:, end, :]
+                    for i in range(1, frequency):
+                        w1 = i / frequency
+                        w0 = 1 - w1
+                        ensemble[:, start + i, :, :] = w0 * fcst_start + w1 * fcst_end
+
+            # Trim the last leadtimes that cannot be interpolated
+            fcst[:, ends[-1]+1:, :] = np.nan
+            if has_ens:
+                ensemble[:, ends[-1]+1:, :, :] = np.nan
+
+        elif args.interpolator == "spline":
+            y = fcst[:, 0::frequency, :]
+            x = leadtimes[0::frequency]
+
+            interpolator = scipy.interpolate.CubicSpline(x, y, axis=1)
+            xx = set([i for i in leadtimes]) - set([i for i in x])
+            xx = np.array(list(xx))
+            temp = interpolator(xx)
+            for i, index in enumerate(xx):
+                fcst[:, int(index), :] = temp[:, i, :]
+
+            if has_ens:
+                y = ensemble[:, 0::frequency, :, :]
+                x = leadtimes[0::frequency]
+
+                # Spline can't deal with NaNs in y, so set NaNs to a 0 and then put the Nans back in
+                # later. This means that if a single timestep is missing for a member, the whole
+                # member is NaN for that time/location.
+                is_nan = np.sum(np.isnan(y), axis=1) > 0
+                y[np.isnan(y)] = 0
+
+                interpolator = scipy.interpolate.CubicSpline(x, y, axis=1)
+                xx = set([i for i in leadtimes]) - set([i for i in x])
+                xx = np.array(list(xx))
+                interp = interpolator(xx)
+                for i, index in enumerate(xx):
+                    temp = interp[:, i, :, :]
+                    temp[is_nan] = np.nan
+                    ensemble[:, int(index), :, :] = temp
+
+        file.variables["fcst"][:] = fcst
+        if has_ens:
+            file.variables["ensemble"][:] = ensemble
+
+        # These variables are no longer valid, since we cannot derive them for interpolated
+        # timeseries
+        invalid_variables = ["x", "cdf", "pdf", "pit", "ensemble_mean", "ensemble_variance"]
+        for variable in invalid_variables:
+            if variable in file.variables:
+                file.variables[variable] = np.nan
+
+        # Compute CRPS
+        if has_ens:
+            print("Computing CRPS...")
+            if not "ensemble_crps" in file.variables:
+                file.createVariable("ensemble_crps", "f4", ("time", "leadtime", "location"))
+
+            num_members = ensemble.shape[-1]
+            crps = compute_crps(ensemble, file.variables["obs"][:], num_members)
+
+            num_valid_members = np.sum(~np.isnan(ensemble), axis=3)
+            crps[num_valid_members != num_members] = np.nan
+
+            file.variables["ensemble_crps"][:] = crps
+
+            if "ensemble_mean" in file.variables:
+                ens_mean = np.nanmean(ensemble, axis=3)
+                file.variables["ensemble_mean"][:] = ens_mean
+            if "ensemble_variance" in file.variables:
+                file.variables["ensemble_mean"][:] = np.nanvar(ensemble, axis=3)
+
+
+def compute_crps(preds, targets, num_members, fair=True) -> np.ndarray:
+    """Continuous Ranked Probability Score (CRPS).
+
+    Taken from bris-inference
+
+    Args:
+        preds: numpy.ndarray
+            Predictions, shape (time, leadtime, location, ens_size)
+        targets: numpy.ndarray
+            Targets, shape (time, leadtime, location)
+        fair: bool
+            Defaults to true
+
+    Returns:
+        crps: numpy.ndarray
+            Shape (time, leadtime, location)
+    """
+
+    coef = (
+        -1.0 / (num_members * (num_members - 1))
+        if fair
+        else -1.0 / (num_members**2)
+    )
+
+    mae = np.mean(np.abs(targets[..., None] - preds), axis=-1)
+
+    # var = np.abs(preds[..., None] - preds[..., None, :])
+    var = np.zeros(preds.shape[:-1])
+    for i in range(num_members):  # loop version to reduce memory usage
+        var += np.sum(np.abs(preds[..., i, None] - preds[..., i + 1 :]), axis=-1)
+    var *= coef
+    return mae + var
+
+
+if __name__ == "__main__":
+    main()

verif/__main__.py

@@ -0,0 +1,6 @@
+import sys
+import verif.driver
+
+
+if __name__ == '__main__':
+    verif.driver.run(sys.argv)

verif/aggregator.py

@@ -168,6 +168,7 @@ def __call__(self, array, axis=None):
         else:
             raise NotImplementedError(f"This function not implemented for axis {axis}")
 
+
 class AbsChange(Aggregator):
     """Absolute value of difference between the last and the first element. Is most useful
     when used with -Tagg since it implies an array representing a sequence

verif/axis.py

@@ -38,6 +38,7 @@ def get_leadtime_axes():
     return [axis[1]() for axis in get_all() if hasattr(axis[1], "compute_from_leadtimes")]
     return [verif.axis.Leadtime(), verif.axis.Leadtimeday()]
 
+
 def get_spatial_axes():
     return [axis[1]() for axis in get_all() if axis[1]().is_location_like]