Skip to content

Commit

Permalink
GRIDEDIT-757 Reduced code duplication
Browse files Browse the repository at this point in the history
  • Loading branch information
@BillSenior
BillSenior committed Oct 27, 2023
1 parent 1984a58 commit 9da956b
Show file tree
Hide file tree
Showing 2 changed files with 153 additions and 12 deletions.
58 changes: 53 additions & 5 deletions libs/MeshKernel/include/MeshKernel/ProjectionConversions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,13 +29,20 @@

#include <concepts>

#include <boost/geometry.hpp>
#include <boost/geometry/core/coordinate_system.hpp>
#include <boost/geometry/geometries/point_xy.hpp>

#include "MeshKernel/Operations.hpp"
#include "MeshKernel/Point.hpp"
#include "MeshKernel/Vector.hpp"

namespace meshkernel
{

/// @brief Namespace alias for boost::geometry
namespace bg = boost::geometry;

/// @brief Ensure any instantiation of the MeshTransformation Compute function is with the correct operation
template <typename Operation>
concept HasTransformationOperation = requires(Operation op, Point p) {{ op(p)} -> std::same_as<Point>; };
Expand All @@ -50,11 +57,17 @@ namespace meshkernel
concept ConversionFunctor = HasTransformationOperation<Operation> && HasConversionProjection<Operation>;

/// @brief Converts points from spherical to Cartesian coordinate system.
class ConvertSphericalToCartesian
template <typename ProjectionConversion>
class ConvertSphericalToCartesianBase
{
public:
/// @brief Construct ConvertSphericalToCartesian with the origin of the target mesh
ConvertSphericalToCartesian(const Point& o) : origin(o) {}
/// @brief
using LongLat = bg::model::d2::point_xy<double, bg::cs::geographic<bg::degree>>;

/// @brief
using UTM = bg::model::d2::point_xy<double, Projection>;

ConvertSphericalToCartesianBase(const ProjectionConversion& proj) : projection(proj) {}

/// @brief The coordinate system of the point parameter to the conversion operation.
Projection SourceProjection() const
Expand All @@ -71,13 +84,48 @@ namespace meshkernel
/// @brief Apply the conversion of a point in Spherical coordinate system to Cartesian
Point operator()(const Point& pnt) const
{
Point result = origin + Vector(GetDx(origin, pnt, Projection::spherical), GetDy(origin, pnt, Projection::spherical));
LongLat longLat(pnt.x, pnt.y);
UTM utm{0.0, 0.0};
projection.forward(longLat, utm);

Point result(utm.x(), utm.y());
return result;
}

/// @brief Apply the reverse conversion of a point in Cartesian coordinate system to Spherical
Point reverse(const Point& pnt) const
{
UTM utm{pnt.x, pnt.y};
LongLat longLat{0.0, 0.0};
projection.inverse(utm, longLat);

Point result(longLat.x(), longLat.y());
return result;
}

private:
/// @brief The origin of the target mesh.
Point origin;
ProjectionConversion projection;
};

/// @brief Converts points from spherical to Cartesian coordinate system.
template <const int EpsgCode>
class ConvertSphericalToCartesianEPSG : public ConvertSphericalToCartesianBase<bg::srs::projection<bg::srs::static_epsg<EpsgCode>>>
{
public:
/// @brief The EPSG projection
using EpsgProjection = bg::srs::projection<bg::srs::static_epsg<EpsgCode>>;

/// @brief Construct spherical to Cartesian with an EPSG code
ConvertSphericalToCartesianEPSG() : ConvertSphericalToCartesianBase<bg::srs::projection<bg::srs::static_epsg<EpsgCode>>>(EpsgProjection()) {}
};

/// @brief Converts points from spherical to Cartesian coordinate system.
class ConvertSphericalToCartesian : public ConvertSphericalToCartesianBase<bg::srs::projection<>>
{
public:
/// @brief Construct spherical to Cartesian with an zone string
ConvertSphericalToCartesian(const std::string& zone) : ConvertSphericalToCartesianBase<bg::srs::projection<>>(bg::srs::proj4(zone)) {}
};

/// @brief Apply a conversion to nodes of a mesh
Expand Down
107 changes: 100 additions & 7 deletions libs/MeshKernel/tests/src/MeshConversionTests.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@

#include <boost/geometry.hpp>
#include <boost/geometry/core/coordinate_system.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/srs/epsg.hpp>

#include <chrono>
#include <gtest/gtest.h>
#include <random>
Expand All @@ -13,19 +19,106 @@

namespace mk = meshkernel;

TEST(MeshConversionTests, SphericalToCartesianClassTest)
TEST(MeshConversionTests, SphericalToCartesianClassTestFromCode)
{
constexpr double tolerance = 1.0e-3;
constexpr double tolerance = 1.0e-6;

constexpr int EpsgCode = 3043;

mk::ConvertSphericalToCartesianEPSG<EpsgCode> conversion;

// Notice that the coordinate are in (longitude, latitude)
mk::Point pnt{4.82898, 52.410695};
mk::Point result = conversion(pnt);
mk::Point expected{624402.8057, 5808292.0905};

EXPECT_EQ(mk::Projection::spherical, conversion.SourceProjection());
EXPECT_EQ(mk::Projection::cartesian, conversion.TargetProjection());
EXPECT_TRUE(mk::IsEqual(expected.x, result.x, tolerance));
EXPECT_TRUE(mk::IsEqual(expected.y, result.y, tolerance));

pnt = mk::Point(2.1734, 41.3851);
result = conversion(pnt);

EXPECT_TRUE(mk::IsEqual(430887.56433058, result.x, tolerance));
EXPECT_TRUE(mk::IsEqual(4581837.85323929693549871, result.y, tolerance));

// Deltares Delft office
pnt = mk::Point(4.382253, 51.986373);
expected = mk::Point(594919.016912, 5760424.72635);
result = conversion(pnt);
EXPECT_TRUE(mk::IsEqual(expected.x, result.x, tolerance));
EXPECT_TRUE(mk::IsEqual(expected.y, result.y, tolerance));

// Brussels
pnt = mk::Point(4.351697, 50.846557);
expected = mk::Point(595158.7545, 5633632.17474);
result = conversion(pnt);
EXPECT_TRUE(mk::IsEqual(expected.x, result.x, tolerance));
EXPECT_TRUE(mk::IsEqual(expected.y, result.y, tolerance));

// Paris
pnt = mk::Point(2.352222, 48.856614);
expected = mk::Point(452484.160, 5411718.719);
result = conversion(pnt);
EXPECT_TRUE(mk::IsEqual(expected.x, result.x, tolerance));
EXPECT_TRUE(mk::IsEqual(expected.y, result.y, tolerance));
}

mk::Point origin{0.0, 0.0};
mk::ConvertSphericalToCartesian conversion(origin);
TEST(MeshConversionTests, SphericalToCartesianClassTestFromStringZone31)
{
constexpr double tolerance = 1.0e-6;

mk::Point pnt{0.001, 0.0001}; // ≃ 111.3 and 11.13 metres
mk::ConvertSphericalToCartesian conversion("+proj=utm +lat_1=0.5 +lat_2=2 +n=0.5 +zone=31");

mk::Point pnt{4.897, 52.371};
mk::Point result = conversion(pnt);
mk::Point expected{624402.8057, 5808292.0905};

EXPECT_EQ(mk::Projection::spherical, conversion.SourceProjection());
EXPECT_EQ(mk::Projection::cartesian, conversion.TargetProjection());
EXPECT_TRUE(mk::IsEqual(111.3, result.x, tolerance)) << "Expected: 111.3, found: " << result.x;
EXPECT_TRUE(mk::IsEqual(11.13, result.y, tolerance)) << "Expected: 11.13, found: " << result.y;
EXPECT_TRUE(mk::IsEqual(629144.771310, result.x, tolerance));
EXPECT_TRUE(mk::IsEqual(5803996.656944, result.y, tolerance));

pnt = mk::Point(2.1734, 41.3851);
result = conversion(pnt);

EXPECT_TRUE(mk::IsEqual(430887.56433058, result.x, tolerance));
EXPECT_TRUE(mk::IsEqual(4581837.85323929693549871, result.y, tolerance));

// Deltares
pnt = mk::Point(4.382253, 51.986373);
expected = mk::Point(594919.016912, 5760424.72635);
result = conversion(pnt);
EXPECT_TRUE(mk::IsEqual(expected.x, result.x, tolerance));
EXPECT_TRUE(mk::IsEqual(expected.y, result.y, tolerance));

// Brussels
pnt = mk::Point(4.351697, 50.846557);
expected = mk::Point(595158.7545, 5633632.17474);
result = conversion(pnt);
EXPECT_TRUE(mk::IsEqual(expected.x, result.x, tolerance));
EXPECT_TRUE(mk::IsEqual(expected.y, result.y, tolerance));

// Paris
pnt = mk::Point(2.352222, 48.856614);
expected = mk::Point(452484.160, 5411718.719);
result = conversion(pnt);
EXPECT_TRUE(mk::IsEqual(expected.x, result.x, tolerance));
EXPECT_TRUE(mk::IsEqual(expected.y, result.y, tolerance));
}

TEST(MeshConversionTests, SphericalToCartesianClassTestFromStringZone30)
{
constexpr double tolerance = 1.0e-6;

mk::ConvertSphericalToCartesian conversion("+proj=utm +lat_1=0.5 +lat_2=2 +n=0.5 +zone=30N");

// Cwmbran
mk::Point pnt{-3.020822, 51.653644};
mk::Point result = conversion(pnt);
mk::Point expected{498559.553, 5722516.858};

EXPECT_TRUE(mk::IsEqual(expected.x, result.x, tolerance));
EXPECT_TRUE(mk::IsEqual(expected.y, result.y, tolerance));
}

0 comments on commit 9da956b

Please sign in to comment.