diff --git a/src/multio/datamod/GribKeys.h b/src/multio/datamod/GribKeys.h index f67e32869..2ccef26eb 100644 --- a/src/multio/datamod/GribKeys.h +++ b/src/multio/datamod/GribKeys.h @@ -62,6 +62,31 @@ constexpr auto TypeOfProcessedDataEntry = .tagOptional() .withAccessor([](auto&& v) { return &v.typeOfProcessedData; }); +// Section 2 - local definition 39 (4D-var model errors, type=eme/me) + +constexpr auto NumberOfComponents = + EntryDef{"misc-numberOfComponents"} + .tagOptional() + .withAccessor([](auto&& v) { return &v.numberOfComponents; }); + +constexpr auto ModelErrorType = + EntryDef{"misc-modelErrorType"} + .tagOptional() + .withAccessor([](auto&& v) { return &v.modelErrorType; }); + +// Section 2 - local definition 38 (4D-var analysis iterations, type=4i) +// +// `iterationNumber` is forwarded directly through MARS as +// `mars.iteration` (see `dm::ITERATION` in MarsKeys.h); the metkit +// encoder reads it via the iteration deduction. Only the +// `totalNumberOfIterations` key, which has no MARS equivalent, is +// forwarded through misc. + +constexpr auto TotalNumberOfIterations = + EntryDef{"misc-totalNumberOfIterations"} + .tagOptional() + .withAccessor([](auto&& v) { return &v.totalNumberOfIterations; }); + // Section3 (more to be moved here from MarsMiscGeo.h) @@ -191,6 +216,10 @@ constexpr auto SatelliteSeries = EntryDef{"misc-satelliteSeries"} .tagOptional() .withAccessor([](auto&& v) { return &v.satelliteSeries; }); +constexpr auto NumberOfFrequencies = + EntryDef{"misc-numberOfFrequencies"} + .tagOptional() + .withAccessor([](auto&& v) { return &v.numberOfFrequencies; }); // Horizontal Keys constexpr auto PressureUnits = EntryDef{"misc-pressureUnits"} diff --git a/src/multio/datamod/MarsMiscGeo.h b/src/multio/datamod/MarsMiscGeo.h index 99f7c195c..d6ddc65d6 100644 --- a/src/multio/datamod/MarsMiscGeo.h +++ b/src/multio/datamod/MarsMiscGeo.h @@ -324,6 +324,7 @@ struct MiscRecord { EntryType_t typeOfEnsembleForecast; EntryType_t numberOfForecastsInEnsemble; EntryType_t satelliteSeries; + EntryType_t numberOfFrequencies; EntryType_t scaleFactorOfCentralWaveNumber; EntryType_t scaledValueOfCentralWaveNumber; EntryType_t pv; @@ -334,15 +335,21 @@ struct MiscRecord { EntryType_t bitsPerValue; EntryType_t laplacianOperator; EntryType_t subCentre; + EntryType_t numberOfComponents; + EntryType_t modelErrorType; + EntryType_t totalNumberOfIterations; + EntryType_t pvPresent; static constexpr std::string_view record_name_ = "misc"; + static constexpr auto record_entries_ = std::make_tuple( TablesVersion, GeneratingProcessIdentifier, TypeOfProcessedDataEntry, InitialStep, TimeIncrementInSeconds, LengthOfTimeWindow, LengthOfTimeWindowInSeconds, BitmapPresent, MissingValue, TypeOfEnsembleForecast, - NumberOfForecastsInEnsemble, SatelliteSeries, ScaleFactorOfCentralWaveNumber, ScaledValueOfCentralWaveNumber, - Pv, ScaleFactorOfWaveDirections, ScaleFactorOfWaveFrequencies, WaveDirections, WaveFrequencies, BitsPerValue, - LaplacianOperator, SubCentre); + NumberOfForecastsInEnsemble, SatelliteSeries, NumberOfFrequencies, ScaleFactorOfCentralWaveNumber, + ScaledValueOfCentralWaveNumber, Pv, ScaleFactorOfWaveDirections, ScaleFactorOfWaveFrequencies, WaveDirections, + WaveFrequencies, BitsPerValue, LaplacianOperator, SubCentre, NumberOfComponents, ModelErrorType, + TotalNumberOfIterations, PVPresent); }; diff --git a/src/multio/tools/grib1-to-grib2.cc b/src/multio/tools/grib1-to-grib2.cc index 0bd1c3ec7..6c3302815 100644 --- a/src/multio/tools/grib1-to-grib2.cc +++ b/src/multio/tools/grib1-to-grib2.cc @@ -10,7 +10,9 @@ #include #include +#include #include +#include #include #include #include @@ -209,6 +211,59 @@ void handleMissingValue(metkit::codes::CodesHandle& h, dm::MiscRecord& misc) { misc.missingValue.set(missingValue); } +// Spectral_complex packing in ecCodes pre-multiplies each spherical harmonic +// coefficient (m,n) by `(n*(n+1))^laplacianOperator` before quantising it. +// The packed referenceValue and the unpacked low-order subset are stored as +// IEEE 32-bit floats; any product that exceeds FLT_MAX (~3.4e38) trips an +// unrecoverable assertion inside `grib_ieee_to_long` and aborts the whole +// process - bypassing the per-message --on-error handling. +// +// IFS occasionally writes type=al SPP random fields with the sentinel value +// `laplacianOperator=9.9` even though the source data carries no +// pre-laplacian scaling. Combined with the truncation=399 of those fields, +// this yields (399*400)^9.9 ~= 1.6e51, easily overflowing FLT_MAX once +// multiplied by the actual coefficient values. Verified by `grib_set +// laplacianOperator=0` on such a message: re-encoding then succeeds and the +// codedValues match the input exactly, confirming the source data is not +// pre-laplacian-scaled. +// +// We can't compute a tight overflow bound generically from the values array +// (the largest coefficient lives at low n while the largest scaling lives at +// high n; without traversing the spectral (m,n) layout we'd get too many +// false positives). Instead, refuse to encode messages whose laplacianOperator +// is the known IFS sentinel `9.9` by throwing eckit::BadValue. The exception +// is recoverable: --on-error=log-and-skip will skip the offending message and +// continue, while --on-error=abort will propagate and stop the conversion. +// +// NOTE (revisit): the underlying issue is an upstream IFS encoding bug; once +// IFS stops emitting `laplacianOperator=9.9` for non-laplacian-scaled SPP +// random fields this guard becomes dead code and can be removed. +void validateSpectralComplexNoOverflow(const dm::FullMarsRecord& mars, const dm::MiscRecord& misc, + const std::vector& /*values*/) { + if (!mars.packing.isSet() || mars.packing.get() != "complex") { + return; + } + if (!misc.laplacianOperator.isSet()) { + return; + } + + constexpr double kIfsSentinel = 9.9; + constexpr double kEpsilon = 1.0e-3; + const double p = misc.laplacianOperator.get(); + if (std::fabs(p - kIfsSentinel) > kEpsilon) { + return; + } + + std::ostringstream oss; + oss << "spectral_complex message carries IFS sentinel laplacianOperator=" << p + << " (truncation=" << (mars.truncation.isSet() ? mars.truncation.get() : -1L) + << ", type=" << (mars.type.isSet() ? mars.type.get() : std::string{""}) + << ", paramId=" << (mars.param.isSet() ? mars.param.get() : -1L) << "). Re-encoding such messages would " + << "trigger an unrecoverable ecCodes IEEE32 overflow in `grib_ieee_to_long` because the " + << "stored data is not actually pre-laplacian-scaled. Refusing to encode."; + throw eckit::BadValue(oss.str(), Here()); +} + std::optional> parseRange(const std::string& s) { static const std::regex re(R"(^\s*([+-]?\d+)\s*-\s*([+-]?\d+)\s*$)"); std::smatch m; @@ -245,7 +300,7 @@ void handleStepRange(metkit::codes::CodesHandle& h, dm::FullMarsRecord& mars, in // Perform grib1ToGrib2 mapping - for a few marskeys we have to rely on eccodes namespace iterator. // E.g. the key "number" may be defined and set, although it has no meaning. void mapGrib1ToGrib2(KeySet& marsKeys, metkit::codes::CodesHandle& h, dm::FullMarsRecord& mars, dm::MiscRecord& misc, - int verbosity) { + int verbosity, long defaultEnsembleSize = 0) { mars.stream = dm::parseEntry(dm::STREAM, h); mars.type = dm::parseEntry(dm::TYPE, h); mars.klass = dm::parseEntry(dm::CLASS, h); @@ -259,6 +314,7 @@ void mapGrib1ToGrib2(KeySet& marsKeys, metkit::codes::CodesHandle& h, dm::FullMa } mars.anoffset = dm::parseEntry(dm::ANOFFSET, h); + mars.iteration = dm::parseEntry(dm::ITERATION, h); mars.ident = dm::parseEntry(dm::IDENT, h); mars.instrument = dm::parseEntry(dm::INSTRUMENT, h); mars.channel = dm::parseEntry(dm::CHANNEL, h); @@ -344,6 +400,19 @@ void mapGrib1ToGrib2(KeySet& marsKeys, metkit::codes::CodesHandle& h, dm::FullMa misc.initialStep = dm::parseEntry(dm::InitialStep.withKey("initialStep"), h); misc.timeIncrementInSeconds = dm::parseEntry(dm::TimeIncrementInSeconds.withKey("timeIncrement"), h); misc.pv = dm::parseEntry(dm::Pv.withKey("pv"), h); + // Signal explicit absence of vertical-coordinate parameters to the + // encoder. Without this, mars2grib's `resolve_PvArray_or_throw` falls back + // to its PV137 default and writes a 1002-double PV array into the output + // for hybrid-level fields whose source GRIB had NV=0 (e.g. lnsp on ml). + // Detect "no PV" via the input handle's NV key and forward as + // `misc.pvPresent=false`; metkit's LevelOp suppresses PV emission when + // it sees this flag. + if (!misc.pv.isSet()) { + long nv = h.has("NV") ? h.getLong("NV") : 0; + if (nv == 0) { + misc.pvPresent.set(false); + } + } misc.bitmapPresent = dm::parseEntry(dm::BitmapPresent.withKey("bitmapPresent"), h); handleMissingValue(h, misc); @@ -351,12 +420,46 @@ void mapGrib1ToGrib2(KeySet& marsKeys, metkit::codes::CodesHandle& h, dm::FullMa // TODO pgeier set it again ?? misc.laplacianOperator = dm::parseEntry(dm::LaplacianOperator.withKey("laplacianOperator"), h); + // For type=eme (ensemble model errors) or type=me (model errors) the + // GRIB1 input carries local definition 39 with componentIndex + // (= mars.number), numberOfComponents and modelErrorType. The metkit + // encoder deduces componentIndex from mars.number, but + // numberOfComponents and modelErrorType have no MARS equivalent and + // must be forwarded from the input handle. + if (mars.type.get() == "eme" || mars.type.get() == "me") { + if (h.has("numberOfComponents")) { + misc.numberOfComponents.set(h.getLong("numberOfComponents")); + } + if (h.has("modelErrorType")) { + misc.modelErrorType.set(h.getLong("modelErrorType")); + } + // mars.number == componentIndex for eme/me; the encoder reads it via + // the analysis ModelErrors path (see metkit analysisEncoding.h). + // Forwarded here unconditionally because the standard ensemble branch + // below requires `numberOfForecastsInEnsemble` (absent for eme/me) + // and would otherwise leave `mars.number` unset. + if (h.has("number")) { + mars.number.set(h.getLong("number")); + } + } + + // For type=4i (4D-var analysis iterations) the GRIB1 input carries + // local definition 38 with iterationNumber and totalNumberOfIterations. + // `iterationNumber` is forwarded directly through MARS as + // `mars.iteration` (set above); only `totalNumberOfIterations` has no + // MARS equivalent and must be forwarded through misc. + if (mars.type.get() == "4i") { + if (h.has("totalNumberOfIterations")) { + misc.totalNumberOfIterations.set(h.getLong("totalNumberOfIterations")); + } + } + // Can not rely on "number" from mars key iterator... for reference data (with hdate) number // can be 0 but is not emitted although numberOfForecastsInEnsemble has a valid value // if (auto searchNumber = marsKeys.find("number"); searchNumber != marsKeys.end()) // Check for derivedEnsembleForecast - if (mars.type.get() == "es" || mars.type.get() == "em") { + if (mars.type.get() == "es" || mars.type.get() == "em" || mars.type.get() == "ses") { long numForecasts = h.getLong("numberOfForecastsInEnsemble"); misc.numberOfForecastsInEnsemble.set(numForecasts); } @@ -366,7 +469,16 @@ void mapGrib1ToGrib2(KeySet& marsKeys, metkit::codes::CodesHandle& h, dm::FullMa // If both are 0 it is likely a control forecast or no ensemble... if (number != 0 && numForecasts == 0) { - throw std::runtime_error("The value for key numberOfForecastsInEnsemble must not be 0"); + if (defaultEnsembleSize > 0) { + std::cout << "Warning: numberOfForecastsInEnsemble is 0 but number=" << number + << "; applying --default-ensemble-size=" << defaultEnsembleSize << std::endl; + numForecasts = defaultEnsembleSize; + } + else { + throw std::runtime_error( + "The value for key numberOfForecastsInEnsemble must not be 0 (use --default-ensemble-size " + "to supply a fallback, or --on-error log-and-skip to skip such messages)"); + } } if (numForecasts != 0) { mars.number.set(number); @@ -389,6 +501,22 @@ void mapGrib1ToGrib2(KeySet& marsKeys, metkit::codes::CodesHandle& h, dm::FullMa misc.lengthOfTimeWindow = dm::parseEntry(dm::LengthOfTimeWindow.withKey("lengthOfTimeWindow"), h); misc.lengthOfTimeWindowInSeconds = dm::parseEntry(dm::LengthOfTimeWindowInSeconds.withKey("lengthOfTimeWindowInSeconds"), h); + + // Fallback: for classic 4D-Var analysis messages, the window length is + // carried in the GRIB1 local-section key `lengthOf4DvarWindow` (in hours) + // rather than in `lengthOfTimeWindow`. If the primary source did not + // populate `misc.lengthOfTimeWindow`, read the 4D-Var key directly. + // + // No derivation from `anoffset` is performed: we only mirror what the + // producer placed in the input handle. The ecCodes "missing" sentinel + // (0xFFFF for a 16-bit unsigned field) is explicitly ignored so that + // metkit's analysis deduction keeps writing GRIB "missing" downstream. + if (!misc.lengthOfTimeWindow.isSet() && h.has("lengthOf4DvarWindow")) { + long lengthOf4DvarWindowHours = h.getLong("lengthOf4DvarWindow"); + if (lengthOf4DvarWindowHours != 0xFFFF) { + misc.lengthOfTimeWindow.set(lengthOf4DvarWindowHours); + } + } misc.bitsPerValue = dm::parseEntry(dm::BitsPerValue.withKey("bitsPerValue").withDefault(24), h); // TODO pgeier readd maybe @@ -419,6 +547,7 @@ void mapGrib1ToGrib2(KeySet& marsKeys, metkit::codes::CodesHandle& h, dm::FullMa } misc.satelliteSeries = dm::parseEntry(dm::SatelliteSeries.withKey("satelliteSeries"), h); + misc.numberOfFrequencies = dm::parseEntry(dm::NumberOfFrequencies.withKey("numberOfFrequencies"), h); misc.scaleFactorOfCentralWaveNumber = dm::parseEntry(dm::ScaleFactorOfCentralWaveNumber.withKey("scaleFactorOfCentralWaveNumber"), h); misc.scaledValueOfCentralWaveNumber @@ -795,6 +924,30 @@ Discipline192Handling parseDiscipline192Handling(const std::string& str) { } +enum class OnErrorHandling : std::size_t +{ + Abort, + LogAndSkip, + Skip, +}; + +const std::unordered_map& onErrorHandlingMap() { + static const std::unordered_map map{ + {"abort", OnErrorHandling::Abort}, + {"log-and-skip", OnErrorHandling::LogAndSkip}, + {"skip", OnErrorHandling::Skip}}; + return map; +} + +OnErrorHandling parseOnErrorHandling(const std::string& str) { + const auto& map = onErrorHandlingMap(); + if (auto search = map.find(str); search != map.end()) { + return search->second; + } + throw std::runtime_error(std::string("Unsupported --on-error value: ") + str); +} + + class Grib1ToGrib2 final : public multio::MultioTool { public: // methods Grib1ToGrib2(int argc, char** argv); @@ -826,6 +979,7 @@ class Grib1ToGrib2 final : public multio::MultioTool { long ncycle_ = 0; std::optional excludeMap_ = {}; + std::optional exceptMap_ = {}; std::optional filterMap_ = {}; std::optional overwritePacking_ = {}; std::optional setModel_ = {}; @@ -835,6 +989,8 @@ class Grib1ToGrib2 final : public multio::MultioTool { std::optional> mappingRules_ = mars2mars::allRulesNoWMOMapping(); Discipline192Handling discipline192Handling_ = Discipline192Handling::LogAndIgnore; + OnErrorHandling onErrorHandling_ = OnErrorHandling::LogAndSkip; + long defaultEnsembleSize_ = 0; }; Grib1ToGrib2::Grib1ToGrib2(int argc, char** argv) : multio::MultioTool{argc, argv} { @@ -858,6 +1014,11 @@ Grib1ToGrib2::Grib1ToGrib2(int argc, char** argv) : multio::MultioTool{argc, arg "exclude", "Keys and values to be excluded. Multiple values are separated by ','. Multiple key-values pairs are separated " "by ';'. Example --exclude paramId=130,131,133;levtype=pl,sfc")); + options_.push_back(new eckit::option::SimpleOption( + "except", + "Keys and values to be copied verbatim (without re-encoding) when --all is active. Same syntax as --exclude. " + "Only applies to GRIB2 messages; matching a GRIB1 message is an error. " + "Example --except paramId=213131")); options_.push_back(new eckit::option::SimpleOption( "filter", "Keys and values to be included. Multiple values are separated by ','. Multiple key-values pairs are separated " @@ -871,6 +1032,13 @@ Grib1ToGrib2::Grib1ToGrib2(int argc, char** argv) : multio::MultioTool{argc, arg "discipline-192", "Options on handling fields with discipline 192 (field that are ill-formed). Values: \"log-and-ignore\" " "(default), \"ignore\", \"try-to-handle\"")); + options_.push_back(new eckit::option::SimpleOption( + "on-error", + "How to handle per-message conversion errors. Values: \"abort\" (stop on first error), \"log-and-skip\" " + "(default, log the error and continue), \"skip\" (silently skip failing messages)")); + options_.push_back(new eckit::option::SimpleOption( + "default-ensemble-size", + "Fallback value used when numberOfForecastsInEnsemble is 0 but number is non-zero. Default: 0 (throw)")); } void Grib1ToGrib2::init(const eckit::option::CmdArgs& args) { @@ -917,6 +1085,16 @@ void Grib1ToGrib2::init(const eckit::option::CmdArgs& args) { excludeMap_ = parseFieldValueMap(std::move(excludeStr), verbosity_); } + std::string exceptStr = ""; + args.get("except", exceptStr); + if (!exceptStr.empty()) { + if (copyGrib2Messages_) { + std::cerr << "Warning: --except has no effect without --all (GRIB2 messages are already copied verbatim)" + << std::endl; + } + exceptMap_ = parseFieldValueMap(std::move(exceptStr), verbosity_); + } + std::string filterStr = ""; args.get("filter", filterStr); if (!filterStr.empty()) { @@ -928,6 +1106,14 @@ void Grib1ToGrib2::init(const eckit::option::CmdArgs& args) { if (!discipline192.empty()) { discipline192Handling_ = parseDiscipline192Handling(discipline192); } + + std::string onError; + args.get("on-error", onError); + if (!onError.empty()) { + onErrorHandling_ = parseOnErrorHandling(onError); + } + + args.get("default-ensemble-size", defaultEnsembleSize_); } void Grib1ToGrib2::finish(const eckit::option::CmdArgs&) {} @@ -975,7 +1161,11 @@ void Grib1ToGrib2::execute(const eckit::option::CmdArgs& args) { metkit::mars2grib::Mars2Grib encoder{}; eckit::message::Message msg; + std::size_t msgIndex = 0; + std::size_t skippedCount = 0; while ((msg = reader.next())) { + ++msgIndex; + try { // Extract message from datahandle... we expect it to be a memory handle // TODO pgeier: Alternative would be to explicitly create a eckit::MemoryHandle and write to it std::unique_ptr dh{msg.readHandle()}; @@ -1014,6 +1204,11 @@ void Grib1ToGrib2::execute(const eckit::option::CmdArgs& args) { bool isDiscipline192 = (edition == "1") ? isDiscipline192Param(paramId) : (inputHandle->getLong("discipline") == 192); if (isDiscipline192) { + if (exceptMap_ && matches(msg, *exceptMap_, verbosity_)) { + std::cerr << "Warning: --except matched a discipline-192 message (paramId=" << paramId + << ") but it was skipped by --discipline-192 policy. " + << "Use --discipline-192 try-to-handle to allow --except to take effect." << std::endl; + } if (discipline192Handling_ == Discipline192Handling::LogAndIgnore) { std::cout << "Excluding message with discipline 192 (paramId: " << paramId << ")" << std::endl; } @@ -1021,6 +1216,24 @@ void Grib1ToGrib2::execute(const eckit::option::CmdArgs& args) { } } + // --except: copy matching GRIB2 messages verbatim instead of re-encoding. + // Only meaningful with --all; GRIB1 matches are an error. + if (exceptMap_ && matches(msg, *exceptMap_, verbosity_)) { + if (edition == "1") { + throw eckit::BadValue(std::string("--except matched a GRIB1 message (paramId=") + + std::to_string(inputHandle->getLong("paramId")) + + "). --except may only match GRIB2 messages.", + Here()); + } + if (verbosity_ >= 1) { + std::cout << "except map matched — copying GRIB2 message verbatim" << std::endl; + } + if (outputFileHandle) { + write(*inputHandle.get(), *outputFileHandle); + } + continue; + } + if (edition == "2" && copyGrib2Messages_) { // Write the message directly if (verbosity_ > 2) { @@ -1041,7 +1254,7 @@ void Grib1ToGrib2::execute(const eckit::option::CmdArgs& args) { std::cout << "Extracting metadata..." << std::endl; } - extract::mapGrib1ToGrib2(marsKeys, *inputHandle.get(), mars, misc, verbosity_); + extract::mapGrib1ToGrib2(marsKeys, *inputHandle.get(), mars, misc, verbosity_, defaultEnsembleSize_); if (overwritePacking_) { @@ -1121,6 +1334,11 @@ void Grib1ToGrib2::execute(const eckit::option::CmdArgs& args) { const auto marsConfig = dm::dumpRecord(mars); const auto miscConfig = dm::dumpUnscopedRecord(misc); + // Pre-encode validation: catch spectral_complex laplacian-scaling + // overflows before they reach ecCodes (where they would trip an + // unrecoverable assertion in `grib_ieee_to_long`). + extract::validateSpectralComplexNoOverflow(mars, misc, values); + // Call the GRIB2 encoder in metkit auto preparedHandle = encoder.encode(values, marsConfig, miscConfig); @@ -1143,6 +1361,23 @@ void Grib1ToGrib2::execute(const eckit::option::CmdArgs& args) { write(*preparedHandle.get(), *outputFileHandle); } } + } + catch (const std::exception& e) { + if (onErrorHandling_ == OnErrorHandling::Abort) { + throw; + } + ++skippedCount; + if (onErrorHandling_ == OnErrorHandling::LogAndSkip) { + std::cerr << "Error converting message #" << msgIndex << ": " << e.what() + << " -- skipping" << std::endl; + } + continue; + } + } + + if (skippedCount > 0) { + std::cerr << "grib1-to-grib2: skipped " << skippedCount << " message(s) due to conversion errors" + << std::endl; } if (outputFileHandle) {