Enable convertSphericalToCartesian to handle latitudes across poles

src/eckit/geometry/CMakeLists.txt

@@ -1,4 +1,6 @@
 list( APPEND eckit_geometry_srcs
+CoordinateHelpers.cc
+CoordinateHelpers.h
 EllipsoidOfRevolution.cc
 EllipsoidOfRevolution.h
 GreatCircle.cc

src/eckit/geometry/CoordinateHelpers.cc

@@ -0,0 +1,42 @@
+/*
+ * (C) Copyright 2023 UCAR
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ */
+
+#include "eckit/geometry/CoordinateHelpers.h"
+#include "eckit/geometry/Point2.h"
+
+namespace eckit {
+namespace geometry {
+
+//----------------------------------------------------------------------------------------------------------------------
+
+double normalise_angle(double a, const double minimum) {
+    while (a < minimum) {
+        a += 360.;
+    }
+    while (a >= minimum + 360.) {
+        a -= 360.;
+    }
+    return a;
+}
+
+//----------------------------------------------------------------------------------------------------------------------
+
+Point2 canonicaliseOnSphere(const Point2& lonlat, const double minimum_lon) {
+    const double lat = normalise_angle(lonlat[1], -90.);
+    const bool across_pole = (lat > 90.);
+
+    if (!across_pole) {
+        return Point2(normalise_angle(lonlat[0], minimum_lon), lat);
+    } else {
+        return Point2(normalise_angle(lonlat[0] + 180., minimum_lon), 180. - lat);
+    }
+}
+
+//----------------------------------------------------------------------------------------------------------------------
+
+}  // namespace geometry
+}  // namespace eckit

src/eckit/geometry/CoordinateHelpers.h

@@ -0,0 +1,43 @@
+/*
+ * (C) Copyright 2023 UCAR.
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ */
+
+#ifndef CoordinateHelpers_H
+#define CoordinateHelpers_H
+
+namespace eckit {
+namespace geometry {
+
+class Point2;
+
+//------------------------------------------------------------------------------------------------------
+
+/// Shift angle in increments of 360° until it lies in [minimum, minimum+360°).
+///
+/// Inputs angle and minimum are in degrees, returned angle is in degrees.
+double normalise_angle(double angle, double minimum);
+
+//------------------------------------------------------------------------------------------------------
+
+/// Shift input point on sphere so its longitude lies in [minimum_lon, minimum_lon+360°)
+/// and its latitude lies in [-90°, 90°].
+///
+/// Latitudes outside the canonical interval [-90°,90°] are first shifted into the interval
+/// [-90°,270°], then any points with latitudes in [90°,270°] are flagged as "across the pole".
+/// Such points are re-labeled with equivalent coordinates that lie within the canonical coordinate
+/// patch by the transformation: (λ, ϕ) -> (λ+180°, 180°-ϕ).
+///
+/// Finally, the longitude is shifted into [minimum_lon, minimum_lon+360°).
+///
+/// Inputs lonlat and minimum_lon are in degrees, returned angles are in degrees.
+Point2 canonicaliseOnSphere(const Point2& lonlat, double minimum_lon = 0.);
+
+//------------------------------------------------------------------------------------------------------
+
+}  // namespace geometry
+}  // namespace eckit
+
+#endif  // CoordinateHelpers_H

src/eckit/geometry/EllipsoidOfRevolution.cc

@@ -16,6 +16,7 @@
 #include <sstream>
 
 #include "eckit/exception/Exceptions.h"
+#include "eckit/geometry/CoordinateHelpers.h"
 #include "eckit/geometry/Point2.h"
 #include "eckit/geometry/Point3.h"
 
@@ -27,16 +28,6 @@ namespace eckit::geometry {
 
 namespace {
 
-static double normalise_longitude(double a, const double& minimum) {
-    while (a < minimum) {
-        a += 360;
-    }
-    while (a >= minimum + 360) {
-        a -= 360;
-    }
-    return a;
-}
-
 static const double degrees_to_radians = M_PI / 180.;
 
 static std::streamsize max_digits10 = 15 + 3;
@@ -52,19 +43,14 @@ void EllipsoidOfRevolution::convertSphericalToCartesian(const double& a, const d
     ASSERT(a > 0.);
     ASSERT(b > 0.);
 
-    if (!(-90. <= Alonlat[1] && Alonlat[1] <= 90.)) {
-        std::ostringstream oss;
-        oss.precision(max_digits10);
-        oss << "Invalid latitude " << Alonlat[1];
-        throw BadValue(oss.str(), Here());
-    }
-
     // See https://en.wikipedia.org/wiki/Reference_ellipsoid#Coordinates
     // numerical conditioning for both ϕ (poles) and λ (Greenwich/Date Line)
 
-    const double lambda_deg = normalise_longitude(Alonlat[0], -180.);
+    const Point2 alonlat = canonicaliseOnSphere(Alonlat, -180.);
+
+    const double lambda_deg = alonlat[0];
     const double lambda     = degrees_to_radians * lambda_deg;
-    const double phi        = degrees_to_radians * Alonlat[1];
+    const double phi        = degrees_to_radians * alonlat[1];
 
     const double sin_phi    = std::sin(phi);
     const double cos_phi    = std::sqrt(1. - sin_phi * sin_phi);

src/eckit/geometry/Sphere.cc

@@ -17,6 +17,7 @@
 #include <sstream>
 
 #include "eckit/exception/Exceptions.h"
+#include "eckit/geometry/CoordinateHelpers.h"
 #include "eckit/geometry/GreatCircle.h"
 #include "eckit/geometry/Point2.h"
 #include "eckit/geometry/Point3.h"
@@ -28,16 +29,6 @@ namespace eckit::geometry {
 
 //----------------------------------------------------------------------------------------------------------------------
 
-static double normalise_longitude(double a, const double& minimum) {
-    while (a < minimum) {
-        a += 360;
-    }
-    while (a >= minimum + 360) {
-        a -= 360;
-    }
-    return a;
-}
-
 static const double radians_to_degrees = 180. * M_1_PI;
 
 static const double degrees_to_radians = M_PI / 180.;
@@ -69,23 +60,12 @@ double Sphere::centralAngle(const Point2& Alonlat, const Point2& Blonlat) {
      * }
      */
 
-    if (!(-90. <= Alonlat[1] && Alonlat[1] <= 90.)) {
-        ostringstream oss;
-        oss.precision(max_digits10);
-        oss << "Invalid latitude " << Alonlat[1];
-        throw BadValue(oss.str(), Here());
-    }
-
-    if (!(-90. <= Blonlat[1] && Blonlat[1] <= 90.)) {
-        ostringstream oss;
-        oss.precision(max_digits10);
-        oss << "Invalid latitude " << Blonlat[1];
-        throw BadValue(oss.str(), Here());
-    }
+    const Point2 alonlat = canonicaliseOnSphere(Alonlat);
+    const Point2 blonlat = canonicaliseOnSphere(Blonlat);
 
-    const double phi1   = degrees_to_radians * Alonlat[1];
-    const double phi2   = degrees_to_radians * Blonlat[1];
-    const double lambda = degrees_to_radians * (Blonlat[0] - Alonlat[0]);
+    const double phi1   = degrees_to_radians * alonlat[1];
+    const double phi2   = degrees_to_radians * blonlat[1];
+    const double lambda = degrees_to_radians * (blonlat[0] - alonlat[0]);
 
     const double cos_phi1   = cos(phi1);
     const double sin_phi1   = sin(phi1);
@@ -139,7 +119,7 @@ double Sphere::area(const double& radius, const Point2& WestNorth, const Point2&
 
     // Set longitude fraction
     double W = WestNorth[0];
-    double E = normalise_longitude(EastSouth[0], W);
+    double E = normalise_angle(EastSouth[0], W);
     double longitude_range(
         types::is_approximately_equal(W, E) && !types::is_approximately_equal(EastSouth[0], WestNorth[0]) ? 360.
                                                                                                           : E - W);
@@ -178,13 +158,6 @@ void Sphere::greatCircleLongitudeGivenLatitude(const Point2& Alonlat, const Poin
 void Sphere::convertSphericalToCartesian(const double& radius, const Point2& Alonlat, Point3& B, double height) {
     ASSERT(radius > 0.);
 
-    if (!(-90. <= Alonlat[1] && Alonlat[1] <= 90.)) {
-        std::ostringstream oss;
-        oss.precision(max_digits10);
-        oss << "Invalid latitude " << Alonlat[1];
-        throw BadValue(oss.str(), Here());
-    }
-
     /*
      * See https://en.wikipedia.org/wiki/Reference_ellipsoid#Coordinates
      * numerical conditioning for both ϕ (poles) and λ (Greenwich/Date Line).
@@ -198,9 +171,11 @@ void Sphere::convertSphericalToCartesian(const double& radius, const Point2& Alo
      * poles and quadrants.
      */
 
-    const double lambda_deg = normalise_longitude(Alonlat[0], -180.);
+    const Point2 alonlat = canonicaliseOnSphere(Alonlat, -180.);
+
+    const double lambda_deg = alonlat[0];
     const double lambda     = degrees_to_radians * lambda_deg;
-    const double phi        = degrees_to_radians * Alonlat[1];
+    const double phi        = degrees_to_radians * alonlat[1];
 
     const double sin_phi    = std::sin(phi);
     const double cos_phi    = std::sqrt(1. - sin_phi * sin_phi);

src/eckit/geometry/polygon/LonLatPolygon.cc

@@ -14,6 +14,7 @@
 #include <ostream>
 
 #include "eckit/exception/Exceptions.h"
+#include "eckit/geometry/CoordinateHelpers.h"
 #include "eckit/types/FloatCompare.h"
 
 //----------------------------------------------------------------------------------------------------------------------
@@ -106,16 +107,10 @@ std::ostream& operator<<(std::ostream& out, const LonLatPolygon& pc) {
 }
 
 bool LonLatPolygon::contains(const Point2& P) const {
-    auto lat = P[LAT];
-    ASSERT(-90 <= lat && lat <= 90);
 
-    auto lon = P[LON];
-    while (lon >= min_[LON] + 360) {
-        lon -= 360;
-    }
-    while (lon < min_[LON]) {
-        lon += 360;
-    }
+    const Point2 p = canonicaliseOnSphere(P, min_[LON]);
+    auto lat = p[LAT];
+    auto lon = p[LON];
 
     // check poles
     if (includeNorthPole_ && is_approximately_equal(lat, 90)) {

tests/geometry/CMakeLists.txt

@@ -1,4 +1,4 @@
-foreach( _test great_circle kdtree sphere kpoint points polygon )
+foreach( _test coordinate_helpers great_circle kdtree sphere kpoint points polygon )
     ecbuild_add_test( TARGET  eckit_test_geometry_${_test}
                       SOURCES test_${_test}.cc
                       LIBS    eckit_geometry )

tests/geometry/test_coordinate_helpers.cc

@@ -0,0 +1,70 @@
+/*
+ * (C) Copyright 2023 UCAR
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ */
+
+#include "eckit/geometry/CoordinateHelpers.h"
+#include "eckit/geometry/Point2.h"
+#include "eckit/testing/Test.h"
+
+
+namespace eckit {
+namespace test {
+
+using namespace geometry;
+
+// -----------------------------------------------------------------------------
+
+CASE("normalise angles") {
+    EXPECT(0. == normalise_angle(360., 0.));
+    EXPECT(14. == normalise_angle(374., 0.));
+    EXPECT(14. == normalise_angle(374., -90.));
+    EXPECT(374. == normalise_angle(374., 90.));
+}
+
+CASE("canonicalise on sphere") {
+    using types::is_approximately_equal;
+
+    const Point2 p1(108., 32.);
+
+    // Worse coordinates for the same point:
+    const Point2 p2(-252., 32.);
+    const Point2 p3(288., 148.);
+    const Point2 p4(108., -328.);
+
+    // Check each of these is correctly shifted back to original point:
+    const Point2 q2 = canonicaliseOnSphere(p2);
+    const Point2 q3 = canonicaliseOnSphere(p3);
+    const Point2 q4 = canonicaliseOnSphere(p4);
+
+    EXPECT(p1.x() == q2.x());
+    EXPECT(p1.y() == q2.y());
+    EXPECT(p1.x() == q3.x());
+    EXPECT(p1.y() == q3.y());
+    EXPECT(p1.x() == q4.x());
+    EXPECT(p1.y() == q4.y());
+
+    // Check with longitude offset
+    const Point2 q2b = canonicaliseOnSphere(p2, -90.);
+    EXPECT(q2b.x() == 108.);
+    EXPECT(q2b.y() == 32.);
+
+    const Point2 q2c = canonicaliseOnSphere(p2, 90.);
+    EXPECT(q2c.x() == 108.);
+    EXPECT(q2c.y() == 32.);
+
+    const Point2 q2d = canonicaliseOnSphere(p2, 180.);
+    EXPECT(q2d.x() == 468.);
+    EXPECT(q2d.y() == 32.);
+}
+
+// -----------------------------------------------------------------------------
+
+}  // namespace test
+}  // namespace eckit
+
+int main(int argc, char** argv) {
+    return eckit::testing::run_tests(argc, argv);
+}

tests/geometry/test_polygon.cc

@@ -264,6 +264,10 @@ CASE("LonLatPolygon") {
                 EXPECT_NOT(poly.contains({lonmax + eps, lat}));
             }
         }
+
+        // Test points at non-canonical coordinates
+        EXPECT(poly.contains({lonmid + 360., latmid}));
+        EXPECT(poly.contains({lonmid, 180. - latmid}));
     }
 
     SECTION("Parallelogram") {