eckit::geo reduced_gg

ecmwf · Oct 11, 2024 · ef71ccc · ef71ccc
1 parent c6ec373
commit ef71ccc
Show file tree

Hide file tree

Showing 5 changed files with 14 additions and 18 deletions.
diff --git a/src/eckit/geo/grid/ReducedGaussian.cc b/src/eckit/geo/grid/ReducedGaussian.cc
@@ -25,31 +25,27 @@
 namespace eckit::geo::grid {
 
 
-static size_t calculate_n(const pl_type& pl) {
-    ASSERT(!pl.empty() && pl.size() % 2 == 0);
-    return pl.size() / 2;
-}
-
-
 ReducedGaussian::ReducedGaussian(const Spec& spec) :
-    ReducedGaussian(spec.get_long_vector("pl"), area::BoundingBox(spec), projection::Rotation::make_from_spec(spec)) {}
+    ReducedGaussian(spec.get_long("N"), spec.get_long_vector("pl"), area::BoundingBox(spec),
+                    projection::Rotation::make_from_spec(spec)) {}
 
 
-ReducedGaussian::ReducedGaussian(const pl_type& pl, const area::BoundingBox& bbox, projection::Rotation* rotation) :
+ReducedGaussian::ReducedGaussian(size_t N, const pl_type& pl, const area::BoundingBox& bbox,
+                                 projection::Rotation* rotation) :
     Reduced(bbox, rotation),
-    N_(calculate_n(pl)),
+    N_(N),
     pl_(pl),
     j_(0),
-    Nj_(N_ * 2),
+    Nj_(pl.size()),
     x_(Nj_),
     y_(range::GaussianLatitude(N_, false).make_range_cropped(bbox.north, bbox.south)) {
-    ASSERT(Nj_ == pl_.size());
+    ASSERT(0 < N_ && Nj_ <= 2 * N_);
     ASSERT(y_);
 }
 
 
 ReducedGaussian::ReducedGaussian(size_t N, const area::BoundingBox& bbox, projection::Rotation* rotation) :
-    ReducedGaussian(util::reduced_octahedral_pl(N), bbox, rotation) {}
+    ReducedGaussian(N, util::reduced_octahedral_pl(N), bbox, rotation) {}
 
 
 Grid::iterator ReducedGaussian::cbegin() const {
@@ -124,7 +120,7 @@ void ReducedGaussian::fill_spec(spec::Custom& custom) const {
 
 Grid* ReducedGaussian::make_grid_cropped(const Area& crop) const {
     if (auto cropped(boundingBox()); crop.intersects(cropped)) {
-        return new ReducedGaussian(pl_, cropped);
+        return new ReducedGaussian(N_, pl_, cropped);
     }
 
     throw UserError("ReducedGaussian: cannot crop grid (empty intersection)", Here());

diff --git a/src/eckit/geo/grid/ReducedGaussian.h b/src/eckit/geo/grid/ReducedGaussian.h
@@ -28,7 +28,7 @@ class ReducedGaussian : public Reduced {
     // -- Constructors
 
     explicit ReducedGaussian(const Spec&);
-    explicit ReducedGaussian(const pl_type&, const area::BoundingBox& = {}, projection::Rotation* = nullptr);
+    explicit ReducedGaussian(size_t N, const pl_type&, const area::BoundingBox& = {}, projection::Rotation* = nullptr);
     explicit ReducedGaussian(size_t N, const area::BoundingBox& = {}, projection::Rotation* = nullptr);
 
     // -- Methods

diff --git a/src/eckit/geo/iterator/Reduced.cc b/src/eckit/geo/iterator/Reduced.cc
@@ -26,8 +26,8 @@ Reduced::Reduced(const Grid& grid, size_t index) :
     grid_(dynamic_cast<const grid::Reduced&>(grid)),
     latitudes_(grid_.latitudes()),
     niacc_(grid_.niacc()),
-    index_(index),
-    size_(grid.size()) {
+    index_(index) {
+    size_ = grid.size();
     if (index_ < size_) {
         longitudes_j_ = grid_.longitudes(j_ = j(index_));
         ASSERT(niacc_[j_] <= index && index_ < niacc_[j_ + 1]);

diff --git a/src/eckit/geo/iterator/Reduced.h b/src/eckit/geo/iterator/Reduced.h
@@ -40,7 +40,7 @@ class Reduced final : public geo::Iterator {
     const std::vector<size_t>& niacc_;
     size_t j_;
     size_t index_;
-    const size_t size_;
+    size_t size_;
 
     // -- Overridden operators
 

diff --git a/tests/geo/grid_reduced_gg.cc b/tests/geo/grid_reduced_gg.cc
@@ -165,7 +165,7 @@ CASE("equals") {
     std::unique_ptr<const Grid> grid1(GridFactory::build(spec::Custom({{"grid", "o3"}})));
     std::unique_ptr<const Grid> grid2(GridFactory::make_from_string("N: 3"));
     std::unique_ptr<const Grid> grid3(new ReducedGaussian(3));
-    std::unique_ptr<const Grid> grid4(new ReducedGaussian(pl_type{20, 24, 28, 28, 24, 20}));
+    std::unique_ptr<const Grid> grid4(new ReducedGaussian(3, pl_type{20, 24, 28, 28, 24, 20}));
 
     EXPECT(*grid1 == *grid2);
     EXPECT(*grid2 == *grid3);