diff --git a/CMakeLists.txt b/CMakeLists.txt index 4498386df..5470681c7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -21,7 +21,7 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON) ecbuild_add_option(FEATURE BUILD_TOOLS DEFAULT OFF DESCRIPTION "build the command line tools (install)") ecbuild_add_option(FEATURE MIR_DOWNLOAD_MASKS DEFAULT ON DESCRIPTION "download LSM files") -ecbuild_add_option(FEATURE NETCDF DEFAULT ON DESCRIPTION "support for netCDF" REQUIRED_PACKAGES "NetCDF COMPONENTS C") +ecbuild_add_option(FEATURE NETCDF DEFAULT ON DESCRIPTION "support for netCDF" REQUIRED_PACKAGES "NetCDF COMPONENTS C CXX") ecbuild_add_option(FEATURE PNG DEFAULT OFF DESCRIPTION "support for PNG encoding" REQUIRED_PACKAGES PNG) ecbuild_add_option(FEATURE OMP DEFAULT OFF DESCRIPTION "support for OpenMP shared memory parallelism" REQUIRED_PACKAGES "OpenMP COMPONENTS CXX") diff --git a/VERSION b/VERSION index a6c2798a4..53cc1a6f9 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -1.23.0 +1.24.0 diff --git a/etc/mir/grib-input.yaml b/etc/mir/grib-input.yaml index 1d0302bb7..fe6a7bb31 100644 --- a/etc/mir/grib-input.yaml +++ b/etc/mir/grib-input.yaml @@ -19,6 +19,10 @@ md5GridSection=f5dc74ec36353f4c83f7de3bf46e1aef, centre:i=98, gridName=O640, edi md5GridSection=6ae0076681a2541abe053ab52fa60f30, centre:i=84/85, gridType=lambert, edition=1: - first_point_bottom_left: true +# gridType=polar_stereographic does not support the default area method +gridType=polar_stereographic: + - area-mode: mask + # gridType=space_view does not support the default interpolation method gridType=space_view: - interpolation: k-nearest diff --git a/etc/mir/parameter-class.yaml b/etc/mir/parameter-class.yaml index 40b516700..d4196082e 100644 --- a/etc/mir/parameter-class.yaml +++ b/etc/mir/parameter-class.yaml @@ -1952,6 +1952,7 @@ intensity: - 261020 # Physiological equivalent temperature - 261021 # Saturation water vapor pressure - 261022 # Wet-bulb potential temperature + - 261023 # Wet-bulb temperature - 262000 # Sea ice thickness - 262002 # Snow thickness over sea ice - 262003 # Eastward sea ice velocity diff --git a/src/mir/CMakeLists.txt b/src/mir/CMakeLists.txt index eff7e1629..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 @@ -246,6 +248,8 @@ list(APPEND mir_srcs input/GribMemoryInput.h input/GribStreamInput.cc input/GribStreamInput.h + input/GriddefInput.cc + input/GriddefInput.h input/MIRInput.cc input/MIRInput.h input/MultiDimensionalGribFileInput.cc @@ -457,6 +461,8 @@ list(APPEND mir_srcs output/GribOutput.h output/GribStreamOutput.cc output/GribStreamOutput.h + output/GriddefOutput.cc + output/GriddefOutput.h output/MIROutput.cc output/MIROutput.h output/MultiDimensionalOutput.cc diff --git a/src/mir/action/area/Area.h b/src/mir/action/area/Area.h index 7b36b270c..a5062e869 100644 --- a/src/mir/action/area/Area.h +++ b/src/mir/action/area/Area.h @@ -38,7 +38,7 @@ class Area : public Action { // -- Destructor - virtual ~Area() override; + ~Area() override; // -- Convertors // None @@ -72,7 +72,7 @@ class Area : public Action { // -- Overridden methods - virtual bool sameAs(const Action&) const override; + bool sameAs(const Action&) const override; // -- Class members // None diff --git a/src/mir/action/filter/AdjustWindsAtPoles.cc b/src/mir/action/filter/AdjustWindsAtPoles.cc new file mode 100644 index 000000000..6fcdecc7a --- /dev/null +++ b/src/mir/action/filter/AdjustWindsAtPoles.cc @@ -0,0 +1,114 @@ +/* + * (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() > 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()); + } + } + + if (field.dimensions() % 2 == 0) { + // two-component mode: 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; + } + } + } + } + } + else { + // single-component mode: set u, v to zero at pole(s) + for (size_t i = 0; i < field.dimensions(); ++i) { + auto& v = field.direct(i); + ASSERT(v.size() == N); + + for (const auto& indices : {north, south}) { + for (const auto& idx : indices) { + if (v[idx] != missingValue) { + 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/action/filter/AdjustWindsScaleCosLatitude.cc b/src/mir/action/filter/AdjustWindsScaleCosLatitude.cc index 2d2d5ec58..e9fd10d9c 100644 --- a/src/mir/action/filter/AdjustWindsScaleCosLatitude.cc +++ b/src/mir/action/filter/AdjustWindsScaleCosLatitude.cc @@ -14,7 +14,6 @@ #include #include -#include #include #include @@ -45,24 +44,24 @@ void AdjustWindsScaleCosLatitude::print(std::ostream& out) const { void AdjustWindsScaleCosLatitude::execute(context::Context& ctx) const { - data::MIRField& field = ctx.field(); - ASSERT(field.dimensions() > 0); + auto& field = ctx.field(); + ASSERT(field.dimensions() > 0); + ASSERT(!field.hasMissing()); // determine scaling vector (all fields share the same representation) const size_t N = field.values(0).size(); ASSERT(N > 0); - const repres::Representation* representation(field.representation()); + const auto* representation(field.representation()); ASSERT(representation); - std::vector scale(N, std::numeric_limits::quiet_NaN()); + std::vector scale(N); for (const std::unique_ptr it(representation->iterator()); it->next();) { - const auto& p = it->pointUnrotated(); - scale[it->index()] = (p.lat() == Latitude::SOUTH_POLE) ? 0. - : (p.lat() == Latitude::NORTH_POLE) + const auto lat = it->pointUnrotated().lat(); + scale[it->index()] = (lat == Latitude::SOUTH_POLE || lat == Latitude::NORTH_POLE) ? 0. - : 1. / std::cos(util::degree_to_radian(p.lat().value())); + : 1. / std::cos(util::degree_to_radian(lat.value())); } // apply scaling to each field component directly @@ -70,10 +69,8 @@ void AdjustWindsScaleCosLatitude::execute(context::Context& ctx) const { auto& values = field.direct(i); ASSERT(values.size() == N); - std::transform( - values.begin(), values.end(), scale.begin(), values.begin(), [](const double& a, const double& b) { - return std::isfinite(a) && std::isfinite(b) ? a * b : std::numeric_limits::quiet_NaN(); - }); + std::transform(values.begin(), values.end(), scale.begin(), values.begin(), + [](const double& a, const double& b) { return a * b; }); } } diff --git a/src/mir/action/io/EndAction.h b/src/mir/action/io/EndAction.h index 690386a7d..a8cc0fb17 100644 --- a/src/mir/action/io/EndAction.h +++ b/src/mir/action/io/EndAction.h @@ -35,7 +35,7 @@ class EndAction : public Action { // -- Destructor - virtual ~EndAction() override; + ~EndAction() override; // -- Convertors // None @@ -87,8 +87,8 @@ class EndAction : public Action { bool sameAs(const Action&) const override; bool isEndAction() const override; - virtual void print(std::ostream&) const override; - virtual void custom(std::ostream&) const override; + void print(std::ostream&) const override; + void custom(std::ostream&) const override; // -- Class members // None diff --git a/src/mir/api/mir_config.h.in b/src/mir/api/mir_config.h.in index d5555e8aa..2bfaed5d4 100644 --- a/src/mir/api/mir_config.h.in +++ b/src/mir/api/mir_config.h.in @@ -1,11 +1,36 @@ +/* + * (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/mir_ecbuild_config.h" #include "mir/api/mir_version.h" + #cmakedefine01 mir_HAVE_ATLAS #cmakedefine01 mir_HAVE_NETCDF #cmakedefine01 mir_HAVE_PNG #cmakedefine01 mir_HAVE_OMP + +constexpr bool _to_bool(int x = 0) { + return x != 0; +} + +constexpr bool HAVE_ATLAS = _to_bool(@mir_HAVE_ATLAS@); +constexpr bool HAVE_NETCDF = _to_bool(@mir_HAVE_NETCDF@); +constexpr bool HAVE_PNG = _to_bool(@mir_HAVE_PNG@); +constexpr bool HAVE_OMP = _to_bool(@mir_HAVE_OMP@); + +constexpr bool HAVE_TESSELATION = _to_bool(@atlas_HAVE_TESSELATION@); +constexpr bool HAVE_PROJ = _to_bool(@atlas_HAVE_PROJ@); diff --git a/src/mir/caching/LegendreCache.cc b/src/mir/caching/LegendreCache.cc index 1f9aa6bc1..37ac4278e 100644 --- a/src/mir/caching/LegendreCache.cc +++ b/src/mir/caching/LegendreCache.cc @@ -12,6 +12,8 @@ #include "mir/caching/LegendreCache.h" +#include "eckit/config/Resource.h" + #include "mir/config/LibMir.h" diff --git a/src/mir/caching/MeshCache.cc b/src/mir/caching/MeshCache.cc index 9101d88b7..5ec66b851 100644 --- a/src/mir/caching/MeshCache.cc +++ b/src/mir/caching/MeshCache.cc @@ -12,6 +12,8 @@ #include "mir/caching/MeshCache.h" +#include "eckit/config/Resource.h" + #include "mir/config/LibMir.h" diff --git a/src/mir/caching/WeightCache.cc b/src/mir/caching/WeightCache.cc index 4c3a4c2ec..a67508214 100644 --- a/src/mir/caching/WeightCache.cc +++ b/src/mir/caching/WeightCache.cc @@ -14,6 +14,8 @@ #include +#include "eckit/config/Resource.h" + #include "mir/caching/matrix/MatrixLoader.h" #include "mir/config/LibMir.h" #include "mir/method/WeightMatrix.h" @@ -81,10 +83,7 @@ void WeightCacheTraits::load(const eckit::CacheManagerBase& manager, value_type& value_type tmp(matrix::MatrixLoaderFactory::build(manager.loader(), path)); w.swap(tmp); - static bool matrixValidate = eckit::Resource("$MIR_MATRIX_VALIDATE", false); - if (matrixValidate) { - w.validate("fromCache"); - } + w.validate("fromCache"); // check matrix structure (only) } diff --git a/src/mir/caching/legendre/MappedMemoryLoader.h b/src/mir/caching/legendre/MappedMemoryLoader.h index 7d66fdfca..7f9aae6e7 100644 --- a/src/mir/caching/legendre/MappedMemoryLoader.h +++ b/src/mir/caching/legendre/MappedMemoryLoader.h @@ -30,9 +30,9 @@ class MappedMemoryLoader : public LegendreLoader { void print(std::ostream&) const override; private: - virtual const void* address() const override; - virtual size_t size() const override; - virtual bool inSharedMemory() const override; + const void* address() const override; + size_t size() const override; + bool inSharedMemory() const override; private: int fd_; diff --git a/src/mir/caching/legendre/NoLoader.h b/src/mir/caching/legendre/NoLoader.h index 4c259d3dd..8dad2106c 100644 --- a/src/mir/caching/legendre/NoLoader.h +++ b/src/mir/caching/legendre/NoLoader.h @@ -30,9 +30,9 @@ class NoLoader : public LegendreLoader { void print(std::ostream&) const override; private: - virtual const void* address() const override; - virtual size_t size() const override; - virtual bool inSharedMemory() const override; + const void* address() const override; + size_t size() const override; + bool inSharedMemory() const override; }; diff --git a/src/mir/input/ArtificialInput.cc b/src/mir/input/ArtificialInput.cc index f8d261ea6..80d34654f 100644 --- a/src/mir/input/ArtificialInput.cc +++ b/src/mir/input/ArtificialInput.cc @@ -12,9 +12,13 @@ #include "mir/input/ArtificialInput.h" -#include +#include +#include #include +#include "mir/data/MIRField.h" +#include "mir/input/GriddefInput.h" +#include "mir/param/CombinedParametrisation.h" #include "mir/param/SimpleParametrisation.h" #include "mir/util/Exceptions.h" #include "mir/util/Log.h" @@ -47,6 +51,12 @@ size_t ArtificialInput::dimensions() const { } +ArtificialInput::ArtificialInput() : + inputParametrisation_(new param::CombinedParametrisation(parametrisation_, *this, parametrisation_)), calls_(0) { + ASSERT(inputParametrisation_); +} + + param::SimpleParametrisation& ArtificialInput::parametrisation(size_t which) { ASSERT(which == 0); return parametrisation_; @@ -55,7 +65,7 @@ param::SimpleParametrisation& ArtificialInput::parametrisation(size_t which) { const param::MIRParametrisation& ArtificialInput::parametrisation(size_t which) const { ASSERT(which == 0); - return parametrisation_; + return *inputParametrisation_; } @@ -70,6 +80,7 @@ void ArtificialInput::setAuxiliaryInformation(const util::ValueMap& map) { parametrisation_.set("grid", grid); } } + if (!parametrisation_.has("rotation")) { std::vector rotation(2, 0); if (parametrisation_.get("south_pole_latitude", rotation[0]) && @@ -77,77 +88,64 @@ void ArtificialInput::setAuxiliaryInformation(const util::ValueMap& map) { parametrisation_.set("rotation", rotation); } } -} + auto load = [](const eckit::PathName& path, std::vector& values) { + Log::info() << "ArtificialInput::setAuxiliaryInformation: '" << path << "'" << std::endl; -void ArtificialInput::print(std::ostream& out) const { - out << "ArtificialInput[parametrisation=" << parametrisation_ << "]"; -} + static const param::SimpleParametrisation empty; + std::unique_ptr input(input::MIRInputFactory::build(path.asString(), empty)); + ASSERT(input->next()); + auto field = input->field(); + ASSERT(field.dimensions() == 1); -bool ArtificialInput::sameAs(const MIRInput& other) const { - const auto* o = dynamic_cast(&other); - return (o != nullptr) && parametrisation_.matchAll(o->parametrisation_); -} + values = field.values(0); + }; + if (std::string path; parametrisation_.get("griddef", path)) { + GriddefInput::load(path, latitudes_, longitudes_); + } -bool ArtificialInput::has(const std::string& name) const { - return parametrisation_.has(name) || FieldParametrisation::has(name); -} + if (std::string lats, lons; parametrisation_.get("latitudes", lats) && parametrisation_.get("longitudes", lons)) { + load(lats, latitudes_); + load(lons, longitudes_); + } + ASSERT(latitudes_.size() == longitudes_.size()); + if (!parametrisation_.has("numberOfDataPoints")) { + parametrisation_.set("numberOfDataPoints", latitudes_.size()); + } -bool ArtificialInput::get(const std::string& name, std::string& value) const { - return parametrisation_.get(name, value) || FieldParametrisation::get(name, value); + if (!parametrisation_.has("paramId")) { + parametrisation_.set("paramId", 255L /*missing*/); + } } -bool ArtificialInput::get(const std::string& name, bool& value) const { - return parametrisation_.get(name, value) || FieldParametrisation::get(name, value); +void ArtificialInput::print(std::ostream& out) const { + out << "ArtificialInput[parametrisation=" << parametrisation_ << "]"; } -bool ArtificialInput::get(const std::string& name, int& value) const { - return parametrisation_.get(name, value) || FieldParametrisation::get(name, value); +bool ArtificialInput::sameAs(const MIRInput& other) const { + const auto* o = dynamic_cast(&other); + return (o != nullptr) && parametrisation_.matchAll(o->parametrisation_); } bool ArtificialInput::get(const std::string& name, long& value) const { - return parametrisation_.get(name, value) || FieldParametrisation::get(name, value); -} - - -bool ArtificialInput::get(const std::string& name, float& value) const { - return parametrisation_.get(name, value) || FieldParametrisation::get(name, value); + static const std::string PARAM_ID("paramId"); + return name == PARAM_ID ? parametrisation_.get(name, value) : FieldParametrisation::get(name, value); } -bool ArtificialInput::get(const std::string& name, double& value) const { - return parametrisation_.get(name, value) || FieldParametrisation::get(name, value); +void ArtificialInput::latitudes(std::vector& values) const { + values = latitudes_; } -bool ArtificialInput::get(const std::string& name, std::vector& value) const { - return parametrisation_.get(name, value) || FieldParametrisation::get(name, value); -} - - -bool ArtificialInput::get(const std::string& name, std::vector& value) const { - return parametrisation_.get(name, value) || FieldParametrisation::get(name, value); -} - - -bool ArtificialInput::get(const std::string& name, std::vector& value) const { - return parametrisation_.get(name, value) || FieldParametrisation::get(name, value); -} - - -bool ArtificialInput::get(const std::string& name, std::vector& value) const { - return parametrisation_.get(name, value) || FieldParametrisation::get(name, value); -} - - -bool ArtificialInput::get(const std::string& name, std::vector& value) const { - return parametrisation_.get(name, value) || FieldParametrisation::get(name, value); +void ArtificialInput::longitudes(std::vector& values) const { + values = longitudes_; } @@ -175,8 +173,6 @@ ArtificialInput* ArtificialInputFactory::build(const std::string& name, const pa util::call_once(once, init); util::lock_guard lock(*local_mutex); - // Log::debug() << "ArtificialInputFactory: looking for '" << name << "'" << std::endl; - auto j = m->find(name); if (j == m->end()) { list(Log::error() << "ArtificialInputFactory: unknown '" << name << "', choices are: "); diff --git a/src/mir/input/ArtificialInput.h b/src/mir/input/ArtificialInput.h index 16e404091..2960198bd 100644 --- a/src/mir/input/ArtificialInput.h +++ b/src/mir/input/ArtificialInput.h @@ -12,6 +12,8 @@ #pragma once +#include + #include "mir/input/MIRInput.h" #include "mir/param/FieldParametrisation.h" #include "mir/param/SimpleParametrisation.h" @@ -54,7 +56,7 @@ class ArtificialInput : public MIRInput, protected param::FieldParametrisation { protected: // -- Constructors - ArtificialInput() : calls_(0) {} + ArtificialInput(); // -- Members // None @@ -74,20 +76,7 @@ class ArtificialInput : public MIRInput, protected param::FieldParametrisation { bool sameAs(const MIRInput&) const override; // From FieldParametrisation - bool has(const std::string& name) const override; - - bool get(const std::string& name, std::string& value) const override; - bool get(const std::string& name, bool& value) const override; - bool get(const std::string& name, int& value) const override; bool get(const std::string& name, long& value) const override; - bool get(const std::string& name, float& value) const override; - bool get(const std::string& name, double& value) const override; - - bool get(const std::string& name, std::vector& value) const override; - bool get(const std::string& name, std::vector& value) const override; - bool get(const std::string& name, std::vector& value) const override; - bool get(const std::string& name, std::vector& value) const override; - bool get(const std::string& name, std::vector& value) const override; // -- Class members // None @@ -99,13 +88,21 @@ class ArtificialInput : public MIRInput, protected param::FieldParametrisation { // -- Members param::SimpleParametrisation parametrisation_; + std::unique_ptr inputParametrisation_; size_t calls_; + // For unstructured grids + std::vector latitudes_; + std::vector longitudes_; + // -- Methods // None // -- Overridden methods - // None + + // From FieldParametrisation + void latitudes(std::vector&) const override; + void longitudes(std::vector&) const override; // -- Class members // None diff --git a/src/mir/input/GribInput.cc b/src/mir/input/GribInput.cc index a238aa7a6..c90a6e6a5 100644 --- a/src/mir/input/GribInput.cc +++ b/src/mir/input/GribInput.cc @@ -26,7 +26,6 @@ #include "eckit/config/Resource.h" #include "eckit/io/Buffer.h" #include "eckit/io/MemoryHandle.h" -#include "eckit/io/StdFile.h" #include "eckit/serialisation/HandleStream.h" #include "eckit/types/FloatCompare.h" #include "eckit/types/Fraction.h" @@ -34,6 +33,8 @@ #include "mir/config/LibMir.h" #include "mir/data/MIRField.h" #include "mir/grib/Config.h" +#include "mir/input/GriddefInput.h" +#include "mir/param/DefaultParametrisation.h" #include "mir/repres/Representation.h" #include "mir/util/Exceptions.h" #include "mir/util/Grib.h" @@ -376,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; }); } @@ -392,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; }); } @@ -728,6 +733,11 @@ data::MIRField GribInput::field() const { GRIB_CALL(codes_get_long_array(grib_, "pl", pl.data(), &size)); ASSERT(count_pl == size); + if (auto pl_sum = static_cast(std::accumulate(pl.begin(), pl.end(), 0L)); pl_sum != values.size()) { + wrongly_encoded_grib("GribInput: sum of pl array (" + std::to_string(pl_sum) + + ") does not match the size of values array (" + std::to_string(values.size()) + ")"); + } + // NOTE: this fix ties with the method get(const std::string &name, std::vector &value) if (std::find(pl.rbegin(), pl.rend(), 0) != pl.rend()) { @@ -755,11 +765,11 @@ data::MIRField GribInput::field() const { for (auto p1 = pl.begin(), p2 = pl_fixed.begin(); p1 != pl.end(); ++p1, ++p2) { if (*p1 == 0) { ASSERT(*p2 > 0); - auto Ni = size_t(*p2); + auto Ni = static_cast(*p2); values_extended.insert(values_extended.end(), Ni, missingValue); } else { - auto Ni = size_t(*p1); + auto Ni = static_cast(*p1); ASSERT(i + Ni <= count); values_extended.insert(values_extended.end(), &values[i], &values[i + Ni]); @@ -772,7 +782,7 @@ data::MIRField GribInput::field() const { values.swap(values_extended); ASSERT(get("pl", pl)); - size_t pl_sum = size_t(std::accumulate(pl.begin(), pl.end(), 0L)); + auto pl_sum = static_cast(std::accumulate(pl.begin(), pl.end(), 0L)); ASSERT(pl_sum == values.size()); } } @@ -1141,52 +1151,32 @@ bool GribInput::handle(grib_handle* h) { } -void GribInput::auxilaryValues(const std::string& path, std::vector& values) const { +void GribInput::setAuxiliaryInformation(const util::ValueMap& map) { util::lock_guard lock(mutex_); - eckit::AutoStdFile f(path); - - int e; - grib_handle* h = nullptr; - - // We cannot use GribFileInput to read these files, because lat/lon files are also - // has grid_type = triangular_grid, and we will create a loop - - try { - h = codes_grib_handle_new_from_file(nullptr, f, &e); - grib_call(e, path.c_str()); - size_t count; - GRIB_CALL(codes_get_size(h, "values", &count)); - - size_t size = count; - values.resize(count); - GRIB_CALL(codes_get_double_array(h, "values", values.data(), &size)); - ASSERT(count == size); - - long missingValuesPresent; - GRIB_CALL(codes_get_long(h, "missingValuesPresent", &missingValuesPresent)); - ASSERT(!missingValuesPresent); + auto load = [](const eckit::PathName& path, std::vector& values) { + Log::info() << "GribInput::setAuxiliaryInformation: '" << path << "'" << std::endl; - codes_handle_delete(h); - } - catch (...) { - codes_handle_delete(h); - throw; - } -} + static const param::DefaultParametrisation param; + std::unique_ptr input(input::MIRInputFactory::build(path.asString(), param)); + ASSERT(input->next()); + auto field = input->field(); + ASSERT(field.dimensions() == 1); -void GribInput::setAuxiliaryInformation(const util::ValueMap& map) { - util::lock_guard lock(mutex_); + values = field.values(0); + }; for (const auto& kv : map) { - if (kv.first == "latitudes") { - Log::debug() << "Loading auxilary file '" << kv.second << "'" << std::endl; - auxilaryValues(kv.second, latitudes_); + if (kv.first == "griddef") { + GriddefInput::load(kv.second, latitudes_, longitudes_); + ASSERT(latitudes_.size() == longitudes_.size()); + } + else if (kv.first == "latitudes") { + load(kv.second, latitudes_); } else if (kv.first == "longitudes") { - Log::debug() << "Loading auxilary file '" << kv.second << "'" << std::endl; - auxilaryValues(kv.second, longitudes_); + load(kv.second, longitudes_); } } } diff --git a/src/mir/input/GribInput.h b/src/mir/input/GribInput.h index bcc19ffd5..709c1795f 100644 --- a/src/mir/input/GribInput.h +++ b/src/mir/input/GribInput.h @@ -46,10 +46,6 @@ class GribInput : public MIRInput, public param::FieldParametrisation { std::vector latitudes_; std::vector longitudes_; - // -- Methods - - void auxilaryValues(const std::string&, std::vector& values) const; - // -- Overridden methods // From MIRInput diff --git a/src/mir/input/GriddefInput.cc b/src/mir/input/GriddefInput.cc new file mode 100644 index 000000000..e5d764a0e --- /dev/null +++ b/src/mir/input/GriddefInput.cc @@ -0,0 +1,106 @@ +/* + * (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/input/GriddefInput.h" + +#include +#include +#include +#include "eckit/filesystem/PathName.h" +#include "eckit/serialisation/IfstreamStream.h" + +#include "mir/data/MIRField.h" +#include "mir/repres/other/UnstructuredGrid.h" +#include "mir/util/Exceptions.h" +#include "mir/util/Log.h" + + +namespace mir::input { + + +GriddefInput::GriddefInput(const eckit::PathName& path) : calls_(0) { + load(path, latitudes_, longitudes_); + ASSERT(latitudes_.size() == longitudes_.size()); + ASSERT(!latitudes_.empty()); +} + + +void GriddefInput::load(const eckit::PathName& path, std::vector& latitudes, std::vector& longitudes) { + Log::info() << "GriddefInput::load '" << path << "'" << std::endl; + + std::ifstream in(path.asString().c_str()); + if (!in) { + throw exception::CantOpenFile(path); + } + + if (std::isprint(in.peek()) == 0) { + // binary + eckit::IfstreamStream s(in); + + size_t version = 0; + s >> version; + ASSERT(version == 1); + + size_t count = 0; + s >> count; + + latitudes.resize(count); + longitudes.resize(count); + + for (size_t i = 0; i < count; ++i) { + s >> latitudes[i]; + s >> longitudes[i]; + } + } + else { + double lat = 0; + double lon = 0; + while (in >> lat >> lon) { + latitudes.push_back(lat); + longitudes.push_back(lon); + } + } +} + + +bool GriddefInput::next() { + return calls_++ == 0; +} + + +size_t GriddefInput::dimensions() const { + return 1; +} + + +const param::MIRParametrisation& GriddefInput::parametrisation(size_t which) const { + ASSERT(which == 0); + return empty_; +} + + +data::MIRField GriddefInput::field() const { + return {new repres::other::UnstructuredGrid{latitudes_, longitudes_}}; +} + + +bool GriddefInput::sameAs(const MIRInput&) const { + return false; +} + + +void GriddefInput::print(std::ostream& out) const { + out << "GriddefInput[size=" << latitudes_.size() << "]"; +} + + +} // namespace mir::input diff --git a/src/mir/input/GriddefInput.h b/src/mir/input/GriddefInput.h new file mode 100644 index 000000000..3c54253bf --- /dev/null +++ b/src/mir/input/GriddefInput.h @@ -0,0 +1,56 @@ +/* + * (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 + +#include "mir/input/MIRInput.h" +#include "mir/param/SimpleParametrisation.h" + + +namespace eckit { +class PathName; +} + + +namespace mir::input { + + +class GriddefInput final : public MIRInput { +public: + GriddefInput(const eckit::PathName&); + + static void load(const eckit::PathName&, std::vector& latitudes, std::vector& longitudes); + +private: + std::vector latitudes_; + std::vector longitudes_; + size_t calls_; + + const param::SimpleParametrisation empty_; + + bool next() override; + size_t dimensions() const override; + const param::MIRParametrisation& parametrisation(size_t which = 0) const override; + data::MIRField field() const override; + bool sameAs(const MIRInput&) const override; + void print(std::ostream&) const override; + + friend std::ostream& operator<<(std::ostream& s, const GriddefInput& p) { + p.print(s); + return s; + } +}; + + +} // namespace mir::input diff --git a/src/mir/input/NetcdfFileInput.h b/src/mir/input/NetcdfFileInput.h index 9db74f7a8..f445b5136 100644 --- a/src/mir/input/NetcdfFileInput.h +++ b/src/mir/input/NetcdfFileInput.h @@ -13,7 +13,6 @@ #pragma once #include "mir/api/mir_config.h" -#if mir_HAVE_NETCDF #include @@ -134,6 +133,3 @@ class NetcdfFileInput : public MIRInput, public param::FieldParametrisation, pub } // namespace mir::input - - -#endif diff --git a/src/mir/key/grid/TypedGrid.h b/src/mir/key/grid/TypedGrid.h index 98d1e9e66..f1b471061 100644 --- a/src/mir/key/grid/TypedGrid.h +++ b/src/mir/key/grid/TypedGrid.h @@ -29,7 +29,7 @@ class TypedGrid : public Grid { // -- Destructor - virtual ~TypedGrid() override; + ~TypedGrid() override; // -- Convertors // None diff --git a/src/mir/key/style/ECMWFStyle.cc b/src/mir/key/style/ECMWFStyle.cc index be9ae1b51..3b68b9fc1 100644 --- a/src/mir/key/style/ECMWFStyle.cc +++ b/src/mir/key/style/ECMWFStyle.cc @@ -24,7 +24,10 @@ #include "mir/key/grid/Grid.h" #include "mir/key/resol/Resol.h" #include "mir/output/MIROutput.h" +#include "mir/param/CombinedParametrisation.h" +#include "mir/param/DefaultParametrisation.h" #include "mir/param/MIRParametrisation.h" +#include "mir/param/RuntimeParametrisation.h" #include "mir/param/SameParametrisation.h" #include "mir/repres/latlon/LatLon.h" #include "mir/util/BoundingBox.h" @@ -270,6 +273,10 @@ void ECMWFStyle::sh2grid(action::ActionPlan& plan) const { bool vod2uv = option(user, "vod2uv", false); bool uv2uv = option(user, "uv2uv", false) || uv_input; // where "MIR knowledge of winds" is hardcoded + if (vod2uv && uv2uv) { + throw exception::UserError("ECMWFStyle: option 'vod2uv' is incompatible with option 'uv2uv'"); + } + if (vod2uv && uv_input) { throw exception::UserError("ECMWFStyle: option 'vod2uv' is incompatible with input U/V"); } @@ -281,13 +288,21 @@ void ECMWFStyle::sh2grid(action::ActionPlan& plan) const { auto target = target_gridded_from_parametrisation(parametrisation_, false); if (!target.empty()) { if (resol.resultIsSpectral()) { - plan.add("transform." + std::string(vod2uv ? "sh-vod-to-uv-" : "sh-scalar-to-") + target); + + if (uv2uv) { + plan.add("filter.adjust-winds-scale-cos-latitude"); + } } else { - resol.prepare(plan); + // apply u = U / cos(theta) first, as the associated supporting grid is typically Gaussian + // hence avoiding bad conditioning at the poles + if (uv2uv) { + plan.add("filter.adjust-winds-scale-cos-latitude"); + } + // if the intermediate grid is the same as the target grid, the interpolation to the // intermediate grid is not followed by an additional interpolation std::string grid; @@ -296,13 +311,11 @@ void ECMWFStyle::sh2grid(action::ActionPlan& plan) const { } } - if (vod2uv || uv2uv) { - ASSERT(vod2uv != uv2uv); - - if (uv2uv) { - plan.add("filter.adjust-winds-scale-cos-latitude"); - } + if (uv2uv) { + plan.add("filter.adjust-winds-at-poles"); + } + if (vod2uv || uv2uv) { if (rotation) { plan.add("filter.adjust-winds-directions"); } @@ -333,7 +346,8 @@ void ECMWFStyle::sh2sh(action::ActionPlan& plan) const { void ECMWFStyle::grid2grid(action::ActionPlan& plan) const { - const auto& user = parametrisation_.userParametrisation(); + const auto& user = parametrisation_.userParametrisation(); + const auto& field = parametrisation_.fieldParametrisation(); bool rotation = user.has("rotation"); bool vod2uv = option(user, "vod2uv", false); @@ -346,6 +360,21 @@ void ECMWFStyle::grid2grid(action::ActionPlan& plan) const { add_formula(plan, user, {"gridded", "raw"}); + if (std::string intint; user.get("intermediate-interpolation", intint) && !intint.empty()) { + if (std::string intgrid; user.get("intgrid", intgrid) && !intgrid.empty()) { + auto runtime = new param::RuntimeParametrisation{parametrisation_}; + runtime->set("interpolation", intint); + runtime->set("grid", intgrid); + runtime->unset("rotation"); + + static param::DefaultParametrisation defaults; + auto recombined = std::make_unique(*runtime, field, defaults); + auto target = target_gridded_from_parametrisation(*recombined, false); + + plan.add("interpolate.grid2" + target, runtime); + } + } + auto target = target_gridded_from_parametrisation(parametrisation_, rotation); if (!target.empty()) { plan.add("interpolate.grid2" + target); diff --git a/src/mir/method/Method.cc b/src/mir/method/Method.cc index 161d197e2..1a58ae528 100644 --- a/src/mir/method/Method.cc +++ b/src/mir/method/Method.cc @@ -22,23 +22,21 @@ #include "mir/util/Log.h" #include "mir/util/Mutex.h" -#if mir_HAVE_OMP extern "C" { void omp_set_num_threads(int); } -#endif namespace mir::method { Method::Method(const param::MIRParametrisation& params) : parametrisation_(params) { -#if mir_HAVE_OMP - int num_threads = 1; - if (params.get("parallel-omp-num-threads", num_threads)) { - omp_set_num_threads(num_threads); + if constexpr (HAVE_OMP) { + int num_threads = 1; + if (params.get("parallel-omp-num-threads", num_threads)) { + omp_set_num_threads(num_threads); + } } -#endif } diff --git a/src/mir/method/MethodWeighted.cc b/src/mir/method/MethodWeighted.cc index 443946e66..2f90c3761 100644 --- a/src/mir/method/MethodWeighted.cc +++ b/src/mir/method/MethodWeighted.cc @@ -19,6 +19,7 @@ #include #include +#include "eckit/config/Resource.h" #include "eckit/log/JSON.h" #include "eckit/types/FloatCompare.h" #include "eckit/utils/MD5.h" @@ -64,7 +65,7 @@ static caching::WeightCache& matrix_cache_disk(const param::MIRParametrisation& const static std::map KNOWN_METHOD{ {"73e1dd539879ffbbbb22d6bc789c2262", "linear"}, {"7738675c7e2c64d463718049ebef6563", "nearest-neighbour"}, - {"a81efab621096650c20a978062cdd169", "grid-box-average"}, + {"0346db910681bffd0d518b49923879dc", "grid-box-average"}, }; @@ -80,7 +81,6 @@ MethodWeighted::MethodWeighted(const param::MIRParametrisation& parametrisation) parametrisation_.get("pole-displacement-in-degree", poleDisplacement_); ASSERT(poleDisplacement_ >= 0); - matrixValidate_ = eckit::Resource("$MIR_MATRIX_VALIDATE", false); matrixAssemble_ = parametrisation_.userParametrisation().has("filter"); std::string nonLinear = "missing-if-heaviest-missing"; @@ -164,14 +164,15 @@ void MethodWeighted::createMatrix(context::Context& ctx, const repres::Represent const repres::Representation& out, WeightMatrix& W, const lsm::LandSeaMasks& masks, const Cropping& /*cropping*/) const { trace::ResourceUsage usage(std::string("MethodWeighted::createMatrix [") + name() + "]"); + const auto checks = validateMatrixWeights(); - computeMatrixWeights(ctx, in, out, W, validateMatrixWeights()); + // matrix validation always happens after creation, because the matrix can/will be cached + computeMatrixWeights(ctx, in, out, W); + W.validate("computeMatrixWeights", checks); if (masks.active() && masks.cacheable()) { applyMasks(W, masks); - if (matrixValidate_) { - W.validate("applyMasks"); - } + W.validate("applyMasks", checks); } } @@ -264,9 +265,7 @@ const WeightMatrix& MethodWeighted::getMatrix(context::Context& ctx, const repre // it will be cached in memory nevertheless if (masks.active() && !masks.cacheable()) { applyMasks(W, masks); - if (matrixValidate_) { - W.validate("applyMasks"); - } + W.validate("applyMasks", validateMatrixWeights()); } @@ -403,8 +402,8 @@ lsm::LandSeaMasks MethodWeighted::getMasks(const repres::Representation& in, con } -bool MethodWeighted::validateMatrixWeights() const { - return true; +WeightMatrix::Check MethodWeighted::validateMatrixWeights() const { + return {}; } @@ -493,9 +492,7 @@ void MethodWeighted::execute(context::Context& ctx, const repres::Representation trace::Timer t(str.str()); if (n->treatment(A, M, B, field.values(i), missingValue)) { - if (matrixValidate_) { - M.validate(str.str().c_str()); - } + M.validate(str.str().c_str(), validateMatrixWeights()); } } @@ -527,7 +524,7 @@ void MethodWeighted::execute(context::Context& ctx, const repres::Representation void MethodWeighted::computeMatrixWeights(context::Context& ctx, const repres::Representation& in, - const repres::Representation& out, WeightMatrix& W, bool validate) const { + const repres::Representation& out, WeightMatrix& W) const { auto timing(ctx.statistics().computeMatrixTimer()); if (in.sameAs(out) && !matrixAssemble_) { @@ -564,11 +561,6 @@ void MethodWeighted::computeMatrixWeights(context::Context& ctx, const repres::R W.swap(w); } } - - // matrix validation always happens after creation, because the matrix can/will be cached - if (validate) { - W.validate("computeMatrixWeights"); - } } diff --git a/src/mir/method/MethodWeighted.h b/src/mir/method/MethodWeighted.h index 8ff6e1556..0c81775fe 100644 --- a/src/mir/method/MethodWeighted.h +++ b/src/mir/method/MethodWeighted.h @@ -114,7 +114,6 @@ class MethodWeighted : public Method { std::unique_ptr reorderRows_; std::unique_ptr reorderCols_; - bool matrixValidate_; bool matrixAssemble_; // -- Methods @@ -127,10 +126,10 @@ class MethodWeighted : public Method { virtual void applyMasks(WeightMatrix&, const lsm::LandSeaMasks&) const; virtual lsm::LandSeaMasks getMasks(const repres::Representation& in, const repres::Representation& out) const; - virtual bool validateMatrixWeights() const; + virtual WeightMatrix::Check validateMatrixWeights() const; void computeMatrixWeights(context::Context&, const repres::Representation& in, const repres::Representation& out, - WeightMatrix&, bool validate) const; + WeightMatrix&) const; void createMatrix(context::Context&, const repres::Representation& in, const repres::Representation& out, WeightMatrix&, const lsm::LandSeaMasks&, const Cropping&) const; diff --git a/src/mir/method/ProxyMethod.h b/src/mir/method/ProxyMethod.h index e0d2399df..62cae2349 100644 --- a/src/mir/method/ProxyMethod.h +++ b/src/mir/method/ProxyMethod.h @@ -61,7 +61,7 @@ class ProxyMethod : public Method { // -- Destructor - virtual ~ProxyMethod() override; + ~ProxyMethod() override; // -- Members // None diff --git a/src/mir/method/WeightMatrix.cc b/src/mir/method/WeightMatrix.cc index 7dc933851..ac425ea44 100644 --- a/src/mir/method/WeightMatrix.cc +++ b/src/mir/method/WeightMatrix.cc @@ -13,7 +13,11 @@ #include "mir/method/WeightMatrix.h" #include +#include +#include +#include +#include "eckit/config/Resource.h" #include "eckit/types/FloatCompare.h" #include "mir/util/Exceptions.h" @@ -100,55 +104,81 @@ void WeightMatrix::cleanup(const double& pruneEpsilon) { } -void WeightMatrix::validate(const char* when) const { - constexpr size_t Nerrors = 50; - constexpr size_t Nvalues = 10; - - size_t errors = 0; +void WeightMatrix::validate(const char* when, Check check) const { + static bool matrixValidate = eckit::Resource("$MIR_MATRIX_VALIDATE", true); + if (!matrixValidate || (!check.duplicates && !check.bounds && !check.sum)) { + return; + } - for (Size r = 0; r < rows(); r++) { + constexpr size_t Nerrors = 10; + size_t errors = 0; - // check for W(i,j)<0, or W(i,j)>1, or sum(W(i,:))!=(0,1) - double sum = 0.; - bool ok = true; + std::ostringstream what; + const char* sep = ""; - for (const_iterator it = begin(r); it != end(r); ++it) { - double a = *it; - ok &= eckit::types::is_approximately_greater_or_equal(a, 0.) && - eckit::types::is_approximately_greater_or_equal(1., a); + // check for weights out of bounds 0 <= W(i,j) <= 1, sum(W(i,:))=(0,1), and dulplicate column indices + for (Size r = 0; r < rows(); r++) { + Scalar sum = 0.; + std::unordered_set cols; + + auto check_bounds = true; + auto check_duplicates = true; + for (auto it = begin(r); it != end(r); ++it) { + auto a = *it; + check_bounds &= eckit::types::is_approximately_greater_or_equal(a, 0.) && + eckit::types::is_approximately_greater_or_equal(1., a); sum += a; - } - ok &= (eckit::types::is_approximately_equal(sum, 0.) || eckit::types::is_approximately_equal(sum, 1.)); - if (ok) { - continue; + check_duplicates &= cols.insert(it.col()).second; } - // log issues, per row - if (Log::debug_active()) { + constexpr auto EPS = 1.e2 * std::numeric_limits::epsilon(); + auto check_sum = + eckit::types::is_approximately_equal(sum, 0., EPS) || eckit::types::is_approximately_equal(sum, 1., EPS); + + // ignore checks as required + check_duplicates |= !check.duplicates; + check_bounds |= !check.bounds; + check_sum |= !check.sum; + + if (!check_bounds || !check_sum || !check_duplicates) { if (errors < Nerrors) { - if (errors == 0) { - Log::debug() << "WeightMatrix::validate(" << when << ") failed " << std::endl; + what << sep << "row " << r << ": "; + const char* s = ""; + + if (!check_bounds) { + what << s << "weights out-of-bounds"; + s = ", "; } - Log::debug() << "Row: " << r; - size_t n = 0; - for (const_iterator it = begin(r); it != end(r); ++it, ++n) { - if (n > Nvalues) { - Log::debug() << " ..."; - break; - } - Log::debug() << " [" << *it << "]"; + if (!check_sum) { + what << s << "weights sum not 0 or 1 (sum=" << sum << ", 1-sum=" << (1 - sum) << ")"; + s = ", "; } - Log::debug() << " sum=" << sum << ", 1-sum " << (1 - sum) << std::endl; - } - else if (errors == Nerrors) { - Log::debug() << "..." << std::endl; + if (!check_duplicates) { + what << s << "duplicate indices"; + s = ", "; + } + + what << s << "contents: "; + s = ""; + for (auto it = begin(r); it != end(r); ++it) { + what << s << '(' << it.row() << ',' << it.col() << ',' << *it << ')'; + s = ", "; + } } + errors++; + sep = ", "; } } + + if (errors > 0) { + std::ostringstream errors_str; + errors_str << Log::Pretty{errors, {"row error"}}; + throw exception::InvalidWeightMatrix{when, errors_str.str() + ", " + what.str()}; + } } diff --git a/src/mir/method/WeightMatrix.h b/src/mir/method/WeightMatrix.h index 7695c879d..908ba0f8b 100644 --- a/src/mir/method/WeightMatrix.h +++ b/src/mir/method/WeightMatrix.h @@ -22,12 +22,18 @@ namespace mir::method { class WeightMatrix final : public eckit::linalg::SparseMatrix { -public: // types +public: using Triplet = eckit::linalg::Triplet; using Scalar = eckit::linalg::Scalar; using Size = eckit::linalg::Size; -public: // methods + struct Check { + bool duplicates = true; + bool bounds = true; + bool sum = true; + }; + +public: WeightMatrix(SparseMatrix::Allocator* = nullptr); WeightMatrix(const eckit::PathName&); @@ -38,9 +44,10 @@ class WeightMatrix final : public eckit::linalg::SparseMatrix { void cleanup(const double& pruneEpsilon = 0); - void validate(const char* when) const; + // Validate interpolation weights (default check matrix structure only) + void validate(const char* when, Check = {true, false, false}) const; -private: // members +private: void print(std::ostream&) const; friend std::ostream& operator<<(std::ostream& out, const WeightMatrix& m) { diff --git a/src/mir/method/fe/FiniteElement.cc b/src/mir/method/fe/FiniteElement.cc index 0afcff281..72b201c68 100644 --- a/src/mir/method/fe/FiniteElement.cc +++ b/src/mir/method/fe/FiniteElement.cc @@ -246,6 +246,11 @@ void FiniteElement::print(std::ostream& out) const { } +int FiniteElement::version() const { + return 1; +} + + bool FiniteElement::sameAs(const Method& other) const { const auto* o = dynamic_cast(&other); return (o != nullptr) && meshGeneratorParams_.sameAs(o->meshGeneratorParams_) && diff --git a/src/mir/method/fe/FiniteElement.h b/src/mir/method/fe/FiniteElement.h index 3eb0da542..a1c77d79e 100644 --- a/src/mir/method/fe/FiniteElement.h +++ b/src/mir/method/fe/FiniteElement.h @@ -105,6 +105,8 @@ class FiniteElement : public MethodWeighted { void json(eckit::JSON&) const override; void print(std::ostream&) const override; + int version() const override; + // -- Overridden methods // None diff --git a/src/mir/method/gridbox/GridBoxAverage.cc b/src/mir/method/gridbox/GridBoxAverage.cc index 72799a97c..9d698dbc6 100644 --- a/src/mir/method/gridbox/GridBoxAverage.cc +++ b/src/mir/method/gridbox/GridBoxAverage.cc @@ -49,12 +49,13 @@ GridBoxAverage::GridBoxAverage(const param::MIRParametrisation& parametrisation) void GridBoxAverage::assemble(util::MIRStatistics& /*unused*/, WeightMatrix& W, const repres::Representation& in, const repres::Representation& out) const { auto& log = Log::debug(); - log << "GridBoxMethod::assemble (input: " << in << ", output: " << out << ")" << std::endl; - + log << "GridBoxMethod::assemble" + "\ninput: " + << in << "\ndomain: " << in.domain() << "\noutput: " << out << "\ndomain: " << out.domain() << std::endl; if (!in.domain().contains(out.domain())) { std::ostringstream msg; - msg << "GridBoxMethod: input must contain output (input:" << in.domain() << ", output:" << out.domain() << ")"; + msg << "GridBoxMethod: input domain must contain output domain"; throw exception::UserError(msg.str()); } @@ -185,7 +186,7 @@ const char* GridBoxAverage::name() const { int GridBoxAverage::version() const { - return 3; + return 5; } diff --git a/src/mir/method/gridbox/GridBoxMethod.cc b/src/mir/method/gridbox/GridBoxMethod.cc index 4e917e3ba..938007914 100644 --- a/src/mir/method/gridbox/GridBoxMethod.cc +++ b/src/mir/method/gridbox/GridBoxMethod.cc @@ -45,8 +45,8 @@ void GridBoxMethod::print(std::ostream& out) const { } -bool GridBoxMethod::validateMatrixWeights() const { - return false; +WeightMatrix::Check GridBoxMethod::validateMatrixWeights() const { + return {true, true, false}; } diff --git a/src/mir/method/gridbox/GridBoxMethod.h b/src/mir/method/gridbox/GridBoxMethod.h index e36210ac2..5708a321c 100644 --- a/src/mir/method/gridbox/GridBoxMethod.h +++ b/src/mir/method/gridbox/GridBoxMethod.h @@ -27,7 +27,7 @@ class GridBoxMethod : public MethodWeighted { bool sameAs(const Method&) const override; void json(eckit::JSON&) const override; void print(std::ostream&) const override; - bool validateMatrixWeights() const override; + WeightMatrix::Check validateMatrixWeights() const override; }; diff --git a/src/mir/method/gridbox/GridBoxStatistics.cc b/src/mir/method/gridbox/GridBoxStatistics.cc index 6401a8119..3dba3ac19 100644 --- a/src/mir/method/gridbox/GridBoxStatistics.cc +++ b/src/mir/method/gridbox/GridBoxStatistics.cc @@ -181,7 +181,7 @@ const char* GridBoxStatistics::name() const { int GridBoxStatistics::version() const { - return 3; + return 5; } diff --git a/src/mir/method/knn/KNearestNeighbours.cc b/src/mir/method/knn/KNearestNeighbours.cc index 2b945b4f5..ca6f82ffa 100644 --- a/src/mir/method/knn/KNearestNeighbours.cc +++ b/src/mir/method/knn/KNearestNeighbours.cc @@ -162,4 +162,14 @@ void KNearestNeighbours::print(std::ostream& out) const { } +int KNearestNeighbours::version() const { + return 1; +} + + +WeightMatrix::Check KNearestNeighbours::validateMatrixWeights() const { + return distanceWeighting().validateMatrixWeights(); +} + + } // namespace mir::method::knn diff --git a/src/mir/method/knn/KNearestNeighbours.h b/src/mir/method/knn/KNearestNeighbours.h index d222bb565..38447102f 100644 --- a/src/mir/method/knn/KNearestNeighbours.h +++ b/src/mir/method/knn/KNearestNeighbours.h @@ -52,9 +52,12 @@ class KNearestNeighbours : public MethodWeighted { void print(std::ostream&) const override; const char* name() const override = 0; + int version() const override; virtual const pick::Pick& pick() const = 0; virtual const distance::DistanceWeighting& distanceWeighting() const = 0; + + WeightMatrix::Check validateMatrixWeights() const override; }; diff --git a/src/mir/method/knn/distance/DistanceWeighting.cc b/src/mir/method/knn/distance/DistanceWeighting.cc index 03f7e5eae..7dd3d2884 100644 --- a/src/mir/method/knn/distance/DistanceWeighting.cc +++ b/src/mir/method/knn/distance/DistanceWeighting.cc @@ -85,4 +85,9 @@ void DistanceWeightingFactory::list(std::ostream& out) { } +WeightMatrix::Check DistanceWeighting::validateMatrixWeights() const { + return {}; +} + + } // namespace mir::method::knn::distance diff --git a/src/mir/method/knn/distance/DistanceWeighting.h b/src/mir/method/knn/distance/DistanceWeighting.h index 5970a8c6b..2a3c85f88 100644 --- a/src/mir/method/knn/distance/DistanceWeighting.h +++ b/src/mir/method/knn/distance/DistanceWeighting.h @@ -47,6 +47,8 @@ class DistanceWeighting { virtual void hash(eckit::MD5&) const = 0; + virtual WeightMatrix::Check validateMatrixWeights() const; + private: virtual void json(eckit::JSON&) const = 0; virtual void print(std::ostream&) const = 0; diff --git a/src/mir/method/knn/distance/PseudoLaplace.cc b/src/mir/method/knn/distance/PseudoLaplace.cc index f0f577ccd..b0f59bc02 100644 --- a/src/mir/method/knn/distance/PseudoLaplace.cc +++ b/src/mir/method/knn/distance/PseudoLaplace.cc @@ -127,6 +127,12 @@ void PseudoLaplace::hash(eckit::MD5& h) const { } +WeightMatrix::Check PseudoLaplace::validateMatrixWeights() const { + // this method does not produce bounded interpolation weights + return {true, false, false}; +} + + static const DistanceWeightingBuilder __distance("pseudo-laplace"); diff --git a/src/mir/method/knn/distance/PseudoLaplace.h b/src/mir/method/knn/distance/PseudoLaplace.h index b8343cc7a..4c87b4860 100644 --- a/src/mir/method/knn/distance/PseudoLaplace.h +++ b/src/mir/method/knn/distance/PseudoLaplace.h @@ -28,6 +28,8 @@ struct PseudoLaplace : DistanceWeighting { void json(eckit::JSON&) const override; void print(std::ostream&) const override; void hash(eckit::MD5&) const override; + + WeightMatrix::Check validateMatrixWeights() const override; }; diff --git a/src/mir/method/voronoi/VoronoiMethod.cc b/src/mir/method/voronoi/VoronoiMethod.cc index 4672cf89f..965559425 100644 --- a/src/mir/method/voronoi/VoronoiMethod.cc +++ b/src/mir/method/voronoi/VoronoiMethod.cc @@ -168,8 +168,8 @@ void VoronoiMethod::print(std::ostream& out) const { } -bool VoronoiMethod::validateMatrixWeights() const { - return false; +WeightMatrix::Check VoronoiMethod::validateMatrixWeights() const { + return {true, true, false}; } @@ -178,4 +178,9 @@ const char* VoronoiMethod::name() const { } +int VoronoiMethod::version() const { + return 1; +} + + } // namespace mir::method::voronoi diff --git a/src/mir/method/voronoi/VoronoiMethod.h b/src/mir/method/voronoi/VoronoiMethod.h index 7a2dbd76c..bfaf22cf5 100644 --- a/src/mir/method/voronoi/VoronoiMethod.h +++ b/src/mir/method/voronoi/VoronoiMethod.h @@ -30,8 +30,9 @@ class VoronoiMethod : public MethodWeighted { bool sameAs(const Method&) const override; void json(eckit::JSON&) const override; void print(std::ostream&) const override; - bool validateMatrixWeights() const override; + WeightMatrix::Check validateMatrixWeights() const override; const char* name() const override; + int version() const override; knn::pick::NClosestOrNearest pick_; }; diff --git a/src/mir/netcdf/Codec.cc b/src/mir/netcdf/Codec.cc index 12a30b35b..a6f01fde1 100644 --- a/src/mir/netcdf/Codec.cc +++ b/src/mir/netcdf/Codec.cc @@ -146,12 +146,12 @@ CodecFactory::CodecFactory(const std::string& name) : name_(name) { util::call_once(once, init); util::lock_guard lock(*local_mutex); - if (m->find(name) != m->end()) { - throw exception::SeriousBug("CodecFactory: duplicate '" + name + "'"); + if (auto j = m->find(name); j == m->end()) { + (*m)[name] = this; + return; } - ASSERT(m->find(name) == m->end()); - (*m)[name] = this; + throw exception::SeriousBug("CodecFactory: duplicate '" + name + "'"); } @@ -168,13 +168,12 @@ Codec* CodecFactory::build(const std::string& name, const Variable& variable) { Log::debug() << "CodecFactory: looking for '" << name << "'" << std::endl; - auto j = m->find(name); - if (j == m->end()) { - list(Log::error() << "CodecFactory: unknown '" << name << "', choices are: "); - throw exception::SeriousBug("CodecFactory: unknown '" + name + "'"); + if (auto j = m->find(name); j != m->end()) { + return j->second->make(variable); } - return j->second->make(variable); + list(Log::error() << "CodecFactory: unknown '" << name << "', choices are: "); + throw exception::SeriousBug("CodecFactory: unknown '" + name + "'"); } @@ -184,7 +183,7 @@ void CodecFactory::list(std::ostream& out) { const char* sep = ""; for (const auto& j : *m) { - out << sep << j.first; + out << sep << "'" << j.first << "'"; sep = ", "; } } diff --git a/src/mir/netcdf/DataInputVariable.cc b/src/mir/netcdf/DataInputVariable.cc index 3abbf42a4..4a068a40c 100644 --- a/src/mir/netcdf/DataInputVariable.cc +++ b/src/mir/netcdf/DataInputVariable.cc @@ -129,8 +129,8 @@ void DataInputVariable::get2DValues(MIRValuesVector& values, size_t index) const dims.pop_back(); std::vector coords(dims.size()); - for (int i = int(dims.size()) - 1; i >= 0; i--) { - auto j = size_t(i); + for (int i = static_cast(dims.size()) - 1; i >= 0; i--) { + auto j = static_cast(i); coords[j] = (index % dims[j]); index /= dims[j]; } @@ -139,11 +139,8 @@ void DataInputVariable::get2DValues(MIRValuesVector& values, size_t index) const start[j] = coords[j]; } - std::vector data(nx * ny); - matrix_->read(data, start, count); - - // FIXME: this copies, so is slow, but we could 'move' data instead - values.assign(data.begin(), data.end()); + values.resize(nx * ny); + matrix_->read(values, start, count); } diff --git a/src/mir/netcdf/Type.cc b/src/mir/netcdf/Type.cc index 5006232da..fa3580191 100644 --- a/src/mir/netcdf/Type.cc +++ b/src/mir/netcdf/Type.cc @@ -17,7 +17,6 @@ #include "mir/netcdf/Type.h" #include -#include #include #include @@ -416,7 +415,7 @@ Value* TypeT::attributeValue(int nc, int id, const char* name, size std::string value(len + 1, '\0'); NC_CALL(nc_get_att_text(nc, id, name, value.data()), path); - auto* v = new ValueT(*this, value); + auto* v = new ValueT(*this, value.data()); return v; } diff --git a/src/mir/output/GeoPointsFileOutputXYV.h b/src/mir/output/GeoPointsFileOutputXYV.h index a09f30be7..769167c9e 100644 --- a/src/mir/output/GeoPointsFileOutputXYV.h +++ b/src/mir/output/GeoPointsFileOutputXYV.h @@ -58,7 +58,7 @@ class GeoPointsFileOutputXYV : public GeoPointsFileOutput { // -- Overridden methods // From MIROutput - virtual size_t save(const param::MIRParametrisation&, context::Context&) override; + size_t save(const param::MIRParametrisation&, context::Context&) override; // -- Class members // None diff --git a/src/mir/output/GeoPointsFileOutputXYVector.h b/src/mir/output/GeoPointsFileOutputXYVector.h index 59e6f52b8..9f6fa9c4c 100644 --- a/src/mir/output/GeoPointsFileOutputXYVector.h +++ b/src/mir/output/GeoPointsFileOutputXYVector.h @@ -58,7 +58,7 @@ class GeoPointsFileOutputXYVector : public GeoPointsFileOutput { // -- Overridden methods // From MIROutput - virtual size_t save(const param::MIRParametrisation&, context::Context&) override; + size_t save(const param::MIRParametrisation&, context::Context&) override; // -- Class members // None diff --git a/src/mir/output/GribOutput.cc b/src/mir/output/GribOutput.cc index e9b15059b..b5946f72d 100644 --- a/src/mir/output/GribOutput.cc +++ b/src/mir/output/GribOutput.cc @@ -97,27 +97,21 @@ size_t GribOutput::copy(const param::MIRParametrisation& /*unused*/, context::Co void GribOutput::estimate(const param::MIRParametrisation& param, api::MIREstimation& estimator, context::Context& ctx) const { - - const data::MIRField& field = ctx.field(); + const auto& field = ctx.field(); ASSERT(field.dimensions() == 1); field.representation()->estimate(estimator); - long bits = 0; - if (param.get("accuracy", bits)) { - estimator.accuracy(size_t(bits)); + if (long accuracy = 0; param.get("accuracy", accuracy)) { + estimator.accuracy(static_cast(accuracy)); } - std::string packing; - if (param.get(packing, packing)) { + if (std::string packing; param.get(packing, packing)) { estimator.packing(packing); - // const grib::Packing &packer = grib::Packing::lookup(packing); - // packer.estimate(estimator, *field.representation()); } - long edition; - if (param.get("edition", edition)) { - estimator.edition(size_t(edition)); + if (long edition; param.get("edition", edition)) { + estimator.edition(static_cast(edition)); } } diff --git a/src/mir/output/GriddefOutput.cc b/src/mir/output/GriddefOutput.cc new file mode 100644 index 000000000..55b0dccb5 --- /dev/null +++ b/src/mir/output/GriddefOutput.cc @@ -0,0 +1,119 @@ +/* + * (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/output/GriddefOutput.h" + +#include +#include +#include + +#include "eckit/filesystem/PathName.h" +#include "eckit/serialisation/FileStream.h" + +#include "mir/action/context/Context.h" +#include "mir/data/MIRField.h" +#include "mir/repres/Iterator.h" +#include "mir/repres/Representation.h" + + +namespace mir::output { + + +GriddefOutput::GriddefOutput(std::string path, bool binary) : path_(std::move(path)), binary_(binary) {} + + +void GriddefOutput::save(const eckit::PathName& path, const std::vector& latitudes, + const std::vector& longitudes, bool binary) { + ASSERT(latitudes.size() == longitudes.size()); + const auto count = latitudes.size(); + + if (binary) { + eckit::FileStream s(path, "w"); + s << static_cast(1); // version + s << static_cast(count); + + for (size_t i = 0; i < count; ++i) { + s << latitudes[i]; + s << longitudes[i]; + } + + s.close(); + } + else { + std::ofstream out(path.asString().c_str()); + std::fixed(out).precision(6); // dot + 6 digits + + for (size_t i = 0; i < count; ++i) { + out << latitudes[i] << ' ' << longitudes[i] << '\n'; + } + } +} + + +size_t GriddefOutput::save(const param::MIRParametrisation& param, context::Context& ctx) { + const auto& field = ctx.field(); + repres::RepresentationHandle repres{field.representation()}; + + if (binary_) { + eckit::FileStream s(path_, "w"); + s << static_cast(1); // version + s << static_cast(repres->numberOfPoints()); + + for (std::unique_ptr it{repres->iterator()}; it->next();) { + const auto& p = **it; + s << p[0]; // latitude + s << p[1]; // longitude + } + + s.close(); + } + else { + std::ofstream out(path_.c_str()); + std::fixed(out).precision(6); // dot + 6 digits + + for (std::unique_ptr it{repres->iterator()}; it->next();) { + const auto& p = **it; + out << p[0] << ' ' << p[1] << '\n'; + } + } + + return 1; +} + + +bool GriddefOutput::sameAs(const MIROutput& other) const { + const auto* o = dynamic_cast(&other); + return (o != nullptr) && eckit::PathName(path_) == eckit::PathName(o->path_); +} + + +bool GriddefOutput::sameParametrisation(const param::MIRParametrisation& /*unused*/, + const param::MIRParametrisation& /*unused*/) const { + return false; +} + + +bool GriddefOutput::printParametrisation(std::ostream& /*unused*/, const param::MIRParametrisation& /*unused*/) const { + return false; +} + + +void GriddefOutput::print(std::ostream& out) const { + out << "GriddefOutput[path=" << path_ << ",binary=" << binary_ << "]"; +} + + +static const MIROutputBuilder _output1("griddef"); +static const MIROutputBuilder _output2("griddef-text"); + + +} // namespace mir::output diff --git a/src/mir/output/GriddefOutput.h b/src/mir/output/GriddefOutput.h new file mode 100644 index 000000000..78f24b7c4 --- /dev/null +++ b/src/mir/output/GriddefOutput.h @@ -0,0 +1,58 @@ +/* + * (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/output/MIROutput.h" + + +namespace eckit { +class PathName; +} + + +namespace mir::output { + + +class GriddefOutput : public MIROutput { +public: + static void save(const eckit::PathName&, const std::vector& latitudes, + const std::vector& longitudes, bool binary = true); + +protected: + GriddefOutput(std::string path, bool binary); + +private: + const std::string path_; + const bool binary_; + + size_t save(const param::MIRParametrisation&, context::Context&) override; + bool sameAs(const MIROutput&) const override; + bool sameParametrisation(const param::MIRParametrisation&, const param::MIRParametrisation&) const override; + bool printParametrisation(std::ostream&, const param::MIRParametrisation&) const override; + void print(std::ostream&) const override; +}; + + +class GriddefTextOutput : public GriddefOutput { +public: + GriddefTextOutput(std::string path) : GriddefOutput(path, false) {} +}; + + +class GriddefBinaryOutput : public GriddefOutput { +public: + GriddefBinaryOutput(std::string path) : GriddefOutput(path, true) {} +}; + + +} // namespace mir::output diff --git a/src/mir/repres/Representation.h b/src/mir/repres/Representation.h index f768b1096..6a1c4f199 100644 --- a/src/mir/repres/Representation.h +++ b/src/mir/repres/Representation.h @@ -14,6 +14,7 @@ #include #include +#include #include "eckit/memory/Counted.h" diff --git a/src/mir/repres/gauss/reduced/Classic.h b/src/mir/repres/gauss/reduced/Classic.h index 56c8a0346..721f14f26 100644 --- a/src/mir/repres/gauss/reduced/Classic.h +++ b/src/mir/repres/gauss/reduced/Classic.h @@ -30,7 +30,7 @@ class Classic : public Reduced { // -- Destructor - virtual ~Classic() override = default; + ~Classic() override = default; // -- Convertors // None @@ -58,8 +58,8 @@ class Classic : public Reduced { // None // -- Overridden methods - virtual void fillGrib(grib_info&) const override; - virtual void fillJob(api::MIRJob&) const override; + void fillGrib(grib_info&) const override; + void fillJob(api::MIRJob&) const override; atlas::Grid atlasGrid() const override; void makeName(std::ostream&) const override; bool sameAs(const Representation&) const override; diff --git a/src/mir/repres/gauss/reduced/FromPL.h b/src/mir/repres/gauss/reduced/FromPL.h index 84d4add1d..a462061b9 100644 --- a/src/mir/repres/gauss/reduced/FromPL.h +++ b/src/mir/repres/gauss/reduced/FromPL.h @@ -32,7 +32,7 @@ class FromPL : public Reduced { // -- Destructor - virtual ~FromPL() override = default; + ~FromPL() override = default; // -- Convertors // None @@ -60,8 +60,8 @@ class FromPL : public Reduced { // None // -- Overridden methods - virtual void fillGrib(grib_info&) const override; - virtual void fillJob(api::MIRJob&) const override; + void fillGrib(grib_info&) const override; + void fillJob(api::MIRJob&) const override; atlas::Grid atlasGrid() const override; void makeName(std::ostream&) const override; bool sameAs(const Representation&) const override; diff --git a/src/mir/repres/gauss/reduced/Octahedral.h b/src/mir/repres/gauss/reduced/Octahedral.h index c5a68c98f..7e473f8d5 100644 --- a/src/mir/repres/gauss/reduced/Octahedral.h +++ b/src/mir/repres/gauss/reduced/Octahedral.h @@ -30,7 +30,7 @@ class Octahedral : public Reduced { // -- Destructor - virtual ~Octahedral() override = default; + ~Octahedral() override = default; // -- Convertors // None @@ -58,8 +58,8 @@ class Octahedral : public Reduced { // None // -- Overridden methods - virtual void fillGrib(grib_info&) const override; - virtual void fillJob(api::MIRJob&) const override; + void fillGrib(grib_info&) const override; + void fillJob(api::MIRJob&) const override; atlas::Grid atlasGrid() const override; void makeName(std::ostream&) const override; bool sameAs(const Representation&) const override; diff --git a/src/mir/repres/gauss/reduced/Reduced.cc b/src/mir/repres/gauss/reduced/Reduced.cc index 16ac8f0fd..51c743bdc 100644 --- a/src/mir/repres/gauss/reduced/Reduced.cc +++ b/src/mir/repres/gauss/reduced/Reduced.cc @@ -29,6 +29,7 @@ #include "mir/repres/gauss/GaussianIterator.h" #include "mir/util/Atlas.h" #include "mir/util/BoundingBox.h" +#include "mir/util/Domain.h" #include "mir/util/Exceptions.h" #include "mir/util/Grib.h" #include "mir/util/GridBox.h" @@ -294,12 +295,19 @@ std::vector Reduced::gridBoxes() const { // latitude edges std::vector latEdges = calculateUnrotatedGridBoxLatitudeEdges(); + const auto dom = domain(); + const auto north = dom.north().value(); + const auto west = dom.west(); + const auto south = dom.south().value(); + const auto east = dom.east(); + + const auto periodic = isPeriodicWestEast(); + // grid boxes std::vector r; r.reserve(numberOfPoints()); - bool periodic = isPeriodicWestEast(); const auto& pl = pls(); for (size_t j = k_; j < k_ + Nj_; ++j) { @@ -307,20 +315,20 @@ std::vector Reduced::gridBoxes() const { const auto inc = eckit::Fraction(360, pl[j]); // latitude edges - const auto n = includesNorthPole() ? latEdges[j] : std::min(bbox_.north().value(), latEdges[j]); - const auto s = includesSouthPole() ? latEdges[j + 1] : std::max(bbox_.south().value(), latEdges[j + 1]); + const auto n = includesNorthPole() ? latEdges[j] : std::min(north, latEdges[j]); + const auto s = includesSouthPole() ? latEdges[j + 1] : std::max(south, latEdges[j + 1]); ASSERT(n >= s); // longitude edges - const auto west = bbox_.west().fraction(); - auto Nw = (west / inc).integralPart(); - if (Nw * inc < west) { + const auto w = west.fraction(); + auto Nw = (w / inc).integralPart(); + if (Nw * inc < w) { Nw += 1; } - const auto east = bbox_.east().fraction(); - auto Ne = (east / inc).integralPart(); - if (Ne * inc > east) { + const auto e = east.fraction(); + auto Ne = (e / inc).integralPart(); + if (Ne * inc > e) { Ne -= 1; } @@ -340,9 +348,9 @@ std::vector Reduced::gridBoxes() const { } else { for (size_t i = 0; i < N; ++i) { - auto w = std::max(bbox_.west().value(), lon1.value()); + auto w = std::max(west.value(), lon1.value()); lon1 += inc; - r.emplace_back(n, w, s, std::min(bbox_.east().value(), lon1.value())); + r.emplace_back(n, w, s, std::min(east.value(), lon1.value())); } ASSERT(lon0 <= lon1.normalise(lon0)); diff --git a/src/mir/repres/gauss/reduced/Reduced.h b/src/mir/repres/gauss/reduced/Reduced.h index 6e6e2ace2..9020f48f8 100644 --- a/src/mir/repres/gauss/reduced/Reduced.h +++ b/src/mir/repres/gauss/reduced/Reduced.h @@ -79,8 +79,8 @@ class Reduced : public Gaussian { // -- Overridden methods - virtual void fillGrib(grib_info&) const override; - virtual void fillJob(api::MIRJob&) const override; + void fillGrib(grib_info&) const override; + void fillJob(api::MIRJob&) const override; bool sameAs(const Representation&) const override; util::BoundingBox extendBoundingBox(const util::BoundingBox&) const override; diff --git a/src/mir/repres/gauss/regular/Regular.h b/src/mir/repres/gauss/regular/Regular.h index a570f2fb5..145fb0149 100644 --- a/src/mir/repres/gauss/regular/Regular.h +++ b/src/mir/repres/gauss/regular/Regular.h @@ -30,7 +30,7 @@ class Regular : public Gaussian { // -- Destructor - virtual ~Regular() override = default; + ~Regular() override = default; // -- Convertors // None @@ -67,8 +67,8 @@ class Regular : public Gaussian { // -- Overridden methods - virtual void fillGrib(grib_info&) const override; - virtual void fillJob(api::MIRJob&) const override; + void fillGrib(grib_info&) const override; + void fillJob(api::MIRJob&) const override; bool sameAs(const Representation&) const override; atlas::Grid atlasGrid() const override; void makeName(std::ostream&) const override; diff --git a/src/mir/repres/gauss/regular/RegularGG.cc b/src/mir/repres/gauss/regular/RegularGG.cc index e49d347a6..f172907cb 100644 --- a/src/mir/repres/gauss/regular/RegularGG.cc +++ b/src/mir/repres/gauss/regular/RegularGG.cc @@ -17,6 +17,7 @@ #include #include "mir/repres/gauss/GaussianIterator.h" +#include "mir/util/Domain.h" #include "mir/util/Exceptions.h" #include "mir/util/GridBox.h" @@ -68,22 +69,29 @@ std::vector RegularGG::gridBoxes() const { // latitude edges std::vector latEdges = calculateUnrotatedGridBoxLatitudeEdges(); + const auto dom = domain(); + const auto north = dom.north().value(); + const auto west = dom.west(); + const auto south = dom.south().value(); + const auto east = dom.east(); + + const auto periodic = isPeriodicWestEast(); + // longitude edges - bool periodic = isPeriodicWestEast(); - auto lon0 = bbox_.west(); - auto inc = (bbox_.east() - bbox_.west()).fraction() / (Ni_ - 1); + auto lon0 = west; + auto we = (east - west).fraction() / (periodic ? Ni_ : Ni_ - 1); eckit::Fraction half(1, 2); - std::vector lonEdges(Ni_ + 1, 0.); - lonEdges[0] = (lon0 - inc / 2).value(); + std::vector lonEdges(Ni_ + 1); + lonEdges[0] = (lon0 - we / 2).value(); for (size_t i = 0; i < Ni_; ++i) { - lonEdges[i + 1] = (lon0 + (i + half) * inc).value(); + lonEdges[i + 1] = (lon0 + (i + half) * we).value(); } if (!periodic) { - lonEdges.front() = std::max(bbox_.west().value(), lonEdges.front()); - lonEdges.back() = std::min(bbox_.east().value(), lonEdges.back()); + lonEdges.front() = std::max(west.value(), lonEdges.front()); + lonEdges.back() = std::min(east.value(), lonEdges.back()); } @@ -96,9 +104,9 @@ std::vector RegularGG::gridBoxes() const { for (size_t i = 0; i < Ni_; ++i) { auto l = lon1; - lon1 = l + Longitude(inc * (i + half)); - r.emplace_back(std::min(bbox_.north().value(), latEdges[j]), lonEdges[i], - std::max(bbox_.south().value(), latEdges[j + 1]), lonEdges[i + 1]); + lon1 = l + Longitude(we * (i + half)); + r.emplace_back(std::min(north, latEdges[j]), lonEdges[i], std::max(south, latEdges[j + 1]), + lonEdges[i + 1]); } ASSERT(periodic ? lon0 == lon1.normalise(lon0) : lon0 <= lon1.normalise(lon0)); diff --git a/src/mir/repres/latlon/ReducedLL.cc b/src/mir/repres/latlon/ReducedLL.cc index 8b5d0bd82..b39b06271 100644 --- a/src/mir/repres/latlon/ReducedLL.cc +++ b/src/mir/repres/latlon/ReducedLL.cc @@ -22,6 +22,7 @@ #include "eckit/utils/MD5.h" #include "mir/api/MIRJob.h" +#include "mir/api/mir_config.h" #include "mir/param/MIRParametrisation.h" #include "mir/repres/Iterator.h" #include "mir/util/Domain.h" @@ -272,11 +273,16 @@ Iterator* ReducedLL::iterator() const { } std::vector ReducedLL::gridBoxes() const { - auto dom = domain(); - bool periodic = isPeriodicWestEast(); + const auto dom = domain(); + const auto north = dom.north(); + const auto west = dom.west(); + const auto south = dom.south(); + const auto east = dom.east(); + + const auto periodic = isPeriodicWestEast(); auto Nj = pl_.size(); - const eckit::Fraction sn = (dom.north() - dom.south()).fraction() / eckit::Fraction(Nj - 1); + const eckit::Fraction sn = (north - south).fraction() / eckit::Fraction(Nj - 1); eckit::Fraction half(1, 2); @@ -294,8 +300,15 @@ std::vector ReducedLL::gridBoxes() const { latEdges[j + 1] = (lat0 - (j + half) * sn).value(); } - latEdges.front() = std::min(dom.north().value(), latEdges.front()); - latEdges.back() = std::max(dom.south().value(), latEdges.back()); + latEdges.front() = std::min(north.value(), std::max(south.value(), latEdges.front())); + if (dom.includesPoleNorth()) { + latEdges.front() = Latitude::NORTH_POLE.value(); + } + + latEdges.back() = std::min(north.value(), std::max(south.value(), latEdges.back())); + if (dom.includesPoleSouth()) { + latEdges.back() = Latitude::SOUTH_POLE.value(); + } for (size_t j = 0; j < Nj; ++j) { @@ -303,9 +316,9 @@ std::vector ReducedLL::gridBoxes() const { // longitude edges auto Ni = pl_[j]; ASSERT(Ni > 1); - const eckit::Fraction we = (dom.east() - dom.west()).fraction() / (Ni - (periodic ? 0 : 1)); + const eckit::Fraction we = (east - west).fraction() / (periodic ? Ni : Ni - 1); - auto lon0 = bbox_.west() - we / 2; + auto lon0 = west - we / 2; auto lon1 = lon0; if (periodic) { @@ -319,9 +332,9 @@ std::vector ReducedLL::gridBoxes() const { } else { for (long i = 0; i < Ni; ++i) { - auto w = std::max(bbox_.west().value(), lon1.value()); + auto w = std::max(west.value(), lon1.value()); lon1 += we; - auto e = std::min(bbox_.east().value(), lon1.value()); + auto e = std::min(east.value(), lon1.value()); r.emplace_back(latEdges[j], w, latEdges[j + 1], e); } diff --git a/src/mir/repres/latlon/RegularLL.cc b/src/mir/repres/latlon/RegularLL.cc index b2aac5a7d..7796505b0 100644 --- a/src/mir/repres/latlon/RegularLL.cc +++ b/src/mir/repres/latlon/RegularLL.cc @@ -145,10 +145,13 @@ util::BoundingBox RegularLL::extendBoundingBox(const util::BoundingBox& bbox) co } std::vector RegularLL::gridBoxes() const { + const auto dom = domain(); + const auto north = dom.north(); + const auto west = dom.west(); + const auto south = dom.south(); + const auto east = dom.east(); - auto dom = domain(); - auto north = dom.north().value(); - auto south = dom.south().value(); + const auto periodic = isPeriodicWestEast(); auto lat0 = bbox_.north(); auto lon0 = bbox_.west(); @@ -166,8 +169,16 @@ std::vector RegularLL::gridBoxes() const { latEdges[j + 1] = (lat0 - (j + half) * sn.fraction()).value(); } - latEdges.front() = std::min(north, std::max(south, latEdges.front())); - latEdges.back() = std::min(north, std::max(south, latEdges.back())); + latEdges.front() = std::min(north.value(), std::max(south.value(), latEdges.front())); + if (dom.includesPoleNorth()) { + latEdges.front() = Latitude::NORTH_POLE.value(); + } + + latEdges.back() = std::min(north.value(), std::max(south.value(), latEdges.back())); + latEdges.back() = std::min(north.value(), std::max(south.value(), latEdges.back())); + if (dom.includesPoleSouth()) { + latEdges.back() = Latitude::SOUTH_POLE.value(); + } // longitude edges @@ -177,10 +188,9 @@ std::vector RegularLL::gridBoxes() const { lonEdges[i + 1] = (lon0 + (i + half) * we).value(); } - bool periodic = isPeriodicWestEast(); if (!periodic) { - lonEdges.front() = std::max(lonEdges.front(), lon0.value()); - lonEdges.back() = std::min(lonEdges.back(), bbox_.east().value()); + lonEdges.front() = std::max(west.value(), lonEdges.front()); + lonEdges.back() = std::min(east.value(), lonEdges.back()); } diff --git a/src/mir/repres/other/UnstructuredGrid.cc b/src/mir/repres/other/UnstructuredGrid.cc index ac5d064ea..9384c81f0 100644 --- a/src/mir/repres/other/UnstructuredGrid.cc +++ b/src/mir/repres/other/UnstructuredGrid.cc @@ -12,21 +12,19 @@ #include "mir/repres/other/UnstructuredGrid.h" -#include #include #include #include #include #include -#include "eckit/filesystem/PathName.h" -#include "eckit/serialisation/FileStream.h" -#include "eckit/serialisation/IfstreamStream.h" #include "eckit/utils/MD5.h" #include "mir/api/MIRJob.h" #include "mir/api/mir_config.h" +#include "mir/input/GriddefInput.h" #include "mir/iterator/UnstructuredIterator.h" +#include "mir/output/GriddefOutput.h" #include "mir/param/MIRParametrisation.h" #include "mir/util/CheckDuplicatePoints.h" #include "mir/util/Domain.h" @@ -68,7 +66,7 @@ UnstructuredGrid::UnstructuredGrid(const param::MIRParametrisation& parametrisat parametrisation.get("longitudes", longitudes_); if (latitudes_.size() != longitudes_.size()) { - throw exception::UserError("UnstructuredGrid: requires 'latitudes' and 'longitudes' with the same size"); + throw exception::UserError("UnstructuredGrid: requires 'latitudes'/'longitudes' with the same size"); } if (latitudes_.empty()) { @@ -76,7 +74,7 @@ UnstructuredGrid::UnstructuredGrid(const param::MIRParametrisation& parametrisat parametrisation.get("numberOfPoints", numberOfPoints); if (numberOfPoints == 0) { - throw exception::UserError("UnstructuredGrid: requires 'latitudes' and 'longitudes', or 'numberOfPoints'"); + throw exception::UserError("UnstructuredGrid: requires 'latitudes'/'longitudes' or 'numberOfPoints'"); } // coordinates are unusable but unique @@ -90,68 +88,16 @@ UnstructuredGrid::UnstructuredGrid(const param::MIRParametrisation& parametrisat UnstructuredGrid::UnstructuredGrid(const eckit::PathName& path) { - std::ifstream in(path.asString().c_str()); - if (!in) { - throw exception::CantOpenFile(path); - } - - if (::isprint(in.peek()) == 0) { - - Log::info() << "UnstructuredGrid::load " << path << std::endl; - - eckit::IfstreamStream s(in); - size_t version = 0; - s >> version; - ASSERT(version == 1); - - size_t count = 0; - s >> count; - - latitudes_.resize(count); - longitudes_.resize(count); - - for (size_t i = 0; i < count; ++i) { - s >> latitudes_[i]; - s >> longitudes_[i]; - // Log::info() << latitudes_[i] << " " << longitudes_[i] << std::endl; - } - } - else { - double lat = 0; - double lon = 0; - while (in >> lat >> lon) { - latitudes_.push_back(lat); - longitudes_.push_back(lon); - } - } - + input::GriddefInput::load(path, latitudes_, longitudes_); util::check_duplicate_points("UnstructuredGrid from " + path.asString(), latitudes_, longitudes_); } void UnstructuredGrid::save(const eckit::PathName& path, const std::vector& latitudes, const std::vector& longitudes, bool binary) { - Log::info() << "UnstructuredGrid::save " << path << std::endl; - util::check_duplicate_points("UnstructuredGrid::save to " + path.asString(), latitudes, longitudes); - - if (binary) { - eckit::FileStream s(path, "w"); - size_t version = 1; - size_t count = latitudes.size(); - s << version; - s << count; - for (size_t i = 0; i < count; ++i) { - s << latitudes[i]; - s << longitudes[i]; - - Log::info() << latitudes[i] << " " << longitudes[i] << std::endl; - } - s.close(); - } - else { - NOTIMP; - } + ASSERT(binary); + output::GriddefOutput::save(path, latitudes, longitudes); } @@ -159,7 +105,7 @@ UnstructuredGrid::UnstructuredGrid(const std::vector& latitudes, const s const util::BoundingBox& bbox) : Gridded(bbox), latitudes_(latitudes), longitudes_(longitudes) { if (latitudes_.size() != longitudes_.size()) { - throw exception::UserError("UnstructuredGrid: requires 'latitudes' and 'longitudes' with the same size"); + throw exception::UserError("UnstructuredGrid: requires 'latitudes'/'longitudes' with the same size"); } util::check_duplicate_points("UnstructuredGrid from arguments", latitudes_, longitudes_); diff --git a/src/mir/repres/regular/PolarStereographic.cc b/src/mir/repres/regular/PolarStereographic.cc index 167d99e00..4caf588c6 100644 --- a/src/mir/repres/regular/PolarStereographic.cc +++ b/src/mir/repres/regular/PolarStereographic.cc @@ -12,7 +12,14 @@ #include "mir/repres/regular/PolarStereographic.h" +#include +#include + +#include "mir/api/MIRJob.h" +#include "mir/param/MIRParametrisation.h" +#include "mir/util/Angles.h" #include "mir/util/Exceptions.h" +#include "mir/util/Grib.h" namespace mir::repres::regular { @@ -22,11 +29,84 @@ static const RepresentationBuilder __builder("polar_stereogr PolarStereographic::PolarStereographic(const param::MIRParametrisation& param) : - RegularGrid(param, make_proj_spec(param)) {} + RegularGrid(param, make_proj_spec(param)), + LaDInDegrees_(0), + orientationOfTheGridInDegrees_(0), + southPoleOnProjectionPlane_(false) { + long edition = 0; + param.get("edition", edition); + + // GRIB1 cannot write LaD + writeLaDInDegrees_ = edition == 2; + param.get("writeLaDInDegrees", writeLaDInDegrees_); + + // GRIB2 cannot write negative longitude values + writeLonPositive_ = edition == 2; + param.get("writeLonPositive", writeLonPositive_); + + // Details + ASSERT(param.get("proj", proj_) && !proj_.empty()); + + param.get("LaDInDegrees", LaDInDegrees_); + param.get("orientationOfTheGridInDegrees", orientationOfTheGridInDegrees_); + + long south = 0; + param.get("southPoleOnProjectionPlane", south); + southPoleOnProjectionPlane_ = south != 0; + + long uvRelativeToGrid = 0; + uvRelativeToGrid_ = param.get("uvRelativeToGrid", uvRelativeToGrid) && uvRelativeToGrid != 0; +} + + +void PolarStereographic::fillGrib(grib_info& info) const { + info.grid.grid_type = CODES_UTIL_GRID_SPEC_POLAR_STEREOGRAPHIC; + + Point2 first = {x().front(), y().front()}; + Point2 firstLL = grid().projection().lonlat(first); + + info.grid.latitudeOfFirstGridPointInDegrees = firstLL[LLCOORDS::LAT]; + info.grid.longitudeOfFirstGridPointInDegrees = + writeLonPositive_ ? util::normalise_longitude(firstLL[LLCOORDS::LON], 0) : firstLL[LLCOORDS::LON]; + + info.grid.Ni = static_cast(x().size()); + info.grid.Nj = static_cast(y().size()); + + info.grid.uvRelativeToGrid = uvRelativeToGrid_ ? 1 : 0; + + info.extra_set("DxInMetres", std::abs(x().step())); + info.extra_set("DyInMetres", std::abs(y().step())); + info.extra_set("orientationOfTheGridInDegrees", util::normalise_longitude(orientationOfTheGridInDegrees_, 0)); + + if (writeLaDInDegrees_) { + info.extra_set("LaDInDegrees", LaDInDegrees_); + } + + // some extra keys are edition-specific, so parent call is here + RegularGrid::fillGrib(info); +} + + +void PolarStereographic::fillJob(api::MIRJob& job) const { + Point2 first = {x().front(), y().front()}; + Point2 firstLL = grid().projection().lonlat(first); + + std::ostringstream grid; + grid << "gridType=polar_stereographic;" + << "proj=" << proj_ << ";" + << "LaDInDegrees=" << LaDInDegrees_ << ";" + << "orientationOfTheGridInDegrees=" << orientationOfTheGridInDegrees_ << ";" + << "southPoleOnProjectionPlane=" << (southPoleOnProjectionPlane_ ? "1" : "0") << ";" + << "Ni=" << x().size() << ";" + << "Nj=" << y().size() << ";" + << "grid=" << std::abs(x().step()) << "/" << std::abs(y().step()) << ";" + << "latitudeOfFirstGridPointInDegrees=" << firstLL[LLCOORDS::LAT] << ";" + << "longitudeOfFirstGridPointInDegrees=" << firstLL[LLCOORDS::LON]; + job.set("grid", grid.str()); -void PolarStereographic::fillGrib(grib_info& /*info*/) const { - NOTIMP; + // some extra keys are edition-specific, so parent call is here + RegularGrid::fillJob(job); } diff --git a/src/mir/repres/regular/PolarStereographic.h b/src/mir/repres/regular/PolarStereographic.h index 35afa7bea..0c1a3b12b 100644 --- a/src/mir/repres/regular/PolarStereographic.h +++ b/src/mir/repres/regular/PolarStereographic.h @@ -12,6 +12,8 @@ #pragma once +#include + #include "mir/repres/regular/RegularGrid.h" @@ -50,7 +52,14 @@ class PolarStereographic final : public RegularGrid { private: // -- Members - // None + + std::string proj_; + double LaDInDegrees_; + double orientationOfTheGridInDegrees_; + bool southPoleOnProjectionPlane_; + bool writeLaDInDegrees_; + bool writeLonPositive_; + bool uvRelativeToGrid_; // -- Methods // None @@ -58,6 +67,7 @@ class PolarStereographic final : public RegularGrid { // -- Overridden methods void fillGrib(grib_info&) const override; + void fillJob(api::MIRJob&) const override; // -- Class members // None diff --git a/src/mir/repres/regular/RegularGrid.cc b/src/mir/repres/regular/RegularGrid.cc index e27f10c5f..e5a3e045e 100644 --- a/src/mir/repres/regular/RegularGrid.cc +++ b/src/mir/repres/regular/RegularGrid.cc @@ -23,6 +23,7 @@ #include "eckit/utils/MD5.h" #include "eckit/utils/StringTools.h" +#include "mir/api/MIRJob.h" #include "mir/param/MIRParametrisation.h" #include "mir/repres/Iterator.h" #include "mir/util/Domain.h" @@ -207,21 +208,31 @@ ::atlas::Grid RegularGrid::atlasGrid() const { } -bool RegularGrid::isPeriodicWestEast() const { - return false; -} - - void RegularGrid::fillGrib(grib_info& info) const { // shape of the reference system shape_.fillGrib(info, grid_.projection().spec()); // scanningMode - info.grid.iScansNegatively = x_.back() < x_.front() ? 1L : 0L; + info.grid.iScansNegatively = x_.front() < x_.back() ? 0L : 1L; info.grid.jScansPositively = y_.front() < y_.back() ? 1L : 0L; } +void RegularGrid::fillJob(api::MIRJob& job) const { + // shape of the reference system + shape_.fillJob(job, grid_.projection().spec()); + + // scanningMode + std::string grid; + ASSERT(job.get("grid", grid) && !grid.empty()); + + grid += ";iScansNegatively=" + std::to_string(x().front() < x().back() ? 0 : 1); + grid += ";jScansPositively=" + std::to_string(y().front() < y().back() ? 1 : 0); + + job.set("grid", grid); +} + + bool RegularGrid::includesNorthPole() const { return bbox_.north() == Latitude::NORTH_POLE; } @@ -232,6 +243,12 @@ bool RegularGrid::includesSouthPole() const { } +bool RegularGrid::isPeriodicWestEast() const { + return includesNorthPole() || includesSouthPole() || + (bbox_.east().value() - bbox_.west().value() >= Longitude::GLOBE.value()); +} + + void RegularGrid::reorder(long /*scanningMode*/, MIRValuesVector& /*unused*/) const { // do not reorder, iterator is doing the right thing // FIXME this function should not be overriding to do nothing diff --git a/src/mir/repres/regular/RegularGrid.h b/src/mir/repres/regular/RegularGrid.h index c8d1000d8..ddde26de0 100644 --- a/src/mir/repres/regular/RegularGrid.h +++ b/src/mir/repres/regular/RegularGrid.h @@ -102,7 +102,9 @@ class RegularGrid : public Gridded { bool isPeriodicWestEast() const override; void fillGrib(grib_info&) const override; + void fillJob(api::MIRJob&) const override; void fillMeshGen(util::MeshGeneratorParameters&) const override; + void reorder(long scanningMode, MIRValuesVector&) const override; void validate(const MIRValuesVector&) const override; void makeName(std::ostream&) const override; diff --git a/src/mir/stats/detail/ModeT.h b/src/mir/stats/detail/ModeT.h index 43035971f..2f90a82cc 100644 --- a/src/mir/stats/detail/ModeT.h +++ b/src/mir/stats/detail/ModeT.h @@ -68,7 +68,7 @@ struct ModeIntegral : Mode { ModeIntegral(bool disambiguateMax); ModeIntegral(const param::MIRParametrisation&); - virtual ~ModeIntegral() override = default; + ~ModeIntegral() override = default; void operator()(const double& value) override { binCount_[static_cast(std::lround(value))]++; } void operator()(const float& value) override { binCount_[static_cast(std::lround(value))]++; } @@ -90,7 +90,7 @@ struct ModeReal : Mode { ModeReal(bool disambiguateMax, const std::vector& values = {0, 1}, const std::vector& mins = {0.5}); ModeReal(const param::MIRParametrisation&); - virtual ~ModeReal() override = default; + ~ModeReal() override = default; void operator()(const double& value) override { size_t bin = 0; diff --git a/src/mir/util/Atlas.cc b/src/mir/util/Atlas.cc index 210d67d73..9f43978b4 100644 --- a/src/mir/util/Atlas.cc +++ b/src/mir/util/Atlas.cc @@ -11,9 +11,8 @@ #include "mir/util/Atlas.h" -#if mir_HAVE_ATLAS -#error "mir/util/Atlas.cc cannot be included with Atlas present" -#endif + +static_assert(!mir_HAVE_ATLAS, "mir/util/Atlas.cc cannot be included with Atlas present"); #include diff --git a/src/mir/util/Atlas.h b/src/mir/util/Atlas.h index 698b721ca..6801e5da8 100644 --- a/src/mir/util/Atlas.h +++ b/src/mir/util/Atlas.h @@ -52,8 +52,6 @@ #include "atlas/util/Point.h" #include "atlas/util/Rotation.h" #else -#define ATLAS_HAVE_PROJ 0 -#define ATLAS_HAVE_TESSELATION 0 namespace eckit { diff --git a/src/mir/util/Exceptions.h b/src/mir/util/Exceptions.h index 0c3eb5f1a..b7bbf8d00 100644 --- a/src/mir/util/Exceptions.h +++ b/src/mir/util/Exceptions.h @@ -55,4 +55,14 @@ class CannotConvert : public eckit::Exception { }; +class InvalidWeightMatrix : public eckit::Exception { +public: + InvalidWeightMatrix(const char* when, const std::string& what) { + std::ostringstream os; + os << "Invalid weight matrix (" << when << "): " << what; + reason(os.str()); + } +}; + + } // namespace mir::exception diff --git a/src/mir/util/MeshGeneratorParameters.h b/src/mir/util/MeshGeneratorParameters.h index 831401e4f..38af100d3 100644 --- a/src/mir/util/MeshGeneratorParameters.h +++ b/src/mir/util/MeshGeneratorParameters.h @@ -42,7 +42,7 @@ class MeshGeneratorParameters : public atlas::MeshGenerator::Parameters { // -- Destructor - virtual ~MeshGeneratorParameters() override = default; + ~MeshGeneratorParameters() override = default; // -- Convertors // None diff --git a/src/mir/util/Point2ToPoint3.cc b/src/mir/util/Point2ToPoint3.cc index 7659e7c7e..8e016b104 100644 --- a/src/mir/util/Point2ToPoint3.cc +++ b/src/mir/util/Point2ToPoint3.cc @@ -37,9 +37,10 @@ Point2ToPoint3::Point3Displaced::Point3Displaced(double poleDisplacement) : eps_ Point3 Point2ToPoint3::Point3Displaced::make_point3(const Point2& p2) const { - const Point2 p2d{p2[0], Latitude::NORTH_POLE == p2[1] ? p2[1] - eps_ - : Latitude::SOUTH_POLE == p2[1] ? p2[1] + eps_ - : p2[1]}; + const Point2 p2d{Latitude::NORTH_POLE <= p2[0] + eps_ ? p2[0] - eps_ + : Latitude::SOUTH_POLE >= p2[0] - eps_ ? p2[0] + eps_ + : p2[0], + p2[1]}; return repres::Iterator::point_3D(p2d); } diff --git a/src/mir/util/Shape.cc b/src/mir/util/Shape.cc index ea16bee55..da39ec82e 100644 --- a/src/mir/util/Shape.cc +++ b/src/mir/util/Shape.cc @@ -12,15 +12,17 @@ #include "mir/util/Shape.h" +#include + +#include "mir/api/MIRJob.h" #include "mir/param/MIRParametrisation.h" +#include "mir/util/Grib.h" namespace mir::util { Shape::Shape(const param::MIRParametrisation& param) { - param.get("edition", edition = 0); - provided = param.get("shapeOfTheEarth", code = 6); bool isOblate = false; @@ -30,7 +32,7 @@ Shape::Shape(const param::MIRParametrisation& param) { } -Shape::Shape(const Projection::Spec& spec) : edition(0) { +Shape::Shape(const Projection::Spec& spec) { if (spec.has("radius")) { code = 1L; a = b = spec.getDouble("radius"); @@ -50,14 +52,20 @@ Shape::Shape(const Projection::Spec& spec) : edition(0) { } -Shape::Shape(const Shape&) = default; - +void Shape::fillGrib(grib_info& info, const Projection::Spec& spec) const { + const static std::string EDITION{"edition"}; -Shape& Shape::operator=(const Shape&) = default; + // GRIB2 encoding of user-provided shape + auto edition = info.packing.editionNumber; + for (long j = 0; j < info.packing.extra_settings_count; ++j) { + const auto& set = info.packing.extra_settings[j]; + if (set.name == EDITION && set.type == CODES_TYPE_LONG) { + edition = set.long_value; + break; + } + } -void Shape::fillGrib(grib_info& info, const Projection::Spec& spec) const { - // GRIB2 encoding of user-provided shape if (edition != 2) { return; } @@ -93,4 +101,44 @@ void Shape::fillGrib(grib_info& info, const Projection::Spec& spec) const { } +void Shape::fillJob(api::MIRJob& job, const Projection::Spec& spec) const { + // MIRJob encoding of user-provided shape + std::ostringstream shape; + + // shape given by radius or semi-major/minor axis + if (provided) { + shape << ";shapeOfTheEarth=" << code; + switch (code) { + case 1: + shape << ";radius", spec.getDouble("radius", a); + break; + case 3: + shape << ";earthMajorAxis=" << spec.getDouble("semi_major_axis", a) / 1000.; + shape << ";earthMinorAxis=" << spec.getDouble("semi_minor_axis", b) / 1000.; + break; + case 7: + shape << ";earthMajorAxis=" << spec.getDouble("semi_major_axis", a); + shape << ";earthMinorAxis=" << spec.getDouble("semi_minor_axis", b); + break; + default: + break; + } + } + else if (spec.has("radius")) { + shape << ";shapeOfTheEarth=" << 1L; + shape << ";radius", spec.getDouble("radius"); + } + else if (spec.has("semi_major_axis") && spec.has("semi_minor_axis")) { + shape << ";shapeOfTheEarth=" << 7L; + shape << ";earthMajorAxis=" << spec.getDouble("semi_major_axis"); + shape << ";earthMinorAxis=" << spec.getDouble("semi_minor_axis"); + } + + std::string grid; + ASSERT(job.get("grid", grid) && !grid.empty()); + + job.set("grid", grid + shape.str()); +} + + } // namespace mir::util diff --git a/src/mir/util/Shape.h b/src/mir/util/Shape.h index a538bd1d6..a9c7e4945 100644 --- a/src/mir/util/Shape.h +++ b/src/mir/util/Shape.h @@ -16,9 +16,14 @@ #include "mir/util/Grib.h" -namespace mir::param { +namespace mir { +namespace param { class MIRParametrisation; -} // namespace mir::param +} +namespace api { +class MIRJob; +} +} // namespace mir namespace mir::util { @@ -30,16 +35,21 @@ struct Shape { Shape(const param::MIRParametrisation&); Shape(const Projection::Spec&); Shape() : code(6L), a(0.), b(0.), provided(false) {} - Shape(const Shape& other); + + Shape(const Shape& other) = default; + Shape(Shape&&) = default; + virtual ~Shape() = default; - Shape& operator=(const Shape& other); - void fillGrib(grib_info& info, const Projection::Spec&) const; + Shape& operator=(const Shape&) = default; + Shape& operator=(Shape&&) = default; + + void fillGrib(grib_info&, const Projection::Spec&) const; + void fillJob(api::MIRJob&, const Projection::Spec&) const; long code; double a; double b; - long edition; bool provided; }; diff --git a/src/tools/CMakeLists.txt b/src/tools/CMakeLists.txt index 3a1419f0f..38472cd14 100644 --- a/src/tools/CMakeLists.txt +++ b/src/tools/CMakeLists.txt @@ -2,6 +2,7 @@ ecbuild_add_executable(TARGET mir-tool OUTPUT_NAME mir SOURCES mir.cc LIBS mir ${mir_INSTALL_TOOLS}) foreach(tool IN ITEMS + mir-bounding-box mir-climate-filter mir-compare # NOTE: set in the testing scripts, however not used mir-compute @@ -30,6 +31,13 @@ if(eckit_HAVE_ECKIT_CODEC) mir-codec-to-weight-matrix) ecbuild_add_executable(TARGET ${tool} SOURCES ${tool}.cc LIBS mir eckit_codec ${mir_INSTALL_TOOLS}) endforeach() + + if(mir_HAVE_NETCDF) + foreach(tool IN ITEMS + mir-scrip-to-weight-matrix) + ecbuild_add_executable(TARGET ${tool} SOURCES ${tool}.cc LIBS mir eckit_codec NetCDF::NetCDF_CXX ${mir_INSTALL_TOOLS}) + endforeach() + endif() endif() if(mir_HAVE_ATLAS) diff --git a/src/tools/mir-bounding-box.cc b/src/tools/mir-bounding-box.cc new file mode 100644 index 000000000..e87eb45ec --- /dev/null +++ b/src/tools/mir-bounding-box.cc @@ -0,0 +1,73 @@ +/* + * (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 "eckit/option/CmdArgs.h" +#include "eckit/option/SimpleOption.h" + +#include "mir/data/MIRField.h" +#include "mir/input/MIRInput.h" +#include "mir/param/ConfigurationWrapper.h" +#include "mir/repres/Representation.h" +#include "mir/tools/MIRTool.h" +#include "mir/util/BoundingBox.h" +#include "mir/util/Domain.h" +#include "mir/util/Log.h" + + +namespace mir::tools { + + +struct MIRBoundingBox : MIRTool { + using MIRTool::MIRTool; + + MIRBoundingBox(int argc, char** argv) : MIRTool(argc, argv) { + options_.push_back(new eckit::option::SimpleOption("domain", "")); + } + + void usage(const std::string& tool) const override { + Log::info() << "\n" + "Usage: " + << tool << "[FILE [FILE [FILE]]]" << std::endl; + } + + void execute(const eckit::option::CmdArgs& args) override { + const param::ConfigurationWrapper args_wrap(args); + + auto domain = false; + args_wrap.get("domain", domain); + + for (const auto& arg : args) { + std::unique_ptr input(input::MIRInputFactory::build(arg, args_wrap)); + + while (input->next()) { + auto field = input->field(); + repres::RepresentationHandle representation(field.representation()); + + if (domain) { + Log::info() << representation->domain() << std::endl; + } + else { + Log::info() << representation->boundingBox() << std::endl; + } + } + } + } +}; + + +} // namespace mir::tools + + +int main(int argc, char** argv) { + mir::tools::MIRBoundingBox tool(argc, argv); + return tool.start(); +} diff --git a/src/tools/mir-matrix-reorder.cc b/src/tools/mir-matrix-reorder.cc index afbfce991..45860774f 100644 --- a/src/tools/mir-matrix-reorder.cc +++ b/src/tools/mir-matrix-reorder.cc @@ -22,7 +22,6 @@ #include "eckit/option/SimpleOption.h" #include "mir/reorder/Reorder.h" -#include "mir/repres/Representation.h" #include "mir/tools/MIRTool.h" #include "mir/util/Exceptions.h" #include "mir/util/Log.h" @@ -65,6 +64,8 @@ void MIRMatrixReorder::execute(const eckit::option::CmdArgs& args) { eckit::linalg::SparseMatrix M; M.load(args(0)); + Log::info() << M << std::endl; + // renumbering maps auto rows = diff --git a/src/tools/mir-scrip-to-weight-matrix.cc b/src/tools/mir-scrip-to-weight-matrix.cc new file mode 100644 index 000000000..91bfd20c1 --- /dev/null +++ b/src/tools/mir-scrip-to-weight-matrix.cc @@ -0,0 +1,221 @@ +/* + * (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 +#include + +#include "netcdf" + +#include "eckit/filesystem/PathName.h" +#include "eckit/linalg/SparseMatrix.h" +#include "eckit/option/CmdArgs.h" +#include "eckit/option/SimpleOption.h" + +#include "mir/tools/MIRTool.h" +#include "mir/util/Exceptions.h" +#include "mir/util/Log.h" + + +namespace mir::tools { + + +using Index = eckit::linalg::Index; +using Scalar = eckit::linalg::Scalar; +using Size = eckit::linalg::Size; + + +static std::string DEFAULT_DIM_NR = "dst_grid_size"; +static std::string DEFAULT_DIM_NC = "src_grid_size"; +static std::string DEFAULT_DIM_NNZ = "num_links"; +static std::string DEFAULT_VAR_IA = "dst_address"; +static std::string DEFAULT_VAR_JA = "src_address"; +static std::string DEFAULT_VAR_A = "remap_matrix"; +static int DEFAULT_BASE = 1; + + +class InPlaceAllocator : public eckit::linalg::SparseMatrix::Allocator { +public: + using Layout = eckit::linalg::SparseMatrix::Layout; + using Shape = eckit::linalg::SparseMatrix::Shape; + + InPlaceAllocator(Size Nr, Size Nc, Size nnz, Index* ia, Index* ja, Scalar* a) : + Nr_(Nr), Nc_(Nc), nnz_(nnz), ia_(ia), ja_(ja), a_(a) { + ASSERT(ia_ != nullptr); + ASSERT(ja_ != nullptr); + ASSERT(a_ != nullptr); + } + + Layout allocate(Shape& shape) override { + shape.size_ = nnz_; + shape.rows_ = Nr_; + shape.cols_ = Nc_; + + Layout layout; + layout.outer_ = reinterpret_cast(ia_); + layout.inner_ = ja_; + layout.data_ = a_; + + return layout; + } + + void deallocate(Layout, Shape) override {} + + void print(std::ostream&) const override { NOTIMP; } + + bool inSharedMemory() const override { return false; } + +private: + const Size Nr_; + const Size Nc_; + const Size nnz_; + Index* ia_; // NOTE: not owned + Index* ja_; + Scalar* a_; +}; + + +struct MIRScripToWeightMatrix : public MIRTool { + MIRScripToWeightMatrix(int argc, char** argv) : MIRTool(argc, argv) { + using eckit::option::SimpleOption; + + options_.push_back( + new SimpleOption("nr", "Dimension number of rows (default '" + DEFAULT_DIM_NR + "')")); + options_.push_back( + new SimpleOption("nc", "Dimension number of columns (default '" + DEFAULT_DIM_NC + "')")); + options_.push_back(new SimpleOption( + "nnz", "Dimension number of non-zeros (for checking) (default '" + DEFAULT_DIM_NNZ + "')")); + options_.push_back(new SimpleOption( + "ia", "Variable Matrix destination/rows indices, not compressed (default '" + DEFAULT_VAR_IA + "')")); + options_.push_back(new SimpleOption( + "ja", "Variable Matrix source/column indices (default '" + DEFAULT_VAR_JA + "')")); + options_.push_back(new SimpleOption( + "a", "Variable Matrix interpolation weights (default '" + DEFAULT_VAR_A + "')")); + + options_.push_back( + new SimpleOption("base", "Indices base (default " + std::to_string(DEFAULT_BASE) + ")")); + } + + int numberOfPositionalArguments() const override { return 2; } + + void usage(const std::string& tool) const override { + Log::info() << "\n" + << "Conversion of specificly formated SCRIP interpolation weights file into " + "eckit::linalg::SparseMatrix file\n" + << "\n" + << "Usage: " << tool << " " << std::endl; + } + + void execute(const eckit::option::CmdArgs& args) override { + ASSERT(args.count() == numberOfPositionalArguments()); + + Size Nr = 0; + Size Nc = 0; + Size nnz = 0; + + std::vector ia; + std::vector ja; + std::vector a; + + + // read SCRIP file + try { + netCDF::NcFile f(args(0), netCDF::NcFile::read); + + Nr = f.getDim(args.getString("nr", DEFAULT_DIM_NR)).getSize(); + Nc = f.getDim(args.getString("nc", DEFAULT_DIM_NC)).getSize(); + nnz = f.getDim(args.getString("nnz", DEFAULT_DIM_NNZ)).getSize(); + + ASSERT(Nr > 0); + ASSERT(Nc > 0); + ASSERT(nnz > 0); + + auto var_ia = f.getVar(args.getString("ia", DEFAULT_VAR_IA)); + ASSERT(var_ia.getDimCount() == 1 && var_ia.getDim(0).getSize() == nnz); + ia.resize(nnz); // NOTE: not compressed + var_ia.getVar(ia.data()); + + auto var_ja = f.getVar(args.getString("ja", DEFAULT_VAR_JA)); + ASSERT(var_ja.getDimCount() == 1 && var_ja.getDim(0).getSize() == nnz); + ja.resize(nnz); + var_ja.getVar(ja.data()); + + auto var_a = f.getVar(args.getString("a", DEFAULT_VAR_A)); + ASSERT(var_a.getDimCount() == 2 && var_a.getDim(0).getSize() == nnz && var_a.getDim(1).getSize() == 1); + a.resize(nnz); + var_a.getVar(a.data()); + } + catch (netCDF::exceptions::NcException& e) { + throw eckit::SeriousBug("SCRIP reading failed", Here()); + } + + + // matrix conversion to 0-based indexing, sorting and compression + { + std::vector sorted(nnz); + for (size_t n = 0; n < nnz; ++n) { + sorted[n] = n; + } + + std::sort(sorted.begin(), sorted.end(), [&](auto i, auto j) { + return ia[i] != ia[j] ? ia[i] < ia[j] : ja[i] != ja[j] ? ja[i] < ja[j] : i < j; + }); + + decltype(ia) ia_unsorted(Nr + 1, 0); + decltype(ja) ja_unsorted(nnz); + ia.swap(ia_unsorted); + ja.swap(ja_unsorted); + + const auto base = static_cast(args.getUnsigned("base", DEFAULT_BASE)); + for (size_t n = 0; n < nnz; ++n) { + auto r = ia_unsorted[sorted[n]] - base; + auto c = ja_unsorted[sorted[n]] - base; + ASSERT(0 <= r && r < Nr); + ASSERT(0 <= c && c < Nc); + + ia[r + 1]++; + ja[n] = c; + } + + ia_unsorted.clear(); + ja_unsorted.clear(); + + for (size_t r = 0; r < Nr; ++r) { + ia[r + 1] += ia[r]; + } + + ASSERT(ia.front() == 0 && ia.back() == nnz); + + decltype(a) a_unsorted(nnz); + a.swap(a_unsorted); + + for (size_t n = 0; n < nnz; ++n) { + a[n] = a_unsorted[sorted[n]]; + ASSERT(0. <= a[n] && a[n] <= 1.); + } + } + + // create matrix + eckit::linalg::SparseMatrix M(new InPlaceAllocator(Nr, Nc, nnz, ia.data(), ja.data(), a.data())); + + // eckit::linalg::SparseMatrix file + eckit::PathName fmat(args(1)); + M.save(fmat); + } +}; + + +} // namespace mir::tools + + +int main(int argc, char** argv) { + return mir::tools::MIRScripToWeightMatrix{argc, argv}.start(); +} diff --git a/src/tools/mir.cc b/src/tools/mir.cc index 2e3e4da72..a97294de5 100644 --- a/src/tools/mir.cc +++ b/src/tools/mir.cc @@ -115,6 +115,9 @@ struct MIR : MIRTool { options_.push_back( new FactoryOption("interpolation", "Grid to grid interpolation method")); + options_.push_back(new FactoryOption("intermediate-interpolation", + "Grid to grid intermediate interpolation method")); + options_.push_back( new FactoryOption("interpolation-statistics", "Statistics interpolation method")); @@ -189,38 +192,41 @@ struct MIR : MIRTool { "Linear algebra sparse backend (default '" + eckit::linalg::LinearAlgebraSparse::backend().name() + "')")); options_.push_back(new FactoryOption("point-search-trees", "k-d tree control")); -#if mir_HAVE_ATLAS - for (const auto& which : std::vector{"input", "output"}) { - options_.push_back( - new SimpleOption(which + "-mesh-generator", "Mesh generator for " + which + " grid")); - options_.push_back( - new SimpleOption(which + "-mesh-cell-centres", "Calculate cell centres for " + which + " mesh")); - options_.push_back(new SimpleOption(which + "-mesh-cell-longest-diagonal", - "Calculate cells longest diagonal for " + which + " mesh")); - options_.push_back(new SimpleOption(which + "-mesh-node-lumped-mass-matrix", - "Calculate node-lumped mass matrix for " + which + " mesh")); - options_.push_back(new SimpleOption(which + "-mesh-node-to-cell-connectivity", - "Calculate node-to-cell connectivity for " + which + " mesh")); - options_.push_back(new SimpleOption( - which + "-mesh-file-ll", - "Output file for " + which + " grid, in lon/lat coordinates (default )")); - options_.push_back(new SimpleOption( - which + "-mesh-file-xy", "Output file for " + which + " grid, in X/Y coordinates (default )")); - options_.push_back(new SimpleOption( - which + "-mesh-file-xyz", - "Output file for " + which + " grid, in X/Y/Z coordinates (default )")); - options_.push_back(new SimpleOption(which + "-mesh-generator-three-dimensional", - "Generate 3-dimensional " + which + " mesh")); - options_.push_back(new SimpleOption(which + "-mesh-generator-triangulate", - "Generate triangulated " + which + " mesh")); - options_.push_back(new SimpleOption( - which + "-mesh-generator-angle", "Generate with quadrilateral tolerance angle on " + which + " mesh")); - options_.push_back(new SimpleOption(which + "-mesh-generator-force-include-north-pole", - "Generate including North pole on " + which + " mesh")); - options_.push_back(new SimpleOption(which + "-mesh-generator-force-include-south-pole", - "Generate including South pole on " + which + " mesh")); + if constexpr (HAVE_ATLAS) { + for (const auto& which : std::vector{"input", "output"}) { + options_.push_back( + new SimpleOption(which + "-mesh-generator", "Mesh generator for " + which + " grid")); + options_.push_back(new SimpleOption(which + "-mesh-cell-centres", + "Calculate cell centres for " + which + " mesh")); + options_.push_back(new SimpleOption(which + "-mesh-cell-longest-diagonal", + "Calculate cells longest diagonal for " + which + " mesh")); + options_.push_back(new SimpleOption(which + "-mesh-node-lumped-mass-matrix", + "Calculate node-lumped mass matrix for " + which + " mesh")); + options_.push_back( + new SimpleOption(which + "-mesh-node-to-cell-connectivity", + "Calculate node-to-cell connectivity for " + which + " mesh")); + options_.push_back(new SimpleOption( + which + "-mesh-file-ll", + "Output file for " + which + " grid, in lon/lat coordinates (default )")); + options_.push_back(new SimpleOption( + which + "-mesh-file-xy", + "Output file for " + which + " grid, in X/Y coordinates (default )")); + options_.push_back(new SimpleOption( + which + "-mesh-file-xyz", + "Output file for " + which + " grid, in X/Y/Z coordinates (default )")); + options_.push_back(new SimpleOption(which + "-mesh-generator-three-dimensional", + "Generate 3-dimensional " + which + " mesh")); + options_.push_back(new SimpleOption(which + "-mesh-generator-triangulate", + "Generate triangulated " + which + " mesh")); + options_.push_back( + new SimpleOption(which + "-mesh-generator-angle", + "Generate with quadrilateral tolerance angle on " + which + " mesh")); + options_.push_back(new SimpleOption(which + "-mesh-generator-force-include-north-pole", + "Generate including North pole on " + which + " mesh")); + options_.push_back(new SimpleOption(which + "-mesh-generator-force-include-south-pole", + "Generate including South pole on " + which + " mesh")); + } } -#endif options_.push_back( new SimpleOption("counter-upper-limit", "Statistics count values below lower limit")); @@ -361,10 +367,10 @@ struct MIR : MIRTool { "legendre-loader", "Select how to load Legendre coefficients in memory")); #endif -#if mir_HAVE_OMP - options_.push_back( - new SimpleOption("parallel-omp-num-threads", "Set number of threads for OMP parallel regions")); -#endif + if constexpr (HAVE_OMP) { + options_.push_back( + new SimpleOption("parallel-omp-num-threads", "Set number of threads for OMP parallel regions")); + } //============================================== // Only show these options if debug channel is active diff --git a/tests/assertions/CMakeLists.txt b/tests/assertions/CMakeLists.txt index 50b11e743..0f7c8ca8a 100644 --- a/tests/assertions/CMakeLists.txt +++ b/tests/assertions/CMakeLists.txt @@ -23,7 +23,8 @@ if(MARS_SCRIPT) ENVIRONMENT ${_testEnvironment}) if((_t MATCHES "[.]fail[.]") OR (_t MATCHES "[.]lsm[.]" AND NOT HAVE_MIR_DOWNLOAD_MASKS) OR - (NOT mir_HAVE_ATLAS AND NOT (_t MATCHES "[.]core[.]"))) + (NOT mir_HAVE_ATLAS AND NOT (_t MATCHES "[.]core[.]")) OR + (NOT atlas_HAVE_PROJ AND (_t MATCHES "[.]proj[.]"))) set_tests_properties(${_t} PROPERTIES WILL_FAIL TRUE) endif() endforeach() diff --git a/tests/assertions/ICON.griddef b/tests/assertions/ICON.griddef new file mode 120000 index 000000000..817c2a817 --- /dev/null +++ b/tests/assertions/ICON.griddef @@ -0,0 +1 @@ +../data/ICON.griddef \ No newline at end of file diff --git a/tests/assertions/MIR-592.003.test b/tests/assertions/MIR-592.003.test index 50064db75..be41e70b2 100644 --- a/tests/assertions/MIR-592.003.test +++ b/tests/assertions/MIR-592.003.test @@ -1,16 +1,16 @@ # input param=u_v,level=1000,resol=20 # mir ---grid=3/3 --interpolation=linear -# grib_get (placeholder until values can be checked) +--grid=3/3 --interpolation=linear --uv2uv +# grib_get numberOfPoints=7320 -values[0]=0.000213623 -values[1]=0.000213623 -values[2]=0.000213623 -values[3]=0.000213623 -values[4]=0.000213623 -values[5]=0.000213623 -values[6]=0.000213623 -values[7]=0.000213623 -values[8]=0.000213623 -values[9]=0.000213623 +values[0]=8.1329 +values[1]=7.73153 +values[2]=7.32626 +values[3]=6.91806 +values[4]=6.51181 +values[5]=6.11434 +values[6]=5.72567 +values[7]=5.35263 +values[8]=5.00399 +values[9]=4.67587 diff --git a/tests/assertions/MIR-673.001.test b/tests/assertions/MIR-673.001.test new file mode 100644 index 000000000..75140c119 --- /dev/null +++ b/tests/assertions/MIR-673.001.test @@ -0,0 +1,11 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--grid=gridType=polar_stereographic +# grib_get +gridType=polar_stereographic +edition=2 +Ni=2869 +Nj=2869 +Dx=2500000 +Dy=2500000 diff --git a/tests/assertions/MIR-673.002.proj.test b/tests/assertions/MIR-673.002.proj.test new file mode 100644 index 000000000..b9aa4e7b7 --- /dev/null +++ b/tests/assertions/MIR-673.002.proj.test @@ -0,0 +1,11 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--area=66.5/-24.5/63.4/-13.5 +# grib_get +gridType=polar_stereographic +edition=2 +Ni=2869 +Nj=2869 +Dx=2500000 +Dy=2500000 diff --git a/tests/assertions/MIR-673.003.proj.fail.test b/tests/assertions/MIR-673.003.proj.fail.test new file mode 100644 index 000000000..d7aef75c2 --- /dev/null +++ b/tests/assertions/MIR-673.003.proj.fail.test @@ -0,0 +1,4 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--area=66.5/-24.5/63.4/-13.5 --area-mode=crop diff --git a/tests/assertions/MIR-673.004.proj.test b/tests/assertions/MIR-673.004.proj.test new file mode 100644 index 000000000..1f3bc5d69 --- /dev/null +++ b/tests/assertions/MIR-673.004.proj.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--grid=1/1 +# grib_get +gridType=regular_ll diff --git a/tests/assertions/MIR-673.005.test b/tests/assertions/MIR-673.005.test new file mode 100644 index 000000000..f057afaca --- /dev/null +++ b/tests/assertions/MIR-673.005.test @@ -0,0 +1,11 @@ +# input +gridType=polar_stereographic.grib1 +# mir +--grid=gridType=polar_stereographic +# grib_get +gridType=polar_stereographic +edition=1 +Ni=1121 +Nj=881 +Dx=4763 +Dy=4763 diff --git a/tests/assertions/MIR-673.006.proj.test b/tests/assertions/MIR-673.006.proj.test new file mode 100644 index 000000000..7255051aa --- /dev/null +++ b/tests/assertions/MIR-673.006.proj.test @@ -0,0 +1,11 @@ +# input +gridType=polar_stereographic.grib1 +# mir +--area=49/-92/41/-76 +# grib_get +gridType=polar_stereographic +edition=1 +Ni=1121 +Nj=881 +Dx=4763 +Dy=4763 diff --git a/tests/assertions/MIR-673.007.proj.fail.test b/tests/assertions/MIR-673.007.proj.fail.test new file mode 100644 index 000000000..a4077c150 --- /dev/null +++ b/tests/assertions/MIR-673.007.proj.fail.test @@ -0,0 +1,4 @@ +# input +gridType=polar_stereographic.grib1 +# mir +--area=49/-92/41/-76 --area-mode=crop diff --git a/tests/assertions/MIR-673.008.proj.test b/tests/assertions/MIR-673.008.proj.test new file mode 100644 index 000000000..1b0713b5b --- /dev/null +++ b/tests/assertions/MIR-673.008.proj.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib1 +# mir +--grid=1/1 +# grib_get +gridType=regular_ll diff --git a/tests/assertions/MIR-679.ICON-to-WP.core.test b/tests/assertions/MIR-679.ICON-to-WP.core.test new file mode 100644 index 000000000..d363a0b6d --- /dev/null +++ b/tests/assertions/MIR-679.ICON-to-WP.core.test @@ -0,0 +1,7 @@ +# input +ICON.shortName=lsm.grib2 +# mir +--input={griddef:$tests_assertions_dir/ICON.griddef} --griddef=$tests_assertions_dir/weather-params.pts --format=grib --interpolation=nn +# grib_get +gridType=unstructured_grid +numberOfDataPoints=7693 diff --git a/tests/assertions/MIR-679.WP-to-ICON.core.test b/tests/assertions/MIR-679.WP-to-ICON.core.test new file mode 100644 index 000000000..1fb372055 --- /dev/null +++ b/tests/assertions/MIR-679.WP-to-ICON.core.test @@ -0,0 +1,7 @@ +# input +/dev/null +# mir +--input={artificialInput:constant,constant:0.,gridded:True,gridType:unstructured_grid,griddef:$tests_assertions_dir/weather-params.pts} --add-random=uniform-int-distribution{a:1,b:3} --griddef=$tests_assertions_dir/ICON.griddef --format=grib --interpolation=nn +# grib_get +gridType=unstructured_grid +numberOfDataPoints=20480 diff --git a/tests/assertions/MIR-681.001.test b/tests/assertions/MIR-681.001.test new file mode 100644 index 000000000..a168b7843 --- /dev/null +++ b/tests/assertions/MIR-681.001.test @@ -0,0 +1,6 @@ +# input +param=u_v,level=1000,resol=20 +# mir +--grid=3/3 --interpolation=grid-box-average --uv2uv +# grib_get (placeholder until values can be checked) +numberOfPoints=7320 diff --git a/tests/assertions/MIR-681.002.test b/tests/assertions/MIR-681.002.test new file mode 100644 index 000000000..a34b4c37c --- /dev/null +++ b/tests/assertions/MIR-681.002.test @@ -0,0 +1,7 @@ +# input +param=u_v,level=1000,resol=20 +# mir +--grid=3/3 --uv2uv +#--interpolation=grid-box-statistics --interpolation-statistics=minimum +# grib_get (placeholder until values can be checked) +numberOfPoints=7320 diff --git a/tests/assertions/gridType=polar_stereographic.grib1 b/tests/assertions/gridType=polar_stereographic.grib1 new file mode 120000 index 000000000..f5a615d69 --- /dev/null +++ b/tests/assertions/gridType=polar_stereographic.grib1 @@ -0,0 +1 @@ +../data/gridType=polar_stereographic.grib1 \ No newline at end of file diff --git a/tests/assertions/gridType=polar_stereographic.grib2 b/tests/assertions/gridType=polar_stereographic.grib2 new file mode 120000 index 000000000..6a4ca1f01 --- /dev/null +++ b/tests/assertions/gridType=polar_stereographic.grib2 @@ -0,0 +1 @@ +../data/gridType=polar_stereographic.grib2 \ No newline at end of file diff --git a/tests/assertions/intermediate-interpolation=nn.core.test b/tests/assertions/intermediate-interpolation=nn.core.test new file mode 100644 index 000000000..9f2461f95 --- /dev/null +++ b/tests/assertions/intermediate-interpolation=nn.core.test @@ -0,0 +1,7 @@ +# input +regular_ll.2-2.grib2 +# mir +--intermediate-interpolation=nn --intgrid=O80 --interpolation=linear --grid=1/1 +# grib_get +Ni=360 +Nj=181 diff --git a/tests/assertions/regrid.001.core.test b/tests/assertions/regrid.001.core.test new file mode 100644 index 000000000..76a0983f0 --- /dev/null +++ b/tests/assertions/regrid.001.core.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=1/1 +# grib_get +# (none) diff --git a/tests/assertions/regrid.002.test b/tests/assertions/regrid.002.test new file mode 100644 index 000000000..2f0a53152 --- /dev/null +++ b/tests/assertions/regrid.002.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=gridType=lambert_azimuthal_equal_area;latitudeOfFirstGridPointInDegrees=66.982143;longitudeOfFirstGridPointInDegrees=-35.034024;grid=5000./5000.;Ni=1000;Nj=950;standardParallelInDegrees=52.;centralLongitudeInDegrees=10. +# grib_get +gridType=lambert_azimuthal_equal_area diff --git a/tests/assertions/regrid.003.test b/tests/assertions/regrid.003.test new file mode 100644 index 000000000..e998a43d8 --- /dev/null +++ b/tests/assertions/regrid.003.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=gridType=lambert;latitudeOfFirstGridPointInDegrees=52.041;longitudeOfFirstGridPointInDegrees=1.639;grid=2500./2500.;Ni=739;Nj=949;LaDInDegrees=63.;LoVInDegrees=15. +# grib_get +gridType=lambert diff --git a/tests/assertions/regrid.004.core.test b/tests/assertions/regrid.004.core.test new file mode 100644 index 000000000..5d4f73f59 --- /dev/null +++ b/tests/assertions/regrid.004.core.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=O32 --interpolation=nearest-neighbour --area=20/-40/-20/40 +# grib_get +gridName=O32 diff --git a/tests/assertions/regrid.005.core.test b/tests/assertions/regrid.005.core.test new file mode 100644 index 000000000..6b8c62d05 --- /dev/null +++ b/tests/assertions/regrid.005.core.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=O64 --interpolation=grid-box-average +# grib_get +gridName=O64 diff --git a/tests/assertions/regrid.006.core.test b/tests/assertions/regrid.006.core.test new file mode 100644 index 000000000..4eda922c0 --- /dev/null +++ b/tests/assertions/regrid.006.core.test @@ -0,0 +1,7 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=1/1 --interpolation=grid-box-statistics --interpolation-statistics=maximum +# grib_get +Ni=360 +Nj=181 diff --git a/tests/assertions/regrid.007.core.test b/tests/assertions/regrid.007.core.test new file mode 100644 index 000000000..bb0295432 --- /dev/null +++ b/tests/assertions/regrid.007.core.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=N160 --interpolation=k-nearest-neighbours --nearest-method=nclosest-or-nearest --nclosest=10 --distance=100000. --distance-weighting=gaussian --distance-weighting-gaussian-stddev=10000. --non-linear=missing-if-all-missing +# grib_get +gridName=N160 diff --git a/tests/assertions/regrid.009.test b/tests/assertions/regrid.009.test new file mode 100644 index 000000000..242b65ad1 --- /dev/null +++ b/tests/assertions/regrid.009.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=O320 --interpolation=structured-bilinear +# grib_get +gridName=O320 diff --git a/tests/assertions/regrid.010.test b/tests/assertions/regrid.010.test new file mode 100644 index 000000000..9ebafc223 --- /dev/null +++ b/tests/assertions/regrid.010.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=O160 --interpolation=structured-bicubic +# grib_get +gridName=O160 diff --git a/tests/assertions/regrid.012.test b/tests/assertions/regrid.012.test new file mode 100644 index 000000000..e43a4f517 --- /dev/null +++ b/tests/assertions/regrid.012.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--nabla=scalar-gradient +# grib_get +# (none) diff --git a/tests/assertions/regrid.013.test b/tests/assertions/regrid.013.test new file mode 100644 index 000000000..e74c6fc35 --- /dev/null +++ b/tests/assertions/regrid.013.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--nabla=scalar-laplacian +# grib_get +# (none) diff --git a/tests/assertions/regrid.014.test b/tests/assertions/regrid.014.test new file mode 100644 index 000000000..726657052 --- /dev/null +++ b/tests/assertions/regrid.014.test @@ -0,0 +1,6 @@ +# input +param=u_v,level=1000,resol=20 +# mir +--nabla=uv-gradient --grid=3/3 --uv2uv +# grib_get +# (none) diff --git a/tests/assertions/regrid.015.test b/tests/assertions/regrid.015.test new file mode 100644 index 000000000..de8892a1b --- /dev/null +++ b/tests/assertions/regrid.015.test @@ -0,0 +1,6 @@ +# input +param=u_v,level=1000,resol=20 +# mir +--nabla=uv-vorticity --grid=3/3 --uv2uv +# grib_get +# (none) diff --git a/tests/assertions/regrid.016.test b/tests/assertions/regrid.016.test new file mode 100644 index 000000000..56b81664f --- /dev/null +++ b/tests/assertions/regrid.016.test @@ -0,0 +1,6 @@ +# input +param=u_v,level=1000,resol=20 +# mir +--nabla=uv-divergence --grid=3/3 --uv2uv +# grib_get +# (none) diff --git a/tests/assertions/regrid.017.core.test b/tests/assertions/regrid.017.core.test new file mode 100644 index 000000000..2e54712a1 --- /dev/null +++ b/tests/assertions/regrid.017.core.test @@ -0,0 +1,6 @@ +# input +param=z,level=1,levtype=ml,resol=20.grib2 +# mir +--truncation=19 +# grib_get +J=19 diff --git a/tests/assertions/regrid.018.test b/tests/assertions/regrid.018.test new file mode 100644 index 000000000..608904ee0 --- /dev/null +++ b/tests/assertions/regrid.018.test @@ -0,0 +1,7 @@ +# input +param=z,level=1,levtype=ml,resol=20.grib2 +# mir +--grid=1/1 --intgrid=O20 --spectral-order=cubic +# grib_get +Ni=360 +Nj=181 diff --git a/tests/assertions/regrid.019.test b/tests/assertions/regrid.019.test new file mode 100644 index 000000000..9bf3ac11f --- /dev/null +++ b/tests/assertions/regrid.019.test @@ -0,0 +1,7 @@ +# input +param=z,level=1,levtype=ml,resol=20.grib2 +# mir +--grid=1/1 --intgrid=O20 --truncation=19 +# grib_get +Ni=360 +Nj=181 diff --git a/tests/assertions/regrid.020.test b/tests/assertions/regrid.020.test new file mode 100644 index 000000000..797c87c19 --- /dev/null +++ b/tests/assertions/regrid.020.test @@ -0,0 +1,7 @@ +# input +param=vo_d,level=1,levtype=ml,resol=20.grib2 +# mir +--vod2uv +# grib_get +gridType=sh +shortName=u diff --git a/tests/assertions/regrid.021.test b/tests/assertions/regrid.021.test new file mode 100644 index 000000000..6c9b0a729 --- /dev/null +++ b/tests/assertions/regrid.021.test @@ -0,0 +1,7 @@ +# input +param=vo_d,level=1,levtype=ml,resol=20.grib2 +# mir +--grid=3/3 --vod2uv +# grib_get +gridType=regular_ll +shortName=u diff --git a/tests/assertions/regrid.022.core.test b/tests/assertions/regrid.022.core.test new file mode 100644 index 000000000..71e28377c --- /dev/null +++ b/tests/assertions/regrid.022.core.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=5/5 +# grib_get +# (none) diff --git a/tests/assertions/regrid.023.core.test b/tests/assertions/regrid.023.core.test new file mode 100644 index 000000000..014450821 --- /dev/null +++ b/tests/assertions/regrid.023.core.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=5/5 --interpolation=grid-box-average +# grib_get +# (none) diff --git a/tests/assertions/regrid.024.core.test b/tests/assertions/regrid.024.core.test new file mode 100644 index 000000000..00f325333 --- /dev/null +++ b/tests/assertions/regrid.024.core.test @@ -0,0 +1,11 @@ +# input +regular_ll.5-5.area.grib1 +# mir +--grid=3/3 --globalise=1 +# grib_get +numberOfCodedValues=1200 +numberOfMissing=6120 +latitudeOfFirstGridPointInDegrees=90 +latitudeOfLastGridPointInDegrees=-90 +longitudeOfFirstGridPointInDegrees=0 +longitudeOfLastGridPointInDegrees=357 diff --git a/tests/assertions/regrid.025.core.test b/tests/assertions/regrid.025.core.test new file mode 100644 index 000000000..ab31467d5 --- /dev/null +++ b/tests/assertions/regrid.025.core.test @@ -0,0 +1,11 @@ +# input +regular_ll.5-5.area.grib1 +# mir +--grid=3/3 --pre-globalise=1 +# grib_get +numberOfCodedValues=1416 +numberOfMissing=5904 +latitudeOfFirstGridPointInDegrees=90 +latitudeOfLastGridPointInDegrees=-90 +longitudeOfFirstGridPointInDegrees=0 +longitudeOfLastGridPointInDegrees=357 diff --git a/tests/assertions/regrid.026.core.test b/tests/assertions/regrid.026.core.test new file mode 100644 index 000000000..452db2248 --- /dev/null +++ b/tests/assertions/regrid.026.core.test @@ -0,0 +1,9 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=0.75/0.75 --basic-angle=decimal +# grib_get +latitudeOfFirstGridPointInDegrees=90 +iDirectionIncrement=750000 +basicAngleOfTheInitialProductionDomain=0 +subdivisionsOfBasicAngle=MISSING diff --git a/tests/assertions/regrid.027.core.test b/tests/assertions/regrid.027.core.test new file mode 100644 index 000000000..6d25006fa --- /dev/null +++ b/tests/assertions/regrid.027.core.test @@ -0,0 +1,9 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=0.75/0.75 --basic-angle=fraction +# grib_get +latitudeOfFirstGridPoint=120 +iDirectionIncrement=1 +basicAngleOfTheInitialProductionDomain=3 +subdivisionsOfBasicAngle=4 diff --git a/tests/assertions/regrid.028.test b/tests/assertions/regrid.028.test new file mode 100644 index 000000000..b01025c5b --- /dev/null +++ b/tests/assertions/regrid.028.test @@ -0,0 +1,7 @@ +# input +param=u_v,level=1000,resol=20 +# mir +--grid=3/3 --uv2uv +# grib_get +gridType=regular_ll +shortName=u diff --git a/tests/assertions/regular_ll.5-5.area.grib1 b/tests/assertions/regular_ll.5-5.area.grib1 new file mode 120000 index 000000000..b1695f091 --- /dev/null +++ b/tests/assertions/regular_ll.5-5.area.grib1 @@ -0,0 +1 @@ +../data/regular_ll.5-5.area.grib1 \ No newline at end of file diff --git a/tests/data/ICON.griddef b/tests/data/ICON.griddef new file mode 100644 index 000000000..ed88e80ae Binary files /dev/null and b/tests/data/ICON.griddef differ diff --git a/tests/data/gridType=polar_stereographic.grib1 b/tests/data/gridType=polar_stereographic.grib1 new file mode 100644 index 000000000..802dd9f4f Binary files /dev/null and b/tests/data/gridType=polar_stereographic.grib1 differ diff --git a/tests/data/gridType=polar_stereographic.grib2 b/tests/data/gridType=polar_stereographic.grib2 new file mode 100644 index 000000000..73c9f77e3 Binary files /dev/null and b/tests/data/gridType=polar_stereographic.grib2 differ diff --git a/tests/data/regular_ll.5-5.area.grib1 b/tests/data/regular_ll.5-5.area.grib1 new file mode 100644 index 000000000..94aec53da Binary files /dev/null and b/tests/data/regular_ll.5-5.area.grib1 differ diff --git a/tests/plans/CMakeLists.txt b/tests/plans/CMakeLists.txt index 9cbd48a9d..295969061 100644 --- a/tests/plans/CMakeLists.txt +++ b/tests/plans/CMakeLists.txt @@ -3,15 +3,24 @@ find_program(MARS_SCRIPT NAMES "mars" "false") ecbuild_configure_file(mir-test.sh.in mir-test.sh @ONLY) -file(GLOB_RECURSE test_files LIST_DIRECTORIES false *.test) -foreach(_t ${test_files}) - ecbuild_add_test( - TARGET ${_t} - CONDITION MARS_SCRIPT - COMMAND mir-test.sh - ARGS ${_t} - ENVIRONMENT ${_testEnvironment}) -endforeach() +if(MARS_SCRIPT) + file(GLOB_RECURSE test_files LIST_DIRECTORIES false *.test) + + foreach(_t ${test_files}) + ecbuild_add_test( + TARGET ${_t} + CONDITION MARS_SCRIPT + COMMAND mir-test.sh + ARGS ${_t} + ENVIRONMENT ${_testEnvironment}) + if((_t MATCHES "[.]fail[.]") OR + (_t MATCHES "[.]lsm[.]" AND NOT HAVE_MIR_DOWNLOAD_MASKS) OR + (NOT mir_HAVE_ATLAS AND NOT (_t MATCHES "[.]core[.]")) OR + (NOT atlas_HAVE_PROJ AND (_t MATCHES "[.]proj[.]"))) + set_tests_properties(${_t} PROPERTIES WILL_FAIL TRUE) + endif() + endforeach() +endif() if(EXISTS ${CMAKE_BINARY_DIR}/share/plugins/atlas-orca.yml) file(GLOB_RECURSE test_files LIST_DIRECTORIES false *.test.atlas-orca) diff --git a/tests/plans/ICON.griddef b/tests/plans/ICON.griddef new file mode 120000 index 000000000..817c2a817 --- /dev/null +++ b/tests/plans/ICON.griddef @@ -0,0 +1 @@ +../data/ICON.griddef \ No newline at end of file diff --git a/tests/plans/ICON.shortName=lsm.grib2 b/tests/plans/ICON.shortName=lsm.grib2 new file mode 120000 index 000000000..dc3cb1c9b --- /dev/null +++ b/tests/plans/ICON.shortName=lsm.grib2 @@ -0,0 +1 @@ +../data/ICON.shortName=lsm.grib2 \ No newline at end of file diff --git a/tests/plans/MIR-219-gridded-uv-from-spectral-uv-not-wind.test b/tests/plans/MIR-219-gridded-uv-from-spectral-uv-not-wind.test index 16ca9ee58..1ec9e7791 100644 --- a/tests/plans/MIR-219-gridded-uv-from-spectral-uv-not-wind.test +++ b/tests/plans/MIR-219-gridded-uv-from-spectral-uv-not-wind.test @@ -3,5 +3,5 @@ param=u_v,level=1000,resol=20 # mir --grid=1/1 --truncation=19 --intgrid=none # plan -ShTruncate[truncation=19]|ShToRegularLL[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,bbox=BoundingBox[north=90,west=0,south=-90,east=360],increments=Increments[west_east=1,south_north=1]]|AdjustWindsScaleCosLatitude[]|Save[packing=simple,output=...] -ShTruncate[truncation=19]|ShToRegularLL[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,bbox=BoundingBox[north=90,west=0,south=-90,east=360],increments=Increments[west_east=1,south_north=1]]|AdjustWindsScaleCosLatitude[]|Save[packing=simple,output=...] +ShTruncate[truncation=19]|ShToRegularLL[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,bbox=BoundingBox[north=90,west=0,south=-90,east=360],increments=Increments[west_east=1,south_north=1]]|AdjustWindsScaleCosLatitude[]|AdjustWindsAtPoles[]|Save[packing=simple,output=...] +ShTruncate[truncation=19]|ShToRegularLL[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,bbox=BoundingBox[north=90,west=0,south=-90,east=360],increments=Increments[west_east=1,south_north=1]]|AdjustWindsScaleCosLatitude[]|AdjustWindsAtPoles[]|Save[packing=simple,output=...] diff --git a/tests/plans/MIR-219-gridded-uv-from-spectral-uv.test b/tests/plans/MIR-219-gridded-uv-from-spectral-uv.test index 516841418..52ab2e8af 100644 --- a/tests/plans/MIR-219-gridded-uv-from-spectral-uv.test +++ b/tests/plans/MIR-219-gridded-uv-from-spectral-uv.test @@ -3,4 +3,4 @@ param=u_v,level=1000,resol=20 # mir --grid=1/1 --truncation=19 --intgrid=none --uv2uv # plan -ShTruncate[truncation=19]|ShToRegularLL[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,bbox=BoundingBox[north=90,west=0,south=-90,east=360],increments=Increments[west_east=1,south_north=1]]|AdjustWindsScaleCosLatitude[]|Save[packing=simple,output=...] +ShTruncate[truncation=19]|ShToRegularLL[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,bbox=BoundingBox[north=90,west=0,south=-90,east=360],increments=Increments[west_east=1,south_north=1]]|AdjustWindsScaleCosLatitude[]|AdjustWindsAtPoles[]|Save[packing=simple,output=...] diff --git a/tests/plans/MIR-268-u-only.test b/tests/plans/MIR-268-u-only.test index 626494a08..bc53c7c94 100644 --- a/tests/plans/MIR-268-u-only.test +++ b/tests/plans/MIR-268-u-only.test @@ -3,4 +3,4 @@ param=u_v,level=1000,resol=20 # mir --truncation=none --intgrid=none --uv2uv --grid=1/1 --rotation=-89/1 --u-only # plan -ShToRotatedRegularLL[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,bbox=BoundingBox[north=90,west=0,south=-90,east=360],increments=Increments[west_east=1,south_north=1],rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0]]|AdjustWindsScaleCosLatitude[]|AdjustWindsDirections[rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0]]|SelectField[0]|Save[packing=simple,output=...] +ShToRotatedRegularLL[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,bbox=BoundingBox[north=90,west=0,south=-90,east=360],increments=Increments[west_east=1,south_north=1],rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0]]|AdjustWindsScaleCosLatitude[]|AdjustWindsAtPoles[]|AdjustWindsDirections[rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0]]|SelectField[0]|Save[packing=simple,output=...] diff --git a/tests/plans/MIR-268-v-only.test b/tests/plans/MIR-268-v-only.test index f17108167..8df70483d 100644 --- a/tests/plans/MIR-268-v-only.test +++ b/tests/plans/MIR-268-v-only.test @@ -3,4 +3,4 @@ param=u_v,level=1000,resol=20 # mir --truncation=none --intgrid=none --uv2uv --grid=1/1 --rotation=-89/1 --v-only # plan -ShToRotatedRegularLL[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,bbox=BoundingBox[north=90,west=0,south=-90,east=360],increments=Increments[west_east=1,south_north=1],rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0]]|AdjustWindsScaleCosLatitude[]|AdjustWindsDirections[rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0]]|SelectField[1]|Save[packing=simple,output=...] +ShToRotatedRegularLL[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,bbox=BoundingBox[north=90,west=0,south=-90,east=360],increments=Increments[west_east=1,south_north=1],rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0]]|AdjustWindsScaleCosLatitude[]|AdjustWindsAtPoles[]|AdjustWindsDirections[rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0]]|SelectField[1]|Save[packing=simple,output=...] diff --git a/tests/plans/MIR-294.uv2uv=false.test b/tests/plans/MIR-294.uv2uv=false.test index e47d73a67..ec8bf9490 100644 --- a/tests/plans/MIR-294.uv2uv=false.test +++ b/tests/plans/MIR-294.uv2uv=false.test @@ -3,5 +3,5 @@ param=u_v,level=1000,resol=20 # mir --grid=1/1 --rotation=-89/1 --uv2uv=false # plan -ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|Gridded2RotatedLL[increments=Increments[west_east=1,south_north=1],bbox=BoundingBox[north=90,west=0,south=-90,east=359],rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06,validateMesh=0,projectionFail=missing-value]]|AdjustWindsScaleCosLatitude[]|AdjustWindsDirections[rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0]]|Save[packing=simple,output=...] -ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|Gridded2RotatedLL[increments=Increments[west_east=1,south_north=1],bbox=BoundingBox[north=90,west=0,south=-90,east=359],rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06,validateMesh=0,projectionFail=missing-value]]|AdjustWindsScaleCosLatitude[]|AdjustWindsDirections[rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0]]|Save[packing=simple,output=...] +ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|AdjustWindsScaleCosLatitude[]|Gridded2RotatedLL[increments=Increments[west_east=1,south_north=1],bbox=BoundingBox[north=90,west=0,south=-90,east=359],rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06,validateMesh=0,projectionFail=missing-value]]|AdjustWindsAtPoles[]|AdjustWindsDirections[rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0]]|Save[packing=simple,output=...] +ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|AdjustWindsScaleCosLatitude[]|Gridded2RotatedLL[increments=Increments[west_east=1,south_north=1],bbox=BoundingBox[north=90,west=0,south=-90,east=359],rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06,validateMesh=0,projectionFail=missing-value]]|AdjustWindsAtPoles[]|AdjustWindsDirections[rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0]]|Save[packing=simple,output=...] diff --git a/tests/plans/MIR-294.uv2uv=true.test b/tests/plans/MIR-294.uv2uv=true.test index a8bfabcc5..ef56bc774 100644 --- a/tests/plans/MIR-294.uv2uv=true.test +++ b/tests/plans/MIR-294.uv2uv=true.test @@ -3,4 +3,4 @@ param=u_v,level=1000,resol=20 # mir --grid=1/1 --rotation=-89/1 --uv2uv=true # plan -ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|Gridded2RotatedLL[increments=Increments[west_east=1,south_north=1],bbox=BoundingBox[north=90,west=0,south=-90,east=359],rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06,validateMesh=0,projectionFail=missing-value]]|AdjustWindsScaleCosLatitude[]|AdjustWindsDirections[rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0]]|Save[packing=simple,output=...] +ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|AdjustWindsScaleCosLatitude[]|Gridded2RotatedLL[increments=Increments[west_east=1,south_north=1],bbox=BoundingBox[north=90,west=0,south=-90,east=359],rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06,validateMesh=0,projectionFail=missing-value]]|AdjustWindsAtPoles[]|AdjustWindsDirections[rotation=Rotation[south_pole_latitude=-89,south_pole_longitude=1,south_pole_rotation_angle=0]]|Save[packing=simple,output=...] diff --git a/tests/plans/MIR-592.003.test b/tests/plans/MIR-592.003.test index 459e446ba..358d9796a 100644 --- a/tests/plans/MIR-592.003.test +++ b/tests/plans/MIR-592.003.test @@ -1,7 +1,6 @@ # input param=u_v,level=1000,resol=20 # mir ---grid=3/3 --interpolation=linear +--grid=3/3 --interpolation=linear --uv2uv # plan -ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|Gridded2RegularLL[increments=Increments[west_east=3,south_north=3],bbox=BoundingBox[north=90,west=0,south=-90,east=357],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06,validateMesh=0,projectionFail=missing-value]]|AdjustWindsScaleCosLatitude[]|Save[packing=simple,output=...] -ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|Gridded2RegularLL[increments=Increments[west_east=3,south_north=3],bbox=BoundingBox[north=90,west=0,south=-90,east=357],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06,validateMesh=0,projectionFail=missing-value]]|AdjustWindsScaleCosLatitude[]|Save[packing=simple,output=...] +ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|AdjustWindsScaleCosLatitude[]|Gridded2RegularLL[increments=Increments[west_east=3,south_north=3],bbox=BoundingBox[north=90,west=0,south=-90,east=357],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06,validateMesh=0,projectionFail=missing-value]]|AdjustWindsAtPoles[]|Save[packing=simple,output=...] diff --git a/tests/plans/MIR-673.001.test b/tests/plans/MIR-673.001.test new file mode 100644 index 000000000..abdf907a1 --- /dev/null +++ b/tests/plans/MIR-673.001.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--grid=gridType=polar_stereographic +# plan +Copy[] diff --git a/tests/plans/MIR-673.002.test b/tests/plans/MIR-673.002.test new file mode 100644 index 000000000..dce943e52 --- /dev/null +++ b/tests/plans/MIR-673.002.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--area=66.5/-24.5/63.4/-13.5 +# plan +AreaMasker[bbox=BoundingBox[north=66.5,west=-24.5,south=63.4,east=-13.5]]|Save[output=...] diff --git a/tests/plans/MIR-673.003.test b/tests/plans/MIR-673.003.test new file mode 100644 index 000000000..a8f82df91 --- /dev/null +++ b/tests/plans/MIR-673.003.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--area=66.5/-24.5/63.4/-13.5 --area-mode=crop +# plan +AreaCropper[bbox=BoundingBox[north=66.5,west=-24.5,south=63.4,east=-13.5]]|Save[output=...] diff --git a/tests/plans/MIR-673.004.test b/tests/plans/MIR-673.004.test new file mode 100644 index 000000000..7afeb35fd --- /dev/null +++ b/tests/plans/MIR-673.004.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--grid=1/1 +# plan +Gridded2RegularLL[increments=Increments[west_east=1,south_north=1],bbox=BoundingBox[north=90,west=0,south=-90,east=359],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[output=...] diff --git a/tests/plans/MIR-673.005.test b/tests/plans/MIR-673.005.test new file mode 100644 index 000000000..6be0f2dda --- /dev/null +++ b/tests/plans/MIR-673.005.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib1 +# mir +--grid=gridType=polar_stereographic +# plan +Copy[] diff --git a/tests/plans/MIR-673.006.test b/tests/plans/MIR-673.006.test new file mode 100644 index 000000000..585764dcd --- /dev/null +++ b/tests/plans/MIR-673.006.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib1 +# mir +--area=49/-92/41/-76 +# plan +AreaMasker[bbox=BoundingBox[north=49,west=-92,south=41,east=-76]]|Save[output=...] diff --git a/tests/plans/MIR-673.007.test b/tests/plans/MIR-673.007.test new file mode 100644 index 000000000..389ff6deb --- /dev/null +++ b/tests/plans/MIR-673.007.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib1 +# mir +--area=49/-92/41/-76 --area-mode=crop +# plan +AreaCropper[bbox=BoundingBox[north=49,west=-92,south=41,east=-76]]|Save[output=...] diff --git a/tests/plans/MIR-673.008.test b/tests/plans/MIR-673.008.test new file mode 100644 index 000000000..fa52e1420 --- /dev/null +++ b/tests/plans/MIR-673.008.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib1 +# mir +--grid=1/1 +# plan +Gridded2RegularLL[increments=Increments[west_east=1,south_north=1],bbox=BoundingBox[north=90,west=0,south=-90,east=359],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[output=...] diff --git a/tests/plans/MIR-673.009.proj.test b/tests/plans/MIR-673.009.proj.test new file mode 100644 index 000000000..3d4fe05be --- /dev/null +++ b/tests/plans/MIR-673.009.proj.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib2 +# mir +--same=data.in.MIR-673.009.proj.test +# plan +Gridded2TypedGrid[grid={"LaDInDegrees":"90","Ni":"2869","Nj":"2869","grid":"2500/2500","gridType":"polar_stereographic","iScansNegatively":"0","jScansPositively":"1","latitudeOfFirstGridPointInDegrees":"42.1138","longitudeOfFirstGridPointInDegrees":"-71.6168","orientationOfTheGridInDegrees":"-30","proj":"+proj=stere +lat_ts=90.000003 +lat_0=90 +lon_0=-30.000001 +k_0=1 +x_0=0 +y_0=0 +R=6371229.000000","shapeOfTheEarth":"6","southPoleOnProjectionPlane":"0"},interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[output=...] diff --git a/tests/plans/MIR-673.010.proj.test b/tests/plans/MIR-673.010.proj.test new file mode 100644 index 000000000..a31427bcd --- /dev/null +++ b/tests/plans/MIR-673.010.proj.test @@ -0,0 +1,6 @@ +# input +gridType=polar_stereographic.grib1 +# mir +--same=data.in.MIR-673.010.proj.test +# plan +Gridded2TypedGrid[grid={"LaDInDegrees":"60","Ni":"1121","Nj":"881","grid":"4763/4763","gridType":"polar_stereographic","iScansNegatively":"0","jScansPositively":"1","latitudeOfFirstGridPointInDegrees":"23.117","longitudeOfFirstGridPointInDegrees":"-119.023","orientationOfTheGridInDegrees":"-105","proj":"+proj=stere +lat_ts=60.000000 +lat_0=90 +lon_0=-105.000000 +k_0=1 +x_0=0 +y_0=0 +R=6367470.000000","shapeOfTheEarth":"0","southPoleOnProjectionPlane":"0"},interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[output=...] diff --git a/tests/plans/MIR-679.ICON-to-WP.test b/tests/plans/MIR-679.ICON-to-WP.test new file mode 100644 index 000000000..db24a8105 --- /dev/null +++ b/tests/plans/MIR-679.ICON-to-WP.test @@ -0,0 +1,6 @@ +# input +ICON.shortName=lsm.grib2 +# mir +--input={griddef:$tests_plans_dir/ICON.griddef} --griddef=$tests_plans_dir/weather-params.pts --format=grib --interpolation=nn +# plan +Gridded2GridDef[griddef=.../weather-params.pts]|Save[output=...] diff --git a/tests/plans/MIR-679.WP-to-ICON.test b/tests/plans/MIR-679.WP-to-ICON.test new file mode 100644 index 000000000..113f1278b --- /dev/null +++ b/tests/plans/MIR-679.WP-to-ICON.test @@ -0,0 +1,6 @@ +# input +/dev/null +# mir +--input={artificialInput:constant,constant:0.,gridded:True,gridType:unstructured_grid,griddef:$tests_plans_dir/weather-params.pts} --add-random=uniform-int-distribution{a:1,b:3} --griddef=$tests_plans_dir/ICON.griddef --format=grib --interpolation=nn +# plan +Gridded2GridDef[griddef=.../ICON.griddef]|AddRandomFilter[distribution=Distribution[uniform-int-distribution{a:1,b:3},min=1,max=3]]|Save[packing=ccsds,edition=2,output=...] diff --git a/tests/plans/MIR-681.001.test b/tests/plans/MIR-681.001.test new file mode 100644 index 000000000..7e6abaf6f --- /dev/null +++ b/tests/plans/MIR-681.001.test @@ -0,0 +1,7 @@ +# input +param=u_v,level=1000,resol=20 +# mir +--grid=3/3 --interpolation=grid-box-average +# plan +ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|AdjustWindsScaleCosLatitude[]|Gridded2RegularLL[increments=Increments[west_east=3,south_north=3],bbox=BoundingBox[north=90,west=0,south=-90,east=357],interpolation=grid-box-average,method=GridBoxMethod[name=grid-box-average,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06]]|AdjustWindsAtPoles[]|Save[packing=simple,output=...] +ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|AdjustWindsScaleCosLatitude[]|Gridded2RegularLL[increments=Increments[west_east=3,south_north=3],bbox=BoundingBox[north=90,west=0,south=-90,east=357],interpolation=grid-box-average,method=GridBoxMethod[name=grid-box-average,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06]]|AdjustWindsAtPoles[]|Save[packing=simple,output=...] diff --git a/tests/plans/MIR-681.002.test b/tests/plans/MIR-681.002.test new file mode 100644 index 000000000..7da5a2052 --- /dev/null +++ b/tests/plans/MIR-681.002.test @@ -0,0 +1,7 @@ +# input +param=u_v,level=1000,resol=20 +# mir +--grid=3/3 --interpolation=grid-box-statistics --interpolation-statistics=maximum +# plan +ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|AdjustWindsScaleCosLatitude[]|Gridded2RegularLL[increments=Increments[west_east=3,south_north=3],bbox=BoundingBox[north=90,west=0,south=-90,east=357],interpolation=grid-box-statistics,method=GridBoxMethod[name=grid-box-statistics,nonLinear[MissingIfHeaviestMissing[]],solver=Statistics[stats=Maximum[nan]],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06]]|AdjustWindsAtPoles[]|Save[packing=simple,output=...] +ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|AdjustWindsScaleCosLatitude[]|Gridded2RegularLL[increments=Increments[west_east=3,south_north=3],bbox=BoundingBox[north=90,west=0,south=-90,east=357],interpolation=grid-box-statistics,method=GridBoxMethod[name=grid-box-statistics,nonLinear[MissingIfHeaviestMissing[]],solver=Statistics[stats=Maximum[nan]],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06]]|AdjustWindsAtPoles[]|Save[packing=simple,output=...] diff --git a/tests/plans/gridType=polar_stereographic.grib1 b/tests/plans/gridType=polar_stereographic.grib1 new file mode 120000 index 000000000..f5a615d69 --- /dev/null +++ b/tests/plans/gridType=polar_stereographic.grib1 @@ -0,0 +1 @@ +../data/gridType=polar_stereographic.grib1 \ No newline at end of file diff --git a/tests/plans/gridType=polar_stereographic.grib2 b/tests/plans/gridType=polar_stereographic.grib2 new file mode 120000 index 000000000..6a4ca1f01 --- /dev/null +++ b/tests/plans/gridType=polar_stereographic.grib2 @@ -0,0 +1 @@ +../data/gridType=polar_stereographic.grib2 \ No newline at end of file diff --git a/tests/plans/intermediate-interpolation=nn.test b/tests/plans/intermediate-interpolation=nn.test new file mode 100644 index 000000000..e77104861 --- /dev/null +++ b/tests/plans/intermediate-interpolation=nn.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--intermediate-interpolation=nn --intgrid=O80 --interpolation=linear --grid=1/1 +# plan +Gridded2NamedGrid[grid=O80,interpolation=nn,method=KNearestNeighbours[nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,nearestMethod=NearestNeighbourWithLowestIndex[nclosest=4],distanceWeighting=InverseDistanceWeightingSquared[]]]|Gridded2RegularLL[increments=Increments[west_east=1,south_north=1],bbox=BoundingBox[north=90,west=0,south=-90,east=359],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[output=...] diff --git a/tests/plans/mir-test.sh.in b/tests/plans/mir-test.sh.in index 180c29929..055a1e5f9 100755 --- a/tests/plans/mir-test.sh.in +++ b/tests/plans/mir-test.sh.in @@ -40,17 +40,23 @@ do case $state in mars) - m=$(echo $n | tr '/' '_') - if [[ -f "${tests_plans_dir}/$m" ]] + m=$(echo $n) + if [[ -a "$m" ]] then - ln -sf "${tests_plans_dir}/$m" data.in.$t - elif [[ -f "${tmp_plans_dir}/$m" ]] - then - ln -sf "${tmp_plans_dir}/$m" data.in.$t + ln -sf "$m" data.in.$t else - echo "retrieve,$n,target=data.in.$t" > mars.r.$t - $mars_script mars.r.$t || rm data.in.$t - cp data.in.$t "${tmp_plans_dir}/$m" + m=$(echo $n | tr '/' '_') + if [[ -f "${tests_plans_dir}/$m" ]] + then + ln -sf "${tests_plans_dir}/$m" data.in.$t + elif [[ -f "${tmp_plans_dir}/$m" ]] + then + ln -sf "${tmp_plans_dir}/$m" data.in.$t + else + echo "retrieve,$n,target=data.in.$t" > mars.r.$t + $mars_script mars.r.$t || rm data.in.$t + cp data.in.$t "${tmp_plans_dir}/$m" + fi fi state=mir ;; diff --git a/tests/plans/regrid.001.test b/tests/plans/regrid.001.test new file mode 100644 index 000000000..5e3bbca9a --- /dev/null +++ b/tests/plans/regrid.001.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=1/1 +# plan +Gridded2RegularLL[increments=Increments[west_east=1,south_north=1],bbox=BoundingBox[north=90,west=0,south=-90,east=359],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[output=...] diff --git a/tests/plans/regrid.002.test b/tests/plans/regrid.002.test new file mode 100644 index 000000000..05572ae6b --- /dev/null +++ b/tests/plans/regrid.002.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=gridType=lambert_azimuthal_equal_area;latitudeOfFirstGridPointInDegrees=66.982143;longitudeOfFirstGridPointInDegrees=-35.034024;grid=5000./5000.;Ni=1000;Nj=950;standardParallelInDegrees=52.;centralLongitudeInDegrees=10. +# plan +Gridded2TypedGrid[grid={"Ni":"1000","Nj":"950","centralLongitudeInDegrees":"10.","grid":"5000./5000.","gridType":"lambert_azimuthal_equal_area","latitudeOfFirstGridPointInDegrees":"66.982143","longitudeOfFirstGridPointInDegrees":"-35.034024","standardParallelInDegrees":"52."},interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[output=...] diff --git a/tests/plans/regrid.003.test b/tests/plans/regrid.003.test new file mode 100644 index 000000000..ba411c4c3 --- /dev/null +++ b/tests/plans/regrid.003.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=gridType=lambert;latitudeOfFirstGridPointInDegrees=52.041;longitudeOfFirstGridPointInDegrees=1.639;grid=2500./2500.;Ni=739;Nj=949;LaDInDegrees=63.;LoVInDegrees=15. +# plan +Gridded2TypedGrid[grid={"LaDInDegrees":"63.","LoVInDegrees":"15.","Ni":"739","Nj":"949","grid":"2500./2500.","gridType":"lambert","latitudeOfFirstGridPointInDegrees":"52.041","longitudeOfFirstGridPointInDegrees":"1.639"},interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[output=...] diff --git a/tests/plans/regrid.004.test b/tests/plans/regrid.004.test new file mode 100644 index 000000000..a1b441f09 --- /dev/null +++ b/tests/plans/regrid.004.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=O32 --interpolation=nearest-neighbour --area=20/-40/-20/40 +# plan +Gridded2NamedGrid[grid=O32,interpolation=nearest-neighbour,method=KNearestNeighbours[nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=BoundingBox[north=20,west=-40,south=-20,east=40],lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,nearestMethod=NearestNeighbourWithLowestIndex[nclosest=4],distanceWeighting=InverseDistanceWeightingSquared[]]]|Save[output=...] diff --git a/tests/plans/regrid.005.test b/tests/plans/regrid.005.test new file mode 100644 index 000000000..b0f2ea6e1 --- /dev/null +++ b/tests/plans/regrid.005.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=O64 --interpolation=grid-box-average +# plan +Gridded2NamedGrid[grid=O64,interpolation=grid-box-average,method=GridBoxMethod[name=grid-box-average,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0]]|Save[output=...] diff --git a/tests/plans/regrid.006.test b/tests/plans/regrid.006.test new file mode 100644 index 000000000..c3da8eb7b --- /dev/null +++ b/tests/plans/regrid.006.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=1/1 --interpolation=grid-box-statistics --interpolation-statistics=maximum +# plan +Gridded2RegularLL[increments=Increments[west_east=1,south_north=1],bbox=BoundingBox[north=90,west=0,south=-90,east=359],interpolation=grid-box-statistics,method=GridBoxMethod[name=grid-box-statistics,nonLinear[MissingIfHeaviestMissing[]],solver=Statistics[stats=Maximum[nan]],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0]]|Save[output=...] diff --git a/tests/plans/regrid.007.test b/tests/plans/regrid.007.test new file mode 100644 index 000000000..cbb371a24 --- /dev/null +++ b/tests/plans/regrid.007.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=N160 --interpolation=k-nearest-neighbours --nearest-method=nclosest-or-nearest --nclosest=10 --distance=100000. --distance-weighting=gaussian --distance-weighting-gaussian-stddev=10000. --non-linear=missing-if-all-missing +# plan +Gridded2NamedGrid[grid=N160,interpolation=k-nearest-neighbours,method=KNearestNeighbours[nonLinear[MissingIfAllMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,nearestMethod=NClosestOrNearest[nclosest=10,distanceTolerance=1],distanceWeighting=GaussianDistanceWeighting[stddev=10000]]]|Save[output=...] diff --git a/tests/plans/regrid.009.test b/tests/plans/regrid.009.test new file mode 100644 index 000000000..4475c1b2e --- /dev/null +++ b/tests/plans/regrid.009.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=O320 --interpolation=structured-bilinear +# plan +Gridded2NamedGrid[grid=O320,interpolation=structured-bilinear,method=ProxyMethod[options=LocalConfiguration[root={type => structured-bilinear , matrix_free => true}],cropping=none]]|Save[output=...] diff --git a/tests/plans/regrid.010.test b/tests/plans/regrid.010.test new file mode 100644 index 000000000..37ae384b7 --- /dev/null +++ b/tests/plans/regrid.010.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=O160 --interpolation=structured-bicubic +# plan +Gridded2NamedGrid[grid=O160,interpolation=structured-bicubic,method=ProxyMethod[options=LocalConfiguration[root={type => structured-bicubic , matrix_free => true}],cropping=none]]|Save[output=...] diff --git a/tests/plans/regrid.012.test b/tests/plans/regrid.012.test new file mode 100644 index 000000000..fa52d9bdd --- /dev/null +++ b/tests/plans/regrid.012.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--nabla=scalar-gradient +# plan +ScalarGradient[meshGeneratorParameters=...]|Save[output=...] diff --git a/tests/plans/regrid.013.test b/tests/plans/regrid.013.test new file mode 100644 index 000000000..e1472a457 --- /dev/null +++ b/tests/plans/regrid.013.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--nabla=scalar-laplacian +# plan +ScalarLaplacian[meshGeneratorParameters=...]|Save[output=...] diff --git a/tests/plans/regrid.014.test b/tests/plans/regrid.014.test new file mode 100644 index 000000000..bc5b05f8f --- /dev/null +++ b/tests/plans/regrid.014.test @@ -0,0 +1,6 @@ +# input +param=u_v,level=1000,resol=20 +# mir +--nabla=uv-gradient --grid=3/3 --uv2uv +# plan +ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|AdjustWindsScaleCosLatitude[]|Gridded2RegularLL[increments=Increments[west_east=3,south_north=3],bbox=BoundingBox[north=90,west=0,south=-90,east=357],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06,validateMesh=0,projectionFail=missing-value]]|AdjustWindsAtPoles[]|UVGradient[meshGeneratorParameters=...]|Save[packing=simple,output=...] diff --git a/tests/plans/regrid.015.test b/tests/plans/regrid.015.test new file mode 100644 index 000000000..9d4d5fead --- /dev/null +++ b/tests/plans/regrid.015.test @@ -0,0 +1,6 @@ +# input +param=u_v,level=1000,resol=20 +# mir +--nabla=uv-vorticity --grid=3/3 --uv2uv +# plan +ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|AdjustWindsScaleCosLatitude[]|Gridded2RegularLL[increments=Increments[west_east=3,south_north=3],bbox=BoundingBox[north=90,west=0,south=-90,east=357],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06,validateMesh=0,projectionFail=missing-value]]|AdjustWindsAtPoles[]|UVVorticity[meshGeneratorParameters=...]|Save[packing=simple,output=...] diff --git a/tests/plans/regrid.016.test b/tests/plans/regrid.016.test new file mode 100644 index 000000000..22a6cb214 --- /dev/null +++ b/tests/plans/regrid.016.test @@ -0,0 +1,7 @@ +# input +param=u_v,level=1000,resol=20 +# mir +--nabla=uv-divergence --grid=3/3 --uv2uv +# plan +ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|AdjustWindsScaleCosLatitude[]|Gridded2RegularLL[increments=Increments[west_east=3,south_north=3],bbox=BoundingBox[north=90,west=0,south=-90,east=357],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06,validateMesh=0,projectionFail=missing-value]]|AdjustWindsAtPoles[]|UVDivergence[meshGeneratorParameters=...]|Save[packing=simple,output=...] + diff --git a/tests/plans/regrid.017.test b/tests/plans/regrid.017.test new file mode 100644 index 000000000..fdacf7652 --- /dev/null +++ b/tests/plans/regrid.017.test @@ -0,0 +1,6 @@ +# input +param=z,level=1,levtype=ml,resol=20.grib2 +# mir +--truncation=19 +# plan +ShTruncate[truncation=19]|Save[output=...] diff --git a/tests/plans/regrid.018.test b/tests/plans/regrid.018.test new file mode 100644 index 000000000..2a612987b --- /dev/null +++ b/tests/plans/regrid.018.test @@ -0,0 +1,6 @@ +# input +param=z,level=1,levtype=ml,resol=20.grib2 +# mir +--grid=1/1 --intgrid=O20 --spectral-order=cubic +# plan +ShTruncate[truncation=19]|ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=O20]|Gridded2RegularLL[increments=Increments[west_east=1,south_north=1],bbox=BoundingBox[north=90,west=0,south=-90,east=359],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[packing=simple,output=...] diff --git a/tests/plans/regrid.019.test b/tests/plans/regrid.019.test new file mode 100644 index 000000000..40f5443d3 --- /dev/null +++ b/tests/plans/regrid.019.test @@ -0,0 +1,6 @@ +# input +param=z,level=1,levtype=ml,resol=20.grib2 +# mir +--grid=1/1 --intgrid=O20 --truncation=19 +# plan +ShTruncate[truncation=19]|ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=O20]|Gridded2RegularLL[increments=Increments[west_east=1,south_north=1],bbox=BoundingBox[north=90,west=0,south=-90,east=359],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[packing=simple,output=...] diff --git a/tests/plans/regrid.020.test b/tests/plans/regrid.020.test new file mode 100644 index 000000000..d2c6382ad --- /dev/null +++ b/tests/plans/regrid.020.test @@ -0,0 +1,6 @@ +# input +param=vo_d,level=1,levtype=ml,resol=20.grib2 +# mir +--vod2uv +# plan +ShVodToUV[type=local]|Save[output=...] diff --git a/tests/plans/regrid.021.test b/tests/plans/regrid.021.test new file mode 100644 index 000000000..a81fb87c6 --- /dev/null +++ b/tests/plans/regrid.021.test @@ -0,0 +1,6 @@ +# input +param=vo_d,level=1,levtype=ml,resol=20.grib2 +# mir +--grid=3/3 --vod2uv +# plan +ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|Gridded2RegularLL[increments=Increments[west_east=3,south_north=3],bbox=BoundingBox[north=90,west=0,south=-90,east=357],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06,validateMesh=0,projectionFail=missing-value]]|Save[packing=simple,output=...] diff --git a/tests/plans/regrid.022.test b/tests/plans/regrid.022.test new file mode 100644 index 000000000..c00834734 --- /dev/null +++ b/tests/plans/regrid.022.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=5/5 +# plan +Gridded2RegularLL[increments=Increments[west_east=5,south_north=5],bbox=BoundingBox[north=90,west=0,south=-90,east=355],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[output=...] diff --git a/tests/plans/regrid.023.test b/tests/plans/regrid.023.test new file mode 100644 index 000000000..f05f9b6e5 --- /dev/null +++ b/tests/plans/regrid.023.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=5/5 --interpolation=grid-box-average +# plan +Gridded2RegularLL[increments=Increments[west_east=5,south_north=5],bbox=BoundingBox[north=90,west=0,south=-90,east=355],interpolation=grid-box-average,method=GridBoxMethod[name=grid-box-average,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0]]|Save[output=...] diff --git a/tests/plans/regrid.024.test b/tests/plans/regrid.024.test new file mode 100644 index 000000000..ea721c074 --- /dev/null +++ b/tests/plans/regrid.024.test @@ -0,0 +1,6 @@ +# input +regular_ll.5-5.area.grib1 +# mir +--grid=3/3 --globalise=1 +# plan +Gridded2RegularLL[increments=Increments[west_east=3,south_north=3],bbox=BoundingBox[north=90,west=0,south=-90,east=357],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|GlobaliseFilter[]|Save[output=...] diff --git a/tests/plans/regrid.025.test b/tests/plans/regrid.025.test new file mode 100644 index 000000000..4ce51b463 --- /dev/null +++ b/tests/plans/regrid.025.test @@ -0,0 +1,6 @@ +# input +regular_ll.5-5.area.grib1 +# mir +--grid=3/3 --pre-globalise=1 +# plan +GlobaliseFilter[]|Gridded2RegularLL[increments=Increments[west_east=3,south_north=3],bbox=BoundingBox[north=90,west=0,south=-90,east=357],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[output=...] diff --git a/tests/plans/regrid.026.test b/tests/plans/regrid.026.test new file mode 100644 index 000000000..e16221e15 --- /dev/null +++ b/tests/plans/regrid.026.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=0.75/0.75 --basic-angle=decimal +# plan +Gridded2RegularLL[increments=Increments[west_east=0.75,south_north=0.75],bbox=BoundingBox[north=90,west=0,south=-90,east=359.25],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[output=...] diff --git a/tests/plans/regrid.027.test b/tests/plans/regrid.027.test new file mode 100644 index 000000000..3415cce91 --- /dev/null +++ b/tests/plans/regrid.027.test @@ -0,0 +1,6 @@ +# input +regular_ll.2-2.grib2 +# mir +--grid=0.75/0.75 --basic-angle=fraction +# plan +Gridded2RegularLL[increments=Increments[west_east=0.75,south_north=0.75],bbox=BoundingBox[north=90,west=0,south=-90,east=359.25],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[output=...] diff --git a/tests/plans/regrid.028.test b/tests/plans/regrid.028.test new file mode 100644 index 000000000..ea558e4ad --- /dev/null +++ b/tests/plans/regrid.028.test @@ -0,0 +1,6 @@ +# input +param=u_v,level=1000,resol=20 +# mir +--grid=3/3 --uv2uv +# plan +ShToNamedGrid[type=local,cropping=none,options=[4886829aa871a5c38f41fff09ed0970e],invtrans=,grid=F21]|AdjustWindsScaleCosLatitude[]|Gridded2RegularLL[increments=Increments[west_east=3,south_north=3],bbox=BoundingBox[north=90,west=0,south=-90,east=357],interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=1e-06,validateMesh=0,projectionFail=missing-value]]|AdjustWindsAtPoles[]|Save[packing=simple,output=...] diff --git a/tests/plans/regular_ll.5-5.area.grib1 b/tests/plans/regular_ll.5-5.area.grib1 new file mode 120000 index 000000000..b1695f091 --- /dev/null +++ b/tests/plans/regular_ll.5-5.area.grib1 @@ -0,0 +1 @@ +../data/regular_ll.5-5.area.grib1 \ No newline at end of file diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 892b6ca23..44fa4bd96 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -13,6 +13,7 @@ foreach(_t action_graph area bounding_box + earthkit-regrid formula gaussian_grid grib_basic_angle @@ -30,6 +31,7 @@ foreach(_t statistics style vector-space + weight_matrix wind) ecbuild_add_test( TARGET mir_tests_unit_${_t} @@ -48,6 +50,15 @@ if(mir_HAVE_ATLAS) WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}) endif() +if(mir_HAVE_NETCDF) + ecbuild_add_test( + TARGET mir_tests_unit_interpolation_netcdf + SOURCES interpolation_netcdf.cc + LIBS mir + ENVIRONMENT ${_testEnvironment} + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}) +endif() + if(ENABLE_BUILD_TOOLS) foreach(_t MIR-637.grib2) ecbuild_add_test( diff --git a/tests/unit/bounding_box.cc b/tests/unit/bounding_box.cc index b26e5539c..8c35197e2 100644 --- a/tests/unit/bounding_box.cc +++ b/tests/unit/bounding_box.cc @@ -301,9 +301,9 @@ CASE("Representation::extendBoundingBox") { SECTION("Gaussian") { std::vector _name{"F16", "O16", "F21", "O21"}; -#if mir_HAVE_ATLAS - _name.emplace_back("N16"); -#endif + if constexpr (HAVE_ATLAS) { + _name.emplace_back("N16"); + } std::vector _bbox{{90, 10, -10, 9}, {90, -350, -10, 9}, {-70, -10, -90, 10}, {0, -350, 0, 9}, {60, -350, 20, 9}, {90, -10, 70, 10}}; diff --git a/tests/unit/earthkit-regrid.cc b/tests/unit/earthkit-regrid.cc new file mode 100644 index 000000000..ea1ec324d --- /dev/null +++ b/tests/unit/earthkit-regrid.cc @@ -0,0 +1,171 @@ +/* + * (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 +#include +#include +#include + +#include "eckit/filesystem/PathName.h" +#include "eckit/parser/YAMLParser.h" +#include "eckit/testing/Test.h" + +#include "mir/api/MIRJob.h" +#include "mir/input/RawInput.h" +#include "mir/output/EmptyOutput.h" +#include "mir/param/SimpleParametrisation.h" +#include "mir/util/Log.h" + + +namespace mir::tests::unit { + + +CASE("InterpolationSpec integration test") { + auto& log = Log::info(); + + + // temporary file for dumping weights info_path (should be local and unique) + eckit::PathName info_path{ + "mir_tests_unit" + "_" + "earthkit-regrid" + "_" + "dump-weights-info"}; + + + param::SimpleParametrisation meta; + meta.set("gridded", true); + meta.set("gridType", "regular_ll"); + meta.set("north", 10.); + meta.set("west", 0.); + meta.set("south", 0.); + meta.set("east", 10.); + meta.set("south_north_increment", 5.); + meta.set("west_east_increment", 5.); + meta.set("Ni", 3); + meta.set("Nj", 3); + + std::vector values(9, 0.); + + // extract "interpolation" node from the dumped info + auto info_interpolation = [](const auto& path) { + auto value = eckit::YAMLParser::decodeFile(path); + std::ostringstream ss; + ss << value["interpolation"]; + return ss.str(); + }; + + + static const std::string EXPECTED_INTERPOLATION_LINEAR{"{engine => mir , version => 16 , method => linear}"}; + static const std::string EXPECTED_INTERPOLATION_NN{"{engine => mir , version => 16 , method => nearest-neighbour}"}; + static const std::string EXPECTED_INTERPOLATION_GBA{"{engine => mir , version => 16 , method => grid-box-average}"}; + + + SECTION("interpolation=linear (default)") { + std::unique_ptr input(new input::RawInput(values.data(), values.size(), meta)); + std::unique_ptr output(new output::EmptyOutput); + + api::MIRJob job; + job.set("grid", std::vector{2., 2.}); + job.set("caching", false); + job.set("interpolation", "linear"); + job.set("dump-weights-info", info_path); + job.execute(*input, *output); + + EXPECT(info_interpolation(info_path) == EXPECTED_INTERPOLATION_LINEAR); + } + + + SECTION("interpolation=linear with non-default option(s)") { + std::unique_ptr input(new input::RawInput(values.data(), values.size(), meta)); + std::unique_ptr output(new output::EmptyOutput); + + api::MIRJob job; + job.set("grid", std::vector{2., 2.}); + job.set("caching", false); + job.set("interpolation", "linear"); + job.set("dump-weights-info", info_path); + job.set("finite-element-validate-mesh", true); // non-default option + job.execute(*input, *output); + + EXPECT(info_interpolation(info_path) != EXPECTED_INTERPOLATION_LINEAR); + } + + + SECTION("interpolation=nn (default)") { + std::unique_ptr input(new input::RawInput(values.data(), values.size(), meta)); + std::unique_ptr output(new output::EmptyOutput); + + api::MIRJob job; + job.set("grid", std::vector{2., 2.}); + job.set("caching", false); + job.set("interpolation", "nn"); + job.set("dump-weights-info", info_path); + job.execute(*input, *output); + + EXPECT(info_interpolation(info_path) == EXPECTED_INTERPOLATION_NN); + } + + + SECTION("interpolation=nn with non-default option(s)") { + std::unique_ptr input(new input::RawInput(values.data(), values.size(), meta)); + std::unique_ptr output(new output::EmptyOutput); + + api::MIRJob job; + job.set("grid", std::vector{2., 2.}); + job.set("caching", false); + job.set("interpolation", "nn"); + job.set("dump-weights-info", info_path); + job.set("nclosest", 2); // non-default option + job.execute(*input, *output); + + EXPECT(info_interpolation(info_path) != EXPECTED_INTERPOLATION_NN); + } + + + SECTION("interpolation=grid-box-average (default)") { + std::unique_ptr input(new input::RawInput(values.data(), values.size(), meta)); + std::unique_ptr output(new output::EmptyOutput); + + api::MIRJob job; + job.set("grid", std::vector{2., 2.}); + job.set("caching", false); + job.set("interpolation", "grid-box-average"); + job.set("dump-weights-info", info_path); + job.execute(*input, *output); + EXPECT(info_interpolation(info_path) == EXPECTED_INTERPOLATION_GBA); + } + + + SECTION("interpolation=grid-box-average with non-default option(s)") { + std::unique_ptr input(new input::RawInput(values.data(), values.size(), meta)); + std::unique_ptr output(new output::EmptyOutput); + + api::MIRJob job; + job.set("grid", std::vector{2., 2.}); + job.set("caching", false); + job.set("interpolation", "grid-box-average"); + job.set("dump-weights-info", info_path); + job.set("prune-epsilon", 1e-2); // non-default option + job.execute(*input, *output); + + EXPECT(info_interpolation(info_path) != EXPECTED_INTERPOLATION_GBA); + } +} + + +} // namespace mir::tests::unit + + +int main(int argc, char** argv) { + return eckit::testing::run_tests(argc, argv); +} diff --git a/tests/unit/gaussian_grid.cc b/tests/unit/gaussian_grid.cc index 15329896b..63cf0e03c 100644 --- a/tests/unit/gaussian_grid.cc +++ b/tests/unit/gaussian_grid.cc @@ -84,25 +84,25 @@ CASE("NamedGrid") { } } -#if mir_HAVE_ATLAS - SECTION("numberOfPoints (N)") { - std::vector tests{ - {"N128", 88838, 88838, {90, 0, -90, 359.297}}, - {"N320", 542080, 13437, {-60, 50, -90, 180}}, - {"N640", 2140702, 52654, {-60, 50, -90, 180}}, - }; - - for (const auto& t : tests) { - log << "Test " + t.grid + " (global)" << std::endl; - repres::RepresentationHandle global(key::grid::Grid::lookup(t.grid).representation()); - EXPECT(global->numberOfPoints() == t.numberOfPoints); - - log << "Test " + t.grid + " (cropped)" << std::endl; - repres::RepresentationHandle local(global->croppedRepresentation(t.bbox)); - EXPECT(local->numberOfPoints() == t.numberOfCroppedPoints); + if constexpr (HAVE_ATLAS) { + SECTION("numberOfPoints (N)") { + std::vector tests{ + {"N128", 88838, 88838, {90, 0, -90, 359.297}}, + {"N320", 542080, 13437, {-60, 50, -90, 180}}, + {"N640", 2140702, 52654, {-60, 50, -90, 180}}, + }; + + for (const auto& t : tests) { + log << "Test " + t.grid + " (global)" << std::endl; + repres::RepresentationHandle global(key::grid::Grid::lookup(t.grid).representation()); + EXPECT(global->numberOfPoints() == t.numberOfPoints); + + log << "Test " + t.grid + " (cropped)" << std::endl; + repres::RepresentationHandle local(global->croppedRepresentation(t.bbox)); + EXPECT(local->numberOfPoints() == t.numberOfCroppedPoints); + } } } -#endif log.precision(old); } diff --git a/tests/unit/grid_box_method.cc b/tests/unit/grid_box_method.cc index b5580554f..98c19a04a 100644 --- a/tests/unit/grid_box_method.cc +++ b/tests/unit/grid_box_method.cc @@ -12,6 +12,8 @@ #include "eckit/testing/Test.h" +#include "eckit/types/FloatCompare.h" + #include "mir/param/SimpleParametrisation.h" #include "mir/repres/gauss/reduced/ReducedFromPL.h" #include "mir/repres/gauss/regular/RegularGG.h" @@ -21,102 +23,169 @@ #include "mir/util/Domain.h" #include "mir/util/GridBox.h" #include "mir/util/Increments.h" +#include "mir/util/Log.h" #include "mir/util/Types.h" -// define EXPECTV(a) log << "\tEXPECT(" << #a <<")" << std::endl; EXPECT(a) + +constexpr double EPS = 1e-6; + +#define EXPECT_APPROX(a, b) \ + Log::info() << "EXPECT(" << #a << " ~= " << #b << ")" << std::endl; \ + EXPECT(eckit::types::is_approximately_equal(static_cast(a), static_cast(b), EPS)) + +#define EXPECT_BOX(box, n, w, s, e) \ + EXPECT_APPROX(box.north(), n); \ + EXPECT_APPROX(box.west(), w); \ + EXPECT_APPROX(box.south(), s); \ + EXPECT_APPROX(box.east(), e) namespace mir::tests::unit { -CASE("grid boxes: West-East periodicity") { +CASE("grid boxes: West-East periodicity, reduced Gaussian grids") { using Longitude = LongitudeDouble; - SECTION("reduced Gaussian grid (periodic)") { + SECTION("periodic") { const std::vector pl{20, 20}; repres::RepresentationHandle r(new repres::gauss::reduced::ReducedFromPL(pl.size() / 2, pl)); ASSERT(r->domain().isPeriodicWestEast()); - auto boxes = r->gridBoxes(); - auto a = boxes.front(); - auto b = boxes.back(); - auto inc = 360. / static_cast(pl[0]); + const auto boxes = r->gridBoxes(); + EXPECT(boxes.size() == r->numberOfPoints()); + EXPECT(boxes.size() == 40); - EXPECT_EQUAL(Longitude(a.west()), Longitude(-inc / 2.)); - EXPECT_EQUAL(Longitude(b.east()), Longitude(-inc / 2. + 360.)); + const auto inc = 360. / static_cast(pl[0]); + + EXPECT_EQUAL(Longitude(boxes.front().west()), Longitude(-inc / 2.)); + EXPECT_EQUAL(Longitude(boxes.back().east()), Longitude(-inc / 2. + 360.)); + + EXPECT_BOX(boxes[0], 90., -9., 0., 9.); + EXPECT_BOX(boxes[1], 90., 9., 0., 27.); + // ... + EXPECT_BOX(boxes[38], 0., 315., -90., 333.); + EXPECT_BOX(boxes[39], 0., 333., -90., 351.); } - SECTION("reduced Gaussian grid (non-periodic)") { + SECTION("non-periodic") { const std::vector pl{20, 20}; util::BoundingBox bbox{90., 0., 0., 180.}; repres::RepresentationHandle r(new repres::gauss::reduced::ReducedFromPL(pl.size() / 2, pl, bbox)); ASSERT(!r->domain().isPeriodicWestEast()); - for (const auto& box : r->gridBoxes()) { + const auto boxes = r->gridBoxes(); + EXPECT(boxes.size() == r->numberOfPoints()); + EXPECT(boxes.size() == 11); + + for (const auto& box : boxes) { EXPECT(bbox.contains({box.north(), box.west(), box.south(), box.east()})); } + + EXPECT_BOX(boxes[0], 90., 0., 0., 9.); + EXPECT_BOX(boxes[1], 90., 9., 0., 27.); + // ... + EXPECT_BOX(boxes[9], 90., 153., 0., 171.); + EXPECT_BOX(boxes[10], 90., 171., 0., 180.); } +} + +CASE("grid boxes: West-East periodicity, regular Gaussian grids") { - SECTION("regular Gaussian grid (periodic)") { + SECTION("periodic") { repres::RepresentationHandle r(new repres::gauss::regular::RegularGG(1)); ASSERT(r->domain().isPeriodicWestEast()); - auto boxes = r->gridBoxes(); - auto a = boxes.front(); - auto b = boxes.back(); - auto inc = 360. / 4.; // assumes Ni == 4 * N + const auto boxes = r->gridBoxes(); + EXPECT(boxes.size() == r->numberOfPoints()); + EXPECT(boxes.size() == 8); + + const auto inc = 360. / 4.; // assumes Ni == 4 * N + + EXPECT_EQUAL(Longitude(boxes.front().west()), Longitude(-inc / 2.)); + EXPECT_EQUAL(Longitude(boxes.back().east()), Longitude(-inc / 2. + 360.)); - EXPECT_EQUAL(Longitude(a.west()), Longitude(-inc / 2.)); - EXPECT_EQUAL(Longitude(b.east()), Longitude(-inc / 2. + 360.)); + EXPECT_BOX(boxes[0], 90., -45., 0., 45.); + EXPECT_BOX(boxes[1], 90., 45., 0., 135.); + // ... + EXPECT_BOX(boxes[6], 0., 135., -90., 225.); + EXPECT_BOX(boxes[7], 0., 225., -90., 315.); } - SECTION("regular Gaussian grid (non-periodic)") { + SECTION("non-periodic") { util::BoundingBox bbox{90., 0., 0., 180.}; repres::RepresentationHandle r(new repres::gauss::regular::RegularGG(1, bbox)); ASSERT(!r->domain().isPeriodicWestEast()); - for (const auto& box : r->gridBoxes()) { + const auto boxes = r->gridBoxes(); + EXPECT(boxes.size() == r->numberOfPoints()); + EXPECT(boxes.size() == 3); + + for (const auto& box : boxes) { EXPECT(bbox.contains({box.north(), box.west(), box.south(), box.east()})); } + + EXPECT_BOX(boxes[0], 90., 0., 35.26439, 45.); + EXPECT_BOX(boxes[1], 90., 45., 35.26439, 135.); + EXPECT_BOX(boxes[2], 90., 135., 35.26439, 180.); } +} - SECTION("regular lat/lon grid (periodic)") { +CASE("grid boxes: West-East periodicity, regular lat/lon grids") { + SECTION("periodic") { double inc = 18.; repres::RepresentationHandle r(new repres::latlon::RegularLL(util::Increments{inc, 90.})); ASSERT(r->domain().isPeriodicWestEast()); - auto boxes = r->gridBoxes(); - auto a = boxes.front(); - auto b = boxes.back(); + const auto boxes = r->gridBoxes(); + EXPECT(boxes.size() == r->numberOfPoints()); + EXPECT(boxes.size() == 60); + + EXPECT_EQUAL(Longitude(boxes.front().west()), Longitude(-inc / 2.)); + EXPECT_EQUAL(Longitude(boxes.back().east()), Longitude(-inc / 2. + 360.)); - EXPECT_EQUAL(Longitude(a.west()), Longitude(-inc / 2.)); - EXPECT_EQUAL(Longitude(b.east()), Longitude(-inc / 2. + 360.)); + EXPECT_BOX(boxes[0], 90., -9., 45., 9.); + EXPECT_BOX(boxes[1], 90., 9., 45., 27.); + // ... + EXPECT_BOX(boxes[58], -45., 315., -90., 333.); + EXPECT_BOX(boxes[59], -45., 333., -90., 351.); } - SECTION("regular lat/lon grid (non-periodic)") { + SECTION("non-periodic") { util::BoundingBox bbox{90., 0., 0., 180.}; double inc = 18.; repres::RepresentationHandle r(new repres::latlon::RegularLL(util::Increments{inc, 90.}, bbox)); ASSERT(!r->domain().isPeriodicWestEast()); - for (const auto& box : r->gridBoxes()) { + const auto boxes = r->gridBoxes(); + EXPECT(boxes.size() == r->numberOfPoints()); + EXPECT(boxes.size() == 22); + + for (const auto& box : boxes) { EXPECT(bbox.contains({box.north(), box.west(), box.south(), box.east()})); } + + EXPECT_BOX(boxes[0], 90., 0., 45., 9.); + EXPECT_BOX(boxes[1], 90., 9., 45., 27.); + // ... + EXPECT_BOX(boxes[20], 45., 153., 0., 171.); + EXPECT_BOX(boxes[21], 45., 171., 0., 180.); } +} - { +CASE("grid boxes: West-East periodicity, reduced lat/lon grids") { + SECTION("periodic") { const std::vector pl{20, 20}; param::SimpleParametrisation param; @@ -127,31 +196,54 @@ CASE("grid boxes: West-East periodicity") { param.set("south", 0.); param.set("east", 360.); + repres::RepresentationHandle r(new repres::latlon::ReducedLL(param)); + ASSERT(r->domain().isPeriodicWestEast()); + + const auto boxes = r->gridBoxes(); + EXPECT(boxes.size() == r->numberOfPoints()); + EXPECT(boxes.size() == 40); - SECTION("reduced lat/lon grid (periodic)") { - repres::RepresentationHandle r(new repres::latlon::ReducedLL(param)); - ASSERT(r->domain().isPeriodicWestEast()); + const auto inc = 360. / static_cast(pl[0]); - auto boxes = r->gridBoxes(); - auto a = boxes.front(); - auto b = boxes.back(); - auto inc = 360. / static_cast(pl[0]); + EXPECT_EQUAL(Longitude(boxes.front().west()), Longitude(-inc / 2.)); + EXPECT_EQUAL(Longitude(boxes.back().east()), Longitude(-inc / 2. + 360.)); + + EXPECT_BOX(boxes[0], 90., -9., 45., 9.); + EXPECT_BOX(boxes[1], 90., 9., 45., 27.); + // ... + EXPECT_BOX(boxes[38], 45., 315., 0., 333.); + EXPECT_BOX(boxes[39], 45., 333., 0., 351.); + } - EXPECT_EQUAL(Longitude(a.west()), Longitude(-inc / 2.)); - EXPECT_EQUAL(Longitude(b.east()), Longitude(-inc / 2. + 360.)); - } + SECTION("non-periodic") { + const std::vector pl{11, 11}; + const util::BoundingBox bbox{90., 0., 0., 180.}; - SECTION("reduced lat/lon grid (non-periodic)") { - util::BoundingBox bbox{90., 0., 0., 180.}; - param.set("east", bbox.east().value()); + param::SimpleParametrisation param; + param.set("pl", pl); + param.set("Nj", pl.size()); + param.set("north", bbox.north().value()); + param.set("west", bbox.west().value()); + param.set("south", bbox.south().value()); + param.set("east", bbox.east().value()); + + repres::RepresentationHandle r(new repres::latlon::ReducedLL(param)); + ASSERT(!r->domain().isPeriodicWestEast()); + + const auto boxes = r->gridBoxes(); + EXPECT(boxes.size() == r->numberOfPoints()); + EXPECT(boxes.size() == 22); - repres::RepresentationHandle r(new repres::latlon::ReducedLL(param)); - ASSERT(!r->domain().isPeriodicWestEast()); - for (const auto& box : r->gridBoxes()) { - EXPECT(bbox.contains({box.north(), box.west(), box.south(), box.east()})); - } + for (const auto& box : boxes) { + EXPECT(bbox.contains({box.north(), box.west(), box.south(), box.east()})); } + + EXPECT_BOX(boxes[0], 90., 0., 45., 9.); + EXPECT_BOX(boxes[1], 90., 9., 45., 27.); + // ... + EXPECT_BOX(boxes[20], 45., 153., 0., 171.); + EXPECT_BOX(boxes[21], 45., 171., 0., 180.); } } @@ -160,7 +252,7 @@ CASE("MIR-666") { // shifted grid, with finely-adjusted bounding box repres::RepresentationHandle r(new repres::latlon::RegularLL( util::Increments{0.05, 0.05}, util::BoundingBox{0.025, 27.875, -0.025, 57.825}, {-0.025, 27.875})); - auto boxes = r->gridBoxes(); + const auto boxes = r->gridBoxes(); } diff --git a/tests/unit/interpolation_netcdf.cc b/tests/unit/interpolation_netcdf.cc new file mode 100644 index 000000000..4a5cd233b --- /dev/null +++ b/tests/unit/interpolation_netcdf.cc @@ -0,0 +1,61 @@ +/* + * (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 + +#include "eckit/testing/Test.h" + +#include "mir/api/MIRJob.h" +#include "mir/api/mir_config.h" +#include "mir/input/MIRInput.h" +#include "mir/output/EmptyOutput.h" +#include "mir/param/SimpleParametrisation.h" + + +namespace mir::tests::unit { + + +CASE("interpolations") { + api::MIRJob jobs[3]; // jobs[0]: no post-processing + + jobs[0].set("caching", false); + jobs[0].set("grid", "3/3"); + + jobs[1].set("caching", false); + jobs[1].set("grid", "1/1"); + jobs[1].set("area", "40/20/20/40"); + jobs[1].set("frame", 2); + + if constexpr (HAVE_TESSELATION) { + SECTION("interpolation gridded to gridded (netCDF)") { + param::SimpleParametrisation args; + args.set("input", "checkDuplicatePoints: False"); + + for (const auto& job : jobs) { + std::unique_ptr input(input::MIRInputFactory::build("../data/nemo.nc", args)); + output::EmptyOutput output; + + while (input->next()) { + job.execute(*input, output); + } + } + } + } +} + + +} // namespace mir::tests::unit + + +int main(int argc, char** argv) { + return eckit::testing::run_tests(argc, argv); +} diff --git a/tests/unit/interpolations.cc b/tests/unit/interpolations.cc index c51956277..bf528a1ca 100644 --- a/tests/unit/interpolations.cc +++ b/tests/unit/interpolations.cc @@ -37,10 +37,10 @@ CASE("interpolations") { jobs[1].set("area", "40/20/20/40"); jobs[1].set("frame", 2); -#if mir_HAVE_ATLAS - jobs[2].set("caching", false); - jobs[2].set("rotation", "-89/10"); -#endif + if constexpr (HAVE_ATLAS) { + jobs[2].set("caching", false); + jobs[2].set("rotation", "-89/10"); + } SECTION("gridded to gridded (GRIB)") { @@ -62,58 +62,39 @@ CASE("interpolations") { } -#if mir_HAVE_NETCDF -#if atlas_HAVE_TESSELATION - SECTION("gridded to gridded (netCDF)") { - param::SimpleParametrisation args; - args.set("input", "checkDuplicatePoints: False"); - - for (const auto& job : jobs) { - std::unique_ptr input(input::MIRInputFactory::build("../data/nemo.nc", args)); - output::EmptyOutput output; + if constexpr (HAVE_ATLAS) { + SECTION("spectral to gridded (scalar)") { + param::SimpleParametrisation args; - while (input->next()) { - job.execute(*input, output); - } - } - } -#endif -#endif - - -#if mir_HAVE_ATLAS - SECTION("spectral to gridded (scalar)") { - param::SimpleParametrisation args; - - for (const auto& job : jobs) { - std::unique_ptr input( - input::MIRInputFactory::build("../data/param=t,level=1000,resol=20", args)); - output::EmptyOutput output; + for (const auto& job : jobs) { + std::unique_ptr input( + input::MIRInputFactory::build("../data/param=t,level=1000,resol=20", args)); + output::EmptyOutput output; - while (input->next()) { - job.execute(*input, output); + while (input->next()) { + job.execute(*input, output); + } } } - } - SECTION("spectral to gridded (vod2uv)") { - param::SimpleParametrisation args; - args.set("vod2uv", true); + SECTION("spectral to gridded (vod2uv)") { + param::SimpleParametrisation args; + args.set("vod2uv", true); - for (auto& job : jobs) { - std::unique_ptr input( - input::MIRInputFactory::build("../data/param=vo_d,level=1000,resol=20", args)); - output::EmptyOutput output; + for (auto& job : jobs) { + std::unique_ptr input( + input::MIRInputFactory::build("../data/param=vo_d,level=1000,resol=20", args)); + output::EmptyOutput output; - job.set("vod2uv", true); + job.set("vod2uv", true); - while (input->next()) { - job.execute(*input, output); + while (input->next()) { + job.execute(*input, output); + } } } } -#endif } diff --git a/tests/unit/iterator.cc b/tests/unit/iterator.cc index ccdfedff6..3a7fabdfe 100644 --- a/tests/unit/iterator.cc +++ b/tests/unit/iterator.cc @@ -241,32 +241,32 @@ CASE("GaussianIterator") { {-25.061510912726, 135.0969827586207}}, }; -#if mir_HAVE_ATLAS - tests.emplace_back(test_t{"N128", - 88838, - 88838, - {90, 0, -90, 359.297}, - {89.462821568577, 0}, - {89.462821568577, 20}, - {89.462821568577, 0}, - {89.462821568577, 20}}); - tests.emplace_back(test_t{"N320", - 542080, - 13437, - {-60, 50, -90, 180}, - {89.784876907219, 0}, - {89.784876907219, 20}, - {-60.280999861572, 50.13333333333333}, - {-60.280999861572, 50.66666666666666}}); - tests.emplace_back(test_t{"N640", - 2140702, - 52654, - {-60, 50, -90, 180}, - {89.89239644558999, 0}, - {89.89239644558999, 20}, - {-60.093705799122, 50.13333333333333}, - {-60.093705799122, 50.4}}); -#endif + if constexpr (HAVE_ATLAS) { + tests.emplace_back(test_t{"N128", + 88838, + 88838, + {90, 0, -90, 359.297}, + {89.462821568577, 0}, + {89.462821568577, 20}, + {89.462821568577, 0}, + {89.462821568577, 20}}); + tests.emplace_back(test_t{"N320", + 542080, + 13437, + {-60, 50, -90, 180}, + {89.784876907219, 0}, + {89.784876907219, 20}, + {-60.280999861572, 50.13333333333333}, + {-60.280999861572, 50.66666666666666}}); + tests.emplace_back(test_t{"N640", + 2140702, + 52654, + {-60, 50, -90, 180}, + {89.89239644558999, 0}, + {89.89239644558999, 20}, + {-60.093705799122, 50.13333333333333}, + {-60.093705799122, 50.4}}); + } for (const auto& t : tests) { log << "Test " + t.grid + " (global)" << std::endl; diff --git a/tests/unit/style.cc b/tests/unit/style.cc index e2104eb23..b4ad01e30 100644 --- a/tests/unit/style.cc +++ b/tests/unit/style.cc @@ -82,8 +82,6 @@ CASE("ECMWFStyle") { Log::info() << std::boolalpha; static const param::DefaultParametrisation defaults; - std::vector _no_yes{false, true}; - std::vector _yes_no{true, false}; std::vector _when{"prologue", "gridded", "spectral", "raw", "epilogue"}; @@ -110,14 +108,15 @@ CASE("ECMWFStyle") { const action::FormulaAction CORRECT_ACTION(p4); -#if mir_HAVE_ATLAS - for (bool input_gridded : _yes_no) { -#else - for (bool input_gridded : {true}) { -#endif + std::vector inputs_gridded{true}; + if constexpr (HAVE_ATLAS) { + inputs_gridded.push_back(false); + } + + for (bool input_gridded : inputs_gridded) { TestingInput in(input_gridded); - for (bool output_gridded : _yes_no) { + for (bool output_gridded : {true, false}) { if (input_gridded && !output_gridded) { // this combination isn't supported continue; @@ -135,7 +134,7 @@ CASE("ECMWFStyle") { : when == "spectral" ? (!input_gridded || !output_gridded) : true; - for (bool addWrongArguments : _no_yes) { + for (bool addWrongArguments : {false, true}) { if (addWrongArguments) { for (const std::string& wrongWhen : _when) { if (wrongWhen != when) { diff --git a/tests/unit/weight_matrix.cc b/tests/unit/weight_matrix.cc new file mode 100644 index 000000000..9d252ec4f --- /dev/null +++ b/tests/unit/weight_matrix.cc @@ -0,0 +1,107 @@ +/* + * (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 "eckit/testing/Test.h" + +#include "mir/method/WeightMatrix.h" +#include "mir/util/Exceptions.h" + + +namespace mir::tests::unit { + + +CASE("WeightMatrix::validate") { + SECTION("out-of-bounds") { + const auto* when{"out-of-bounds"}; + const std::string what{ + "Invalid weight matrix (out-of-bounds): 1 row error, " + "row 2: weights out-of-bounds, contents: (2,2,-0.1)"}; + + method::WeightMatrix W(3, 3); + W.setFromTriplets({{0, 0, 1.}, {1, 1, 0.}, {2, 2, -0.1}}); + + try { + W.validate(when, {false, true, false}); + ASSERT(false); + } + catch (exception::InvalidWeightMatrix& e) { + EXPECT(e.what() == what); + } + } + + + SECTION("weights sum") { + const auto* when{"weights sum"}; + const std::string what{ + "Invalid weight matrix (weights sum): 2 row errors, " + "row 1: weights sum not 0 or 1 (sum=0.5, 1-sum=0.5), contents: (1,1,0.5), " + "row 2: weights sum not 0 or 1 (sum=0.1, 1-sum=0.9), contents: (2,2,0.1)"}; + + method::WeightMatrix W(3, 3); + W.setFromTriplets({{0, 0, 1.}, {1, 1, 0.5}, {2, 2, 0.1}}); + + try { + W.validate(when, {false, false, true}); + ASSERT(false); + } + catch (exception::InvalidWeightMatrix& e) { + EXPECT(e.what() == what); + } + } + + + SECTION("duplicate indices") { + const auto* when{"duplicate indices"}; + const std::string what{ + "Invalid weight matrix (duplicate indices): 1 row error, " + "row 1: duplicate indices, contents: (1,1,0.5), (1,1,0.5)"}; + + method::WeightMatrix W(3, 3); + W.setFromTriplets({{0, 0, 1.}, {1, 1, 0.5}, {1, 1, 0.5}}); + + try { + W.validate(when, {true, false, false}); + ASSERT(false); + } + catch (exception::InvalidWeightMatrix& e) { + EXPECT(e.what() == what); + } + } + + + SECTION("mixed") { + const auto* when{"mixed"}; + const std::string what{ + "Invalid weight matrix (mixed): 2 row errors, " + "row 0: duplicate indices, contents: (0,0,0.5), (0,0,0.5), " + "row 1: weights sum not 0 or 1 (sum=0.5, 1-sum=0.5), contents: (1,1,0.5)"}; + + method::WeightMatrix W(3, 3); + W.setFromTriplets({{0, 0, 0.5}, {0, 0, 0.5}, {1, 1, 0.5}}); + + try { + W.validate(when, {true, false, true}); + ASSERT(false); + } + catch (exception::InvalidWeightMatrix& e) { + EXPECT(e.what() == what); + } + } +} + + +} // namespace mir::tests::unit + + +int main(int argc, char** argv) { + return eckit::testing::run_tests(argc, argv); +}