From b06e7aae539ad5b6491ea999794e6ab6a5e62b17 Mon Sep 17 00:00:00 2001 From: Thomas Nipen Date: Wed, 10 Jan 2024 09:18:23 +0100 Subject: [PATCH] Fixed bug in optimal_interpolation across dateline Due to a bug in gridpp::KDTree::calc_distance_fast, going across the dateline (giong from -180 to +180 longitude) resulted in correlations of 0, since the distance was estimated to be large. Thanks to @joewkr for reporting the problem. --- include/gridpp.h | 1 + src/api/kdtree.cpp | 12 ++++++-- tests/kdtree_test.py | 72 ++++++++++++++++++++++++++++++++++++++------ 3 files changed, 72 insertions(+), 13 deletions(-) diff --git a/include/gridpp.h b/include/gridpp.h index c20d8da9..66158f63 100644 --- a/include/gridpp.h +++ b/include/gridpp.h @@ -1584,6 +1584,7 @@ namespace gridpp { static float rad2deg(float deg); static float calc_distance(float lat1, float lon1, float lat2, float lon2, CoordinateType type=Geodetic); static float calc_distance(float x0, float y0, float z0, float x1, float y1, float z1); + static float calc_distance(const Point& p1, const Point& p2); static float calc_distance_fast(float lat1, float lon1, float lat2, float lon2, CoordinateType type=Geodetic); static float calc_distance_fast(const Point& p1, const Point& p2); vec get_lats() const; diff --git a/src/api/kdtree.cpp b/src/api/kdtree.cpp index 916d60b3..4e6bbccb 100644 --- a/src/api/kdtree.cpp +++ b/src/api/kdtree.cpp @@ -176,14 +176,20 @@ float gridpp::KDTree::calc_distance_fast(float lat1, float lon1, float lat2, flo double lat2r = deg2rad(lat2); double lon1r = deg2rad(lon1); double lon2r = deg2rad(lon2); - float dx2 = pow(cos((lat1r+lat2r)/2),2)*(lon1r-lon2r)*(lon1r-lon2r); - float dy2 = (lat1r-lat2r)*(lat1r-lat2r); - return gridpp::radius_earth*sqrt(dx2+dy2); + double dlon = fmod(fabs(lon1r - lon2r), 2 * M_PI); + double mean_lat = (lat1r + lat2r) / 2; + float dx2 = pow(cos(mean_lat), 2) * dlon * dlon; + float dy2 = (lat1r - lat2r) * (lat1r - lat2r); + return gridpp::radius_earth * sqrt(dx2 + dy2); } else { throw std::runtime_error("Unknown coordinate type"); } } +float gridpp::KDTree::calc_distance(const Point& p1, const Point& p2) { + assert(p1.type == p2.type); + return calc_distance(p1.lat, p1.lon, p2.lat, p2.lon, p1.type); +} float gridpp::KDTree::calc_distance_fast(const Point& p1, const Point& p2) { assert(p1.type == p2.type); return calc_distance_fast(p1.lat, p1.lon, p2.lat, p2.lon, p1.type); diff --git a/tests/kdtree_test.py b/tests/kdtree_test.py index 7e83512e..ab5fe570 100644 --- a/tests/kdtree_test.py +++ b/tests/kdtree_test.py @@ -58,18 +58,70 @@ def test_rad2deg(self): self.assertAlmostEqual(gridpp.KDTree_rad2deg(0), 0, 5) def test_calc_distance(self): - self.assertAlmostEqual(0, gridpp.KDTree_calc_distance(60,10,60,10)); - self.assertAlmostEqual(20037508, gridpp.KDTree_calc_distance(90,10,-90,10)); - self.assertAlmostEqual(20037508, gridpp.KDTree_calc_distance(0,0,0,180)); - self.assertAlmostEqual(16879114, gridpp.KDTree_calc_distance(60.5,5.25,-84.75,-101.75)); - self.assertAlmostEqual(124080.79, gridpp.KDTree_calc_distance(60,10,61,11), delta=0.1); + config = list() + # lat0, lon0, lat1, lon1, delta, expected + config += [[60,10,61,11,0.1,124080.79]] + config += [[60, 10, 60, 10, 0, 0]] + config += [[90, 10, -90, 10, 0, 20037508]] + config += [[0, 0, 0, 180, 0, 20037508]] + config += [[60.5, 5.25, -84.75, -101.75, 0, 16879114]] + config += [[60, 10, 61, 11, 0.1, 124080.79]] + + for c in config: + self.assertEqual(len(c), 6) + p0 = gridpp.Point(c[0], c[1]) + p1 = gridpp.Point(c[2], c[3]) + delta = c[4] + expected = c[5] + + # Point version + with self.subTest(config=c, type="Point"): + self.assertAlmostEqual(expected, gridpp.KDTree.calc_distance(p0, p1), delta=delta) + + # Scalar version + with self.subTest(config=c, type="Scalar"): + self.assertAlmostEqual(expected, gridpp.KDTree.calc_distance(c[0], c[1], c[2], c[3]), delta=delta) def test_calc_distance_fast(self): - self.assertEqual(0, gridpp.KDTree_calc_distance_fast(60,10,60,10)); - self.assertAlmostEqual(20037508, gridpp.KDTree_calc_distance_fast(90,10,-90,10), delta=100); - self.assertAlmostEqual(20037508, gridpp.KDTree_calc_distance_fast(0,0,0,180), delta=100); - # self.assertAlmostEqual(16879114, gridpp.KDTree_calc_distance_fast(60.5,5.25,-84.75,-101.75), delta=100); - self.assertAlmostEqual(124080.79, gridpp.KDTree_calc_distance_fast(60,10,61,11), delta=100); + config = list() + # lat0,lon0,lat1,lon1, delta, dist + config += [[60,10, 60,10, 10, 0]] + config += [[90,10, -90,10, 10, 20037508]] + config += [[0,0, 0,180, 10, 20037508]] + config += [[60,10, 61,11, 10, 124080.79]] + + for c in config: + self.assertEqual(len(c), 6) + p0 = gridpp.Point(c[0], c[1]) + p1 = gridpp.Point(c[2], c[3]) + delta = c[4] + expected = c[5] + + # Point version + self.assertAlmostEqual(expected, gridpp.KDTree.calc_distance_fast(p0, p1), delta=delta) + + # Scalar version + self.assertAlmostEqual(expected, gridpp.KDTree.calc_distance_fast(c[0], c[1], c[2], c[3]), delta=delta) + + def test_calc_distance_fast_across_date_line(self): + config = list() + for lat in [-90, -89, 0, 89, 90]: + config += [[lat, 180, lat, -180, 10, 0]] + for lat in [-90, -89, 0, 89]: + config += [[lat, 180, lat+1, -180, 10, 111319.4921875]] + + for c in config: + self.assertEqual(len(c), 6) + p0 = gridpp.Point(c[0], c[1]) + p1 = gridpp.Point(c[2], c[3]) + delta = c[4] + expected = c[5] + + # Point version + self.assertAlmostEqual(expected, gridpp.KDTree.calc_distance_fast(p0, p1), delta=delta) + + # Scalar version + self.assertAlmostEqual(expected, gridpp.KDTree.calc_distance_fast(c[0], c[1], c[2], c[3]), delta=delta) def test_radius_match(self): """Check that points right on the radius edge count as a match"""