diff --git a/src/atlas/grid/CubedSphereGrid2.cc b/src/atlas/grid/CubedSphereGrid2.cc index 339e422df..eec7fbce2 100644 --- a/src/atlas/grid/CubedSphereGrid2.cc +++ b/src/atlas/grid/CubedSphereGrid2.cc @@ -13,4 +13,7 @@ CubedSphereGrid2::CubedSphereGrid2(idx_t resolution): CubedSphereGrid2::CubedSphereGrid2(const Grid& grid): Grid(grid), grid_(cubedsphere_grid2(get())) {} +CubedSphereGrid2::CubedSphereGrid2(idx_t resolution, Projection projection): + Grid(new grid::detail::grid::CubedSphere2(resolution, projection)), grid_(cubedsphere_grid2(get())) {} + } // namespace atlas diff --git a/src/atlas/grid/CubedSphereGrid2.h b/src/atlas/grid/CubedSphereGrid2.h index 5cf036c71..a2abab114 100644 --- a/src/atlas/grid/CubedSphereGrid2.h +++ b/src/atlas/grid/CubedSphereGrid2.h @@ -12,6 +12,7 @@ class CubedSphereGrid2 : public atlas::Grid { public: CubedSphereGrid2(idx_t resolution); CubedSphereGrid2(const Grid& grid); + CubedSphereGrid2(idx_t resolution, Projection projection); bool valid() const { return grid_; }; operator bool() const { return valid(); } diff --git a/src/atlas/grid/detail/grid/CubedSphere2.cc b/src/atlas/grid/detail/grid/CubedSphere2.cc index b51de450a..aab885aea 100644 --- a/src/atlas/grid/detail/grid/CubedSphere2.cc +++ b/src/atlas/grid/detail/grid/CubedSphere2.cc @@ -1,25 +1,49 @@ #include "atlas/grid/detail/grid/CubedSphere2.h" #include +#include +#include "atlas/grid/CubedSphereGrid2.h" +#include "atlas/grid/detail/grid/GridBuilder.h" +#include "atlas/grid/detail/grid/GridFactory.h" +#include "atlas/runtime/Exception.h" #include "eckit/geometry/Sphere.h" #include "eckit/utils/Hash.h" +#include "eckit/utils/Translator.h" namespace atlas { namespace grid { namespace detail { namespace grid { +static eckit::Translator to_int; + // Public methods CubedSphere2::CubedSphere2(idx_t resolution) : N_(resolution) {} +CubedSphere2::CubedSphere2(idx_t resolution, Projection projection) : + Grid(), N_(resolution) { + + // Copy members + util::Config defaultProjConfig; + defaultProjConfig.set("type", "lonlat"); + projection_ = projection ? projection : Projection(defaultProjConfig); + + // Domain + domain_ = GlobalDomain(); +} + std::string CubedSphere2::name() const { return "CS-LFR-" + std::to_string(N_) + "-2"; } std::string CubedSphere2::type() const { - return type_; + return static_type(); +} + +std::string CubedSphere2::static_type() { + return "cubedsphere2"; } // Provide a unique identification hash for the grid and the projection. @@ -125,6 +149,57 @@ PointXYZ CubedSphere2::tangent_to_xyz_coord(const PointXY& tan_coord, idx_t tile return PointXYZ::normalize(xyz); } +namespace { +GridFactoryBuilder __register_CubedSphere2(CubedSphere2::static_type()); +} + +static class cubedsphere2_lfric : public GridBuilder { +public: + cubedsphere2_lfric(): + GridBuilder("cubedsphere2_lfric", {"^[Cc][Ss][_-][Ll][Ff][Rr][-_]([1-9][0-9]*)[_-][2]$"}, + {"CS-LFR--2"}) {} + + void print(std::ostream& os) const override { + os << std::left << std::setw(20) << "CS-LFR--2" + << "Cubed sphere for LFRic"; + } + + // Factory constructor + const atlas::Grid::Implementation* create(const std::string& name, const Grid::Config& config) const override { + int id; + std::vector matches; + if (match(name, matches, id)) { + util::Config gridconf(config); + int N = to_int(matches[0]); + gridconf.set("type", type()); + gridconf.set("N", N); + return create(gridconf); + } + return nullptr; + } + + // Factory constructor + const atlas::Grid::Implementation* create(const Grid::Config& config) const override { + int N = 0; + if (not config.get("N", N)) { + throw_AssertionFailed("Could not find \"N\" in configuration of cubed sphere grid 2", Here()); + } + + std::string name = "CS-LFR-" + std::to_string(N) + "-2"; + util::Config projconf; + projconf.set("type", "lonlat"); + + return new CubedSphereGrid2::grid_t(N, Projection(projconf)); + } + + void force_link() {} + +} cubedsphere2_lfric_; + +void force_link_CubedSphere2() { + cubedsphere2_lfric_.force_link(); +} + } // namespace grid } // namespace detail } // namespace grid diff --git a/src/atlas/grid/detail/grid/CubedSphere2.h b/src/atlas/grid/detail/grid/CubedSphere2.h index 988bfc4a2..fd752e09b 100644 --- a/src/atlas/grid/detail/grid/CubedSphere2.h +++ b/src/atlas/grid/detail/grid/CubedSphere2.h @@ -109,9 +109,11 @@ class CubedSphere2 : public Grid { using Spec = atlas::util::Config; CubedSphere2(idx_t resolution); + CubedSphere2(idx_t resolution, Projection projection); std::string name() const override; std::string type() const override; + static std::string static_type(); idx_t N() const {return N_;} void hash(eckit::Hash&) const override; @@ -157,8 +159,6 @@ class CubedSphere2 : public Grid { static constexpr idx_t nTiles_ = 6; private: - std::string type_ = {"cubedsphere2"}; - using Matrix = std::array, 3>; /* @@ -181,7 +181,6 @@ class CubedSphere2 : public Grid { {{ {-1, 0, 0}, {0, 1, 0}, {0, 0, 1} }}, {{ {-1, 0, 0}, {0, -1, 0}, {0, 0, -1} }} }}; - }; } // namespace grid diff --git a/src/atlas/grid/detail/grid/GridBuilder.cc b/src/atlas/grid/detail/grid/GridBuilder.cc index adfdda5f3..cfd63dfc4 100644 --- a/src/atlas/grid/detail/grid/GridBuilder.cc +++ b/src/atlas/grid/detail/grid/GridBuilder.cc @@ -121,6 +121,7 @@ static void init() { namespace detail { namespace grid { void force_link_CubedSphere(); +void force_link_CubedSphere2(); void force_link_Gaussian(); void force_link_LonLat(); void force_link_Regional(); @@ -130,6 +131,7 @@ void force_link_Regional_var_resolution(); const GridBuilder::Registry& GridBuilder::nameRegistry() { detail::grid::force_link_CubedSphere(); + detail::grid::force_link_CubedSphere2(); detail::grid::force_link_Gaussian(); detail::grid::force_link_LonLat(); detail::grid::force_link_Regional(); diff --git a/src/tests/grid/test_cubedsphere_2.cc b/src/tests/grid/test_cubedsphere_2.cc index 873371ec6..0939aa310 100644 --- a/src/tests/grid/test_cubedsphere_2.cc +++ b/src/tests/grid/test_cubedsphere_2.cc @@ -8,16 +8,24 @@ namespace test { namespace { using Point2 = eckit::geometry::Point2; +// XY/lonlat for a Cubed Sphere, with N = 2 +std::vector kgo_lonlat { + {-22.5,20.941}, {22.5,20.941}, {-22.5,-20.941}, {22.5,-20.941}, {67.5,20.941}, {112.5,20.941}, + {67.5,-20.941}, {112.5,-20.941}, {157.5,20.941}, {-157.5,20.941}, {157.5,-20.941}, {-157.5,-20.941}, + {-112.5,20.941}, {-67.5,20.941}, {-112.5,-20.941}, {-67.5,-20.941}, {-45,59.6388}, {-135,59.6388}, + {45,59.6388}, {135,59.6388}, {45,-59.6388}, {135,-59.6388}, {-45,-59.6388}, {-135,-59.6388} +}; + template -bool compare_2D_points(std::vector a, std::vector b, double tolerance = 1e-4) { +bool compare_2D_points(std::vector a, std::vector b, bool print_diff = true, double tolerance = 1e-4) { // Uses a tolerance as the values are stored more precisely than they are printed ATLAS_ASSERT(a.size() == b.size()); bool equal = true; for (int i = 0; i < b.size(); ++i) { for (int j = 0; j < 2; ++j) { if (std::abs(b[i][j] - a[i][j]) > tolerance) { - std::cout << "[" << i << ", " << j << "]\n\t" << a[i][j] << " != " << b[i][j] - << "\n\tdiff = " << b[i][j] - a[i][j] << std::endl; + if (print_diff) std::cout << "[" << i << ", " << j << "]\n\t" << a[i][j] << " != " << b[i][j] + << "\n\tdiff = " << b[i][j] - a[i][j] << std::endl; equal = false; } } @@ -27,17 +35,17 @@ bool compare_2D_points(std::vector a, std::vector b, dou CASE("cubed_sphere_instantiation") { const int n = 2; - const Grid grid = CubedSphereGrid2(n); - - EXPECT(grid.name() == "CS-LFR-" + std::to_string(n) + "-2"); + const std::string name = "CS-LFR-" + std::to_string(n) + "-2"; + + const Grid grid = Grid(name); + EXPECT(grid.name() == name); EXPECT(grid.type() == "cubedsphere2"); EXPECT(grid.size() == n * n * 6); } CASE("constructor_with_grid") { auto grid_og = Grid("O32"); - // auto grid_cs = Grid("CS-LFR-4-2"); // The grid factory is implemented in the next PR - auto grid_cs = CubedSphereGrid2(4); + auto grid_cs = Grid("CS-LFR-4-2"); EXPECT( CubedSphereGrid2( grid_og ).valid() == false ); EXPECT( bool(CubedSphereGrid2( grid_og )) == false ); EXPECT( CubedSphereGrid2( grid_cs ).valid() == true ); @@ -45,14 +53,6 @@ CASE("constructor_with_grid") { } CASE("cubed_sphere_grid_kgo") { - // Lonlat and XY are both currently lonlat positions - std::vector kgo_lonlat { // N = 2 - {-22.5,20.941}, {22.5,20.941}, {-22.5,-20.941}, {22.5,-20.941}, {67.5,20.941}, {112.5,20.941}, - {67.5,-20.941}, {112.5,-20.941}, {157.5,20.941}, {-157.5,20.941}, {157.5,-20.941}, {-157.5,-20.941}, - {-112.5,20.941}, {-67.5,20.941}, {-112.5,-20.941}, {-67.5,-20.941}, {-45,59.6388}, {-135,59.6388}, - {45,59.6388}, {135,59.6388}, {45,-59.6388}, {135,-59.6388}, {-45,-59.6388}, {-135,-59.6388} - }; - const Grid grid = CubedSphereGrid2(2); // LonLat @@ -74,6 +74,24 @@ CASE("cubed_sphere_grid_kgo") { EXPECT(compare_2D_points(points_xy, kgo_lonlat)); } +CASE("cubed_sphere_rotated_lonlat") { + const auto grid_rotated = CubedSphereGrid2(2, Projection(util::Config("type", "rotated_lonlat")("north_pole", std::vector{4., 54.}))); + + // Expect XY points to still match + std::vector points_xy; + for (const auto &xy : grid_rotated.xy()) { + points_xy.push_back(xy); + } + EXPECT(compare_2D_points(points_xy, kgo_lonlat) == true); + + // Expect lonlats to be different + std::vector points_lonlat; + for (const auto &lonlat : grid_rotated.lonlat()) { + points_lonlat.push_back(lonlat); + } + EXPECT(compare_2D_points(points_lonlat, kgo_lonlat, false) == false); +} + } // namespace } // namespace test } // namespace atlas