Skip to content

Fix ocean encoding for CY49r2/CY50r1 #234

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
dsarmany merged 10 commits into from
Jan 28, 2026
Merged

Fix ocean encoding for CY49r2/CY50r1 #234

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
99086a8
Fix significanceOfReferenceTime for ocean in ERA6 production
MircoValentiniECMWF Jan 16, 2026
e950151
Fix ocean level
MircoValentiniECMWF Jan 16, 2026
245536a
Fix productionStatusOfProcessedData for Re-analysis
MircoValentiniECMWF Jan 21, 2026
7ae83cf
Fix subCentre
MircoValentiniECMWF Jan 21, 2026
8d614fd
Fix/Hack for float values encoding (Span not defined for float array)
MircoValentiniECMWF Jan 21, 2026
834073f
Add logic for significanceOfReferenceTime active only for type=e6
MircoValentiniECMWF Jan 22, 2026
741dce6
Fix *AfterDataCutOff and backgroundProcess
MircoValentiniECMWF Jan 22, 2026
bbe51b3
Fix typeOfProcessedData for type=4v
MircoValentiniECMWF Jan 23, 2026
e878a9a
Make CI happy until the reference data is properly updated
dsarmany Jan 28, 2026
b5982e6
Bump version number
dsarmany Jan 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion VERSION
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1 +1 @@
2.8.0
2.8.2
2 changes: 1 addition & 1 deletion src/multio/action/encode/Encode.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -154,7 +154,7 @@ std::unique_ptr<GribEncoder> makeEncoder(const eckit::LocalConfiguration& conf,
}

std::string encodingExceptionReason(const std::string& r) {
std::string s("Enocding exception: ");
std::string s("Encoding exception: ");
s.append(r);
return s;
}
Expand Down
134 changes: 111 additions & 23 deletions src/multio/action/encode/GribEncoder.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,10 +14,12 @@

#include "GribEncoder.h"

#include <algorithm>
#include <cstring>
#include <functional>
#include <iomanip>
#include <iostream>
#include <optional>
#include <unordered_map>
#include <unordered_set>

Expand Down Expand Up @@ -59,6 +61,7 @@ const std::map<const std::string, const std::int64_t> type_of_statistical_proces
{"average", 0}, {"accumulate", 1}, {"maximum", 2}, {"minimum", 3},
{"difference", 4}, {"stddev", 6}, {"inverse-difference", 8}};


const std::map<const std::string, const std::string> category_to_levtype{
{"ocean-grid-coordinate", "oceanSurface"}, {"ocean-2d", "oceanSurface"}, {"ocean-3d", "oceanModelLevel"}};

Expand Down Expand Up @@ -285,18 +288,24 @@ void setMissingFixedSurface(GribEncoder& g, const std::string& typeOfLevel, long
g.setMissing(dm::legacy::ScaledValueOfSecondFixedSurface);
}

void setOceanSurface(GribEncoder& g, const std::string& typeOfLevel, long level) {
g.setValue("typeOfLevel", typeOfLevel);

g.setMissing(dm::legacy::ScaleFactorOfSecondFixedSurface);
g.setMissing(dm::legacy::ScaledValueOfSecondFixedSurface);
}

using TypeOfLevelSetter = std::function<void(GribEncoder&, const std::string&, long)>;

const std::map<std::string, TypeOfLevelSetter> typeOfLevelSetters{
{"snowLayer", &setLayerTypeOfLevel},
{"soilLayer", &setSoilLayerTypeOfLevel},
{"seaIceLayer", &setLayerTypeOfLevel},
{"mediumCloudLayer", &setLevelUnrelatedTypeOfLevel},
{"lowCloudLayer", &setLevelUnrelatedTypeOfLevel},
{"highCloudLayer", &setLevelUnrelatedTypeOfLevel},
{"meanSea", &setLevelUnrelatedTypeOfLevel},
{"iceLayerOnWater", &setMissingFixedSurface},
};
const std::map<std::string, TypeOfLevelSetter> typeOfLevelSetters{{"snowLayer", &setLayerTypeOfLevel},
{"soilLayer", &setSoilLayerTypeOfLevel},
{"seaIceLayer", &setLayerTypeOfLevel},
{"mediumCloudLayer", &setLevelUnrelatedTypeOfLevel},
{"lowCloudLayer", &setLevelUnrelatedTypeOfLevel},
{"highCloudLayer", &setLevelUnrelatedTypeOfLevel},
{"meanSea", &setLevelUnrelatedTypeOfLevel},
{"iceLayerOnWater", &setMissingFixedSurface},
{"oceanSurface", &setOceanSurface}};

template <typename Dict>
QueriedMarsKeys setMarsKeys(GribEncoder& g, const Dict& md) {
Expand All @@ -309,6 +318,9 @@ QueriedMarsKeys setMarsKeys(GribEncoder& g, const Dict& md) {
const auto levtype = lookUp<std::string>(md, dm::legacy::Levtype)();
const auto gribEdition = lookUp<std::string>(md, dm::legacy::GribEdition)().value_or("2");

// Initialize subCenter to 0 -> ERA6 request
g.setValue(dm::legacy::SubCentre, 0);

if ((gribEdition == "2") && (gridType != "sh")) {
withFirstOf(valueSetter(g, dm::legacy::SetPackingType), lookUp<std::string>(md, dm::legacy::SetPackingType));
}
Expand Down Expand Up @@ -363,12 +375,25 @@ QueriedMarsKeys setMarsKeys(GribEncoder& g, const Dict& md) {
}
}

const auto productionStatusOfProcessedData
= lookUp<std::int64_t>(md, dm::legacy::ProductionStatusOfProcessedData)();
auto productionStatusOfProcessedData = lookUp<std::int64_t>(md, dm::legacy::ProductionStatusOfProcessedData)();

// ERA6 hack
const auto marsClass = lookUp<std::string>(md, dm::legacy::ClassKey)();
if (marsClass == "e6") { // Re-analysis

auto type = firstOf(lookUp<std::string>(md, dm::legacy::Type), lookUp<std::string>(md, dm::legacy::MarsType));
g.setValue("backgroundProcess", 255);
productionStatusOfProcessedData = 3;
if ( type && (*type == "an" || *type == "4v")) {
g.setValue("typeOfProcessedData", 0);
}
}

if (productionStatusOfProcessedData) {
g.setValue(dm::legacy::ProductionStatusOfProcessedData, *productionStatusOfProcessedData);
}


ret.paramId = firstOf(
lookUp<std::int64_t>(md, dm::legacy::ParamId),
lookUpTranslate<std::int64_t>(md, dm::legacy::Param)); // param might be a string, separated by . for GRIB1.
Expand Down Expand Up @@ -668,6 +693,33 @@ std::string getTimeReference(GribEncoder& g, const message::Metadata& md, const
return isReferingToStart ? "start" : (isTimeRange ? "previous" : "current");
}

// Hack introduced for ERA6 logic is copied from the new encoder
std::optional<std::int64_t> significanceOfReferenceTimeFromType(const std::string& marsType) {
const std::array<std::string_view, 17> analysisTypes
= {{"an", "ia", "oi", "3v", "3g", "4g", "ea", "pa", "tpa", "ga", "gai", "ai", "af", "ab", "oai", "ga", "gai"}};

const std::array<std::string_view, 4> startOfDataAssimilationTypes = {{"4i", "4v", "me", "eme"}};

const std::array<std::string_view, 32> forecastTypes
= {{"fc", "cf", "pf", "cm", "fp", "em", "es", "fa", "efi", "efic", "bf",
"cd", "wem", "wes", "cr", "ses", "taem", "taes", "sg", "sf", "if", "fcmean",
"fcmax", "fcmin", "fcstdev", "ssd", "tf", "bf", "cd", "hcmean", "s3", "si"}};

if (std::any_of(analysisTypes.begin(), analysisTypes.end(), [&marsType](auto v) { return marsType == v; })) {
return std::optional{static_cast<std::int64_t>(0)};
}
else if (std::any_of(forecastTypes.begin(), forecastTypes.end(), [&marsType](auto v) { return marsType == v; })) {
return std::optional{static_cast<std::int64_t>(1)};
}
else if (std ::any_of(startOfDataAssimilationTypes.begin(), startOfDataAssimilationTypes.end(),
[&marsType](auto v) { return marsType == v; })) {
return std::optional{static_cast<std::int64_t>(6)};
}

return std::nullopt;
}


void setDateAndStatisticalFields(GribEncoder& g, const message::Metadata& in,
const QueriedMarsKeys& queriedMarsFields) {
message::Metadata md = in; // Copy to allow modification
Expand All @@ -676,6 +728,13 @@ void setDateAndStatisticalFields(GribEncoder& g, const message::Metadata& in,
// std::string forecastTimeKey = gribEdition == "2" ? "forecastTime" : "startStep";


std::optional<std::string> marsType
= firstOf(lookUp<std::string>(md, dm::legacy::Type), lookUp<std::string>(md, dm::legacy::MarsType),
lookUp<std::string>(md, "marsType"));
std::optional<std::string> marsClass
= firstOf(lookUp<std::string>(md, dm::legacy::ClassKey), lookUp<std::string>(md, "marsClass"));


auto operation = lookUp<std::string>(md, "operation")();
auto startStep = lookUp<std::int64_t>(md, "startStep")();
auto endStep = lookUp<std::int64_t>(md, "endStep")();
Expand All @@ -691,6 +750,11 @@ void setDateAndStatisticalFields(GribEncoder& g, const message::Metadata& in,
significanceOfReferenceTime = lookUp<std::int64_t>(overwrites, "significanceOfReferenceTime")();
}
}
if (marsType.has_value() && marsClass.has_value() && marsClass.value() == "e6") {
significanceOfReferenceTime = significanceOfReferenceTimeFromType(marsType.value());
}


if ((gribEdition == "2") && significanceOfReferenceTime) {
g.setValue("significanceOfReferenceTime", *significanceOfReferenceTime);
}
Expand All @@ -714,6 +778,11 @@ void setDateAndStatisticalFields(GribEncoder& g, const message::Metadata& in,
auto currentDateTime = util::wrapDateTime({util::toDateInts(md.get<std::int64_t>(dm::legacy::CurrentDate)),
util::toTimeInts(md.get<std::int64_t>(dm::legacy::CurrentTime))});

if (g.hasKey("hoursAfterDataCutoff") && g.hasKey("minutesAfterDataCutoff")) {
g.setMissing("hoursAfterDataCutoff");
g.setMissing("minutesAfterDataCutoff");
}

if (!isTimeRange) {
if (timeRef == std::string("start")) {
// Compute diff to current time in some appropriate unit
Expand Down Expand Up @@ -778,7 +847,8 @@ void setDateAndStatisticalFields(GribEncoder& g, const message::Metadata& in,
if (searchStat == std::end(type_of_statistical_processing)) {
std::ostringstream oss;
oss << "setDateAndStatisticalFields - Cannot map value \"" << *operation
<< "\"for key \"operation\" (statistical output) to a valid grib2 type of statistical processing.";
<< "\"for key \"operation\" (statistical output) to a valid grib2 type of statistical "
"processing.";
throw eckit::UserError(oss.str(), Here());
}
g.setValue("typeOfStatisticalProcessing", searchStat->second);
Expand All @@ -788,21 +858,21 @@ void setDateAndStatisticalFields(GribEncoder& g, const message::Metadata& in,
// # CODE TABLE 4.11, Type of time intervals
// 1 1 Successive times processed have same forecast time, start time of forecast is incremented
// 2 2 Successive times processed have same start time of forecast, forecast time is incremented
// 3 3 Successive times processed have start time of forecast incremented and forecast time decremented so that
// valid time remains constant 4 4 Successive times processed have start time of forecast decremented and
// forecast time incremented so that valid time remains constant 5 5 Floating subinterval of time between
// forecast time and end of overall time interval
// g.setValue("typeOfTimeIncrement", (timeRef == "start" ? 2 : 1));
// 3 3 Successive times processed have start time of forecast incremented and forecast time decremented so
// that valid time remains constant 4 4 Successive times processed have start time of forecast decremented
// and forecast time incremented so that valid time remains constant 5 5 Floating subinterval of time
// between forecast time and end of overall time interval g.setValue("typeOfTimeIncrement", (timeRef ==
// "start" ? 2 : 1));
//
// #### Work around ####
// Eccodes has problems showing stepRange correctly for averaging fields with typeOfTimeIncrement == 1.
// It seems that from this combination eccodes is infering a `stepKey` of avgd (daily average).
// For daily average the stepRange is shown as 0 instead of 0-24 (desired). Hence with DGov we decided to put
// 255 (MISSING) as typeOfTimeIncrement
// For daily average the stepRange is shown as 0 instead of 0-24 (desired). Hence with DGov we decided to
// put 255 (MISSING) as typeOfTimeIncrement
//
// TO BE DISCUSSED - obviously there is some confusion about typeOfTimeIncrement=1. For analysis I read that it
// should be set to 1. However eccodes thinks different and will not consider it as time range then... hence I
// explicily set it to 255 now g.setValue(
// TO BE DISCUSSED - obviously there is some confusion about typeOfTimeIncrement=1. For analysis I read that
// it should be set to 1. However eccodes thinks different and will not consider it as time range then...
// hence I explicily set it to 255 now g.setValue(
// "typeOfTimeIncrement",
// (timeRef == "start" ? 2
// : ((gribEdition == "2") && (significanceOfReferenceTime &&
Expand Down Expand Up @@ -1029,6 +1099,24 @@ message::Message GribEncoder::setFieldValues(message::Message&& msg) {
}


template <>
inline void GribEncoder::setDataValues<double>(const double* data, size_t count) {
encoder_->set("values", metkit::codes::Span<const double>(data, count));
}

template <>
inline void GribEncoder::setDataValues<float>(const float* data, size_t count) {
std::vector<double> tmp;
tmp.reserve(count);

for (size_t i = 0; i < count; ++i) {
tmp.push_back(static_cast<double>(data[i]));
}

setDataValues<double>(tmp.data(), tmp.size());
}


void GribEncoder::print(std::ostream& os) const {
os << "GribEncoder(config=" << config_ << ")";
};
Expand Down
8 changes: 4 additions & 4 deletions src/multio/action/encode/GribEncoder.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ class GribEncoder {
void setValue(const std::string& key, const std::vector<T>& v) {
encoder_->set(key, v);
};

// Dirty implementation - before metkit::codes::CodesHandle
// was used, the former implementation explicitly ignored key setting errors
// when the key was read-only.
Expand All @@ -80,10 +80,10 @@ class GribEncoder {
<< ") failed: " << std::endl;
}
}

template <typename T>
void setDataValues(const T* data, size_t count) {
encoder_->set("values", metkit::codes::Span<const T>(reinterpret_cast<const T*>(data), count));
};
void setDataValues(const T* data, size_t count) = delete;

bool hasKey(const char* key);

message::Message encodeOceanCoordinates(message::Message&& msg, const message::Metadata& additionalMetadata);
Expand Down
4 changes: 2 additions & 2 deletions tests/multio/action/statistics/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,8 +80,8 @@ foreach(_test ${_statistics_tests_names} )
ecbuild_add_test(
TARGET ${PREFIX}_check_metadata_${_kind}_${_start}_${_test}
TEST_REQUIRES ${PREFIX}_run_${_kind}_${_start}_${_test}
ARGS -H ${_reference} ${_result}
COMMAND grib_compare
ARGS -H -b hoursAfterDataCutoff,minutesAfterDataCutoff ${_reference} ${_result}
)

ecbuild_add_test(
Expand All @@ -102,8 +102,8 @@ foreach(_test ${_statistics_tests_names} )
ecbuild_add_test(
TARGET ${PREFIX}_check_metadata_${_kind}_${_start}_${_test}
TEST_REQUIRES ${PREFIX}_run_${_kind}_${_start}_${_test}
ARGS -H ${_reference} ${_result}
COMMAND grib_compare
ARGS -H -b hoursAfterDataCutoff,minutesAfterDataCutoff ${_reference} ${_result}
)

ecbuild_add_test(
Expand Down
Loading