diff --git a/src/mir/CMakeLists.txt b/src/mir/CMakeLists.txt index 58777e0ca..5c9a15535 100644 --- a/src/mir/CMakeLists.txt +++ b/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 diff --git a/src/mir/action/filter/AdjustWindsAtPoles.cc b/src/mir/action/filter/AdjustWindsAtPoles.cc new file mode 100644 index 000000000..d23dce342 --- /dev/null +++ b/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 +#include +#include + +#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(&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::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 north; + std::vector south; + + for (const std::unique_ptr 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 __action("filter.adjust-winds-at-poles"); + + +} // namespace mir::action diff --git a/src/mir/action/filter/AdjustWindsAtPoles.h b/src/mir/action/filter/AdjustWindsAtPoles.h new file mode 100644 index 000000000..96552dd34 --- /dev/null +++ b/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 diff --git a/src/mir/input/GribInput.cc b/src/mir/input/GribInput.cc index d58742c6a..6918326f5 100644 --- a/src/mir/input/GribInput.cc +++ b/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* is_wind_component_uv() { return new ProcessingT([](grib_handle* h, long& value) { long paramId = 0; GRIB_CALL(codes_get_long(h, "paramId", ¶mId)); - 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* is_wind_component_vod() { return new ProcessingT([](grib_handle* h, long& value) { long paramId = 0; GRIB_CALL(codes_get_long(h, "paramId", ¶mId)); - 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; }); } diff --git a/src/mir/key/style/ECMWFStyle.cc b/src/mir/key/style/ECMWFStyle.cc index 043b3698b..3b68b9fc1 100644 --- a/src/mir/key/style/ECMWFStyle.cc +++ b/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");