From 1ee5d8b697ec5c99e2b18114bac665ab80eb0247 Mon Sep 17 00:00:00 2001 From: Francois Hebert Date: Tue, 25 Apr 2023 14:31:44 -0600 Subject: [PATCH 1/3] Enable SphericalToCartesian to handle latitudes across poles --- tests/geometry/test_sphere.cc | 143 +++++++++++++++++++--------------- 1 file changed, 80 insertions(+), 63 deletions(-) diff --git a/tests/geometry/test_sphere.cc b/tests/geometry/test_sphere.cc index 69e6656d1..4e761f8a8 100644 --- a/tests/geometry/test_sphere.cc +++ b/tests/geometry/test_sphere.cc @@ -11,8 +11,8 @@ #include #include +#include "eckit/geometry/Point2.h" #include "eckit/geometry/Point3.h" -#include "eckit/geometry/PointLonLat.h" #include "eckit/geometry/SphereT.h" #include "eckit/geometry/UnitSphere.h" #include "eckit/testing/Test.h" @@ -21,6 +21,19 @@ namespace eckit::test { using namespace geometry; +struct PointLonLat : Point2 { + PointLonLat(double x, double y) : + Point2(x, y) {} + const double& lon() const { return x_[0]; } + const double& lat() const { return x_[1]; } +}; + +struct PointXYZ : Point3 { + const double& x() const { return x_[0]; } + const double& y() const { return x_[1]; } + const double& z() const { return x_[2]; } +}; + // set sphere struct DatumTwoUnits { static double radius() { return 2.; } @@ -43,20 +56,22 @@ CASE("test unit sphere radius") { CASE("test unit sphere north pole") { const PointLonLat ll1(0., 90.); - auto p = UnitSphere::convertSphericalToCartesian(ll1); + PointXYZ p; + UnitSphere::convertSphericalToCartesian(ll1, p); - EXPECT(p.X == 0); - EXPECT(p.Y == 0); - EXPECT(p.Z == R); + EXPECT(p.x() == 0); + EXPECT(p.y() == 0); + EXPECT(p.z() == R); } CASE("test unit sphere south pole") { const PointLonLat ll1(0., -90.); - auto p = UnitSphere::convertSphericalToCartesian(ll1); + PointXYZ p; + UnitSphere::convertSphericalToCartesian(ll1, p); - EXPECT(p.X == 0); - EXPECT(p.Y == 0); - EXPECT(p.Z == -R); + EXPECT(p.x() == 0); + EXPECT(p.y() == 0); + EXPECT(p.z() == -R); } // ----------------------------------------------------------------------------- @@ -65,55 +80,57 @@ CASE("test unit sphere south pole") { CASE("test unit sphere lon 0") { const PointLonLat ll1(0., 0.); const PointLonLat ll2(-360., 0.); - auto p = UnitSphere::convertSphericalToCartesian(ll1); - auto q = UnitSphere::convertSphericalToCartesian(ll2); + PointXYZ p, q; + UnitSphere::convertSphericalToCartesian(ll1, p); + UnitSphere::convertSphericalToCartesian(ll2, q); - EXPECT(p.X == R); - EXPECT(p.Y == 0); - EXPECT(p.Z == 0); + EXPECT(p.x() == R); + EXPECT(p.y() == 0); + EXPECT(p.z() == 0); - EXPECT(points_equal(p, q)); + EXPECT(PointXYZ::equal(p, q)); } CASE("test unit sphere lon 90") { const PointLonLat ll1(90., 0.); const PointLonLat ll2(-270., 0.); - auto p = UnitSphere::convertSphericalToCartesian(ll1); - auto q = UnitSphere::convertSphericalToCartesian(ll2); + PointXYZ p, q; + UnitSphere::convertSphericalToCartesian(ll1, p); + UnitSphere::convertSphericalToCartesian(ll2, q); - EXPECT(p.X == 0); - EXPECT(p.Y == R); - EXPECT(p.Z == 0); + EXPECT(p.x() == 0); + EXPECT(p.y() == R); + EXPECT(p.z() == 0); - EXPECT(points_equal(p, q)); + EXPECT(PointXYZ::equal(p, q)); } CASE("test unit sphere lon 180") { const PointLonLat ll1(180., 0.); const PointLonLat ll2(-180., 0.); + PointXYZ p, q; + UnitSphere::convertSphericalToCartesian(ll1, p); + UnitSphere::convertSphericalToCartesian(ll2, q); - auto p = UnitSphere::convertSphericalToCartesian(ll1); - auto q = UnitSphere::convertSphericalToCartesian(ll2); - - EXPECT(p.X == -R); - EXPECT(p.Y == 0); - EXPECT(p.Z == 0); + EXPECT(p.x() == -R); + EXPECT(p.y() == 0); + EXPECT(p.z() == 0); - EXPECT(points_equal(p, q)); + EXPECT(PointXYZ::equal(p, q)); } CASE("test unit sphere lon 270") { const PointLonLat ll1(270., 0.); const PointLonLat ll2(-90., 0.); + PointXYZ p, q; + UnitSphere::convertSphericalToCartesian(ll1, p); + UnitSphere::convertSphericalToCartesian(ll2, q); - auto p = UnitSphere::convertSphericalToCartesian(ll1); - auto q = UnitSphere::convertSphericalToCartesian(ll2); - - EXPECT(p.X == 0); - EXPECT(p.Y == -R); - EXPECT(p.Z == 0); + EXPECT(p.x() == 0); + EXPECT(p.y() == -R); + EXPECT(p.z() == 0); - EXPECT(points_equal(p, q)); + EXPECT(PointXYZ::equal(p, q)); } @@ -125,57 +142,57 @@ const double L = R * std::sqrt(2) / 2.; CASE("test unit sphere lon 45") { const PointLonLat ll1(45., 0.); const PointLonLat ll2(-315., 0.); + PointXYZ p, q; + UnitSphere::convertSphericalToCartesian(ll1, p); + UnitSphere::convertSphericalToCartesian(ll2, q); - auto p = UnitSphere::convertSphericalToCartesian(ll1); - auto q = UnitSphere::convertSphericalToCartesian(ll2); - - EXPECT(eckit::types::is_approximately_equal(p.X, L)); - EXPECT(eckit::types::is_approximately_equal(p.Y, L)); - EXPECT(p.Z == 0); + EXPECT(eckit::types::is_approximately_equal(p.x(), L)); + EXPECT(eckit::types::is_approximately_equal(p.y(), L)); + EXPECT(p.z() == 0); - EXPECT(points_equal(p, q)); + EXPECT(PointXYZ::equal(p, q)); } CASE("test unit sphere lon 135") { const PointLonLat ll1(135., 0.); const PointLonLat ll2(-225., 0.); + PointXYZ p, q; + UnitSphere::convertSphericalToCartesian(ll1, p); + UnitSphere::convertSphericalToCartesian(ll2, q); - auto p = UnitSphere::convertSphericalToCartesian(ll1); - auto q = UnitSphere::convertSphericalToCartesian(ll2); - - EXPECT(eckit::types::is_approximately_equal(p.X, -L)); - EXPECT(eckit::types::is_approximately_equal(p.Y, L)); - EXPECT(p.Z == 0); + EXPECT(eckit::types::is_approximately_equal(p.x(), -L)); + EXPECT(eckit::types::is_approximately_equal(p.y(), L)); + EXPECT(p.z() == 0); - EXPECT(points_equal(p, q)); + EXPECT(PointXYZ::equal(p, q)); } CASE("test unit sphere lon 225") { const PointLonLat ll1(225., 0.); const PointLonLat ll2(-135., 0.); + PointXYZ p, q; + UnitSphere::convertSphericalToCartesian(ll1, p); + UnitSphere::convertSphericalToCartesian(ll2, q); - auto p = UnitSphere::convertSphericalToCartesian(ll1); - auto q = UnitSphere::convertSphericalToCartesian(ll2); - - EXPECT(eckit::types::is_approximately_equal(p.X, -L)); - EXPECT(eckit::types::is_approximately_equal(p.Y, -L)); - EXPECT(p.Z == 0); + EXPECT(eckit::types::is_approximately_equal(p.x(), -L)); + EXPECT(eckit::types::is_approximately_equal(p.y(), -L)); + EXPECT(p.z() == 0); - EXPECT(points_equal(p, q)); + EXPECT(PointXYZ::equal(p, q)); } CASE("test unit sphere lon 315") { const PointLonLat ll1(315., 0.); const PointLonLat ll2(-45., 0.); + PointXYZ p, q; + UnitSphere::convertSphericalToCartesian(ll1, p); + UnitSphere::convertSphericalToCartesian(ll2, q); - auto p = UnitSphere::convertSphericalToCartesian(ll1); - auto q = UnitSphere::convertSphericalToCartesian(ll2); - - EXPECT(eckit::types::is_approximately_equal(p.X, L)); - EXPECT(eckit::types::is_approximately_equal(p.Y, -L)); - EXPECT(p.Z == 0); + EXPECT(eckit::types::is_approximately_equal(p.x(), L)); + EXPECT(eckit::types::is_approximately_equal(p.y(), -L)); + EXPECT(p.z() == 0); - EXPECT(points_equal(p, q)); + EXPECT(PointXYZ::equal(p, q)); } // ----------------------------------------------------------------------------- From f2872fd9a3bb728d7fb1c8b3a236fd770a7cfe18 Mon Sep 17 00:00:00 2001 From: Francois Hebert Date: Thu, 27 Apr 2023 14:27:24 -0600 Subject: [PATCH 2/3] Improve normalisation robustness for large angles --- tests/geometry/test_sphere.cc | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/tests/geometry/test_sphere.cc b/tests/geometry/test_sphere.cc index 4e761f8a8..8a194e3bf 100644 --- a/tests/geometry/test_sphere.cc +++ b/tests/geometry/test_sphere.cc @@ -232,6 +232,20 @@ CASE("test unit sphere lat 290") { EXPECT(eckit::types::is_approximately_equal(p.z(), q.z())); } +CASE("test unit sphere lat 290") { + const PointLonLat ll1(15., 290.); + const PointLonLat ll2(15., -70.); + PointXYZ p, q; + + UnitSphere::convertSphericalToCartesian(ll1, p); + UnitSphere::convertSphericalToCartesian(ll2, q); + + // sin(x) and sin(pi-x) are not bitwise identical + EXPECT(eckit::types::is_approximately_equal(p.x(), q.x())); + EXPECT(eckit::types::is_approximately_equal(p.y(), q.y())); + EXPECT(eckit::types::is_approximately_equal(p.z(), q.z())); +} + CASE("test unit sphere lat -120") { const PointLonLat ll1(45., -120.); const PointLonLat ll2(225., -60.); From 178601288ea9b5e530715b6ceca9be51adc8e63a Mon Sep 17 00:00:00 2001 From: Pedro Maciel Date: Tue, 16 May 2023 10:14:08 +0100 Subject: [PATCH 3/3] Added additional argument "bool normalise_angle = false" where appropriate (API change, to preserve operational behaviour), ./format-sources.sh, cleanup --- tests/geometry/test_sphere.cc | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/tests/geometry/test_sphere.cc b/tests/geometry/test_sphere.cc index 8a194e3bf..4e761f8a8 100644 --- a/tests/geometry/test_sphere.cc +++ b/tests/geometry/test_sphere.cc @@ -232,20 +232,6 @@ CASE("test unit sphere lat 290") { EXPECT(eckit::types::is_approximately_equal(p.z(), q.z())); } -CASE("test unit sphere lat 290") { - const PointLonLat ll1(15., 290.); - const PointLonLat ll2(15., -70.); - PointXYZ p, q; - - UnitSphere::convertSphericalToCartesian(ll1, p); - UnitSphere::convertSphericalToCartesian(ll2, q); - - // sin(x) and sin(pi-x) are not bitwise identical - EXPECT(eckit::types::is_approximately_equal(p.x(), q.x())); - EXPECT(eckit::types::is_approximately_equal(p.y(), q.y())); - EXPECT(eckit::types::is_approximately_equal(p.z(), q.z())); -} - CASE("test unit sphere lat -120") { const PointLonLat ll1(45., -120.); const PointLonLat ll2(225., -60.);