Update interpolate action to output HEALPix_nested directly

src/multio/action/interpolate/CMakeLists.txt

@@ -27,6 +27,8 @@ ecbuild_add_executable(
 
     SOURCES
         FesomCacheGenerator.cc
+        HEALPix.cc
+        HEALPix.h
         ../../tools/MultioTool.cc
 
     CONDITION

src/multio/action/interpolate/FesomCacheGenerator.cc

@@ -22,6 +22,8 @@
 #include "mir/method/WeightMatrix.h"
 #include "multio/tools/MultioTool.h"
 
+#include "HEALPix.h"
+
 namespace multio::action::interpolateFESOM2MIR {
 
 namespace {
@@ -96,7 +98,7 @@ class Fesom2mirCacheGenerator final : public multio::MultioTool {
     void usage(const std::string& tool) const override {
         eckit::Log::info() << std::endl << "Usage: " << tool << " [options]" << std::endl;
         eckit::Log::info() << "EXAMPLE: " << std::endl
-                           << "fesom-cache-generator --inputPath=. --inputFile=CORE2_ngrid_NSIDE32_0_ring.csv "
+                           << tool << " --inputPath=. --inputFile=CORE2_ngrid_NSIDE32_0_ring.csv --outputPath=. --nCols=126858 --inputOrdering=ring --outputOrdering=nested"
                            << std::endl
                            << std::endl;
     }
@@ -122,30 +124,30 @@ class Fesom2mirCacheGenerator final : public multio::MultioTool {
     size_t level_;
     size_t Nrow_;
     size_t Ncol_;
-    orderingConvention_e orderingConvention_;
+    orderingConvention_e inputOrdering_;
+    orderingConvention_e outputOrdering_;
+
 
     void loadTriplets(std::vector<eckit::linalg::Triplet>& triplets) const;
 };
 
 
 Fesom2mirCacheGenerator::Fesom2mirCacheGenerator(int argc, char** argv) :
-    multio::MultioTool{argc, argv},
-    inputPath_{"."},
-    outputPath_{"."},
-    inputFile_{"CORE2_ngrid_NSIDE32_0_ring.csv"},
-    fesomName_{"CORE2"},
-    domain_{"ngrid"},
-    NSide_{0},
-    level_{0},
-    orderingConvention_{orderingConvention_e::RING} {
+    multio::MultioTool{argc, argv}
+    {
 
     options_.push_back(new eckit::option::SimpleOption<std::string>(
-        "inputPath", "Path of the input files with the triplets. Default( \".\" )"));
+        "inputPath", "Path of the input files with the triplets. Default( \".\" )", "."));
+    options_.push_back(new eckit::option::SimpleOption<std::string>(
+        "outputPath", "Path of the output files with the triplets. Default( \".\" )", "."));
+    options_.push_back(new eckit::option::SimpleOption<std::string>(
+        "inputFile", "Name of the input file. Default( \"CORE2_ngrid_NSIDE32_0_ring.csv\" )", "CORE2_ngrid_NSIDE32_0_ring.csv"));
+    options_.push_back(new eckit::option::SimpleOption<size_t>(
+        "nCols", "Size of the fesom grid.", 126858));
     options_.push_back(new eckit::option::SimpleOption<std::string>(
-        "outputPath", "Path of the output files with the triplets. Default( \".\" )"));
+        "inputOrdering", "Ordering of the input files. Options( \"ring\", \"nested\") Default( \"ring\" )", "ring"));
     options_.push_back(new eckit::option::SimpleOption<std::string>(
-        "inputFile", "Name of the input file. Default( \"CORE2_ngrid_NSIDE32_0_ring.csv\" )"));
-    options_.push_back(new eckit::option::SimpleOption<size_t>("nCols", "Size of the fesom grid."));
+        "outputOrdering", "Ordering of the output files. Options( \"ring\", \"nested\", \"input\") Default( \"input\" )", "input"));
 
     return;
 }
@@ -167,7 +169,20 @@ void Fesom2mirCacheGenerator::loadTriplets(std::vector<eckit::linalg::Triplet>&
             "([0-9][0-9]*)\\s+([0-9][0-9]*)\\s*([+]?([0-9]*[.])?[0-9]+([eE][-+][0-9]+)?)");
         std::smatch matchLine;
         if (std::regex_match(line, matchLine, lineGrammar)) {
-            iHEALPix = std::stoi(matchLine[1].str());
+            if ( inputOrdering_ ==  outputOrdering_ ){
+              iHEALPix = std::stoi(matchLine[1].str());
+            }
+            else if (inputOrdering_ == orderingConvention_e::RING && outputOrdering_ == orderingConvention_e::NESTED) {
+              HEALPix Idx(static_cast<int>(NSide_));
+              iHEALPix = Idx.ring_to_nest(static_cast<int>(std::stoi(matchLine[1].str())));
+            }
+            else if (inputOrdering_ == orderingConvention_e::NESTED && outputOrdering_ == orderingConvention_e::RING){
+              HEALPix Idx(static_cast<int>(NSide_));
+              iHEALPix = Idx.nest_to_ring(static_cast<int>(std::stoi(matchLine[1].str())));
+            }
+            else {
+              throw eckit::SeriousBug("Ordering not supported: " + line, Here());
+            }
             iFESOM = std::stoi(matchLine[2].str());
             weight = std::stod(matchLine[3].str());
         }
@@ -185,13 +200,27 @@ void Fesom2mirCacheGenerator::init(const eckit::option::CmdArgs& args) {
     args.get("outputPath", outputPath_);
     args.get("inputFile", inputFile_);
 
+    std::string inputOrderingType;
+    std::string outputOrderingType;
+    args.get("inputOrdering", inputOrderingType);
+    args.get("outputOrdering", outputOrderingType);
+
     eckit::PathName inputPath_tmp{inputPath_};
     ASSERT(inputPath_tmp.exists());
     eckit::PathName inputFile_tmp{inputPath_ + "/" + inputFile_};
     ASSERT(inputFile_tmp.exists());
     eckit::PathName outputPath_tmp{outputPath_};
     outputPath_tmp.mkdir();
-    parseInputFileName(inputFile_, fesomName_, domain_, NSide_, level_, orderingConvention_);
+    parseInputFileName(inputFile_, fesomName_, domain_, NSide_, level_, inputOrdering_);
+
+    if (inputOrderingType == "ring") { inputOrdering_ = orderingConvention_e::RING; }
+    else if (inputOrderingType == "nested") { inputOrdering_ = orderingConvention_e::NESTED; }
+    else { throw eckit::SeriousBug("Unsupported input ordering convention", Here()); }
+
+    if (outputOrderingType == "input") { outputOrdering_ = inputOrdering_; }
+    else if (outputOrderingType == "ring") { outputOrdering_ = orderingConvention_e::RING; }
+    else if (outputOrderingType == "nested") { outputOrdering_ = orderingConvention_e::NESTED; }
+    else { throw eckit::SeriousBug("Unsupported output ordering convention", Here()); }
 
     args.get("nCols", Ncol_);
     Nrow_ = NSide_ * NSide_ * 12;
@@ -203,7 +232,7 @@ void Fesom2mirCacheGenerator::execute(const eckit::option::CmdArgs& args) {
 
     std::sort(begin(triplets), end(triplets), [](const auto& a, const auto& b) { return a.row() < b.row(); });
 
-    const auto orderingConvention = orderingConvention_enum2string(orderingConvention_);
+    const auto orderingConvention = orderingConvention_enum2string(outputOrdering_);
     const auto cacheFileName = fesomCacheName(fesomName_, domain_, "double", NSide_, orderingConvention, level_);
 
     mir::method::WeightMatrix W(Nrow_, Ncol_);

src/multio/action/interpolate/HEALPix.cc

@@ -0,0 +1,198 @@
+/*
+ * (C) Copyright 1996- ECMWF.
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation nor
+ * does it submit to any jurisdiction.
+ */
+
+/// @author Pedro Maciel
+
+/// @date Sep 2023
+
+#include "HEALPix.h"
+
+#include <bitset>
+#include <cassert>
+#include <cmath>
+#include <cstdint>
+#include <iostream>
+#include <tuple>
+
+namespace {
+
+
+struct CodecFijNest {
+    static constexpr uint64_t __masks[] = {0x00000000ffffffff, 0x0000ffff0000ffff, 0x00ff00ff00ff00ff,
+                                           0x0f0f0f0f0f0f0f0f, 0x3333333333333333, 0x5555555555555555};
+
+    inline static int nest_encode_bits(int n) {
+        auto b = static_cast<uint64_t>(n) & __masks[0];
+        b = (b ^ (b << 16)) & __masks[1];
+        b = (b ^ (b << 8)) & __masks[2];
+        b = (b ^ (b << 4)) & __masks[3];
+        b = (b ^ (b << 2)) & __masks[4];
+        b = (b ^ (b << 1)) & __masks[5];
+        return static_cast<int>(b);
+    }
+
+    inline static int nest_decode_bits(int n) {
+        auto b = static_cast<uint64_t>(n) & __masks[5];
+        b = (b ^ (b >> 1)) & __masks[4];
+        b = (b ^ (b >> 2)) & __masks[3];
+        b = (b ^ (b >> 4)) & __masks[2];
+        b = (b ^ (b >> 8)) & __masks[1];
+        b = (b ^ (b >> 16)) & __masks[0];
+        return static_cast<int>(b);
+    }
+
+    static std::tuple<int, int, int> nest_to_fij(int n, int k) {
+        assert(0 <= n);
+        auto f = n >> (2 * k);    // f = n / (Nside * Nside)
+        n &= (1 << (2 * k)) - 1;  // n = n % (Nside * Nside)
+        auto i = nest_decode_bits(n);
+        auto j = nest_decode_bits(n >> 1);
+        return {f, i, j};
+    }
+
+    static int fij_to_nest(int f, int i, int j, int k) {
+        return (f << (2 * k)) + nest_encode_bits(i) + (nest_encode_bits(j) << 1);
+    }
+};
+
+
+inline int sqrt(int n) {
+    return static_cast<int>(std::sqrt(static_cast<double>(n) + 0.5));
+}
+
+
+// for division result within [0; 3]
+inline int div_03(int a, int b) {
+    int t = (a >= (b << 1)) ? 1 : 0;
+    a -= t * (b << 1);
+    return (t << 1) + (a >= b ? 1 : 0);
+}
+
+
+inline bool is_power_of_2(int n) {
+    return std::bitset<sizeof(int) * 8>(n).count() == 1;
+}
+
+
+inline int pll(int f) {
+    constexpr int __pll[] = {1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};
+    return __pll[f];
+}
+
+
+}  // unnamed namespace
+
+
+HEALPix::HEALPix(int Nside) :
+    Nside_(Nside),
+    Npix_(size()),
+    Ncap_((Nside * (Nside - 1)) << 1),
+    k_(is_power_of_2(Nside_) ? static_cast<int>(std::log2(Nside)) : -1) {
+    assert(0 <= k_);  // (specific to nested ordering)
+    assert(0 < Nside_);
+}
+
+
+int HEALPix::ring_to_nest(int r) const {
+    auto to_nest = [&](int f,      //!< base pixel index
+                       int ring,   //!< 1-based ring number
+                       int Nring,  //!< number of pixels in ring
+                       int phi,    //!< index in longitude
+                       int shift   //!< if ring's first pixel is not at phi=0
+                       ) -> int {
+        int r = ((2 + (f >> 2)) << k_) - ring - 1;
+        int p = 2 * phi - pll(f) * Nring - shift - 1;
+        if (p >= 2 * Nside_) {
+            p -= 8 * Nside_;
+        }
+
+        int i = (r + p) >> 1;
+        int j = (r - p) >> 1;
+
+        assert(f < 12 && i < Nside_ && j < Nside_);
+        return CodecFijNest::fij_to_nest(f, i, j, k_);
+    };
+
+    if (r < Ncap_) {
+        // North polar cap
+        int Nring = (1 + sqrt(2 * r + 1)) >> 1;
+        int phi = 1 + r - 2 * Nring * (Nring - 1);
+        int r = Nring;
+        int f = div_03(phi - 1, Nring);
+
+        return to_nest(f, r, Nring, phi, 0);
+    }
+
+    if (Npix_ - Ncap_ <= r) {
+        // South polar cap
+        int Nring = (1 + sqrt(2 * Npix_ - 2 * r - 1)) >> 1;
+        int phi = 1 + r + 2 * Nring * (Nring - 1) + 4 * Nring - Npix_;
+        int ring = 4 * Nside_ - Nring;  // (from South pole)
+        int f = div_03(phi - 1, Nring) + 8;
+
+        return to_nest(f, ring, Nring, phi, 0);
+    }
+
+    // Equatorial belt
+    int ip = r - Ncap_;
+    int tmp = ip >> (k_ + 2);
+
+    int Nring = Nside_;
+    int phi = ip - tmp * 4 * Nside_ + 1;
+    int ring = tmp + Nside_;
+
+    int ifm = 1 + ((phi - 1 - ((1 + tmp) >> 1)) >> k_);
+    int ifp = 1 + ((phi - 1 - ((1 - tmp + 2 * Nside_) >> 1)) >> k_);
+    int f = (ifp == ifm) ? (ifp | 4) : ((ifp < ifm) ? ifp : (ifm + 8));
+
+    return to_nest(f, ring, Nring, phi, ring & 1);
+}
+
+
+int HEALPix::nest_to_ring(int n) const {
+    auto [f, i, j] = CodecFijNest::nest_to_fij(n, k_);
+    assert(f < 12 && i < Nside_ && j < Nside_);
+
+    auto to_ring_local = [&](int f, int i, int j,
+                             int Nring,  //!< number of pixels in ring
+                             int shift   //!< if ring's first pixel is/is not at phi=0
+                             ) -> int {
+        Nring >>= 2;
+        int r = (pll(f) * Nring + i - j + 1 + shift) / 2 - 1;
+        assert(r < 4 * Nring);
+
+        return r < 0 ? r + 4 * Nside_ : r;
+    };
+
+    const int ring = ((f >> 2) + 2) * Nside_ - i - j - 1;  // 1-based ring number
+    if (ring < Nside_) {
+        // North polar cap
+        int Nring = 4 * ring;
+        int r0 = 2 * ring * (ring - 1);  // index of first ring pixel (ring numbering)
+
+        return r0 + to_ring_local(f, i, j, Nring, 0);
+    }
+
+    if (ring < 3 * Nside_) {
+        // South polar cap
+        int Nring = 4 * Nside_;
+        int r0 = Ncap_ + (ring - Nside_) * Nring;  // index of first ring pixel (ring numbering)
+        int shift = (ring - Nside_) & 1;
+
+        return r0 + to_ring_local(f, i, j, Nring, shift);
+    }
+
+    // Equatorial belt
+    int N = 4 * Nside_ - ring;
+    int Nring = 4 * N;
+    int r0 = Npix_ - 2 * N * (N + 1);  // index of first ring pixel (ring numbering)
+
+    return r0 + to_ring_local(f, i, j, Nring, 0);
+}

src/multio/action/interpolate/HEALPix.h

@@ -0,0 +1,32 @@
+/*
+ * (C) Copyright 1996- ECMWF.
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation nor
+ * does it submit to any jurisdiction.
+ */
+
+/// @author Pedro Maciel
+
+/// @date Sep 2023
+
+#pragma once
+
+class HEALPix {
+public:
+    explicit HEALPix(int Nside);
+
+    int size() const { return 12 * Nside_ * Nside_; }
+    int nside() const { return Nside_; }
+
+    int nest_to_ring(int) const;
+    int ring_to_nest(int) const;
+
+private:
+    const int Nside_;  // up to 2^13
+    const int Npix_;
+    const int Ncap_;
+    const int k_;
+};