Merge branch 'develop' into feature/fix_holland_2010

CHANGELOG.md

@@ -14,8 +14,12 @@ Code freeze date: YYYY-MM-DD
 
 ### Changed
 
+- 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

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/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

@@ -357,10 +357,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