MIR-268 MIR-294 MIR-592 MIR-681 Correction at the poles for input dat…

…a spectral U/V

src/mir/CMakeLists.txt

@@ -41,6 +41,8 @@ list(APPEND mir_srcs
     action/context/Context.h
     action/filter/AddRandomFilter.cc
     action/filter/AddRandomFilter.h
+    action/filter/AdjustWindsAtPoles.cc
+    action/filter/AdjustWindsAtPoles.h
     action/filter/AdjustWindsDirections.cc
     action/filter/AdjustWindsDirections.h
     action/filter/AdjustWindsScaleCosLatitude.cc

src/mir/action/filter/AdjustWindsAtPoles.cc

@@ -0,0 +1,99 @@
+/*
+ * (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.
+ */
+
+
+#include "mir/action/filter/AdjustWindsAtPoles.h"
+
+#include <limits>
+#include <memory>
+#include <ostream>
+
+#include "mir/action/context/Context.h"
+#include "mir/data/MIRField.h"
+#include "mir/repres/Iterator.h"
+#include "mir/repres/Representation.h"
+#include "mir/util/Exceptions.h"
+
+
+namespace mir::action {
+
+
+bool AdjustWindsAtPoles::sameAs(const Action& other) const {
+    const auto* o = dynamic_cast<const AdjustWindsAtPoles*>(&other);
+    return (o != nullptr);
+}
+
+
+void AdjustWindsAtPoles::print(std::ostream& out) const {
+    out << "AdjustWindsAtPoles[]";
+}
+
+
+void AdjustWindsAtPoles::execute(context::Context& ctx) const {
+    auto& field = ctx.field();
+
+    ASSERT(field.dimensions() % 2 == 0);
+    ASSERT(field.dimensions() > 0);
+
+    const auto missingValue = field.hasMissing() ? field.missingValue() : std::numeric_limits<double>::quiet_NaN();
+
+    const auto N = field.values(0).size();
+    ASSERT(N > 0);
+
+    // find points at pole(s)
+    const auto* representation(field.representation());
+    ASSERT(representation);
+
+    std::vector<size_t> north;
+    std::vector<size_t> south;
+
+    for (const std::unique_ptr<repres::Iterator> it(representation->iterator()); it->next();) {
+        const auto lat = it->pointUnrotated().lat();
+        if (lat == Latitude::NORTH_POLE) {
+            north.push_back(it->index());
+        }
+        else if (lat == Latitude::SOUTH_POLE) {
+            south.push_back(it->index());
+        }
+    }
+
+
+    // set v to zero at pole(s)
+    for (size_t i = 0; i < field.dimensions(); i += 2) {
+        auto& u = field.direct(i);
+        auto& v = field.direct(i + 1);
+        ASSERT(u.size() == N);
+        ASSERT(v.size() == N);
+
+        for (const auto& indices : {north, south}) {
+            for (const auto& idx : indices) {
+                if (u[idx] == missingValue || v[idx] == missingValue) {
+                    u[idx] = missingValue;
+                    v[idx] = missingValue;
+                }
+                else {
+                    v[idx] = 0;
+                }
+            }
+        }
+    }
+}
+
+
+const char* AdjustWindsAtPoles::name() const {
+    return "AdjustWindsAtPoles";
+}
+
+
+static const ActionBuilder<AdjustWindsAtPoles> __action("filter.adjust-winds-at-poles");
+
+
+}  // namespace mir::action

src/mir/action/filter/AdjustWindsAtPoles.h

@@ -0,0 +1,37 @@
+/*
+ * (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.
+ */
+
+
+#pragma once
+
+#include "mir/action/plan/Action.h"
+
+
+namespace mir::action {
+
+
+class AdjustWindsAtPoles : public Action {
+public:
+    // -- Constructors
+
+    using Action::Action;
+
+private:
+    // -- Overridden methods
+
+    void execute(context::Context&) const override;
+    bool sameAs(const Action&) const override;
+    const char* name() const override;
+    void print(std::ostream&) const override;
+};
+
+
+}  // namespace mir::action

src/mir/input/GribInput.cc

@@ -377,14 +377,20 @@ struct ProcessingT {
 };
 
 
+struct WindCache {
+    static const util::Wind::Defaults& defaults() {
+        static const util::Wind::Defaults def;
+        return def;
+    };
+} static const WIND;
+
+
 static ProcessingT<long>* is_wind_component_uv() {
     return new ProcessingT<long>([](grib_handle* h, long& value) {
         long paramId = 0;
         GRIB_CALL(codes_get_long(h, "paramId", &paramId));
-        static const util::Wind::Defaults def;
-        long ind = paramId % 1000;
-        value    = (ind == def.u ? 1 : ind == def.v ? 2 : 0);
-        return value;
+        const auto id = paramId % 1000;
+        return value  = id == WIND.defaults().u ? 1 : id == WIND.defaults().v ? 2 : 0;
     });
 }
 
@@ -393,10 +399,8 @@ static ProcessingT<long>* is_wind_component_vod() {
     return new ProcessingT<long>([](grib_handle* h, long& value) {
         long paramId = 0;
         GRIB_CALL(codes_get_long(h, "paramId", &paramId));
-        static const util::Wind::Defaults def;
-        long ind = paramId % 1000;
-        value    = (ind == def.vo ? 1 : ind == def.d ? 2 : 0);
-        return value;
+        const auto id = paramId % 1000;
+        return value  = id == WIND.defaults().vo ? 1 : id == WIND.defaults().d ? 2 : 0;
     });
 }
 

src/mir/key/style/ECMWFStyle.cc

@@ -311,6 +311,10 @@ void ECMWFStyle::sh2grid(action::ActionPlan& plan) const {
             }
         }
 
+        if (uv2uv) {
+            plan.add("filter.adjust-winds-at-poles");
+        }
+
         if (vod2uv || uv2uv) {
             if (rotation) {
                 plan.add("filter.adjust-winds-directions");