Merge branch 'develop' into feature/pickle_free_centroids_store

Author: emanuel-schmid

CHANGELOG.md

@@ -13,6 +13,8 @@ Code freeze date: YYYY-MM-DD
 Removed:
 
 - `pandas-datareader`
+- `intensity_thres` of `Hazard`, `StormEurope`, and `TropCyclones` object [#1065](https://github.com/CLIMADA-project/climada_python/pull/1065)
+- Deprecated method `climada.hazard.trop_cyclone.trop_cyclone.TropCyclone.set_from_tracks` [#1065](https://github.com/CLIMADA-project/climada_python/pull/1065)
 
 ### Added
 - Added optional parameter to `geo_im_from_array`, `plot_from_gdf`, `plot_rp_imp`, `plot_rp_intensity`,

climada/engine/unsequa/test/test_unsequa.py

@@ -650,7 +650,7 @@ def test_calc_sensitivity_all_pass(self):
         }
 
         def test_sensitivity_method(
-            exp_unc, impf_unc, haz_unc, sensitivity_method, param_dict
+            exp_unc, impf_unc, haz_unc, sensitivity_method, param_dict, places
         ):
             """Function to test each seaprate sensitivity method"""
             unc_calc = CalcImpact(exp_unc, impf_unc, haz_unc)
@@ -681,7 +681,7 @@ def test_sensitivity_method(
             self.assertAlmostEqual(
                 param_dict["test_si_value"][0],
                 unc_data.aai_agg_sens_df["aai_agg"][param_dict["test_si_value"][1]],
-                places=5,
+                places=places,
             )
 
             self.assertEqual(unc_data.aai_agg_unc_df.size, unc_data.n_samples)
@@ -695,7 +695,12 @@ def test_sensitivity_method(
         # loop over each method and do test
         for sensitivity_method, method_params in test_dict.items():
             test_sensitivity_method(
-                exp_unc, impf_unc, haz_unc, sensitivity_method, method_params
+                exp_unc,
+                impf_unc,
+                haz_unc,
+                sensitivity_method,
+                method_params,
+                places=2 if sensitivity_method == "rbd_fast" else 5,
             )
 
 

climada/hazard/base.py

@@ -101,10 +101,6 @@ class Hazard(HazardIO, HazardPlot):
         (i.e., is equivalent to fraction is 1 everywhere).
     """
 
-    intensity_thres = 10
-    """Intensity threshold per hazard used to filter lower intensities. To be
-    set for every hazard type"""
-
     vars_oblig = {
         "units",
         "centroids",
@@ -485,7 +481,7 @@ def local_exceedance_intensity(
         self,
         return_periods=(25, 50, 100, 250),
         method="interpolate",
-        min_intensity=None,
+        min_intensity=0,
         log_frequency=True,
         log_intensity=True,
         bin_decimals=None,
@@ -512,8 +508,7 @@ def local_exceedance_intensity(
             periods and extends the interpolation between these points to the given return period
             (similar for small return periods). Defauls to "interpolate".
         min_intensity : float, optional
-            Minimum threshold to filter the hazard intensity. If set to None, self.intensity_thres
-            will be used. Defaults to None.
+            Minimum threshold to filter the hazard intensity. Defaults to 0.
         log_frequency : bool, optional
             If set to True, (cummulative) frequency values are converted to log scale before
             inter- and extrapolation. Defaults to True.
@@ -554,8 +549,6 @@ def local_exceedance_intensity(
         intensities range from 1e6 to 1e9, you could use bin_decimals=-5, if your intensities
         range from 0.0001 to .01, you could use bin_decimals=5.
         """
-        if not min_intensity and min_intensity != 0:
-            min_intensity = self.intensity_thres
         # check frequency unit
         return_period_unit = u_dt.convert_frequency_unit_to_time_unit(
             self.frequency_unit
@@ -639,7 +632,7 @@ def local_return_period(
         self,
         threshold_intensities=(10.0, 20.0),
         method="interpolate",
-        min_intensity=None,
+        min_intensity=0,
         log_frequency=True,
         log_intensity=True,
         bin_decimals=None,
@@ -667,8 +660,7 @@ def local_return_period(
             points to the given threshold intensity (similar for small threshold intensites).
             Defaults to "interpolate".
         min_intensity : float, optional
-            Minimum threshold to filter the hazard intensity. If set to None, self.intensity_thres
-            will be used. Defaults to None.
+            Minimum threshold to filter the hazard intensity. Defaults to 0.
         log_frequency : bool, optional
             If set to True, (cummulative) frequency values are converted to log scale before
             inter- and extrapolation. Defaults to True.
@@ -710,8 +702,6 @@ def local_return_period(
         intensities range from 1e6 to 1e9, you could use bin_decimals=-5, if your intensities
         range from 0.0001 to .01, you could use bin_decimals=5.
         """
-        if not min_intensity and min_intensity != 0:
-            min_intensity = self.intensity_thres
         # check frequency unit
         return_period_unit = u_dt.convert_frequency_unit_to_time_unit(
             self.frequency_unit

climada/hazard/storm_europe.py

@@ -57,6 +57,12 @@
 N_PROB_EVENTS = 5 * 6
 """Number of events per historic event in probabilistic dataset"""
 
+DEF_INTENSITY_THRES = 14.7
+"""
+Default value for the threshold below which wind speeds (in m/s) are stored as 0.
+Same as used by WISC SSI calculations.
+"""
+
 
 class StormEurope(Hazard):
     """A hazard set containing european winter storm events. Historic storm
@@ -79,10 +85,8 @@ class StormEurope(Hazard):
         values only.
     ssi : np.array, float
         SSI as set by set_ssi; uses the Dawkins definition by default.
-    """
 
-    intensity_thres = 14.7
-    """Intensity threshold for storage in m/s; same as used by WISC SSI calculations."""
+    """
 
     vars_opt = Hazard.vars_opt.union({"ssi_wisc", "ssi", "ssi_full_area"})
     """Name of the variables that aren't need to compute the impact."""
@@ -139,7 +143,7 @@ def from_footprints(
         centroids=None,
         files_omit="fp_era20c_1990012515_701_0.nc",
         combine_threshold=None,
-        intensity_thres=None,
+        intensity_thres=DEF_INTENSITY_THRES,
     ):
         """Create new StormEurope object from WISC footprints.
 
@@ -171,18 +175,15 @@ def from_footprints(
             events are combined into one.
             Default is None, Advised for WISC is 2
         intensity_thres : float, optional
-            Intensity threshold for storage in m/s. Default: class attribute
-            StormEurope.intensity_thres (same as used by WISC SSI calculations)
+            Intensity threshold for storage in m/s. Default: 14.7
+            (same as used by WISC SSI calculations)
 
         Returns
         -------
         haz : StormEurope
             StormEurope object with data from WISC footprints.
         """
         # pylint: disable=protected-access
-        intensity_thres = (
-            cls.intensity_thres if intensity_thres is None else intensity_thres
-        )
         file_names = get_file_names(path)
 
         if ref_raster is not None and centroids is not None:
@@ -300,7 +301,7 @@ def from_cosmoe_file(
         event_date=None,
         model_name="COSMO-2E",
         description=None,
-        intensity_thres=None,
+        intensity_thres=DEF_INTENSITY_THRES,
     ):
         """Create a new StormEurope object with gust footprint from weather forecast.
 
@@ -331,17 +332,14 @@ def from_cosmoe_file(
             description of the events, defaults
             to a combination of model_name and run_datetime
         intensity_thres : float, optional
-            Intensity threshold for storage in m/s. Default: class attribute
-            StormEurope.intensity_thres (same as used by WISC SSI calculations)
+            Intensity threshold for storage in m/s. Default: 14.7
+            (same as used by WISC SSI calculations)
 
         Returns
         -------
         haz : StormEurope
             StormEurope object with data from COSMO ensemble file.
         """
-        intensity_thres = (
-            cls.intensity_thres if intensity_thres is None else intensity_thres
-        )
 
         # read intensity from file
         with xr.open_dataset(fp_file) as ncdf:
@@ -435,7 +433,7 @@ def from_icon_grib(
         description=None,
         grib_dir=None,
         delete_raw_data=True,
-        intensity_thres=None,
+        intensity_thres=DEF_INTENSITY_THRES,
     ):
         """Create new StormEurope object from DWD icon weather forecast footprints.
 
@@ -472,18 +470,15 @@ def from_icon_grib(
             .grib.bz2 file format should be stored on the computer or
             removed
         intensity_thres : float, optional
-            Intensity threshold for storage in m/s. Default: class attribute
-            StormEurope.intensity_thres (same as used by WISC SSI calculations)
+            Intensity threshold for storage in m/s. Default: 14.7
+            (same as used by WISC SSI calculations)
 
         Returns
         -------
         haz : StormEurope
             StormEurope object with data from DWD icon weather forecast footprints.
         """
         # pylint: disable=protected-access
-        intensity_thres = (
-            cls.intensity_thres if intensity_thres is None else intensity_thres
-        )
 
         if not (run_datetime.hour == 0 or run_datetime.hour == 12):
             LOGGER.warning(
@@ -724,8 +719,7 @@ def calc_ssi(
             ignoring the intensities at sea. Defaults to true, whereas
             the MATLAB version did not.
         threshold : float, optional
-            Intensity threshold used in index
-            definition. Cannot be lower than the read-in value.
+            Intensity threshold used in index definition.
         sel_cen : np.array, bool
             A boolean vector selecting centroids.
             Takes precendence over on_land.
@@ -739,12 +733,7 @@ def calc_ssi(
             intensity = self.intensity
 
         if threshold is not None:
-            assert (
-                threshold >= self.intensity_thres
-            ), "threshold cannot be below threshold upon read_footprint"
             intensity = intensity.multiply(intensity > threshold)
-        else:
-            intensity = intensity.multiply(intensity > self.intensity_thres)
 
         cent = self.centroids
         area_pixel = cent.get_area_pixel()

climada/hazard/test/test_base.py

@@ -1156,7 +1156,6 @@ def test_local_exceedance_intensity(self):
         """Test local exceedance frequencies with lin lin interpolation"""
         haz = dummy_hazard()
         haz.intensity = sparse.csr_matrix([[1.0, 3.0, 1.0], [2.0, 3.0, 0.0]])
-        haz.intensity_thres = 0.5
         haz.frequency = np.full(2, 1.0)
         return_period = np.array([0.5, 2.0 / 3.0, 1.0])
         # first centroid has intensities 1,2 with cum frequencies 2,1
@@ -1180,17 +1179,19 @@ def test_local_exceedance_intensity_methods(self):
         haz.intensity = sparse.csr_matrix(
             [[0, 0, 1e1], [0.2, 1e1, 1e2], [1e3, 1e3, 1e3]]
         )
-        haz.intensity_thres = 0.5
+        min_intensity = 0.5
         haz.frequency = np.array([1.0, 0.1, 0.01])
-        return_period = (1000, 30, 0.1)
+        return_periods = (1000, 30, 0.1)
         # first centroid has intensities 1e3 with frequencies .01, cum freq .01
         # second centroid has intensities 1e1, 1e3 with cum frequencies .1, .01, cum freq .11, .01
         # third centroid has intensities 1e1, 1e2, 1e3 with cum frequencies 1., .1, .01, cum freq 1.11, .11, .01
         # testing at frequencies .001, .033, 10.
 
         # test stepfunction
         inten_stats, _, _ = haz.local_exceedance_intensity(
-            return_periods=(1000, 30, 0.1), method="stepfunction"
+            return_periods=return_periods,
+            method="stepfunction",
+            min_intensity=min_intensity,
         )
         np.testing.assert_allclose(
             inten_stats.values[:, 1:].astype(float),
@@ -1199,7 +1200,9 @@ def test_local_exceedance_intensity_methods(self):
 
         # test log log extrapolation
         inten_stats, _, _ = haz.local_exceedance_intensity(
-            return_periods=(1000, 30, 0.1), method="extrapolate"
+            return_periods=return_periods,
+            method="extrapolate",
+            min_intensity=min_intensity,
         )
         np.testing.assert_allclose(
             inten_stats.values[:, 1:].astype(float),
@@ -1209,7 +1212,9 @@ def test_local_exceedance_intensity_methods(self):
 
         # test log log interpolation and extrapolation with constant
         inten_stats, _, _ = haz.local_exceedance_intensity(
-            return_periods=(1000, 30, 0.1), method="extrapolate_constant"
+            return_periods=return_periods,
+            method="extrapolate_constant",
+            min_intensity=min_intensity,
         )
         np.testing.assert_allclose(
             inten_stats.values[:, 1:].astype(float),
@@ -1219,7 +1224,7 @@ def test_local_exceedance_intensity_methods(self):
 
         # test log log interpolation and no extrapolation
         inten_stats, _, _ = haz.local_exceedance_intensity(
-            return_periods=(1000, 30, 0.1)
+            return_periods=return_periods, min_intensity=min_intensity
         )
         np.testing.assert_allclose(
             inten_stats.values[:, 1:].astype(float),
@@ -1235,6 +1240,7 @@ def test_local_exceedance_intensity_methods(self):
             log_frequency=False,
             log_intensity=False,
             method="extrapolate_constant",
+            min_intensity=min_intensity,
         )
         np.testing.assert_allclose(
             inten_stats.values[:, 1:].astype(float),
@@ -1270,7 +1276,6 @@ def test_local_return_period_methods(self):
         haz.intensity = sparse.csr_matrix(
             [[0, 0, 1e1], [0.0, 1e1, 1e2], [1e3, 1e3, 1e3]]
         )
-        haz.intensity_thres = 0.5
         haz.frequency = np.array([1.0, 0.1, 0.01])
         # first centroid has intensities 1e3 with frequencies .01, cum freq .01
         # second centroid has intensities 1e1, 1e3 with cum frequencies .1, .01, cum freq .11, .01

climada/hazard/trop_cyclone/trop_cyclone.py

@@ -88,10 +88,8 @@ class TropCyclone(Hazard):
         For each event, the full velocity vectors at each centroid and track position in a sparse
         matrix of shape (npositions, ncentroids * 2) that can be reshaped to a full ndarray of
         shape (npositions, ncentroids, 2).
-    """
 
-    intensity_thres = DEF_INTENSITY_THRES
-    """intensity threshold for storage in m/s"""
+    """
 
     vars_opt = Hazard.vars_opt.union({"category"})
     """Name of the variables that are not needed to compute the impact."""
@@ -142,19 +140,6 @@ def __init__(
         self.basin = basin if basin is not None else []
         self.windfields = windfields if windfields is not None else []
 
-    def set_from_tracks(self, *args, **kwargs):
-        """This function is deprecated, use TropCyclone.from_tracks instead."""
-        LOGGER.warning(
-            "The use of TropCyclone.set_from_tracks is deprecated."
-            "Use TropCyclone.from_tracks instead."
-        )
-        if "intensity_thres" not in kwargs:
-            # some users modify the threshold attribute before calling `set_from_tracks`
-            kwargs["intensity_thres"] = self.intensity_thres
-        if self.pool is not None and "pool" not in kwargs:
-            kwargs["pool"] = self.pool
-        self.__dict__ = TropCyclone.from_tracks(*args, **kwargs).__dict__
-
     @classmethod
     def from_tracks(
         cls,
@@ -381,7 +366,6 @@ def from_tracks(
         LOGGER.debug("Concatenate events.")
         haz = cls.concat(tc_haz_list)
         haz.pool = pool
-        haz.intensity_thres = intensity_thres
         LOGGER.debug("Compute frequency.")
         haz.frequency_from_tracks(tracks.data)
         return haz
@@ -709,7 +693,6 @@ def from_single_track(
         )
 
         new_haz = cls(haz_type=HAZ_TYPE)
-        new_haz.intensity_thres = intensity_thres
         new_haz.intensity = intensity_sparse
         if store_windfields:
             new_haz.windfields = [windfields_sparse]

doc/misc/citation.rst

@@ -20,7 +20,7 @@ If you use specific tools and modules of CLIMADA, please cite the appropriate pu
    * - *Any*
      - The `Zenodo archive <https://doi.org/10.5281/zenodo.4598943>`_ of the CLIMADA version you are using
    * - :doc:`Impact calculations </user-guide/climada_engine_Impact>`
-     - Aznar-Siguan, G. and Bresch, D. N. (2019): CLIMADA v1: A global weather and climate risk assessment platform, Geosci. Model Dev., 12, 3085–3097, https://doi.org/10.5194/gmd-14-351-2021
+     - Aznar-Siguan, G. and Bresch, D. N. (2019): CLIMADA v1: A global weather and climate risk assessment platform, Geosci. Model Dev., 12, 3085–3097, https://doi.org/10.5194/gmd-12-3085-2019
    * - :doc:`Cost-benefit analysis </user-guide/climada_engine_CostBenefit>`
      - Bresch, D. N. and Aznar-Siguan, G. (2021): CLIMADA v1.4.1: Towards a globally consistent adaptation options appraisal tool, Geosci. Model Dev., 14, 351–363, https://doi.org/10.5194/gmd-14-351-2021
    * - :doc:`Uncertainty and sensitivity analysis </user-guide/climada_engine_unsequa>`

doc/user-guide/climada_util_local_exceedance_values.ipynb

@@ -40,7 +40,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -86,7 +86,6 @@
     "    frequency_unit=\"1/year\",\n",
     "    units=\"m/s\",\n",
     ")\n",
-    "hazard.intensity_thres = 0\n",
     "\n",
     "# plot first event of Hazard object\n",
     "hazard.plot_intensity(event=1, smooth=False, figsize=(4, 4));"