diff --git a/src/eckit/geo/Projection.h b/src/eckit/geo/Projection.h index b5ca1fc66..daeb26948 100644 --- a/src/eckit/geo/Projection.h +++ b/src/eckit/geo/Projection.h @@ -21,7 +21,8 @@ namespace eckit { class Configuration; -} +class MappedConfiguration; +} // namespace eckit namespace eckit::geo { @@ -34,6 +35,8 @@ class Projection { using builder_t = BuilderT1; using ARG1 = const Configuration&; + using Spec = MappedConfiguration; + // -- Exceptions // None @@ -57,10 +60,12 @@ class Projection { // -- Methods + static std::string className() { return "projection"; } + virtual Point fwd(const Point&) const = 0; virtual Point inv(const Point&) const = 0; - static std::string className() { return "projection"; } + virtual Spec spec() const = 0; // -- Overridden methods // None diff --git a/src/eckit/geo/projection/LonLatToXYZ.cc b/src/eckit/geo/projection/LonLatToXYZ.cc index e0db54f83..ac8c30497 100644 --- a/src/eckit/geo/projection/LonLatToXYZ.cc +++ b/src/eckit/geo/projection/LonLatToXYZ.cc @@ -12,7 +12,7 @@ #include "eckit/geo/projection/LonLatToXYZ.h" -#include "eckit/config/Configuration.h" +#include "eckit/config/MappedConfiguration.h" #include "eckit/geo/EllipsoidOfRevolution.h" #include "eckit/geo/Sphere.h" #include "eckit/types/FloatCompare.h" @@ -25,8 +25,7 @@ static ProjectionBuilder __projection("ll_to_xyz"); LonLatToXYZ::LonLatToXYZ(double a, double b) { - ASSERT(0. < a); - ASSERT(0. < b); + ASSERT_MSG(types::is_strictly_greater(b, 0.) && b <= a, "LonLatToXYZ requires 0 < b <= a"); struct LonLatToSphereXYZ final : Implementation { using S = Sphere; @@ -36,6 +35,7 @@ LonLatToXYZ::LonLatToXYZ(double a, double b) { R_(R) {} Point3 operator()(const PointLonLat& p) const override { return S::convertSphericalToCartesian(R_, p, 0.); } PointLonLat operator()(const Point3& q) const override { return S::convertCartesianToSpherical(R_, q); } + Spec spec() const override { return Spec{{{"R", R_}}}; } }; struct LonLatToSpheroidXYZ final : Implementation { @@ -47,6 +47,7 @@ LonLatToXYZ::LonLatToXYZ(double a, double b) { a_(a), b_(b) {} Point3 operator()(const PointLonLat& p) const override { return S::convertSphericalToCartesian(a_, b_, p, 0.); } PointLonLat operator()(const Point3& q) const override { NOTIMP; } + Spec spec() const override { return Spec{{{"a", a_}, {"b", b_}}}; } }; impl_.reset(types::is_approximately_equal(a, b) ? static_cast(new LonLatToSphereXYZ(a)) @@ -62,4 +63,9 @@ LonLatToXYZ::LonLatToXYZ(const Configuration& config) : LonLatToXYZ(config.getDouble("a", config.getDouble("R", 1.)), config.getDouble("b", config.getDouble("R", 1.))) {} +Projection::Spec LonLatToXYZ::spec() const { + return impl_->spec(); +} + + } // namespace eckit::geo::projection diff --git a/src/eckit/geo/projection/LonLatToXYZ.h b/src/eckit/geo/projection/LonLatToXYZ.h index c47061fd0..9b2ea6998 100644 --- a/src/eckit/geo/projection/LonLatToXYZ.h +++ b/src/eckit/geo/projection/LonLatToXYZ.h @@ -47,7 +47,8 @@ class LonLatToXYZ final : public Projection { PointLonLat inv(const Point3& q) const { return (*impl_)(q); } // -- Overridden methods - // None + + Spec spec() const override; // -- Class members // None @@ -69,6 +70,7 @@ class LonLatToXYZ final : public Projection { virtual Point3 operator()(const PointLonLat&) const = 0; virtual PointLonLat operator()(const Point3&) const = 0; + virtual Spec spec() const = 0; }; // -- Members diff --git a/src/eckit/geo/projection/None.cc b/src/eckit/geo/projection/None.cc index c3df3b4ad..88f168cba 100644 --- a/src/eckit/geo/projection/None.cc +++ b/src/eckit/geo/projection/None.cc @@ -12,6 +12,8 @@ #include "eckit/geo/projection/None.h" +#include "eckit/config/MappedConfiguration.h" + namespace eckit::geo::projection { @@ -25,4 +27,9 @@ static ProjectionBuilder __projection4("plate-carree"); None::None(const Configuration&) {} +Projection::Spec None::spec() const { + return Spec{}; +} + + } // namespace eckit::geo::projection diff --git a/src/eckit/geo/projection/None.h b/src/eckit/geo/projection/None.h index 1a9fc3ff5..3aa618f72 100644 --- a/src/eckit/geo/projection/None.h +++ b/src/eckit/geo/projection/None.h @@ -44,7 +44,8 @@ class None final : public Projection { // None // -- Overridden methods - // None + + Spec spec() const override; // -- Class members // None diff --git a/src/eckit/geo/projection/PROJ.cc b/src/eckit/geo/projection/PROJ.cc index f84130307..9b8e531b7 100644 --- a/src/eckit/geo/projection/PROJ.cc +++ b/src/eckit/geo/projection/PROJ.cc @@ -12,7 +12,7 @@ #include "eckit/geo/projection/PROJ.h" -#include "eckit/config/Configuration.h" +#include "eckit/config/MappedConfiguration.h" #include "eckit/exception/Exceptions.h" @@ -123,4 +123,9 @@ Point PROJ::inv(const Point& q) const { } +Projection::Spec PROJ::spec() const { + return Spec{{{"source", source_}, {"target", target_}}}; +} + + } // namespace eckit::geo::projection diff --git a/src/eckit/geo/projection/PROJ.h b/src/eckit/geo/projection/PROJ.h index fe87a74c2..cc38f535d 100644 --- a/src/eckit/geo/projection/PROJ.h +++ b/src/eckit/geo/projection/PROJ.h @@ -83,6 +83,8 @@ class PROJ final : public Projection { Point fwd(const Point&) const override; Point inv(const Point&) const override; + Spec spec() const override; + // -- Class members // None diff --git a/src/eckit/geo/projection/Rotation.cc b/src/eckit/geo/projection/Rotation.cc index 66b4bfc3b..77bb02cde 100644 --- a/src/eckit/geo/projection/Rotation.cc +++ b/src/eckit/geo/projection/Rotation.cc @@ -15,7 +15,7 @@ #include #include -#include "eckit/config/Configuration.h" +#include "eckit/config/MappedConfiguration.h" #include "eckit/geo/UnitSphere.h" #include "eckit/geo/util.h" #include "eckit/maths/Matrix3.h" @@ -38,22 +38,32 @@ Rotation::Rotation(double south_pole_lon, double south_pole_lat, double angle) : struct NonRotated final : Rotate { PointLonLat operator()(const PointLonLat& p) const override { return p; } + Spec spec() const override { return Spec{}; } }; struct RotationAngle final : Rotate { explicit RotationAngle(double angle) : angle_(angle) {} PointLonLat operator()(const PointLonLat& p) const override { return {p.lon + angle_, p.lat}; } + Spec spec() const override { return Spec{{{"angle", angle_}}}; } const double angle_; }; struct RotationMatrix final : Rotate { explicit RotationMatrix(M&& R) : - R_(R) {} + RotationMatrix(std::move(R), 0, 0, 0) {} + RotationMatrix(M&& R, double south_pole_lon, double south_pole_lat, double angle) : + R_(R), south_pole_lon_(south_pole_lon), south_pole_lat_(south_pole_lat), angle_(angle) {} PointLonLat operator()(const PointLonLat& p) const override { return UnitSphere::convertCartesianToSpherical(R_ * UnitSphere::convertSphericalToCartesian(p)); } + Spec spec() const override { + return Spec{{{"south_pole_lon", south_pole_lon_}, {"south_pole_lat", south_pole_lat_}, {"angle", angle_}}}; + } const M R_; + const double south_pole_lon_; + const double south_pole_lat_; + const double angle_; }; const auto alpha = util::degree_to_radian * angle; @@ -114,4 +124,9 @@ Rotation::Rotation(const Configuration& config) : Rotation(config.getDouble("south_pole_lon"), config.getDouble("south_pole_lat"), config.getDouble("angle", 0)) {} +Projection::Spec Rotation::spec() const { + return fwd_->spec(); +} + + } // namespace eckit::geo::projection diff --git a/src/eckit/geo/projection/Rotation.h b/src/eckit/geo/projection/Rotation.h index 2ad36a4da..947ef5c81 100644 --- a/src/eckit/geo/projection/Rotation.h +++ b/src/eckit/geo/projection/Rotation.h @@ -48,7 +48,8 @@ class Rotation final : public Projection { PointLonLat inv(const PointLonLat& q) const { return (*inv_)(q); } // -- Overridden methods - // None + + Spec spec() const override; // -- Class members // None @@ -69,6 +70,7 @@ class Rotation final : public Projection { void operator=(Rotate&&) = delete; virtual PointLonLat operator()(const PointLonLat&) const = 0; + virtual Spec spec() const = 0; }; // -- Members diff --git a/tests/geo/test_projection.cc b/tests/geo/test_projection.cc index e4d69f5a9..e39cb7327 100644 --- a/tests/geo/test_projection.cc +++ b/tests/geo/test_projection.cc @@ -15,12 +15,18 @@ #include "eckit/config/MappedConfiguration.h" #include "eckit/geo/Projection.h" +#include "eckit/geo/figure/Sphere.h" #include "eckit/geo/projection/LonLatToXYZ.h" +#include "eckit/geo/projection/Mercator.h" #include "eckit/geo/projection/Rotation.h" +#include "eckit/log/Log.h" #include "eckit/testing/Test.h" -int main(int argc, char* argv[]) { +namespace eckit::test { + + +CASE("projection") { using eckit::MappedConfiguration; using namespace eckit::geo; @@ -30,13 +36,15 @@ int main(int argc, char* argv[]) { Point p = PointLonLat{1, 1}; - { + + SECTION("projection type: none") { Projection projection(ProjectionFactory::instance().get("none").create(MappedConfiguration{})); EXPECT(points_equal(p, projection->inv(p))); EXPECT(points_equal(p, projection->fwd(p))); } - { + + SECTION("projection type: rotation") { MappedConfiguration param({ {"projection", "rotation"}, {"south_pole_lat", -91.}, @@ -49,7 +57,8 @@ int main(int argc, char* argv[]) { EXPECT(points_equal(p, projection->fwd(projection->inv(p)))); } - { + + SECTION("projection type: ll_to_xyz") { Projection s1(ProjectionFactory::instance().get("ll_to_xyz").create(MappedConfiguration({{"R", 1.}}))); Projection s2( ProjectionFactory::instance().get("ll_to_xyz").create(MappedConfiguration({{"a", 1.}, {"b", 1.}}))); @@ -85,35 +94,31 @@ int main(int argc, char* argv[]) { } - { - + SECTION("projection type: ll_to_xyz") { const PointLonLat p(723., 1.); // <- FIXME - { - projection::LonLatToXYZ to_xyz(1.); + projection::LonLatToXYZ to_xyz_r(1.); - auto q = to_xyz.fwd(p); - auto r = to_xyz.inv(q); - std::cout << "p(lat, lon): " << p << " -> p(x,y,z): " << q << " -> p(lat, lon): " << r << std::endl; + auto q = to_xyz_r.fwd(p); + auto r = to_xyz_r.inv(q); + std::cout << "p(lat, lon): " << p << " -> p(x,y,z): " << q << " -> p(lat, lon): " << r << std::endl; - EXPECT(points_equal(p, r)); - } + EXPECT(points_equal(p, r)); + // oblate spheroid (supported) + projection::LonLatToXYZ to_xyz_ab(3., 2.); - { - projection::LonLatToXYZ to_xyz_ab(3., 2.); // oblate - projection::LonLatToXYZ to_xyz_ba(2., 3.); // problate - - for (const auto& lon : {0., 90., 180., 270.}) { - PointLonLat p{lon, 0.}; - std::cout << "p(lat, lon): " << p << " -> p_ab(x,y,z): " << to_xyz_ab.fwd(p) - << ", p_ba(x,y,z): " << to_xyz_ba.fwd(p) << std::endl; - } + for (const auto& lon : {0., 90., 180., 270.}) { + PointLonLat p{lon, 0.}; + std::cout << "p(lat, lon): " << p << " -> p_ab(x,y,z): " << to_xyz_ab.fwd(p) << std::endl; } + + // problate spheroid (not supported) + EXPECT_THROWS(projection::LonLatToXYZ(2., 3.)); } - { + SECTION("projection type: rotation") { const PointLonLat p(1, 1); int delta[] = {-360, -180, -1, 0, 1, 90, 91, 180}; @@ -133,7 +138,8 @@ int main(int argc, char* argv[]) { } - { +#if 0 + SECTION("projection type: rotation") { const int Ni = 12; const int Nj = 3; @@ -230,14 +236,17 @@ int main(int argc, char* argv[]) { PointLonLat a(static_cast(i) * 360. / static_cast(Ni), static_cast(j - Nj) * 90. / static_cast(Nj)); auto b = rot.fwd(a); + auto c = rot.inv(b); + EXPECT(points_equal(b, ref[k])); - EXPECT(points_equal(a, rot.inv(b))); + EXPECT(points_equal(a, c)); } } } +#endif - { + SECTION("projection type: rotation") { const projection::Rotation non_rotated(0., -90., 0.); const projection::Rotation rotation_angle(0., -90., -180.); const projection::Rotation rotation_matrix(4., -40., 180.); @@ -277,52 +286,87 @@ int main(int argc, char* argv[]) { } - if (ProjectionFactory::instance().exists("proj")) { - std::cout.precision(14); + SECTION("projection type: proj") { + if (ProjectionFactory::instance().exists("proj")) { + std::cout.precision(14); - PointLonLat a{12., 55.}; + PointLonLat a{12., 55.}; - struct { - const Point b; - const std::string target; - } tests[] = { - {Point2{691875.632137542, 6098907.825129169}, "+proj=utm +zone=32 +datum=WGS84"}, - {Point2{691875.632137542, 6098907.825129169}, "EPSG:32632"}, - {a, "EPSG:4326"}, - {a, "EPSG:4979"}, - {Point3{3586469.6567764, 762327.65877826, 5201383.5232023}, "EPSG:4978"}, - {Point3{3574529.7050235, 759789.74368715, 5219005.2599833}, "+proj=cart +R=6371229."}, - {Point3{3574399.5431832, 759762.07693392, 5218815.216709}, "+proj=cart +ellps=sphere"}, - {a, "+proj=latlon +ellps=sphere"}, - }; + struct { + const Point b; + const std::string target; + } tests[] = { + {Point2{691875.632137542, 6098907.825129169}, "+proj=utm +zone=32 +datum=WGS84"}, + {Point2{691875.632137542, 6098907.825129169}, "EPSG:32632"}, + {a, "EPSG:4326"}, + {a, "EPSG:4979"}, + {Point3{3586469.6567764, 762327.65877826, 5201383.5232023}, "EPSG:4978"}, + {Point3{3574529.7050235, 759789.74368715, 5219005.2599833}, "+proj=cart +R=6371229."}, + {Point3{3574399.5431832, 759762.07693392, 5218815.216709}, "+proj=cart +ellps=sphere"}, + {a, "+proj=latlon +ellps=sphere"}, + }; - for (const auto& test : tests) { - Projection projection(ProjectionFactory::instance().get("proj").create( - MappedConfiguration{{{"source", "EPSG:4326"}, {"target", test.target}}})); + for (const auto& test : tests) { + Projection projection(ProjectionFactory::instance().get("proj").create( + MappedConfiguration{{{"source", "EPSG:4326"}, {"target", test.target}}})); #if 0 std::cout << "ellipsoid: '" << PROJ::ellipsoid(projection.target()) << std::endl; #endif - auto b = projection->fwd(a); - auto c = projection->inv(b); + auto b = projection->fwd(a); + auto c = projection->inv(b); - std::cout << "-> a:" << a << " -> fwd(a):" << b << " -> inv(fwd(a)):" << c << std::endl; + std::cout << "-> a:" << a << " -> fwd(a):" << b << " -> inv(fwd(a)):" << c << std::endl; - EXPECT(points_equal(b, test.b)); - EXPECT(points_equal(c, a)); + EXPECT(points_equal(b, test.b)); + EXPECT(points_equal(c, a)); - Projection reverse(ProjectionFactory::instance().get("proj").create( - MappedConfiguration({{"source", test.target}, {"target", "EPSG:4326"}}))); + Projection reverse(ProjectionFactory::instance().get("proj").create( + MappedConfiguration({{"source", test.target}, {"target", "EPSG:4326"}}))); - auto d = reverse->fwd(test.b); - auto e = reverse->inv(d); + auto d = reverse->fwd(test.b); + auto e = reverse->inv(d); - std::cout << "-> b:" << test.b << " -> fwd(b):" << d << " -> inv(fwd(b)):" << e << std::endl; + std::cout << "-> b:" << test.b << " -> fwd(b):" << d << " -> inv(fwd(b)):" << e << std::endl; - EXPECT(points_equal(d, a)); - EXPECT(points_equal(e, test.b)); + EXPECT(points_equal(d, a)); + EXPECT(points_equal(e, test.b)); + } } } + + + SECTION("projection type: mercator") { + PointLonLat first{262.036, 14.7365}; + + { + projection::Mercator projection(0., 14., new figure::Sphere(6371229.), first); + + auto p = projection.inv({0, 0}); + Log::info() << p << std::endl; + + auto q = projection.fwd(p); + Log::info() << q << std::endl; + } + + { + projection::Mercator projection(-180., 0., new figure::Sphere(1.), {0, 0}); + + auto q = projection.fwd({-75, 35}); + Log::info() << q << std::endl; + + auto p = projection.inv(q); + Log::info() << p << std::endl; + } + } +} + + +} // namespace eckit::test + + +int main(int argc, char** argv) { + return eckit::testing::run_tests(argc, argv); }