Implement SAL verification method (#248)

* Add initial structure

* Add SAL method

* Refactor imports

* Add MissingOptionalDependency checks

* Rename methods

* Apply suggestions from code review

Co-authored-by: Andres Perez Hortal <16256571+aperezhortal@users.noreply.github.com>

* Add basic test

* Add interface

* More testing

* Return a tuple instead of pandas dataframe

* Some refactoring

* More refactoring

* Avoid re-computing the same calculations

* Extract sal method into separate module

* Minor refactoring

* Rename tests

* Fix compability issue and rename module

* Remove unused argument

* Add thr_factor argument

* Fix backward compatibility

* Rename variable

* Fix black

* Fix black part 2

* Replace max_precip with quantile95_precip

Replace max_precip with quantile95_precip for calculating the threshold.

* Convert masked arrays to numpy arrays

This avoids problems with nanquantile

* Fix skimage backward compatibility

* Add bibliography

* Compute quantile on wet pixels only

* Fix property name

* Set fill value in all cases

* Make sure that it works with empty inputs

Return nan when a score is undefined

* Aesthetics

* Add no-cache-dir flag

* Pin numpy build version to match the one available on conda-forge

* Update module name

* Pin maximum numpy version

* Revert commit 94ea041

* Set minmax=minref

Following suggestion from @feldmann-m

* Subpress all splits if thr_factor is passed

Following suggestion from @feldmann-m

* New set of parameters for tstorm

* Add option to specifiy the detection quantile

* Improve docstrings

* Set minimum python version to 3.7  (#253)

* Set minimum python version to 3.7

and maximum version to 3.9

* Update github actions

Co-authored-by: Andres Perez Hortal <16256571+aperezhortal@users.noreply.github.com>

* Add docstring references

* Use mamba to setup test environments (#256)

* Try using mamba to speedup installing dependencies

* Use mamba.bat command for windows

Co-authored-by: Esmail Ghaemi <esi@Esmails-Air.fritz.box>
Co-authored-by: Esmail Ghaemi <esi@Esmails-MacBook-Air.local>
Co-authored-by: Daniele Nerini <daniele.nerini@gmail.com>
Co-authored-by: Andres Perez Hortal <16256571+aperezhortal@users.noreply.github.com>
Co-authored-by: Daniele Nerini <daniele.nerini@meteoswiss.ch>

Author: 5 people

doc/source/pysteps_reference/verification.rst

@@ -11,4 +11,5 @@ Methods for verification of deterministic, probabilistic and ensemble forecasts.
 .. automodule:: pysteps.verification.lifetime
 .. automodule:: pysteps.verification.plots
 .. automodule:: pysteps.verification.probscores
+.. automodule:: pysteps.verification.salscores
 .. automodule:: pysteps.verification.spatialscores

doc/source/references.bib

@@ -52,6 +52,16 @@ @ARTICLE{EWWM2013
   DOI = "10.1002/met.1392"
 }
 
+@ARTICLE{Feldmann2021,
+  AUTHOR = "M. Feldmann and U. Germann and M. Gabella and A. Berne",
+  TITLE = "A Characterisation of Alpine Mesocyclone Occurrence",
+  JOURNAL = "Weather and Climate Dynamics Discussions",
+  PAGES = "1--26",
+  URL = "https://wcd.copernicus.org/preprints/wcd-2021-53/",
+  DOI = "10.5194/wcd-2021-53",
+  YEAR = 2021
+}
+
 @ARTICLE{FSNBG2019,
   AUTHOR = "Foresti, L. and Sideris, I.V. and Nerini, D. and Beusch, L. and Germann, U.",
   TITLE = "Using a 10-Year Radar Archive for Nowcasting Precipitation Growth and Decay: A Probabilistic Machine Learning Approach",
@@ -264,6 +274,35 @@ @ARTICLE{SPN2013
   DOI = "10.1002/wrcr.20536"
 }
 
+@ARTICLE{TRT2004,
+  AUTHOR = "A. M. Hering and C. Morel and G. Galli and P. Ambrosetti and M. Boscacci",
+  TITLE = "Nowcasting thunderstorms in the Alpine Region using a radar based adaptive thresholding scheme",
+  JOURNAL = "Proceedings of ERAD Conference 2004",
+  NUMBER = "January",
+  PAGES = "206--211",
+  YEAR = 2004
+}
+
+@ARTICLE{WHZ2009,
+  AUTHOR = "Heini  Wernli and Christiane  Hofmann and Matthias  Zimmer",
+  TITLE = "Spatial Forecast Verification Methods Intercomparison Project: Application of the SAL Technique",
+  JOURNAL = "Weather and Forecasting",
+  NUMBER = "6",
+  VOLUME = "24",
+  PAGES = "1472 - 14847",
+  YEAR = 2009
+}
+
+@ARTICLE{WPHF2008,
+  AUTHOR = "Heini Wernli and Marcus Paulat and Martin Hagen and Christoph Frei",
+  TITLE = "SAL—A Novel Quality Measure for the Verification of Quantitative Precipitation Forecasts",
+  JOURNAL = "Monthly Weather Review",
+  NUMBER = "11",
+  VOLUME = "136",
+  PAGES = "4470 - 4487",
+  YEAR = 2008
+}
+
 @ARTICLE{XWF2005,
   AUTHOR = "K. Xu and C. K Wikle and N. I. Fox",
   TITLE = "A Kernel-Based Spatio-Temporal Dynamical Model for Nowcasting Weather Radar Reflectivities",
@@ -284,23 +323,3 @@ @ARTICLE{ZR2009
   YEAR = 2009,
   DOI = "10.1016/j.atmosres.2009.03.004"
 }
-
-@ARTICLE{TRT2004,
-  AUTHOR = "A. M. Hering and C. Morel and G. Galli and P. Ambrosetti and M. Boscacci",
-  TITLE = "Nowcasting thunderstorms in the Alpine Region using a radar based adaptive thresholding scheme",
-  JOURNAL = "Proceedings of ERAD Conference 2004",
-  NUMBER = "January",
-  PAGES = "206--211",
-  YEAR = 2004
-}
-
-@ARTICLE{Feldmann2021,
-  AUTHOR = "M. Feldmann and U. Germann and M. Gabella and A. Berne",
-  TITLE = "A Characterisation of Alpine Mesocyclone Occurrence",
-  JOURNAL = "Weather and Climate Dynamics Discussions",
-  PAGES = "1--26",
-  URL = "https://wcd.copernicus.org/preprints/wcd-2021-53/",
-  DOI = "10.5194/wcd-2021-53",
-  YEAR = 2021
-}
-

pyproject.toml

@@ -24,4 +24,4 @@ exclude = '''
   | build
   | dist
 )/
-'''
+'''

pysteps/tests/helpers.py

@@ -195,9 +195,9 @@ def get_precipitation_fields(
                 reference_field, ref_metadata, threshold=0.1, zerovalue=-15.0
             )
 
-            # Set missing values with the fill value
-            np.ma.set_fill_value(reference_field, -15.0)
-            reference_field.data[reference_field.mask] = -15.0
+        # Set missing values with the fill value
+        np.ma.set_fill_value(reference_field, ref_metadata["zerovalue"])
+        reference_field.data[reference_field.mask] = ref_metadata["zerovalue"]
 
     if metadata:
         return reference_field, ref_metadata

pysteps/tests/test_ensemblestats.py renamed to pysteps/tests/test_postprocessing_ensemblestats.py

@@ -0,0 +1,56 @@
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import pytest
+
+from pysteps.tests.helpers import get_precipitation_fields
+from pysteps.verification.salscores import sal
+from pysteps.utils import to_rainrate, to_reflectivity
+
+test_data = [
+    (to_rainrate, 1 / 15),
+    (to_reflectivity, None),
+]
+
+
+@pytest.mark.parametrize("converter, thr_factor", test_data)
+class TestSAL:
+    pytest.importorskip("pandas")
+    pytest.importorskip("skimage")
+
+    def test_sal_zeros(self, converter, thr_factor):
+        """Test the SAL verification method."""
+        precip, metadata = get_precipitation_fields(
+            num_prev_files=0, log_transform=False, metadata=True
+        )
+        precip, metadata = converter(precip.filled(np.nan), metadata)
+        result = sal(precip * 0, precip * 0, thr_factor)
+        assert np.isnan(result).all()
+        result = sal(precip * 0, precip, thr_factor)
+        assert result[:2] == (-2, -2)
+        assert np.isnan(result[2])
+        result = sal(precip, precip * 0, thr_factor)
+        assert result[:2] == (2, 2)
+        assert np.isnan(result[2])
+
+    def test_sal_same_image(self, converter, thr_factor):
+        """Test the SAL verification method."""
+        precip, metadata = get_precipitation_fields(
+            num_prev_files=0, log_transform=False, metadata=True
+        )
+        precip, metadata = converter(precip.filled(np.nan), metadata)
+        result = sal(precip, precip, thr_factor)
+        assert isinstance(result, tuple)
+        assert len(result) == 3
+        assert np.allclose(result, [0, 0, 0])
+
+    def test_sal_translation(self, converter, thr_factor):
+        precip, metadata = get_precipitation_fields(
+            num_prev_files=0, log_transform=False, metadata=True
+        )
+        precip, metadata = converter(precip.filled(np.nan), metadata)
+        precip_translated = np.roll(precip, 10, axis=0)
+        result = sal(precip, precip_translated, thr_factor)
+        assert np.allclose(result[0], 0)
+        assert np.allclose(result[1], 0)
+        assert not np.allclose(result[2], 0)