New cubed sphere 2 grid

src/atlas/CMakeLists.txt

@@ -131,6 +131,8 @@ domain/detail/DomainFactory.cc
 grid.h
 grid/CubedSphereGrid.cc
 grid/CubedSphereGrid.h
+grid/CubedSphereGrid2.cc
+grid/CubedSphereGrid2.h
 grid/Grid.cc
 grid/Grid.h
 grid/SpecRegistry.cc
@@ -155,6 +157,8 @@ grid/Tiles.cc
 
 grid/detail/grid/CubedSphere.cc
 grid/detail/grid/CubedSphere.h
+grid/detail/grid/CubedSphere2.cc
+grid/detail/grid/CubedSphere2.h
 grid/detail/grid/GridBuilder.h
 grid/detail/grid/GridBuilder.cc
 grid/detail/grid/GridFactory.h

src/atlas/grid/CubedSphereGrid2.cc

@@ -0,0 +1,9 @@
+#include "atlas/grid/CubedSphereGrid2.h"
+#include "atlas/grid/detail/grid/CubedSphere2.h"
+
+namespace atlas {
+
+CubedSphereGrid2::CubedSphereGrid2(idx_t resolution):
+    Grid(new grid::detail::grid::CubedSphere2(resolution)) {}
+
+}  // namespace atlas

src/atlas/grid/CubedSphereGrid2.h

@@ -0,0 +1,12 @@
+#pragma once
+
+#include "atlas/grid/Grid.h"
+
+namespace atlas {
+
+class CubedSphereGrid2 : public atlas::Grid {
+public:
+    CubedSphereGrid2(idx_t resolution);
+};
+
+}  // namespace atlas

src/atlas/grid/detail/grid/CubedSphere2.cc

@@ -0,0 +1,137 @@
+#include "atlas/grid/detail/grid/CubedSphere2.h"
+
+#include <cmath>
+
+#include "eckit/geometry/Sphere.h"
+#include "eckit/utils/Hash.h"
+
+namespace atlas {
+namespace grid {
+namespace detail {
+namespace grid {
+
+// Public methods
+
+CubedSphere2::CubedSphere2(idx_t resolution) : N_(resolution) {}
+
+std::string CubedSphere2::name() const {
+    return "CS-LFR-" + std::to_string(N_) + "-2";
+}
+
+std::string CubedSphere2::type() const {
+    return type_;
+}
+
+// Provide a unique identification hash for the grid and the projection.
+void CubedSphere2::hash(eckit::Hash& h) const {
+    h.add(name());
+    h.add(int(N_));
+
+    // also add projection information
+    projection().hash(h);
+
+    // also add domain information, even though already encoded in grid.
+    domain().hash(h);
+}
+
+// Return the bounding box for the grid, global
+RectangularLonLatDomain CubedSphere2::lonlatBoundingBox() const {
+    return GlobalDomain();
+}
+
+// Return the total number of points
+idx_t CubedSphere2::size() const {
+    return N_ * N_ * nTiles_;
+}
+
+// Return the specification for the grid.
+Grid::Spec CubedSphere2::spec() const {
+    Grid::Spec grid_spec;
+
+    grid_spec.set("name", name());
+    grid_spec.set("type", type());
+    grid_spec.set("projection", projection().spec());
+    grid_spec.set("domain", domain());
+
+    return grid_spec;
+}
+
+// Get the xy for a given index
+void CubedSphere2::xy(idx_t n, Point2& point) const {
+    auto [t, i, j] = get_cs_indices(n);
+
+    PointAlphaBeta ab = ij_to_curvilinear_coord(i, j);
+    PointXY tangent_xy = curvilinear_to_tangent_coord(ab);
+    PointXYZ xyz = tangent_to_xyz_coord(tangent_xy, t);
+    Matrix xyz_m { {xyz[0], xyz[1], xyz[2]} };
+
+    eckit::geometry::Sphere::convertCartesianToSpherical(xyz_m.norm(), xyz, point);
+}
+
+// Get the xy for a given index
+Point2 CubedSphere2::xy(idx_t n) const {
+    Point2 point;
+    xy(n, point);
+    return point;
+}
+
+// Get the lonlat for a given index
+void CubedSphere2::lonlat(idx_t n, Point2& point) const {
+    xy(n, point);
+    projection_.xy2lonlat(point);
+}
+
+// Get the lonlat for a given index
+Point2 CubedSphere2::lonlat(idx_t n) const {
+    Point2 point;
+    lonlat(n, point);
+    return point;
+}
+
+// Protected methods
+
+// Print the name of the Grid
+void CubedSphere2::print(std::ostream& os) const {
+    os << "CubedSphere2(Name:" << name() << ")";
+}
+
+// Private methods
+
+// Get t, i, and j for a given index
+CubedSphere2::CSIndices CubedSphere2::get_cs_indices(gidx_t n) const {
+    ATLAS_ASSERT(n <= size());
+    const idx_t tile_size = N() * N();
+    const idx_t t = n / tile_size;
+    const idx_t ij = n % tile_size;
+    const idx_t j = ij / N();
+    const idx_t i = ij % N();
+    return {t, i, j};
+}
+
+// Get the curvilinear coordinate for a given ij index
+CubedSphere2::PointAlphaBeta CubedSphere2::ij_to_curvilinear_coord(idx_t i, idx_t j) const {
+    const auto get_coord = [&](idx_t idx) {
+        return M_PI / 2 * (-0.5 + (0.5 + idx) / N());
+    };
+    return {get_coord(i), get_coord(j)};
+}
+
+// Get the point on the tangent plane for a given curvilinear coordinate
+PointXY CubedSphere2::curvilinear_to_tangent_coord(CubedSphere2::PointAlphaBeta& curvi_coord) const {
+    return {std::tan(curvi_coord[0]), std::tan(curvi_coord[1])};
+}
+
+// Transform a point on the tangent plane to a point on a cube
+PointXYZ CubedSphere2::tangent_to_xyz_coord(PointXY& tan_coord, idx_t tile) const {
+    Matrix tan_point { {tan_coord[0]}, {tan_coord[1]}, {1} };
+
+    Matrix xyz(3, 1);
+    xyz = lfric_rotations_[tile].transpose() * tan_point;
+
+    return {xyz(0), xyz(1), xyz(2)};
+}
+
+}  // namespace grid
+}  // namespace detail
+}  // namespace grid
+}  // namespace atlas

src/atlas/grid/detail/grid/CubedSphere2.h

@@ -0,0 +1,182 @@
+#pragma once
+
+#include "atlas/grid/detail/grid/Grid.h"
+#include "atlas/runtime/Exception.h"
+#include "atlas/util/Config.h"
+#include "atlas/util/Point.h"
+#include "eckit/maths/Matrix.h"
+
+namespace atlas {
+namespace grid {
+namespace detail {
+namespace grid {
+
+class CubedSphere2 : public Grid {
+private:
+
+    // Get the lonlat and return as PointLonLat object
+    struct ComputePointXY {
+        ComputePointXY(const CubedSphere2& grid): grid_(grid) {}
+        void operator()(idx_t n, PointXY& point) { grid_.lonlat(n, point); }
+        const CubedSphere2& grid_;
+    };
+
+    // Get the lonlat and return as PointLonLat object
+    struct ComputePointLonLat {
+        ComputePointLonLat(const CubedSphere2& grid): grid_(grid) {}
+        void operator()(idx_t n, PointLonLat& point) { grid_.lonlat(n, point); }
+        const CubedSphere2& grid_;
+    };
+
+    template <typename Base, typename ComputePoint>
+    class CubedSphere2Iterator : public Base {
+    public:
+        CubedSphere2Iterator(const CubedSphere2& grid, bool begin = true):
+            grid_(grid),
+            size_(grid.size()),
+            n_(begin ? 0 : size_),
+            compute_point{grid_}
+            {
+                if (n_ < size_) {
+                    compute_point(n_, point_);
+                }
+            }
+
+        // Return the point and move iterator to the next location
+        bool next(typename Base::value_type& point) {
+            if (n_ < size_) {
+                compute_point(n_, point);
+                ++n_;
+                return true;
+            }
+            return false;
+        }
+
+        // * operator    
+        const typename Base::reference operator*() const {return point_; }
+
+        // ++ operator, move to next element and return point
+        const Base& operator++() {
+            ++n_;
+            if (n_ < size_) {
+                compute_point(n_, point_);
+            }
+            return *this;
+        }
+
+        // +=  operator, move some distance through the iterator and return point
+        const Base& operator+=(typename Base::difference_type distance) {
+            n_ += distance;
+            if (n_ >= 0 && n_ < size_) {
+                compute_point(n_, point_);
+            }
+            return *this;
+        }
+
+        // == operator, check two iterators are on the same index
+        bool operator==(const Base& other) const {
+            return n_ == static_cast<const CubedSphere2Iterator&>(other).n_;
+        }
+
+        // != operator, check two iterators are not on the same index
+        bool operator!=(const Base& other) const {
+            return n_ != static_cast<const CubedSphere2Iterator&>(other).n_;
+        }
+
+        // Return the number of points between two iterators
+        typename Base::difference_type distance(const Base& other) const {
+            return static_cast<const CubedSphere2Iterator&>(other).n_ - n_;
+        }
+
+        // Clone the iterator in its current position
+        std::unique_ptr<Base> clone() const {
+            auto result = new CubedSphere2Iterator(grid_, false);
+            result->n_ = n_;
+            result->compute_point(n_, result->point_);
+            return std::unique_ptr<Base>(result);
+        }
+
+        const CubedSphere2& grid_;
+        idx_t size_;
+        idx_t n_;
+        typename Base::value_type point_;
+        ComputePoint compute_point;
+    };
+
+public:
+    using IteratorXY = CubedSphere2Iterator<Grid::IteratorXY, ComputePointXY>;
+    using IteratorLonLat = CubedSphere2Iterator<Grid::IteratorLonLat, ComputePointLonLat>;
+
+    using Spec = atlas::util::Config;
+
+    CubedSphere2(idx_t resolution);
+
+    std::string name() const override;
+    std::string type() const override;
+    idx_t N() const {return N_;}
+
+    void hash(eckit::Hash&) const override;
+    RectangularLonLatDomain lonlatBoundingBox() const override;
+
+    idx_t size() const override;
+    Spec spec() const override;
+
+    virtual std::unique_ptr<Grid::IteratorXY> xy_begin() const override {
+        return std::make_unique<IteratorXY>(*this);
+    }
+    virtual std::unique_ptr<Grid::IteratorXY> xy_end() const override {
+        return std::make_unique<IteratorXY>(*this, false);
+    }
+    virtual std::unique_ptr<Grid::IteratorLonLat> lonlat_begin() const override {
+        return std::make_unique<IteratorLonLat>(*this);
+    }
+    virtual std::unique_ptr<Grid::IteratorLonLat> lonlat_end() const override {
+        return std::make_unique<IteratorLonLat>(*this, false);
+    }
+
+    void xy(idx_t n, Point2& point) const;
+    Point2 xy(idx_t n) const;
+
+    void lonlat(idx_t n, Point2& point) const;
+    Point2 lonlat(idx_t n) const;
+
+protected:
+    void print(std::ostream&) const override;
+
+private:
+    using CSIndices = std::array<idx_t, 3>;
+    using PointAlphaBeta = Point2;
+
+    CSIndices get_cs_indices(gidx_t n) const;
+
+    PointAlphaBeta ij_to_curvilinear_coord(idx_t i, idx_t j) const;
+    PointXY curvilinear_to_tangent_coord(PointAlphaBeta& curvi_coord) const;
+    PointXYZ tangent_to_xyz_coord(PointXY& tan_coord, idx_t tile) const;
+
+protected:
+    // (N_ * N_) = number of cells on a tile
+    idx_t N_;
+
+    // Number of tiles
+    static constexpr idx_t nTiles_ = 6;
+
+private:
+    std::string type_ = {"cubedsphere2"};
+
+    using Matrix = eckit::maths::Matrix<double>;
+
+    std::array<Matrix, 6> lfric_rotations_ = {
+        Matrix({{0, 1, 0}, {0, 0, -1}, {1, 0, 0}}),
+        Matrix({{-1, 0, 0}, {0, 0, -1}, {0, 1, 0}}),
+        Matrix({{0, -1, 0}, {0, 0, -1}, {-1, 0, 0}}),
+        Matrix({{1, 0, 0}, {0, 0, -1}, {0, -1, 0}}),
+        Matrix({{-1, 0, 0}, {0, 1, 0}, {0, 0, 1}}),
+        Matrix({{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}})
+    };
+
+};
+
+}  // namespace grid
+}  // namespace detail
+}  // namespace grid
+}  // namespace atlas

src/tests/grid/CMakeLists.txt

@@ -41,6 +41,7 @@ foreach(test
         test_grids
         test_state
         test_cubedsphere
+        test_cubedsphere_2
         )
     ecbuild_add_test( TARGET atlas_${test} SOURCES ${test}.cc LIBS atlas ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} )
 endforeach()