Merge branch 'develop' into calibrate-impact-functions

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/coordinates.py

@@ -358,10 +358,13 @@ def dist_approx(lat1, lon1, lat2, lon2, log=False, normalize=True,
         if log:
             vec1, vbasis = latlon_to_geosph_vector(lat1, lon1, rad=True, basis=True)
             vec2 = latlon_to_geosph_vector(lat2, lon2, rad=True)
-            scal = 1 - 2 * hav
-            fact = dist / np.fmax(np.spacing(1), np.sqrt(1 - scal**2))
-            vtan = fact[..., None] * (vec2[:, None, :] - scal[..., None] * vec1[:, :, None])
+            vtan = vec2[:, None, :] - (1 - 2 * hav[..., None]) * vec1[:, :, None]
             vtan = np.einsum('nkli,nkji->nklj', vtan, vbasis)
+            # faster version of `vtan_norm = np.linalg.norm(vtan, axis=-1)`
+            vtan_norm = np.sqrt(np.einsum("...l,...l->...", vtan, vtan))
+            # for consistency, set dist to 0 if vtan is 0
+            dist[vtan_norm < np.spacing(1)] = 0
+            vtan *= dist[..., None] / np.fmax(np.spacing(1), vtan_norm[..., None])
     else:
         raise KeyError("Unknown distance approximation method: %s" % method)
     return (dist, vtan) if log else dist

climada/util/test/test_coordinates.py

@@ -277,13 +277,20 @@ def test_geosph_vector(self):
     def test_dist_approx_pass(self):
         """Test approximate distance functions"""
         data = np.array([
-            # lat1, lon1, lat2, lon2, dist, dist_sph
+            # lat1, lon1, lat2, lon2, dist_equirect, dist_geosphere
             [45.5, -32.1, 14, 56, 7702.88906574, 8750.64119051],
             [45.5, 147.8, 14, -124, 7709.82781473, 8758.34146833],
             [45.5, 507.9, 14, -124, 7702.88906574, 8750.64119051],
             [45.5, -212.2, 14, -124, 7709.82781473, 8758.34146833],
             [-3, -130.1, 4, -30.5, 11079.7217421, 11087.0352544],
         ])
+        # conversion factors from reference data (in km, see above) to other units
+        factors_km_to_x = {
+            "m": 1e3,
+            "radian": np.radians(1.0) / u_coord.ONE_LAT_KM,
+            "degree": 1.0 / u_coord.ONE_LAT_KM,
+            "km": 1.0,
+        }
         compute_dist = np.stack([
             u_coord.dist_approx(data[:, None, 0], data[:, None, 1],
                                 data[:, None, 2], data[:, None, 3],
@@ -297,9 +304,7 @@ def test_dist_approx_pass(self):
             self.assertAlmostEqual(d[0], cd[0])
             self.assertAlmostEqual(d[1], cd[1])
 
-        for units, factor in zip(["radian", "degree", "km"],
-                                 [np.radians(1.0), 1.0, u_coord.ONE_LAT_KM]):
-            factor /= u_coord.ONE_LAT_KM
+        for units, factor in factors_km_to_x.items():
             compute_dist = np.stack([
                 u_coord.dist_approx(data[:, None, 0], data[:, None, 1],
                                     data[:, None, 2], data[:, None, 3],
@@ -309,21 +314,29 @@ def test_dist_approx_pass(self):
                                     method="geosphere", units=units)[:, 0, 0],
             ], axis=-1)
             self.assertEqual(compute_dist.shape[0], data.shape[0])
+            places = 4 if units == "m" else 7
             for d, cd in zip(data[:, 4:], compute_dist):
-                self.assertAlmostEqual(d[0] * factor, cd[0])
-                self.assertAlmostEqual(d[1] * factor, cd[1])
+                self.assertAlmostEqual(d[0] * factor, cd[0], places=places)
+                self.assertAlmostEqual(d[1] * factor, cd[1], places=places)
 
     def test_dist_approx_log_pass(self):
         """Test log-functionality of approximate distance functions"""
         data = np.array([
-            # lat1, lon1, lat2, lon2, dist, dist_sph
+            # lat1, lon1, lat2, lon2, dist_equirect, dist_geosphere
             [0, 0, 0, 1, 111.12, 111.12],
             [-13, 179, 5, -179, 2011.84774049, 2012.30698122],
+            [24., 85., 23.99999967, 85., 3.666960e-5, 3.666960e-5],
+            [24., 85., 24., 85., 0, 0],
         ])
+        # conversion factors from reference data (in km, see above) to other units
+        factors_km_to_x = {
+            "m": 1e3,
+            "radian": np.radians(1.0) / u_coord.ONE_LAT_KM,
+            "degree": 1.0 / u_coord.ONE_LAT_KM,
+            "km": 1.0,
+        }
         for i, method in enumerate(["equirect", "geosphere"]):
-            for units, factor in zip(["radian", "degree", "km"],
-                                     [np.radians(1.0), 1.0, u_coord.ONE_LAT_KM]):
-                factor /= u_coord.ONE_LAT_KM
+            for units, factor in factors_km_to_x.items():
                 dist, vec = u_coord.dist_approx(data[:, None, 0], data[:, None, 1],
                                                 data[:, None, 2], data[:, None, 3],
                                                 log=True, method=method, units=units)
@@ -613,7 +626,7 @@ def test_match_centroids(self):
             u_coord.match_centroids(gdf, centroids)
         self.assertIn('Set hazard and GeoDataFrame to same CRS first!',
                       str(cm.exception))
-        
+
     def test_dist_sqr_approx_pass(self):
         """Test approximate distance helper function."""
         lats1 = 45.5