diff --git a/CMakeLists.txt b/CMakeLists.txt index 05bd33688..571b24f77 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -112,6 +112,14 @@ ecbuild_add_option( FEATURE GEOGRAPHY DESCRIPTION "Support for Geoiterator and nearest neighbour" DEFAULT ON ) +ecbuild_add_option( FEATURE ECKIT_GEO + DESCRIPTION "Support for Geoiterator and nearest neighbour (additional backend)" + CONDITION ENABLE_GEOGRAPHY + DEFAULT OFF ) +if( ECKIT_GEO AND NOT TARGET eckit_geo ) + find_library(ECKIT_GEO_LIB NAMES eckit_geo REQUIRED) +endif() + ecbuild_add_option( FEATURE JPG DESCRIPTION "Support for JPG decoding/encoding" DEFAULT ON ) @@ -485,7 +493,7 @@ ecbuild_pkgconfig( IGNORE_INCLUDE_DIRS ${PYTHON_INCLUDE_DIRS} ${NUMPY_INCLUDE_DIRS} ${NETCDF_INCLUDE_DIRS} VARIABLES HAVE_MEMFS HAVE_GEOGRAPHY HAVE_JPEG HAVE_LIBJASPER HAVE_LIBOPENJPEG HAVE_ECCODES_THREADS HAVE_ECCODES_OMP_THREADS - HAVE_NETCDF HAVE_FORTRAN HAVE_PNG HAVE_AEC + HAVE_NETCDF HAVE_FORTRAN HAVE_PNG HAVE_AEC HAVE_ECKIT_GEO ) if( HAVE_FORTRAN ) ecbuild_pkgconfig( @@ -497,7 +505,7 @@ if( HAVE_FORTRAN ) ${PYTHON_INCLUDE_DIRS} ${NUMPY_INCLUDE_DIRS} ${NETCDF_INCLUDE_DIRS} VARIABLES HAVE_MEMFS HAVE_GEOGRAPHY HAVE_JPEG HAVE_LIBJASPER HAVE_LIBOPENJPEG HAVE_ECCODES_THREADS HAVE_ECCODES_OMP_THREADS - HAVE_NETCDF HAVE_PNG HAVE_AEC + HAVE_NETCDF HAVE_PNG HAVE_AEC HAVE_ECKIT_GEO ) endif() diff --git a/eccodes_config.h.in b/eccodes_config.h.in index 29877a3d9..a003e802d 100644 --- a/eccodes_config.h.in +++ b/eccodes_config.h.in @@ -114,6 +114,7 @@ #cmakedefine HAVE_NETCDF #cmakedefine HAVE_GEOGRAPHY +#cmakedefine HAVE_ECKIT_GEO #cmakedefine HAVE_MEMFS #cmakedefine HAVE_FORTRAN diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 01d497bb1..8a5dcd280 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -411,6 +411,15 @@ if( HAVE_MEMFS ) list(APPEND ECCODES_EXTRA_LIBRARIES eccodes_memfs) endif() +if( eccodes_HAVE_GEOGRAPHY AND eccodes_HAVE_ECKIT_GEO ) + list( APPEND eccodes_src_files + eccodes/geo/GeoIterator.cc + eccodes/geo/GeoIterator.h + eccodes/geo/GribSpec.cc + eccodes/geo/GribSpec.h ) + list( APPEND ECCODES_EXTRA_LIBRARIES eckit_geo ) +endif() + ecbuild_add_library( TARGET eccodes SOURCES ${CMAKE_CURRENT_BINARY_DIR}/grib_api_version.cc # griby.cc gribl.cc diff --git a/src/eccodes/geo/GeoIterator.cc b/src/eccodes/geo/GeoIterator.cc new file mode 100644 index 000000000..f20af7bab --- /dev/null +++ b/src/eccodes/geo/GeoIterator.cc @@ -0,0 +1,104 @@ +/* + * (C) Copyright 2024- 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 "eccodes/geo/GeoIterator.h" + +#include "eckit/exception/Exceptions.h" + +#include "eccodes/geo/GribSpec.h" + + +namespace eccodes::geo +{ + + +GeoIterator::GeoIterator(grib_handle* h, unsigned long flags) : + spec_(new GribSpec(h)), grid_(eckit::geo::GridFactory::build(*spec_)), iter_(grid_->cbegin().release()), end_(grid_->cend().release()) +{ + h_ = h; + class_name_ = "geo_iterator"; + flags_ = flags; + Assert(h_ != nullptr); + + CODES_CHECK(codes_get_size(h_, "values", &nv_), ""); + Assert(nv_ > 0); + + data_ = (flags_ & GRIB_GEOITERATOR_NO_VALUES) ? nullptr : static_cast(grib_context_malloc(h_->context, nv_ * sizeof(double))); + Assert(data_ != nullptr); + + auto size = nv_; + CODES_CHECK(codes_get_double_array(h_, "values", data_, &size), ""); + Assert(nv_ == size); +} + + +int GeoIterator::init(grib_handle*, grib_arguments*) +{ + NOTIMP; +} + + +int GeoIterator::next(double* lat, double* lon, double* val) const +{ + if (iter_ == end_) { + return 0; // (false) + } + + const auto p = *iter_; + const auto& q = std::get(p); + + *lat = q.lat; + *lon = q.lon; + if (val != nullptr && data_ != nullptr) { + *val = data_[iter_->index()]; + } + + ++iter_; + return 1; // (true) +} + + +int GeoIterator::previous(double*, double*, double*) const +{ + return GRIB_NOT_IMPLEMENTED; +} + + +int GeoIterator::reset() +{ + iter_.reset(grid_->cbegin().release()); + return GRIB_SUCCESS; +} + + +int GeoIterator::destroy() +{ + if (data_ != nullptr) { + grib_context_free(h_->context, data_); + } + return Iterator::destroy(); +} + + +bool GeoIterator::has_next() const +{ + return iter_ != end_; +} + + +geo_iterator::Iterator* +GeoIterator::create() const +{ + return new GeoIterator{ h_, flags_ }; +} + + +} // namespace eccodes::geo diff --git a/src/eccodes/geo/GeoIterator.h b/src/eccodes/geo/GeoIterator.h new file mode 100644 index 000000000..312b0357c --- /dev/null +++ b/src/eccodes/geo/GeoIterator.h @@ -0,0 +1,51 @@ +/* + * (C) Copyright 2024- 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 "eckit/geo/Grid.h" + +// eccodes macros conflict with eckit +#ifdef Assert + #undef Assert +#endif + +#include "geo_iterator/grib_iterator.h" + + +namespace eccodes::geo +{ + + +class GeoIterator : public geo_iterator::Iterator +{ +public: + explicit GeoIterator(grib_handle*, unsigned long flags); + +private: + std::unique_ptr spec_; + std::unique_ptr grid_; + mutable eckit::geo::Grid::Iterator iter_; + eckit::geo::Grid::Iterator end_; + + int init(grib_handle*, grib_arguments*) override; + int next(double* lat, double* lon, double* val) const override; + int previous(double* lat, double* lon, double* val) const override; + int reset() override; + int destroy() override; + bool has_next() const override; + geo_iterator::Iterator* create() const override; +}; + + +} // namespace eccodes::geo diff --git a/src/eccodes/geo/GribSpec.cc b/src/eccodes/geo/GribSpec.cc new file mode 100644 index 000000000..3d6a58fe6 --- /dev/null +++ b/src/eccodes/geo/GribSpec.cc @@ -0,0 +1,1104 @@ +/* + * (C) Copyright 2024- 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 "eccodes/geo/GribSpec.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "eckit/config/Resource.h" +#include "eckit/exception/Exceptions.h" +#include "eckit/geo/PointLonLat.h" +#include "eckit/geo/util/mutex.h" +#include "eckit/log/JSON.h" +#include "eckit/types/FloatCompare.h" +#include "eckit/types/Fraction.h" + + +namespace eckit::geo::util +{ +const std::vector& gaussian_latitudes(size_t N, bool increasing); +} + + +namespace eccodes::geo +{ + + +bool codes_check_error(int e, const char* call) +{ + if (e != CODES_SUCCESS) { + std::ostringstream os; + os << call << ": " << codes_get_error_message(e); + throw ::eckit::SeriousBug(os.str()); + } + return true; +} + + +#define CHECK_ERROR(a, b) codes_check_error(a, b) +#define CHECK_CALL(a) codes_check_error(a, #a) + + +namespace +{ + + +using eckit::Log; + + +struct Condition +{ + Condition() = default; + + Condition(const Condition&) = delete; + Condition(Condition&&) = delete; + Condition& operator=(const Condition&) = delete; + Condition& operator=(Condition&&) = delete; + + virtual ~Condition() = default; + virtual bool eval(codes_handle*) const = 0; +}; + + +template +struct ConditionT : Condition +{ + const char* key_; + T value_; + bool eval(codes_handle* /*unused*/) const override; + +public: + ConditionT(const char* key, const T& value) : + key_(key), value_(value) {} +}; + + +template <> +bool ConditionT::eval(codes_handle* h) const +{ + ASSERT(h != nullptr); + + long value = 0; + int err = codes_get_long(h, key_, &value); + + if (err == CODES_NOT_FOUND) { + return false; + } + + CHECK_ERROR(err, key_); + return value_ == value; +} + + +template <> +bool ConditionT::eval(codes_handle* h) const +{ + ASSERT(h != nullptr); + + double value = 0; + int err = codes_get_double(h, key_, &value); + + if (err == CODES_NOT_FOUND) { + return false; + } + + CHECK_ERROR(err, key_); + return value_ == value; // Want an epsilon? +} + + +template <> +bool ConditionT::eval(codes_handle* h) const +{ + ASSERT(h != nullptr); + + char buffer[10240]; + size_t size = sizeof(buffer); + int err = codes_get_string(h, key_, buffer, &size); + + if (err == CODES_NOT_FOUND) { + return false; + } + + CHECK_ERROR(err, key_); + return value_ == buffer; +} + + +struct ConditionOR : Condition +{ + const std::unique_ptr left_; + const std::unique_ptr right_; + bool eval(codes_handle* h) const override { return left_->eval(h) || right_->eval(h); } + + ConditionOR(const Condition* left, const Condition* right) : + left_(left), right_(right) {} +}; + + +struct ConditionAND : Condition +{ + const std::unique_ptr left_; + const std::unique_ptr right_; + bool eval(codes_handle* h) const override { return left_->eval(h) && right_->eval(h); } + + ConditionAND(const Condition* left, const Condition* right) : + left_(left), right_(right) {} +}; + + +/* +class ConditionNOT : public Condition { + const Condition* c_; + bool eval(codes_handle* h) const override { return !c_->eval(h); } + ~ConditionNOT() override { + delete c_; + } + +public: + ConditionNOT(const Condition* c) : c_(c) {} +}; +*/ + + +void wrongly_encoded_grib(const std::string& msg) +{ + static bool abortIfWronglyEncodedGRIB = eckit::Resource("$MIR_ABORT_IF_WRONGLY_ENCODED_GRIB", false); + + if (abortIfWronglyEncodedGRIB) { + Log::error() << msg << std::endl; + throw eckit::UserError(msg); + } + + Log::warning() << msg << std::endl; +} + + +template +Condition* is(const char* key, const T& value) +{ + return new ConditionT(key, value); +} + + +Condition* is(const char* key, const char* value) +{ + return new ConditionT(key, value); +} + + +Condition* _and(const Condition* left, const Condition* right) +{ + return new ConditionAND(left, right); +} + + +Condition* _or(const Condition* left, const Condition* right) +{ + return new ConditionOR(left, right); +} + + +/* + Condition *_not(const Condition *c) { + return new ConditionNOT(c); +} +*/ + + +Condition* is_gaussian() +{ + const auto* key = "gridType"; + return _or(is(key, "reduced_gg"), + _or(is(key, "regular_gg"), + _or(is(key, "reduced_rotated_gg"), + _or(is(key, "regular_rotated_gg"), + _or(is(key, "rotated_gg"), + _or(is(key, "reduced_stretched_gg"), + _or(is(key, "regular_stretched_gg"), + _or(is(key, "reduced_stretched_rotated_gg"), + _or(is(key, "regular_stretched_rotated_gg"), + _or(is(key, "stretched_gg"), is(key, "stretched_rotated_gg"))))))))))); +} + + +const char* get_key(const std::string& name, codes_handle* h) +{ + struct P + { + const std::string name; + const char* key; + const std::unique_ptr condition; + P(const std::string _name, const char* _key, const Condition* _condition = nullptr) : + name(_name), key(_key), condition(_condition) {} + }; + + static const std::initializer_list

mappings{ + { "type", "gridType" }, + + { "west_east_increment", "iDirectionIncrementInDegrees_fix_for_periodic_regular_grids", + is("gridType", "regular_ll") }, + { "west_east_increment", "iDirectionIncrementInDegrees" }, + { "south_north_increment", "jDirectionIncrementInDegrees" }, + + { "west", "longitudeOfFirstGridPointInDegrees" }, + { "east", "longitudeOfLastGridPointInDegrees_fix_for_global_reduced_grids", is("gridType", "reduced_gg") }, + { "east", "longitudeOfLastGridPointInDegrees" }, + + { "north", "latitudeOfFirstGridPointInDegrees_fix_for_gaussian_grids", + _and(is("scanningMode", 0L), is_gaussian()) }, + { "south", "latitudeOfLastGridPointInDegrees_fix_for_gaussian_grids", + _and(is("scanningMode", 0L), is_gaussian()) }, + { "north", "latitudeOfLastGridPointInDegrees_fix_for_gaussian_grids", + _and(is("scanningMode", 1L), is_gaussian()) }, + { "south", "latitudeOfFirstGridPointInDegrees_fix_for_gaussian_grids", + _and(is("scanningMode", 1L), is_gaussian()) }, + + { "north", "latitudeOfLastGridPointInDegrees", is("jScansPositively", 1L) }, + { "south", "latitudeOfFirstGridPointInDegrees", is("jScansPositively", 1L) }, + { "north", "latitudeOfFirstGridPointInDegrees", is("scanningMode", 0L) }, + { "south", "latitudeOfLastGridPointInDegrees", is("scanningMode", 0L) }, + + { "north", "latitudeOfLastGridPointInDegrees", is("jScansPositively", 1L) }, + { "south", "latitudeOfFirstGridPointInDegrees", is("jScansPositively", 1L) }, + { "north", "latitudeOfFirstGridPointInDegrees" }, + { "south", "latitudeOfLastGridPointInDegrees" }, + + { "truncation", "pentagonalResolutionParameterJ" }, // Assumes triangular truncation + { "accuracy", "bitsPerValue" }, + + { "south_pole_latitude", "latitudeOfSouthernPoleInDegrees" }, + { "south_pole_longitude", "longitudeOfSouthernPoleInDegrees" }, + { "south_pole_rotation_angle", "angleOfRotationInDegrees" }, + + { "proj", "projTargetString" }, + { "projSource", "projSourceString" }, + + // This will be just called for has() + { + "gridded", + "Nx", + _or(_or(_or(is("gridType", "polar_stereographic"), is("gridType", "lambert_azimuthal_equal_area")), + is("gridType", "lambert")), + is("gridType", "space_view")), + }, + { + "gridded", + "Ni", + is("gridType", "triangular_grid"), + }, + { + "gridded", + "numberOfGridInReference" /*just a dummy*/, + is("gridType", "unstructured_grid"), + }, + { "gridded", "numberOfPointsAlongAMeridian" }, // Is that always true? + { "gridded_regular_ll", "Ni", _or(is("gridType", "regular_ll"), is("gridType", "rotated_ll")) }, + { "gridded_named", "gridName" }, + + { "grid", "gridName", + _or(_or(_or(_or(is("gridType", "regular_gg"), is("gridType", "reduced_gg")), is("gridType", "rotated_gg")), + is("gridType", "reduced_rotated_gg")), + is("gridType", "unstructured_grid")) }, + + { "spectral", "pentagonalResolutionParameterJ" }, + + { "uid", "uuidOfHGrid", is("gridType", "unstructured_grid") }, + + /// FIXME: Find something that does no clash + { "reduced", "numberOfParallelsBetweenAPoleAndTheEquator", is("isOctahedral", 0L) }, + { "regular", "N", is("gridType", "regular_gg") }, + { "octahedral", "numberOfParallelsBetweenAPoleAndTheEquator", is("isOctahedral", 1L) }, + + /// TODO: is that a good idea? + { "param", "paramId" }, + { "statistics", "" }, // (avoid ecCodes error "statistics: Function not yet implemented") + }; + + for (const auto& m : mappings) { + if (name == m.name) { + if (!m.condition || m.condition->eval(h)) { + return m.key; + } + } + } + + const auto* key = name.c_str(); + return key; +} + + +template +struct ProcessingT +{ + using fun_t = std::function; + fun_t fun_; + explicit ProcessingT(fun_t&& fun) : + fun_(fun) {} + ~ProcessingT() = default; + ProcessingT(const ProcessingT&) = delete; + ProcessingT(ProcessingT&&) = delete; + void operator=(const ProcessingT&) = delete; + void operator=(ProcessingT&&) = delete; + bool eval(codes_handle* h, T& v) const { return fun_(h, v); } +}; + + +ProcessingT* angular_precision() +{ + return new ProcessingT([](codes_handle* h, double& value) { + auto well_defined = [](codes_handle* h, const char* key) -> bool { + long dummy = 0; + int err = 0; + return (codes_is_defined(h, key) != 0) && (codes_is_missing(h, key, &err) == 0) && (err == CODES_SUCCESS) && + (codes_get_long(h, key, &dummy) == CODES_SUCCESS) && (dummy != 0); + }; + + if (well_defined(h, "basicAngleOfTheInitialProductionDomain") && well_defined(h, "subdivisionsOfBasicAngle")) { + value = 0.; + return true; + } + + long angleSubdivisions = 0; + CHECK_CALL(codes_get_long(h, "angleSubdivisions", &angleSubdivisions)); + + value = angleSubdivisions > 0 ? 1. / static_cast(angleSubdivisions) : 0.; + return true; + }); +} + + +ProcessingT* longitudeOfLastGridPointInDegrees_fix_for_global_reduced_grids() +{ + return new ProcessingT([](codes_handle* h, double& Lon2) { + Lon2 = 0; + CHECK_CALL(codes_get_double(h, "longitudeOfLastGridPointInDegrees", &Lon2)); + + if (codes_is_defined(h, "pl") != 0) { + double Lon1 = 0; + CHECK_CALL(codes_get_double(h, "longitudeOfFirstGridPointInDegrees", &Lon1)); + + if (eckit::types::is_approximately_equal(Lon1, 0)) { + // get pl array maximum and sum + // if sum equals values size the grid must be global + size_t plSize = 0; + CHECK_CALL(codes_get_size(h, "pl", &plSize)); + ASSERT(plSize > 0); + + std::vector pl(plSize, 0); + size_t plSizeAsRead = plSize; + CHECK_CALL(codes_get_long_array(h, "pl", pl.data(), &plSizeAsRead)); + ASSERT(plSize == plSizeAsRead); + + long plMax = 0; + long plSum = 0; + for (auto& p : pl) { + plSum += p; + if (plMax < p) { + plMax = p; + } + } + ASSERT(plMax > 0); + + size_t valuesSize = 0; + CHECK_CALL(codes_get_size(h, "values", &valuesSize)); + + if (static_cast(plSum) == valuesSize) { + double eps = 0.; + ASSERT(std::unique_ptr>(angular_precision())->eval(h, eps)); + + eckit::Fraction Lon2_expected(360L * (plMax - 1L), plMax); + + if (!eckit::types::is_approximately_greater_or_equal(Lon2, Lon2_expected, eps)) { + std::ostringstream msgs; + msgs.precision(32); + msgs << "GribParametrisation: wrongly encoded longitudeOfLastGridPointInDegrees:" + << "\n" + "encoded: " + << Lon2 + << "\n" + "expected: " + << static_cast(Lon2_expected) << " (" << Lon2_expected << " +- " << eps << ")"; + + wrongly_encoded_grib(msgs.str()); + } + + Lon2 = 360.; + } + } + } + + return true; + }); +}; + + +ProcessingT* latitudeOfFirstGridPointInDegrees_fix_for_gaussian_grids() +{ + return new ProcessingT([](codes_handle* h, double& lat) { + double eps = 0.; + ASSERT(std::unique_ptr>(angular_precision())->eval(h, eps)); + + long N = 0L; + CHECK_CALL(codes_get_long(h, "N", &N)); + ASSERT(N > 0); + + double& Lat1 = lat; + CHECK_CALL(codes_get_double(h, "latitudeOfFirstGridPointInDegrees", &Lat1)); + + double Lat2 = 0.; + CHECK_CALL(codes_get_double(h, "latitudeOfLastGridPointInDegrees", &Lat2)); + + const auto snap = eckit::geo::util::gaussian_latitudes(N, Lat1 < Lat2).front(); + if (eckit::types::is_approximately_equal(Lat1, snap, eps)) { + lat = snap < 0. ? eckit::geo::SOUTH_POLE.lat : eckit::geo::NORTH_POLE.lat; + } + + return true; + }); +} + + +ProcessingT* latitudeOfLastGridPointInDegrees_fix_for_gaussian_grids() +{ + return new ProcessingT([](codes_handle* h, double& lat) { + double eps = 0.; + ASSERT(std::unique_ptr>(angular_precision())->eval(h, eps)); + + long N = 0L; + CHECK_CALL(codes_get_long(h, "N", &N)); + ASSERT(N > 0); + + double Lat1 = 0.; + CHECK_CALL(codes_get_double(h, "latitudeOfFirstGridPointInDegrees", &Lat1)); + + double& Lat2 = lat; + CHECK_CALL(codes_get_double(h, "latitudeOfLastGridPointInDegrees", &Lat2)); + + const auto snap = eckit::geo::util::gaussian_latitudes(N, Lat1 < Lat2).back(); + if (eckit::types::is_approximately_equal(Lat2, snap, eps)) { + lat = snap < 0. ? eckit::geo::SOUTH_POLE.lat : eckit::geo::NORTH_POLE.lat; + } + + return true; + }); +} + + +ProcessingT* iDirectionIncrementInDegrees_fix_for_periodic_regular_grids() +{ + return new ProcessingT([](codes_handle* h, double& we) { + long iScansPositively = 0L; + CHECK_CALL(codes_get_long(h, "iScansPositively", &iScansPositively)); + ASSERT(iScansPositively == 1L); + + CHECK_CALL(codes_get_double(h, "iDirectionIncrementInDegrees", &we)); + ASSERT(we > 0.); + + double Lon1 = 0.; + double Lon2 = 0.; + long Ni = 0; + CHECK_CALL(codes_get_double(h, "longitudeOfFirstGridPointInDegrees", &Lon1)); + CHECK_CALL(codes_get_double(h, "longitudeOfLastGridPointInDegrees", &Lon2)); + CHECK_CALL(codes_get_long(h, "Ni", &Ni)); + ASSERT(Ni > 0); + + Lon2 = eckit::geo::PointLonLat::normalise_angle_to_minimum(Lon2, Lon1); + ASSERT(Lon2 >= Lon1); + + // angles are within +-1/2 precision, so (Lon2 - Lon1 + we) uses factor 3*1/2 + double eps = 0.; + ASSERT(std::unique_ptr>(angular_precision())->eval(h, eps)); + eps *= 1.5; + + constexpr double GLOBE = 360; + + auto Nid = static_cast(Ni); + if (eckit::types::is_approximately_equal(Lon2 - Lon1 + we, GLOBE, eps)) { + we = GLOBE / Nid; + } + else if (!eckit::types::is_approximately_equal(Lon1 + (Nid - 1) * we, Lon2, eps)) { + // TODO refactor, not really specific to "periodic regular grids", but useful + std::ostringstream msgs; + msgs.precision(32); + msgs << "GribParametrisation: wrongly encoded iDirectionIncrementInDegrees:" + "\n" + "encoded: " + << we + << "\n" + "Ni: " + << Ni + << "\n" + "longitudeOfFirstGridPointInDegree: " + << Lon1 + << "\n" + "longitudeOfLastGridPointInDegrees: " + << Lon2; + + wrongly_encoded_grib(msgs.str()); + } + + return true; + }); +}; + + +ProcessingT>* vector_double(std::initializer_list keys) +{ + const std::vector keys_(keys); + return new ProcessingT>([=](codes_handle* h, std::vector& values) { + ASSERT(keys.size()); + + values.assign(keys_.size(), 0); + size_t i = 0; + for (const auto& key : keys_) { + if (codes_is_defined(h, key.c_str()) == 0) { + return false; + } + CHECK_CALL(codes_get_double(h, key.c_str(), &values[i++])); + } + return true; + }); +} + + +ProcessingT* packing() +{ + return new ProcessingT([](codes_handle* h, std::string& value) { + auto get = [](codes_handle* h, const char* key) -> std::string { + if (codes_is_defined(h, key) != 0) { + char buffer[64]; + size_t size = sizeof(buffer); + + CHECK_CALL(codes_get_string(h, key, buffer, &size)); + ASSERT(size < sizeof(buffer) - 1); + + if (std::strcmp(buffer, "MISSING") != 0) { + return buffer; + } + } + return ""; + }; + + auto packingType = get(h, "packingType"); + for (const auto& prefix : std::vector{ "grid_", "spectral_" }) { + if (packingType.find(prefix) == 0) { + value = packingType.substr(prefix.size()); + std::replace(value.begin(), value.end(), '_', '-'); + return true; + } + } + + return false; + }); +} + + +template +struct ConditionedProcessingT +{ + const std::string name; + const std::unique_ptr processing; + const std::unique_ptr condition; + ConditionedProcessingT(const std::string& _name, const T* _processing, const Condition* _condition = nullptr) : + name(_name), processing(_processing), condition(_condition) {} +}; + + +template +using ProcessingList = std::initializer_list>>; + + +template +bool get_value(const std::string& name, codes_handle* h, T& value, const ProcessingList& process) +{ + for (auto& p : process) { + if (name == p.name) { + if (!p.condition || p.condition->eval(h)) { + ASSERT(p.processing); + return p.processing->eval(h, value); + } + } + } + return false; +} + + +eckit::geo::util::recursive_mutex MUTEX; + + +class lock_type +{ + eckit::geo::util::lock_guard lock_guard_{ MUTEX }; +}; + + +} // namespace + + +GribSpec::GribSpec(codes_handle* h) : + handle_(h) +{ + ASSERT(handle_ != nullptr); +} + + +bool GribSpec::has(const std::string& name) const +{ + lock_type lock; + + if (cache_.has(name)) { + return true; + } + + const auto* key = get_key(name, handle_); + + ASSERT(key != nullptr); + if (std::strlen(key) == 0) { + return false; + } + + return codes_is_defined(handle_, key) != 0; +} + + +bool GribSpec::get(const std::string& name, std::string& value) const +{ + lock_type lock; + + if (cache_.get(name, value)) { + return true; + } + + const auto* key = get_key(name, handle_); + + ASSERT(key != nullptr); + if (std::strlen(key) == 0) { + return false; + } + + char buffer[10240]; + size_t size = sizeof(buffer); + int err = codes_get_string(handle_, key, buffer, &size); + + if (err == CODES_NOT_FOUND) { + static const ProcessingList process{ + { "packing", packing() }, + }; + + return get_value(key, handle_, value, process); + } + + CHECK_ERROR(err, key); + + ASSERT(size < sizeof(buffer) - 1); + + if (std::strcmp(buffer, "MISSING") == 0) { + return false; + } + + cache_.set(name, value = buffer); + return true; +} + + +bool GribSpec::get(const std::string& name, bool& value) const +{ + lock_type lock; + + if (cache_.get(name, value)) { + return true; + } + + const auto* key = get_key(name, handle_); + + ASSERT(key != nullptr); + if (std::strlen(key) == 0) { + return false; + } + + // FIXME: make sure that 'temp' is not set if CODES_MISSING_LONG + long temp = CODES_MISSING_LONG; + int err = codes_get_long(handle_, key, &temp); + CHECK_ERROR(err, key); + + cache_.set(name, value = temp != 0); + return true; +} + + +bool GribSpec::get(const std::string& name, int& value) const +{ + if (long v = 0; get(name, v)) { + ASSERT(static_cast(static_cast(v)) == v); + value = static_cast(v); + return true; + } + + return false; +} + + +bool GribSpec::get(const std::string& name, long& value) const +{ + lock_type lock; + + if (cache_.get(name, value)) { + return true; + } + + const std::string key = get_key(name, handle_); + if (key.empty()) { + return false; + } + + // FIXME: make sure that 'value' is not set if CODES_MISSING_LONG + int err = codes_get_long(handle_, key.c_str(), &value); + if (err == CODES_NOT_FOUND || codes_is_missing(handle_, key.c_str(), &err) != 0) { + return false; + } + + CHECK_ERROR(err, key.c_str()); + + cache_.set(name, value); + return true; +} + + +bool GribSpec::get(const std::string& /*name*/, long long& /*value*/) const +{ + return false; +} + + +bool GribSpec::get(const std::string& /*name*/, std::size_t& /*value*/) const +{ + return false; +} + + +bool GribSpec::get(const std::string& name, float& value) const +{ + if (cache_.get(name, value)) { + return true; + } + + if (double v = 0; get(name, v)) { + cache_.set(name, value = static_cast(v)); + return true; + } + + return false; +} + + +bool GribSpec::get(const std::string& name, double& value) const +{ + lock_type lock; + + if (cache_.get(name, value)) { + return true; + } + + ASSERT(name != "grid"); + + const auto* key = get_key(name, handle_); + + ASSERT(key != nullptr); + if (std::strlen(key) == 0) { + return false; + } + + // FIXME: make sure that 'value' is not set if CODES_MISSING_DOUBLE + int err = codes_get_double(handle_, key, &value); + if (err == CODES_NOT_FOUND || codes_is_missing(handle_, key, &err) != 0) { + static const ProcessingList process{ + { "angular_precision", angular_precision() }, + { "longitudeOfLastGridPointInDegrees_fix_for_global_reduced_grids", + longitudeOfLastGridPointInDegrees_fix_for_global_reduced_grids() }, + { "iDirectionIncrementInDegrees_fix_for_periodic_regular_grids", + iDirectionIncrementInDegrees_fix_for_periodic_regular_grids() }, + { "latitudeOfFirstGridPointInDegrees_fix_for_gaussian_grids", + latitudeOfFirstGridPointInDegrees_fix_for_gaussian_grids() }, + { "latitudeOfLastGridPointInDegrees_fix_for_gaussian_grids", + latitudeOfLastGridPointInDegrees_fix_for_gaussian_grids() }, + }; + + if (get_value(key, handle_, value, process)) { + cache_.set(name, value); + return true; + } + + return false; + } + + CHECK_ERROR(err, key); + + cache_.set(name, value); + return true; +} + + +bool GribSpec::get(const std::string& /*name*/, std::vector& /*value*/) const +{ + return false; +} + + +bool GribSpec::get(const std::string& name, std::vector& value) const +{ + lock_type lock; + + if (cache_.get(name, value)) { + return true; + } + + const auto* key = get_key(name, handle_); + + ASSERT(key != nullptr); + if (std::strlen(key) == 0) { + return false; + } + + size_t count = 0; + int err = codes_get_size(handle_, key, &count); + CHECK_ERROR(err, key); + + size_t size = count; + + value.resize(count); + + CHECK_CALL(codes_get_long_array(handle_, key, value.data(), &size)); + ASSERT(count == size); + + ASSERT(!value.empty()); + + if (name == "pl") { + if (std::find(value.rbegin(), value.rend(), 0) != value.rend()) { + wrongly_encoded_grib("GribParametrisation: pl array contains zeros"); + } + } + + cache_.set(name, value); + return true; +} + + +bool GribSpec::get(const std::string& /*name*/, std::vector& /*value*/) const +{ + return false; +} + + +bool GribSpec::get(const std::string& /*name*/, std::vector& /*value*/) const +{ + return false; +} + + +bool GribSpec::get(const std::string& name, std::vector& value) const +{ + if (cache_.get(name, value)) { + return true; + } + + if (std::vector v; get(name, v)) { + value.clear(); + value.reserve(v.size()); + for (const auto& d : v) { + value.push_back(static_cast(d)); + } + + cache_.set(name, value); + return true; + } + + return false; +} + + +bool GribSpec::get(const std::string& name, std::vector& value) const +{ + lock_type lock; + + if (cache_.get(name, value)) { + return true; + } + + const auto* key = get_key(name, handle_); + + // NOTE: MARS client sets 'grid=vector' (deprecated) which needs to be compared against GRIB gridName + ASSERT(key != nullptr); + if (std::strlen(key) == 0 || std::strncmp(key, "gridName", 8) == 0) { + return false; + } + + static const ProcessingList> process{ + { "grid", vector_double({ "iDirectionIncrementInDegrees", "jDirectionIncrementInDegrees" }), + _or(is("gridType", "regular_ll"), is("gridType", "rotated_ll")) }, + { "grid", vector_double({ "xDirectionGridLengthInMetres", "yDirectionGridLengthInMetres" }), + is("gridType", "lambert_azimuthal_equal_area") }, + { "grid", vector_double({ "DxInMetres", "DyInMetres" }), + _or(is("gridType", "lambert"), is("gridType", "polar_stereographic")) }, + { "grid", vector_double({ "DiInMetres", "DjInMetres" }), is("gridType", "mercator") }, + { "grid", vector_double({ "dx", "dy" }), is("gridType", "space_view") }, + { "rotation", vector_double({ "latitudeOfSouthernPoleInDegrees", "longitudeOfSouthernPoleInDegrees" }), + _or(_or(_or(is("gridType", "rotated_ll"), is("gridType", "rotated_gg")), is("gridType", "rotated_sh")), + is("gridType", "reduced_rotated_gg")) }, + }; + + if (get_value(key, handle_, value, process)) { + cache_.set(name, value); + return true; + } + + // FIXME make logic consistent for ::get(,*) + if (codes_is_defined(handle_, key) == 0) { + return false; + } + + size_t count = 0; + int err = codes_get_size(handle_, key, &count); + CHECK_ERROR(err, key); + + ASSERT(count > 0); + size_t size = count; + + value.resize(count); + + CHECK_CALL(codes_get_double_array(handle_, key, value.data(), &size)); + ASSERT(count == size); + + ASSERT(!value.empty()); + + cache_.set(name, value); + return true; +} + + +bool GribSpec::get(const std::string& /*name*/, std::vector& /*value*/) const +{ + return false; +} + + +void GribSpec::json(eckit::JSON& j) const +{ + j.startObject(); + + auto* kit = codes_keys_iterator_new( + handle_, CODES_KEYS_ITERATOR_SKIP_DUPLICATES | CODES_KEYS_ITERATOR_SKIP_READ_ONLY, "geography"); + ASSERT(kit != nullptr); + + for (; codes_keys_iterator_next(kit) != 0;) { + const auto* name = codes_keys_iterator_get_name(kit); + + int err = 0; + bool is_missing = 0 != codes_is_missing(handle_, name, &err); + CHECK_ERROR(err, name); + + if (is_missing) { + j << name << "MISSING"; + continue; + } + + size_t size = 0; + CHECK_CALL(codes_get_size(handle_, name, &size)); + + int type = 0; + CHECK_CALL(codes_get_native_type(handle_, name, &type)); + + if (size > 1) { + if (type == CODES_TYPE_LONG) { + std::vector array(size); + + auto size_read = size; + CHECK_CALL(codes_get_long_array(handle_, name, array.data(), &size_read)); + ASSERT(size == size_read); + + j << name; + j.startList(); + for (const auto& value : array) { + j << value; + } + j.endList(); + continue; + } + + if (type == CODES_TYPE_DOUBLE) { + std::vector array(size); + + auto size_read = size; + CHECK_CALL(codes_get_double_array(handle_, name, array.data(), &size_read)); + ASSERT(size == size_read); + + j << name; + j.startList(); + for (const auto& value : array) { + j << value; + } + j.endList(); + continue; + } + + (j << name).startList().endList(); + continue; + } + + if (type == CODES_TYPE_LONG) { + long value = 0; + CHECK_CALL(codes_get_long(handle_, name, &value)); + j << name << value; + continue; + } + + if (type == CODES_TYPE_DOUBLE) { + double value = 0; + CHECK_CALL(codes_get_double(handle_, name, &value)); + j << name << value; + continue; + } + + if (type == CODES_TYPE_STRING) { + size_t length = 1024; + char value[length]; + CHECK_CALL(codes_get_string(handle_, name, value, &length)); + j << name << value; + continue; + } + + j << name + << (type == CODES_TYPE_UNDEFINED ? "UNDEFINED" + : type == CODES_TYPE_BYTES ? "BYTES" + : type == CODES_TYPE_SECTION ? "SECTION" + : type == CODES_TYPE_LABEL ? "LABEL" + : type == CODES_TYPE_MISSING ? "MISSING" + : "?"); + } + + codes_keys_iterator_delete(kit); + j.endObject(); +} + + +} // namespace eccodes::geo + + +#undef CHECK_ERROR +#undef CHECK_CALL diff --git a/src/eccodes/geo/GribSpec.h b/src/eccodes/geo/GribSpec.h new file mode 100644 index 000000000..01e2628f2 --- /dev/null +++ b/src/eccodes/geo/GribSpec.h @@ -0,0 +1,63 @@ +/* + * (C) Copyright 2024- 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 "eccodes.h" + +#include +#include + +#include "eckit/geo/Spec.h" +#include "eckit/geo/spec/Custom.h" + + +namespace eccodes::geo +{ + + +bool codes_check_error(int e, const char* call); + + +class GribSpec final : public eckit::geo::Spec +{ +public: + explicit GribSpec(codes_handle*); + + 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, long long& value) const override; + bool get(const std::string& name, size_t& 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; + bool get(const std::string& name, std::vector& value) const override; + bool get(const std::string& name, std::vector& value) const override; + +private: + mutable eckit::geo::spec::Custom cache_; + codes_handle* handle_; + + void json(eckit::JSON&) const final; +}; + + +} // namespace eccodes::geo diff --git a/src/geo_iterator/grib_iterator.cc b/src/geo_iterator/grib_iterator.cc index 56e34c251..4b08c5d93 100644 --- a/src/geo_iterator/grib_iterator.cc +++ b/src/geo_iterator/grib_iterator.cc @@ -12,6 +12,17 @@ * Jean Baptiste Filippi - 01.11.2005 * ***************************************************************************/ +#include "eccodes_config.h" + +#if defined(HAVE_ECKIT_GEO) + #include "eccodes/geo/GeoIterator.h" + + // eccodes macros conflict with eckit + #ifdef Assert + #undef Assert + #endif +#endif + #include "grib_iterator.h" #include "grib_iterator_factory.h" #include "accessor/grib_accessor_class_iterator.h" @@ -95,7 +106,18 @@ int grib_iterator_destroy(grib_context* c, grib_iterator* i) grib_iterator* grib_iterator_new(const grib_handle* ch, unsigned long flags, int* error) { grib_iterator* i = (grib_iterator*)grib_context_malloc_clear(ch->context, sizeof(grib_iterator)); - i->iterator = eccodes::geo_iterator::gribIteratorNew(ch, flags, error); + + #if defined(HAVE_ECKIT_GEO) + static const auto* eckit_geo = codes_getenv("ECCODES_ECKIT_GEO"); + if (eckit_geo != nullptr && strcmp(eckit_geo, "1") == 0) { + i->iterator = new eccodes::geo::GeoIterator(const_cast(ch), flags); + } + else + #endif + { + i->iterator = eccodes::geo_iterator::gribIteratorNew(ch, flags, error); + } + if (!i->iterator) { grib_context_free(ch->context, i); return NULL; diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 9cb1eec4e..0a3ee9486 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -24,6 +24,9 @@ ecbuild_add_library( TARGET ecc_tools NOINSTALL SOURCES ${ecc_tools_sources} PRIVATE_LIBS eccodes ) +if( eccodes_HAVE_ECKIT_GEO ) + target_link_libraries( ecc_tools PRIVATE eckit ) +endif() # tools binaries list( APPEND ecc_tools_binaries diff --git a/tools/grib_tools.cc b/tools/grib_tools.cc index 2c0598b0b..8c05058b2 100644 --- a/tools/grib_tools.cc +++ b/tools/grib_tools.cc @@ -9,9 +9,14 @@ */ #include "grib_tools.h" + #include #include +#ifdef HAVE_ECKIT_GEO + #include "eckit/runtime/Main.h" +#endif + #if HAVE_LIBJASPER /* Remove compiler warnings re macros being redefined */ #undef PACKAGE_BUGREPORT @@ -141,6 +146,10 @@ static grib_handle* grib_handle_new_from_file_x(grib_context* c, FILE* f, int mo int grib_tool(int argc, char** argv) { +#ifdef HAVE_ECKIT_GEO + eckit::Main::initialise(argc, argv); +#endif + int ret = 0; int i = 0; grib_context* c = grib_context_get_default();