diff --git a/src/atlas/util/Geometry.h b/src/atlas/util/Geometry.h index 90b868ce8..71bd219e6 100644 --- a/src/atlas/util/Geometry.h +++ b/src/atlas/util/Geometry.h @@ -71,7 +71,11 @@ template class GeometrySphereT : public GeometryBase { public: void lonlat2xyz(const Point2& lonlat, Point3& xyz) const override { +#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 24, 0) + SphereT::convertSphericalToCartesian(lonlat, xyz, 0., true); +#else SphereT::convertSphericalToCartesian(lonlat, xyz); +#endif } void xyz2lonlat(const Point3& xyz, Point2& lonlat) const override { SphereT::convertCartesianToSpherical(xyz, lonlat); diff --git a/src/tests/AtlasTestEnvironment.h b/src/tests/AtlasTestEnvironment.h index 503fa43e4..43cb0238d 100644 --- a/src/tests/AtlasTestEnvironment.h +++ b/src/tests/AtlasTestEnvironment.h @@ -217,6 +217,12 @@ bool approx_eq(const Point2& v1, const Point2& v2) { bool approx_eq(const Point2& v1, const Point2& v2, const double& t) { return approx_eq(v1[0], v2[0], t) && approx_eq(v1[1], v2[1], t); } +bool approx_eq(const Point3& v1, const Point3& v2) { + return approx_eq(v1[0], v2[0]) && approx_eq(v1[1], v2[1]) && approx_eq(v1[2], v2[2]); +} +bool approx_eq(const Point3& v1, const Point3& v2, const double& t) { + return approx_eq(v1[0], v2[0], t) && approx_eq(v1[1], v2[1], t) && approx_eq(v1[2], v2[2], t); +} template std::string expect_message(const std::string& condition, const T1& lhs, const T2& rhs, const eckit::CodeLocation& loc) { diff --git a/src/tests/util/test_earth.cc b/src/tests/util/test_earth.cc index 0551acdeb..3cb86a5d6 100644 --- a/src/tests/util/test_earth.cc +++ b/src/tests/util/test_earth.cc @@ -13,6 +13,7 @@ #include "atlas/util/Earth.h" #include "atlas/util/Point.h" +#include "atlas/util/Geometry.h" #include "tests/AtlasTestEnvironment.h" @@ -169,6 +170,13 @@ CASE("test_earth_great_circle_latitude_given_longitude") { EXPECT(eckit::types::is_approximately_equal(lat, -6.81, 0.01)); } +CASE("test_across_pole") { + const PointLonLat P1{-16.875,-105.255}; + atlas::geometry::Earth geo; + EXPECT_APPROX_EQ(geo.xyz(P1), PointXYZ(-1.60418e+06,486624,-6.14673e+06), 50.); +} + + } // namespace test } // namespace atlas