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

@@ -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;

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);

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