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,55 @@
+/*
+ * (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 <limits>
+#include <sstream>
+
+#include "eckit/exception/Exceptions.h"
+#include "eckit/geometry/CoordinateHelpers.h"
+#include "eckit/geometry/Point2.h"
+
+namespace eckit::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 {normalise_angle(lonlat[0], minimum_lon), lat};
+    }
+
+    return {normalise_angle(lonlat[0] + 180., minimum_lon), 180. - lat};
+}
+
+//----------------------------------------------------------------------------------------------------------------------
+
+void assert_latitude_range(double lat) {
+    if (!(-90. <= lat && lat <= 90.)) {
+        std::ostringstream oss;
+        oss.precision(std::numeric_limits<double>::max_digits10);
+        oss << "Invalid latitude " << lat;
+        throw BadValue(oss.str(), Here());
+    }
+}
+
+//----------------------------------------------------------------------------------------------------------------------
+
+}  // namespace eckit::geometry

src/eckit/geometry/CoordinateHelpers.h

@@ -0,0 +1,47 @@
+/*
+ * (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::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.);
+
+//----------------------------------------------------------------------------------------------------------------------
+
+/// Assert latitude lies in [-90°, 90°].
+void assert_latitude_range(double lat);
+
+
+//----------------------------------------------------------------------------------------------------------------------
+
+}  // namespace eckit::geometry
+
+#endif  // CoordinateHelpers_H

src/eckit/geometry/EllipsoidOfRevolution.cc

@@ -11,11 +11,9 @@
 #include "eckit/geometry/EllipsoidOfRevolution.h"
 
 #include <cmath>
-#include <ios>
-// #include <limits>  // for std::numeric_limits
-#include <sstream>
 
 #include "eckit/exception/Exceptions.h"
+#include "eckit/geometry/CoordinateHelpers.h"
 #include "eckit/geometry/Point2.h"
 #include "eckit/geometry/Point3.h"
 
@@ -27,44 +25,33 @@ 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;
-
-// C++-11: std::numeric_limits<double>::max_digits10;
-
 }  // namespace
 
 //----------------------------------------------------------------------------------------------------------------------
 
-void EllipsoidOfRevolution::convertSphericalToCartesian(const double& a, const double& b, const Point2& Alonlat,
-                                                        Point3& B, double height) {
+void EllipsoidOfRevolution::convertSphericalToCartesian(const double& a,
+                                                        const double& b,
+                                                        const Point2& Alonlat,
+                                                        Point3& B,
+                                                        double height,
+                                                        bool normalise_angle) {
     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.);
+    if (!normalise_angle) {
+        assert_latitude_range(Alonlat[1]);
+    }
+
+    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/EllipsoidOfRevolution.h

@@ -11,24 +11,28 @@
 #ifndef EllipsoidOfRevolution_H
 #define EllipsoidOfRevolution_H
 
-//------------------------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------------------------------
 
 namespace eckit::geometry {
 
-//------------------------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------------------------------
 
 class Point2;
 class Point3;
 
-//------------------------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------------------------------
 
 struct EllipsoidOfRevolution {
     // Convert elliptic coordinates to Cartesian
-    static void convertSphericalToCartesian(const double& radiusA, const double& radiusB, const Point2& Alonlat,
-                                            Point3& B, double height = 0.);
+    static void convertSphericalToCartesian(const double& radiusA,
+                                            const double& radiusB,
+                                            const Point2& Alonlat,
+                                            Point3& B,
+                                            double height        = 0.,
+                                            bool normalise_angle = false);
 };
 
-//------------------------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------------------------------------
 
 }  // namespace eckit::geometry
 

src/eckit/geometry/Sphere.cc

@@ -12,11 +12,10 @@
 
 #include <algorithm>
 #include <cmath>
-#include <ios>
 #include <limits>
-#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,31 +27,17 @@ 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.;
 
-static constexpr std::streamsize max_digits10 = std::numeric_limits<double>::max_digits10;
-
 inline double squared(double x) {
     return x * x;
 }
 
 //----------------------------------------------------------------------------------------------------------------------
 
-double Sphere::centralAngle(const Point2& Alonlat, const Point2& Blonlat) {
-    using namespace std;
-
+double Sphere::centralAngle(const Point2& Alonlat, const Point2& Blonlat, bool normalise_angle) {
     /*
      * Δσ = atan( ((cos(ϕ2) * sin(Δλ))^2 + (cos(ϕ1) * sin(ϕ2) - sin(ϕ1) * cos(ϕ2) * cos(Δλ))^2) /
      *            (sin(ϕ1) * sin(ϕ2) + cos(ϕ1) * cos(ϕ2) * cos(Δλ)) )
@@ -69,33 +54,29 @@ 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 (!normalise_angle) {
+        assert_latitude_range(Alonlat[1]);
+        assert_latitude_range(Blonlat[1]);
     }
 
-    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);
-    const double cos_phi2   = cos(phi2);
-    const double sin_phi2   = sin(phi2);
-    const double cos_lambda = cos(lambda);
-    const double sin_lambda = sin(lambda);
+    const double cos_phi1   = std::cos(phi1);
+    const double sin_phi1   = std::sin(phi1);
+    const double cos_phi2   = std::cos(phi2);
+    const double sin_phi2   = std::sin(phi2);
+    const double cos_lambda = std::cos(lambda);
+    const double sin_lambda = std::sin(lambda);
 
-    const double angle = atan2(sqrt(squared(cos_phi2 * sin_lambda) + squared(cos_phi1 * sin_phi2 - sin_phi1 * cos_phi2 * cos_lambda)),
-                               sin_phi1 * sin_phi2 + cos_phi1 * cos_phi2 * cos_lambda);
+    const double angle = std::atan2(std::sqrt(squared(cos_phi2 * sin_lambda)
+                                              + squared(cos_phi1 * sin_phi2
+                                                        - sin_phi1 * cos_phi2 * cos_lambda)),
+                                    sin_phi1 * sin_phi2 + cos_phi1 * cos_phi2 * cos_lambda);
 
     if (types::is_approximately_equal(angle, 0.)) {
         return 0.;
@@ -122,7 +103,7 @@ double Sphere::centralAngle(const double& radius, const Point3& A, const Point3&
 }
 
 double Sphere::distance(const double& radius, const Point2& Alonlat, const Point2& Blonlat) {
-    return radius * centralAngle(Alonlat, Blonlat);
+    return radius * centralAngle(Alonlat, Blonlat, true);
 }
 
 double Sphere::distance(const double& radius, const Point3& A, const Point3& B) {
@@ -139,7 +120,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);
@@ -175,16 +156,10 @@ void Sphere::greatCircleLongitudeGivenLatitude(const Point2& Alonlat, const Poin
     Clon2 = lon.size() > 1 ? lon[1] : std::numeric_limits<double>::signaling_NaN();
 }
 
-void Sphere::convertSphericalToCartesian(const double& radius, const Point2& Alonlat, Point3& B, double height) {
+void Sphere::convertSphericalToCartesian(
+    const double& radius, const Point2& Alonlat, Point3& B, double height, bool normalise_angle) {
     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 +173,15 @@ void Sphere::convertSphericalToCartesian(const double& radius, const Point2& Alo
      * poles and quadrants.
      */
 
-    const double lambda_deg = normalise_longitude(Alonlat[0], -180.);
+    if (!normalise_angle) {
+        assert_latitude_range(Alonlat[1]);
+    }
+
+    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);