Skip to content

Commit

Permalink
Bugfix/sla do not assimilate and _FillValue (#56)
Browse files Browse the repository at this point in the history 
* Add tests and do-not-assimilate for finalReject QC
* use enum class for clarity
* Fillvalue normalisation
  • Loading branch information
@twsearle
twsearle authored Jul 13, 2023
1 parent b78e395 commit 7730763
Show file tree
Hide file tree
Showing 16 changed files with 284 additions and 127 deletions.
129 changes: 54 additions & 75 deletions src/nemo-feedback/NemoFeedback.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@

#include "nemo-feedback/NemoFeedback.h"

#include <string>
#include <string_view>
#include <string>

#include <algorithm>
#include <utility>
Expand Down Expand Up @@ -46,35 +46,6 @@ namespace nemo_feedback {
constexpr std::string_view defaultDepthGroup{"MetaData"};
constexpr std::string_view defaultDepthVariable{"depthBelowWaterSurface"};

/// \brief create whole report variables from a diagnostic flag
void wholeReportQCFromDiagnosticFlags(const
NemoFeedbackDataCreator& creator, const std::string& group, const
std::string& name, feedback_io::Data<int32_t>& wholeReportQCData) {
if (wholeReportQCData.n_obs() == 0) {
wholeReportQCData = feedback_io::Data<int32_t>(creator.indexer(),
std::vector<int32_t>(creator.indexer()->n_source_data(), 0));
}
feedback_io::Data<int32_t> QCData =
creator.create(group, name,
ufo::DiagnosticFlag(0), 4, 1);
for (size_t iProfile = 0;
iProfile < QCData.n_obs(); ++iProfile) {
size_t nBadObs = 0;
for (size_t iLevel = 0;
iLevel < QCData.length(iProfile);
++iLevel) {
if (QCData(iProfile, iLevel) == 4) {
++nBadObs;
}
}
if (wholeReportQCData[iProfile] == 0 ||
wholeReportQCData[iProfile] == 4) {
wholeReportQCData[iProfile] =
(nBadObs == QCData.length(iProfile) ? 4 : 1);
}
}
}


NemoFeedback::NemoFeedback(
ioda::ObsSpace & obsdb,
Expand Down Expand Up @@ -222,10 +193,11 @@ template <typename T>
void NemoFeedback::write_all_data(feedback_io::Writer<T>& writer,
const NemoFeedbackDataCreator& creator) const
{
feedback_io::Data<int32_t> wholeReportQCData(creator.indexer(),
std::vector<int32_t>(creator.indexer()->n_source_data(), 0));
feedback_io::Data<int32_t> wholeReportPositionQCData;
feedback_io::Data<int32_t> wholeReportTimeQCData;
feedback_io::Data<feedback_io::QC::Level> wholeReportQCData(creator.indexer(),
std::vector<feedback_io::QC::Level>(creator.indexer()->n_source_data(),
feedback_io::QC::Level::None));
feedback_io::Data<feedback_io::QC::Level> wholeReportPositionQCData;
feedback_io::Data<feedback_io::QC::Level> wholeReportTimeQCData;
std::vector<ufo::DiagnosticFlag> do_not_assimilate;
for (const NemoFeedbackVariableParameters& nemoVariableParams :
parameters_.variables.value()) {
Expand Down Expand Up @@ -270,70 +242,78 @@ void NemoFeedback::write_all_data(feedback_io::Writer<T>& writer,

// Quality control rank variables
if (obsdb_.has("DiagnosticFlags/FinalReject", ufo_name)) {

feedback_io::Data<int32_t> variableFinalQCData =
feedback_io::Data<feedback_io::QC::Level> variableFinalQCData =
creator.create("DiagnosticFlags/FinalReject", ufo_name,
ufo::DiagnosticFlag(0), 4, 1);
ufo::DiagnosticFlag(0), feedback_io::QC::Level::Bad,
feedback_io::QC::Level::Good);

// Add do not assimilate flag if required.
// Add do not assimilate flag if required to the final QC information.
if (obsdb_.has("DiagnosticFlags/DoNotAssimilate", ufo_name)) {
feedback_io::Data<int32_t> variableDoNotAssimilateData(
feedback_io::Data<feedback_io::QC::Level> variableDoNotAssimilateData(
creator.create("DiagnosticFlags/DoNotAssimilate", ufo_name,
ufo::DiagnosticFlag(0), 128, 0));
ufo::DiagnosticFlag(0), feedback_io::QC::Level::DoNotAssimilate,
feedback_io::QC::Level::None));
for (size_t iProfile = 0;
iProfile < variableDoNotAssimilateData.n_obs(); ++iProfile) {
for (size_t iLevel = 0;
iLevel < variableDoNotAssimilateData.length(iProfile);
++iLevel) {
variableFinalQCData(iProfile, iLevel) +=
variableDoNotAssimilateData(iProfile, iLevel);
if (variableDoNotAssimilateData(iProfile, iLevel) ==
feedback_io::QC::Level::DoNotAssimilate) {
variableFinalQCData(iProfile, iLevel) =
static_cast<feedback_io::QC::Level>(
static_cast<int32_t>(variableFinalQCData(iProfile, iLevel)) +
static_cast<int32_t>(feedback_io::QC::Level::DoNotAssimilate));
}
}
}
}

// Per-profile quality rank
feedback_io::Data<int32_t> wholeVariableQCData(creator.indexer(),
std::vector<int32_t>(creator.indexer()->n_source_data(), 0));
for (size_t iProfile = 0;
iProfile < variableFinalQCData.n_obs(); ++iProfile) {
size_t nBadObs = 0;
for (size_t iLevel = 0;
iLevel < variableFinalQCData.length(iProfile);
++iLevel) {
if (variableFinalQCData(iProfile, iLevel) == 4) {
++nBadObs;
}
}
wholeVariableQCData[iProfile] =
(nBadObs == variableFinalQCData.length(iProfile) ? 4 : 1);
}
// set per-profile quality rank for the variable
feedback_io::Data<feedback_io::QC::Level> wholeVariableQCData(
creator.indexer(), std::vector<feedback_io::QC::Level>(
creator.indexer()->n_source_data(), feedback_io::QC::Level::None));
feedback_io::wholeReportFromPerProfile(variableFinalQCData,
wholeVariableQCData);
// update set per-profile quality rank for the whole report
feedback_io::wholeReportFromPerProfile(variableFinalQCData,
wholeReportQCData);

writer.write_variable_surf_qc(nemo_name + "_QC",
wholeVariableQCData);
writer.write_variable_level_qc(nemo_name + "_LEVEL_QC",
variableFinalQCData);

// Whole Observation report QC
// Whole Observation Position report QC
const std::string positionQCGroup("DiagnosticFlags/PositionReject");
if (obsdb_.has(positionQCGroup, ufo_name)) {
wholeReportQCFromDiagnosticFlags(creator, positionQCGroup, ufo_name,
feedback_io::Data<feedback_io::QC::Level> QCData = creator.create(
positionQCGroup, ufo_name, ufo::DiagnosticFlag(0),
feedback_io::QC::Level::Bad, feedback_io::QC::Level::Good);
if (wholeReportPositionQCData.n_obs() == 0) {
wholeReportPositionQCData = feedback_io::Data<feedback_io::QC::Level>(
creator.indexer(), std::vector<feedback_io::QC::Level>(
creator.indexer()->n_source_data(),
feedback_io::QC::Level::None));
}
feedback_io::wholeReportFromPerProfile(QCData,
wholeReportPositionQCData);
}

// Whole Observation time report QC
const std::string timeQCGroup("DiagnosticFlags/TimeReject");
if (obsdb_.has(timeQCGroup, ufo_name)) {
wholeReportQCFromDiagnosticFlags(creator, timeQCGroup, ufo_name,
wholeReportTimeQCData);
}

{
for (size_t iProfile = 0; iProfile < variableData.n_obs(); ++iProfile) {
if (wholeReportQCData[iProfile] == 0 ||
wholeReportQCData[iProfile] == 4) {
wholeReportQCData[iProfile] =
(4 == wholeVariableQCData[iProfile] ? 4 : 1);
}
feedback_io::Data<feedback_io::QC::Level> QCData =
creator.create(timeQCGroup, ufo_name, ufo::DiagnosticFlag(0),
feedback_io::QC::Level::Bad, feedback_io::QC::Level::Good);
if (wholeReportTimeQCData.n_obs() == 0) {
wholeReportTimeQCData = feedback_io::Data<feedback_io::QC::Level>(
creator.indexer(), std::vector<feedback_io::QC::Level>(
creator.indexer()->n_source_data(),
feedback_io::QC::Level::None));
}
feedback_io::wholeReportFromPerProfile(QCData,
wholeReportTimeQCData);
}
}

Expand All @@ -356,9 +336,9 @@ void NemoFeedback::write_all_data(feedback_io::Writer<T>& writer,
const std::string depthVariable = parameters_.depthVariable.value()
.value_or(static_cast<std::string>(defaultDepthVariable));
if (obsdb_.has(depthQCGroup, depthVariable)) {
feedback_io::Data<int32_t> depthQCData(
creator.create(depthQCGroup, depthVariable,
ufo::DiagnosticFlag(0), 4, 1));
feedback_io::Data<feedback_io::QC::Level> depthQCData(
creator.create(depthQCGroup, depthVariable, ufo::DiagnosticFlag(0),
feedback_io::QC::Level::Bad, feedback_io::QC::Level::Good));
writer.write_variable_level_qc("DEPTH_QC", depthQCData);
}
if (wholeReportPositionQCData.n_obs() != 0) {
Expand Down Expand Up @@ -424,7 +404,6 @@ std::tuple<feedback_io::Data<std::string>, feedback_io::Data<std::string>>
stationIdentificationAvailable = true;
}

const int32_t buoyIDMissingValueInt = util::missingValue<int32_t>();
constexpr size_t buoyIDWidth = 8;
const std::string missingStringFeedback =
feedback_io::typeToFill::value<std::string>();
Expand Down
20 changes: 11 additions & 9 deletions src/nemo-feedback/NemoFeedbackDataCreator.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,14 +144,15 @@ feedback_io::Data<std::string> NemoFeedbackDataCreator::create_from_obsdb(
return feedback_io::Data<std::string>(indexer_, std::move(data));
}

feedback_io::Data<int32_t> NemoFeedbackDataCreator::create_from_obsdb(const
std::string& obsGroup, const std::string& ufoName, const
ufo::DiagnosticFlag TypeInstance, int32_t whenTrue,
int32_t whenFalse) const {
feedback_io::Data<feedback_io::QC::Level>
NemoFeedbackDataCreator::create_from_obsdb(
const std::string& obsGroup, const std::string& ufoName, const
ufo::DiagnosticFlag TypeInstance, feedback_io::QC::Level whenTrue,
feedback_io::QC::Level whenFalse) const {
oops::Log::trace() << NemoFeedbackDataCreator::className()
<< ":create_from_obsdb ufo::DiagnosticFlag "
<< obsGroup << "/" << ufoName << std::endl;
std::vector<int32_t> data;
std::vector<feedback_io::QC::Level> data;
std::vector<ufo::DiagnosticFlag> flagData(obsdb_.nlocs(), 0);
obsdb_.get_db(obsGroup, ufoName, flagData);
for (ufo::DiagnosticFlag flag : flagData) {
Expand All @@ -161,7 +162,7 @@ feedback_io::Data<int32_t> NemoFeedbackDataCreator::create_from_obsdb(const
throw eckit::BadValue(NemoFeedbackDataCreator::className()
+ ":create_from_obsdb ufo::DiagnosticFlag no data.", Here());

return feedback_io::Data<int32_t>(indexer_, std::move(data));
return feedback_io::Data<feedback_io::QC::Level>(indexer_, std::move(data));
}

template<typename T>
Expand Down Expand Up @@ -304,9 +305,10 @@ NemoFeedbackDataCreator::create_altimeter_IDs() const {
feedback_io::Data<std::string>(indexer_, std::move(station_types)));
}

feedback_io::Data<int32_t> NemoFeedbackDataCreator::create(const std::string&
obsGroup, const std::string& ufoName, const ufo::DiagnosticFlag
typeInstance, int32_t whenTrue, int32_t whenFalse) const {
feedback_io::Data<feedback_io::QC::Level> NemoFeedbackDataCreator::create(
const std::string& obsGroup, const std::string& ufoName, const
ufo::DiagnosticFlag typeInstance, feedback_io::QC::Level whenTrue,
feedback_io::QC::Level whenFalse) const {
if (obsGroup == "hofx" || obsGroup == "HofX") {
std::ostringstream err_stream;
err_stream << NemoFeedbackDataCreator::className()
Expand Down
13 changes: 7 additions & 6 deletions src/nemo-feedback/NemoFeedbackDataCreator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,11 +63,11 @@ class NemoFeedbackDataCreator {
create_altimeter_IDs() const;

/// \brief create DiagnosticFlags feedbackData
feedback_io::Data<int32_t> create(const std::string& obsGroup,
feedback_io::Data<feedback_io::QC::Level> create(const std::string& obsGroup,
const std::string& ufoName,
const ufo::DiagnosticFlag typeInstance,
const int32_t whenTrue,
const int32_t whenFalse) const;
const feedback_io::QC::Level whenTrue,
const feedback_io::QC::Level whenFalse) const;

/// \brief create feedbackData
template<typename T>
Expand All @@ -91,11 +91,12 @@ std::shared_ptr<feedback_io::DataIndexer> indexer() const {return indexer_;}
size_t width, bool leftJustify = false) const;

/// \brief create DiagnosticFlags feedbackData from obsdb
feedback_io::Data<int32_t> create_from_obsdb(const std::string& obsGroup,
feedback_io::Data<feedback_io::QC::Level> create_from_obsdb(
const std::string& obsGroup,
const std::string& ufoName,
const ufo::DiagnosticFlag TypeInstance,
int32_t whenTrue,
int32_t whenFalse) const;
feedback_io::QC::Level whenTrue,
feedback_io::QC::Level whenFalse) const;

/// \brief create feedbackData from obsdb
template<typename T>
Expand Down
64 changes: 64 additions & 0 deletions src/nemo-feedback/feedback_io/Data.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,71 @@ std::vector<T> Data<T>::raw_profile(size_t iProfile) const {
template class Data<std::string>;
template class Data<size_t>;
template class Data<int32_t>;
template class Data<feedback_io::QC::Level>;
template class Data<float>;
template class Data<double>;

/// \brief update the whole report QC information for a profile based on the
/// current whole variable QC information
feedback_io::QC::Level wholeReportUpdate(feedback_io::QC::Level current, const
size_t length, const size_t nGoodObs, const size_t nBadObs, const size_t
nBadDoNotAssimilate, const size_t nDoNotAssimilate) {
// If good obs were previously found, the report is good
if (current == QC::Level::Good) {
return QC::Level::Good;
}
// Upgrade any reports if they contain any good obs (goodDoNotAssimilate ->
// good)
if (nGoodObs > 0) {
return QC::Level::Good;
}
// Check currently bad, do-not-assimilate or not-yet-checked reports and
// downgrade where necessary
if (nBadObs == length) {
return QC::Level::Bad;
} else if (nBadDoNotAssimilate == length) {
// No good observations, all bad-do-not-assimilate
return QC::Level::BadDoNotAssimilate;
} else if (nDoNotAssimilate == length) {
// No good observations, some mix of good and bad do-not-assimilate types
// -> GoodDoNotAssimilate
return QC::Level::GoodDoNotAssimilate;
} else {
// No good observations, some mix of bad and do-not-assimilate -> Bad
return QC::Level::Bad;
}
}

void wholeReportFromPerProfile(const Data<feedback_io::QC::Level>& QCData,
Data<feedback_io::QC::Level>& wholeReportQCData) {
for (size_t iProfile = 0;
iProfile < QCData.n_obs(); ++iProfile) {
size_t nGoodObs = 0;
size_t nBadObs = 0;
size_t nBadDoNotAssimilate = 0;
size_t nDoNotAssimilate = 0;
for (size_t iLevel = 0;
iLevel < QCData.length(iProfile);
++iLevel) {
if (QCData(iProfile, iLevel) == QC::Level::Bad) {
++nBadObs;
} else if (QCData(iProfile, iLevel) ==
QC::Level::BadDoNotAssimilate) {
++nBadDoNotAssimilate;
}
if (QCData(iProfile, iLevel) == QC::Level::BadDoNotAssimilate ||
QCData(iProfile, iLevel) == QC::Level::GoodDoNotAssimilate ||
QCData(iProfile, iLevel) == QC::Level::DoNotAssimilate) {
++nDoNotAssimilate;
}
if (QCData(iProfile, iLevel) == QC::Level::Good) {
++nGoodObs;
}
}
wholeReportQCData[iProfile] = wholeReportUpdate(
wholeReportQCData[iProfile], QCData.length(iProfile), nGoodObs,
nBadObs, nBadDoNotAssimilate, nDoNotAssimilate);
}
}
} // namespace feedback_io
} // namespace nemo_feedback
8 changes: 7 additions & 1 deletion src/nemo-feedback/feedback_io/Data.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,12 @@

#include <nemo-feedback/feedback_io/DataIndexer.h>

#include <vector>
#include <string>
#include <vector>
#include <memory>

#include <nemo-feedback/feedback_io/Utils.h>

namespace nemo_feedback {
namespace feedback_io {

Expand Down Expand Up @@ -119,5 +121,9 @@ class Data {
std::shared_ptr<DataIndexer> indexer_;
std::shared_ptr<std::vector<T>> data_;
};

/// \brief create whole report data from per-profile QC data
void wholeReportFromPerProfile(const Data<feedback_io::QC::Level>& QCData,
Data<feedback_io::QC::Level>& wholeReportQCData);
} // namespace feedback_io
} // namespace nemo_feedback
9 changes: 9 additions & 0 deletions src/nemo-feedback/feedback_io/Utils.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,12 @@
#include "oops/util/Logger.h"
#include "oops/util/Duration.h"

#define DOUBLE_FILLVALUE 99999.0
#define FLOAT_FILLVALUE 99999.0
#define INT32_FILLVALUE 0
#define STRING_FILLVALUE "MISSING "


namespace nemo_feedback {
namespace feedback_io {

Expand All @@ -41,6 +47,9 @@ template <> const float typeToFill::value<float>() {
template <> const int32_t typeToFill::value<int32_t>() {
return INT32_FILLVALUE;
}
template <> const QC::Level typeToFill::value<QC::Level>() {
return QC::Level::None;
}
template <> const std::string typeToFill::value<std::string>() {
return STRING_FILLVALUE;
}
Expand Down
Loading

0 comments on commit 7730763

Please sign in to comment.