Merge branch 'calibrate-impact-functions' of https://github.com/CLIMADA-project/climada_python into calibrate-impact-functions

Author: peanutfun

CHANGELOG.md

@@ -14,8 +14,13 @@ Code freeze date: YYYY-MM-DD
 
 ### Changed
 
+- Update `CONTRIBUTING.md` to better explain types of contributions to this repository [#797](https://github.com/CLIMADA-project/climada_python/pull/797)
+- The default tile layer in Exposures maps is not Stamen Terrain anymore, but [CartoDB Positron](https://github.com/CartoDB/basemap-styles). Affected methods are `climada.engine.Impact.plot_basemap_eai_exposure`,`climada.engine.Impact.plot_basemap_impact_exposure` and `climada.entity.Exposures.plot_basemap`. [#798](https://github.com/CLIMADA-project/climada_python/pull/798)
+
 ### Fixed
 
+- Fix the dist_approx util function when used with method="geosphere" and log=True and points that are very close. [#792](https://github.com/CLIMADA-project/climada_python/pull/792)
+
 ### Deprecated
 
 ### Removed

CONTRIBUTING.md

@@ -1,8 +1,26 @@
 # CLIMADA Contribution Guide
 
-Thank you for contributing to CLIMADA!
+We welcome any contribution to CLIMADA and want to express our thanks to everybody who contributes.
 
-## Overview
+## What Warrants a Contribution?
+
+Anything!
+For orientation, these are some categories of possible contributions we can think of:
+
+* **Technical problems and bugs:** Did you encounter a problem when using CLIMADA? Raise an [issue](https://github.com/CLIMADA-project/climada_python/issues) in our repository, providing a description or ideally a code replicating the error. Did you already find a solution to the problem? Please raise a pull request to help us resolve the issue!
+* **Documentation and Tutorial Updates:** Found a typo in the documentation? Is a tutorial lacking some information you find important? Simply fix a line, or add a paragraph. We are happy to incorporate you additions! Please raise a pull request!
+* **New Modules and Utility Functions:** Did you create a function or an entire module you find useful for your work? Maybe you are not the only one! Feel free to simply raise a pull request for functions that improve, e.g., plotting or data handling. As an entire module has to be carefully integrated into the framework, it might help if you talk to us first so we can design the module and plan the next steps. You can do that by raising an issue or starting a [discussion](https://github.com/CLIMADA-project/climada_python/discussions) on GitHub.
+
+A good place to start a personal discussion is our monthly CLIMADA developers call.
+Please contact the [lead developers](https://wcr.ethz.ch/research/climada.html) if you want to join.
+
+## Why Should You Contribute?
+
+* You will be listed as author of the CLIMADA repository in the [AUTHORS](AUTHORS.md) file.
+* You will improve the quality of the CLIMADA software for you and for everybody else using it.
+* You will gain insights into scientific software development.
+
+## Minimal Steps to Contribute
 
 Before you start, please have a look at our [Developer Guide][devguide].
 

README.md

@@ -61,7 +61,9 @@ Please use the following logo if you are presenting results obtained with or thr
 
 ## Contributing
 
-See the [Contribution Guide](CONTRIBUTING.md).
+We welcome any contribution to this repository, be it bugfixes and other code changes and additions, documentation improvements, or tutorial updates.
+
+If you would like to contribute, please refer to our [Contribution Guide](CONTRIBUTING.md).
 
 ## Versioning
 

climada/engine/impact.py

@@ -673,7 +673,7 @@ def plot_raster_eai_exposure(self, res=None, raster_res=None, save_tiff=None,
 
     def plot_basemap_eai_exposure(self, mask=None, ignore_zero=False, pop_name=True,
                                   buffer=0.0, extend='neither', zoom=10,
-                                  url=ctx.providers.Stamen.Terrain,
+                                  url=ctx.providers.CartoDB.Positron,
                                   axis=None, **kwargs):
         """Plot basemap expected impact of each exposure within a period of 1/frequency_unit.
 
@@ -694,7 +694,7 @@ def plot_basemap_eai_exposure(self, mask=None, ignore_zero=False, pop_name=True,
         zoom : int, optional
             zoom coefficient used in the satellite image
         url : str, optional
-            image source, e.g. ctx.providers.OpenStreetMap.Mapnik
+            image source, default: ctx.providers.CartoDB.Positron
         axis : matplotlib.axes.Axes, optional
             axis to use
         kwargs : dict, optional
@@ -764,7 +764,7 @@ def plot_hexbin_impact_exposure(self, event_id=1, mask=None, ignore_zero=False,
 
     def plot_basemap_impact_exposure(self, event_id=1, mask=None, ignore_zero=False,
                                      pop_name=True, buffer=0.0, extend='neither', zoom=10,
-                                     url=ctx.providers.Stamen.Terrain,
+                                     url=ctx.providers.CartoDB.Positron,
                                      axis=None, **kwargs):
         """Plot basemap impact of an event at each exposure.
         Requires attribute imp_mat.
@@ -789,7 +789,7 @@ def plot_basemap_impact_exposure(self, event_id=1, mask=None, ignore_zero=False,
         zoom : int, optional
             zoom coefficient used in the satellite image
         url : str, optional
-            image source, e.g. ctx.providers.OpenStreetMap.Mapnik
+            image source, default: ctx.providers.CartoDB.Positron
         axis : matplotlib.axes.Axes, optional
             axis to use
         kwargs : dict, optional

climada/entity/exposures/base.py

@@ -749,7 +749,7 @@ def plot_raster(self, res=None, raster_res=None, save_tiff=None,
 
     def plot_basemap(self, mask=None, ignore_zero=False, pop_name=True,
                      buffer=0.0, extend='neither', zoom=10,
-                     url=None, axis=None, **kwargs):
+                     url=ctx.providers.CartoDB.Positron, axis=None, **kwargs):
         """Scatter points over satellite image using contextily
 
         Parameters
@@ -771,7 +771,7 @@ def plot_basemap(self, mask=None, ignore_zero=False, pop_name=True,
             zoom coefficient used in the satellite image
         url : Any, optional
             image source, e.g., ``ctx.providers.OpenStreetMap.Mapnik``.
-            Default: ``ctx.providers.Stamen.Terrain``
+            Default: ``ctx.providers.CartoDB.Positron``
         axis : matplotlib.axes._subplots.AxesSubplot, optional
             axis to use
         kwargs : optional

climada/hazard/test/test_trop_cyclone.py

@@ -304,22 +304,61 @@ def test_v_max_s_holland_2008_pass(self):
 
     def test_holland_2010_pass(self):
         """Test Holland et al. 2010 wind field model."""
-        # test at centroids within and outside of radius of max wind
+        # The parameter "x" is designed to be exactly 0.5 inside the radius of max wind (RMW) and
+        # to increase or decrease linearly outside of it in radial direction.
+        #
+        # An increase (decrease) of "x" outside of the RMW is for cases where the max wind is very
+        # high (low), but the RMW is still comparably large (small). This means, wind speeds need
+        # to decay very sharply (only moderately) outside of the RMW to reach the low prescribed
+        # peripheral wind speeds.
+        #
+        # The "hol_b" parameter tunes the meaning of a "comparably" large or small RMW.
         si_track = xr.Dataset({
-            "rad": ("time", KM_TO_M * np.array([75, 40])),
-            "vmax": ("time", [35.0, 40.0]),
-            "hol_b": ("time", [1.80, 2.5]),
+            # four test cases:
+            # - low vmax, moderate RMW: x decreases moderately
+            # - large hol_b: x decreases sharply
+            # - very low vmax: x decreases so much, it needs to be clipped at 0
+            # - large vmax, large RMW: x increases
+            "rad": ("time", KM_TO_M * np.array([75, 75, 75, 90])),
+            "vmax": ("time", [35.0, 35.0, 16.0, 90.0]),
+            "hol_b": ("time", [1.75, 2.5, 1.9, 1.6]),
         })
-        d_centr = KM_TO_M * np.array([[35, 75, 220], [30, 1000, 300]], dtype=float)
-        close_centr = np.array([[True, True, True], [True, False, True]], dtype=bool)
+        d_centr = KM_TO_M * np.array([
+            # first column is for locations within the storm eye
+            # second column is for locations at or close to the radius of max wind
+            # third column is for locations outside the storm eye
+            # fourth column is for locations exactly at the peripheral radius
+            # fifth column is for locations outside the peripheral radius
+            [0., 75, 220, 300, 490],
+            [30, 74, 170, 300, 501],
+            [21, 76, 230, 300, 431],
+            [32, 91, 270, 300, 452],
+        ], dtype=float)
+        close_centr = np.array([
+            # note that we set one of these to "False" for testing
+            [True, True, True, True, True],
+            [True, True, True, True, False],
+            [True, True, True, True, True],
+            [True, True, True, True, True],
+        ], dtype=bool)
         hol_x = _x_holland_2010(si_track, d_centr, close_centr)
-        np.testing.assert_array_almost_equal(
-            hol_x, [[0.5, 0.5, 0.47273], [0.5, 0, 0.211602]])
+        np.testing.assert_array_almost_equal(hol_x, [
+            [0.5, 0.500000, 0.485077, 0.476844, 0.457291],
+            [0.5, 0.500000, 0.410997, 0.289203, 0.000000],
+            [0.5, 0.497620, 0.131072, 0.000000, 0.000000],
+            [0.5, 0.505022, 1.403952, 1.554611, 2.317948],
+        ])
 
-        # test exactly at radius of maximum wind (35 m/s) and at peripheral radius (17 m/s)
         v_ang_norm = _stat_holland_2010(si_track, d_centr, close_centr, hol_x)
-        np.testing.assert_array_almost_equal(v_ang_norm,
-            [[15.957853, 35.0, 20.99411], [33.854826, 0, 17.0]])
+        np.testing.assert_array_almost_equal(v_ang_norm, [
+            # first column: converge to 0 when approaching storm eye
+            # second column: vmax at RMW
+            # fourth column: peripheral speed (17 m/s) at peripheral radius (unless x is clipped!)
+            [0.0000000, 35.000000, 21.181497, 17.00000, 12.103461],
+            [1.2964800, 34.990037, 21.593755, 17.00000, 0.0000000],
+            [0.3219518, 15.997500, 13.585498, 16.00000, 16.000000],
+            [24.823469, 89.992938, 24.381965, 17.00000, 1.9292020],
+        ])
 
     def test_stat_holland_1980(self):
         """Test _stat_holland_1980 function. Compare to MATLAB reference."""

climada/hazard/trop_cyclone.py

@@ -1424,7 +1424,13 @@ def _x_holland_2010(
     # linearly interpolate between max exponent and peripheral exponent
     x_max = 0.5
     hol_x[close_centr] = x_max + np.fmax(0, d_centr - r_max) * (x_n - x_max) / (r_n - r_max)
-    hol_x[close_centr] = np.clip(hol_x[close_centr], 0.0, 0.5)
+
+    # Negative hol_x values appear when v_max_s is very close to or even lower than v_n (which
+    # should never happen in theory). In those cases, wind speeds might decrease outside of the eye
+    # wall and increase again towards the peripheral radius (which is actually unphysical).
+    # We clip hol_x to 0, otherwise wind speeds keep increasing indefinitely away from the eye:
+    hol_x[close_centr] = np.fmax(hol_x[close_centr], 0.0)
+
     return hol_x
 
 def _stat_holland_2010(

climada/test/test_plot.py

@@ -162,10 +162,7 @@ def test_ctx_osm_pass(self):
         myexp.gdf['value'] = np.array([1, 1, 1])
         myexp.check()
 
-        try:
-            myexp.plot_basemap(url=ctx.providers.OpenStreetMap.Mapnik)
-        except urllib.error.HTTPError:
-            self.assertEqual(1, 0)
+        myexp.plot_basemap(url=ctx.providers.OpenStreetMap.Mapnik)
 
     def test_disc_rates(self):
         """Test plot function of discount rates."""

climada/util/api_client.py

@@ -341,7 +341,7 @@ def _divide_straight_from_multi(properties):
             elif isinstance(_v, list):
                 multis[k] = _v
             else:
-                raise ValueError("properties must be a string or a list of strings")
+                raise ValueError("the value of a property must be a string or a list of strings")
         return straights, multis
 
     @staticmethod

climada/util/calibrate/base.py

@@ -8,6 +8,7 @@
 import pandas as pd
 import numpy as np
 from scipy.optimize import Bounds, LinearConstraint, NonlinearConstraint
+import matplotlib.pyplot as plt
 import seaborn as sns
 
 from climada.hazard import Hazard
@@ -154,6 +155,131 @@ def plot_impf_set(self, **plot_kwargs):
         """
         return self.impf_set.plot(**plot_kwargs)
 
+    def plot_impf_variability(
+        self,
+        cost_func_diff: float = 0.1,
+        p_space_df: Optional[pd.DataFrame] = None,
+        plot_haz: bool = True,
+        plot_impf_kws: Optional[dict] = None,
+        plot_hist_kws: Optional[dict] = None,
+    ):
+
+        """Plot impact function variability with parameter combinations of
+        almost equal cost function values
+
+        Args:
+            cost_func_diff (float, optional): Max deviation from optimal cost
+                function value (as fraction). Defaults to 0.1 (i.e. 10%).
+            p_space_df (pd.DataFrame, optional): parameter space. Defaults to None.
+            plot_haz (bool, optional): Whether or not to plot hazard intensity
+                distibution. Defaults to False.
+            plot_impf_kws (dict, optional): Keyword arguments for impact
+                function plot. Defaults to None.
+            plot_hist_kws (dict, optional): Keyword arguments for hazard
+                intensity distribution plot. Defaults to None.
+        """
+
+        # Initialize plot keyword arguments
+        if plot_impf_kws is None:
+            plot_impf_kws = {}
+        if plot_hist_kws is None:
+            plot_hist_kws = {}
+
+        # Retrieve hazard type and parameter space
+        haz_type = self.input.hazard.haz_type
+        if p_space_df is None:
+            # Assert that self.output has the p_space_to_dataframe() method,
+            # which is defined for the BayesianOptimizerOutput class
+            if not hasattr(self.output,"p_space_to_dataframe"):
+                raise TypeError(
+                    "To derive the full impact function parameter space, "
+                    "plot_impf_variability() requires BayesianOptimizerOutput "
+                    "as OutputEvaluator.output attribute, which provides the "
+                    "method p_space_to_dataframe()."
+                )
+            p_space_df = self.output.p_space_to_dataframe()
+
+        # Retrieve list of parameters required for creating impact functions
+        # and remove the dimension 'Cost Function'.
+        params = p_space_df.columns.tolist()
+        try:
+            params.remove('Cost Function')
+        except ValueError:
+            pass
+
+        # Retrieve parameters of impact functions with cost function values
+        # within 'cost_func_diff' % of the best estimate
+        params_within_range = p_space_df[params]
+        plot_space_label = 'Parameter space'
+        if cost_func_diff is not None:
+            max_cost_func_val = (p_space_df['Cost Function'].min()*
+                                 (1+cost_func_diff))
+            params_within_range = p_space_df.loc[
+                p_space_df['Cost Function'] <=max_cost_func_val,params
+            ]
+            plot_space_label = (f"within {int(cost_func_diff*100)} percent "
+                                f"of best fit")
+
+        # Set plot defaults
+        color = plot_impf_kws.pop('color','tab:blue')
+        lw = plot_impf_kws.pop('lw',2)
+        zorder = plot_impf_kws.pop('zorder',3)
+        label = plot_impf_kws.pop('label','best fit')
+
+        #get number of impact functions and create a plot for each
+        n_impf = len(self.impf_set.get_func(haz_type=haz_type))
+        axes=[]
+
+        for impf_idx in range(n_impf):
+
+            _,ax = plt.subplots()
+
+            #Plot best-fit impact function
+            best_impf = self.impf_set.get_func(haz_type=haz_type)[impf_idx]
+            ax.plot(best_impf.intensity,best_impf.mdd*best_impf.paa*100,
+                    color=color,lw=lw,zorder=zorder,label=label,**plot_impf_kws)
+
+            #Plot all impact functions within 'cost_func_diff' % of best estimate
+            for row in range(params_within_range.shape[0]):
+                label_temp = plot_space_label if row == 0 else None
+
+                sel_params = params_within_range.iloc[row,:].to_dict()
+                temp_impf_set = self.input.impact_func_creator(**sel_params)
+                temp_impf = temp_impf_set.get_func(haz_type=haz_type)[impf_idx]
+
+                ax.plot(temp_impf.intensity,temp_impf.mdd*temp_impf.paa*100,
+                        color='grey',alpha=0.4,label=label_temp)
+
+            # Plot hazard intensity value distributions
+            if plot_haz:
+                haz_vals = self.input.hazard.intensity[
+                    :, self.input.exposure.gdf[f"centr_{haz_type}"]
+                ]
+
+                #Plot defaults
+                color_hist = plot_hist_kws.pop('color','tab:orange')
+                alpha_hist = plot_hist_kws.pop('alpha',0.3)
+
+                ax2 = ax.twinx()
+                ax2.hist(haz_vals.data,bins=40,color=color_hist,
+                        alpha=alpha_hist,label='Hazard intensity\noccurence')
+                ax2.set(ylabel='Hazard intensity occurence (#Exposure points)')
+                ax.axvline(x=haz_vals.max(),label='Maximum hazard value',
+                        color='tab:orange')
+                ax2.legend(loc='lower right')
+
+            ax.set(xlabel=f"Intensity ({self.input.hazard.units})",
+                ylabel="Mean Damage Ratio (MDR) in %",
+                xlim=(min(best_impf.intensity),max(best_impf.intensity)))
+            ax.legend()
+            axes.append(ax)
+
+        if n_impf > 1:
+            return axes
+
+        return ax
+
+
     def plot_at_event(
         self,
         data_transf: Callable[[pd.DataFrame], pd.DataFrame] = lambda x: x,