From c8a05872c9cd7c67655b269799b123b0489f692f Mon Sep 17 00:00:00 2001 From: emanuel-schmid Date: Fri, 27 Oct 2023 10:49:27 +0200 Subject: [PATCH 1/3] temporary solution for #796: pin xyzservices version --- requirements/env_climada.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/requirements/env_climada.yml b/requirements/env_climada.yml index 95221a38d..3b2885cc8 100644 --- a/requirements/env_climada.yml +++ b/requirements/env_climada.yml @@ -40,3 +40,4 @@ dependencies: - xarray>=2023.8 - xlrd>=2.0 - xlsxwriter>=3.1 + - xyzservices<2023.10.1 From 4be3ab576038340ea40ebe8df8b9582fde6acc2c Mon Sep 17 00:00:00 2001 From: emanuel-schmid Date: Fri, 27 Oct 2023 14:21:16 +0200 Subject: [PATCH 2/3] temporary solution for #796: pin xyzservices version --- requirements/env_climada.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements/env_climada.yml b/requirements/env_climada.yml index 3b2885cc8..1fad27e0b 100644 --- a/requirements/env_climada.yml +++ b/requirements/env_climada.yml @@ -40,4 +40,4 @@ dependencies: - xarray>=2023.8 - xlrd>=2.0 - xlsxwriter>=3.1 - - xyzservices<2023.10.1 + - xyzservices=2023.10.0 From f362c2bf8aa63578424b48cb9e77b298f4ec93b7 Mon Sep 17 00:00:00 2001 From: Thomas Vogt <57705593+tovogt@users.noreply.github.com> Date: Mon, 30 Oct 2023 08:53:41 +0100 Subject: [PATCH 3/3] Fix dist_approx with "geosphere" method and close points (#792) * u_coord.dist_approx: fix geosphere with log * Update CHANGELOG --- CHANGELOG.md | 2 ++ climada/util/coordinates.py | 9 ++++--- climada/util/test/test_coordinates.py | 35 ++++++++++++++++++--------- 3 files changed, 32 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 483fe8bc7..36d511ef2 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/climada/util/coordinates.py b/climada/util/coordinates.py index 2699344da..c719b233b 100644 --- a/climada/util/coordinates.py +++ b/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 diff --git a/climada/util/test/test_coordinates.py b/climada/util/test/test_coordinates.py index 46c39456d..4e2a855f4 100644 --- a/climada/util/test/test_coordinates.py +++ b/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