Fix dist_approx with "geosphere" method and close points (#792)

* u_coord.dist_approx: fix geosphere with log

* Update CHANGELOG

Author: tovogt

CHANGELOG.md

@@ -16,6 +16,8 @@ Code freeze date: YYYY-MM-DD
 
 ### 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/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