Improve profile QC output (#54)

 * OBSERVATION_QC only rejects if all data rejected at a given location
 * fix whole profile QC variables for profiles
 * add test of profiles obs to write behaviour
* Add DEPTH_QC, POSITION_QC, JULD_QC with test coverage

src/nemo-feedback/NemoFeedback.cc

@@ -6,7 +6,8 @@
 
 #include "nemo-feedback/NemoFeedback.h"
 
-#include <string.h>
+#include <string>
+#include <string_view>
 
 #include <algorithm>
 #include <utility>
@@ -37,11 +38,43 @@
 #include "ufo/filters/DiagnosticFlag.h"
 #include "ufo/filters/processWhere.h"
 #include "ufo/filters/ObsAccessor.h"
-#include "ufo/utils/metoffice/MetOfficeQCFlags.h"
 
 
 namespace nemo_feedback {
 
+/// helpers
+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,
@@ -189,6 +222,10 @@ 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;
   std::vector<ufo::DiagnosticFlag> do_not_assimilate;
   for (const NemoFeedbackVariableParameters& nemoVariableParams :
         parameters_.variables.value()) {
@@ -210,12 +247,10 @@ void NemoFeedback::write_all_data(feedback_io::Writer<T>& writer,
 
     writer.write_variable(nemo_name + "_OBS", variableData);
 
-    // Write Met Office QC flag data for this variable if they exist.
+    // Pull out QC flags from UFO
     feedback_io::Data<int32_t> variableQCFlagsData;
-    bool hasMetOfficeQC = false;
     if (obsdb_.has("QCFlags", ufo_name)) {
       variableQCFlagsData = creator.create("QCFlags", ufo_name, int32_t(0));
-      hasMetOfficeQC = true;
     } else {
       const size_t iv = flags_->varnames().find(ufo_name);
       std::vector<int32_t> variable_qcFlags;
@@ -233,35 +268,12 @@ void NemoFeedback::write_all_data(feedback_io::Writer<T>& writer,
           nemo_name + "_LEVEL_QC_FLAGS", variableQCFlagsData, 0);
     }
 
-    // Whole Observation report QC flags
-    feedback_io::Data<int32_t> variableQCData(variableData.indexer(),
-        std::vector<int32_t>(variableData.indexer()->n_source_data(), 1));
-    {
-      for (size_t iProf = 0; iProf < variableData.n_obs(); ++iProf) {
-        size_t badObs = 0;
-        for (size_t iLevel = 0; iLevel < variableData.length(iProf); ++iLevel) {
-          if (hasMetOfficeQC &&
-              (variableQCFlagsData(iProf, iLevel)
-               & ufo::MetOfficeQCFlags::WholeObReport::FinalRejectReport)) {
-            ++badObs;
-          }
-          if (!hasMetOfficeQC &&
-              variableQCFlagsData(iProf, iLevel)) {
-            ++badObs;
-          }
-        }
-        variableQCData[iProf] = (badObs == variableData.length(iProf) ? 4 : 1);
-      }
-      writer.write_variable_surf_qc("OBSERVATION_QC", variableQCData);
-    }
-
-    // Overall quality control flags
-    feedback_io::Data<int32_t> variableFinalQCData;
-
+    // Quality control rank variables
     if (obsdb_.has("DiagnosticFlags/FinalReject", ufo_name)) {
-      std::vector<ufo::DiagnosticFlag> final_qc;
-      variableFinalQCData = creator.create("DiagnosticFlags/FinalReject",
-          ufo_name, ufo::DiagnosticFlag(0), 4, 1);
+
+      feedback_io::Data<int32_t> variableFinalQCData =
+        creator.create("DiagnosticFlags/FinalReject", ufo_name,
+                       ufo::DiagnosticFlag(0), 4, 1);
 
       // Add do not assimilate flag if required.
       if (obsdb_.has("DiagnosticFlags/DoNotAssimilate", ufo_name)) {
@@ -279,10 +291,50 @@ void NemoFeedback::write_all_data(feedback_io::Writer<T>& writer,
         }
       }
 
+      // 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);
+      }
+
       writer.write_variable_surf_qc(nemo_name + "_QC",
-          variableFinalQCData);
+        wholeVariableQCData);
       writer.write_variable_level_qc(nemo_name + "_LEVEL_QC",
           variableFinalQCData);
+
+      // Whole Observation report QC
+      const std::string positionQCGroup("DiagnosticFlags/PositionReject");
+      if (obsdb_.has(positionQCGroup, ufo_name)) {
+        wholeReportQCFromDiagnosticFlags(creator, positionQCGroup, ufo_name,
+            wholeReportPositionQCData);
+      }
+
+      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);
+          }
+        }
+      }
     }
 
     // Write additional variables for this variable
@@ -295,6 +347,25 @@ void NemoFeedback::write_all_data(feedback_io::Writer<T>& writer,
             ufo_name, T(0)));
       writer.write_variable(add_name, variableAdditionalData);
     }
+  }  // loop: variables
+
+  // Write whole report variables
+  writer.write_variable_surf_qc("OBSERVATION_QC", wholeReportQCData);
+
+  const std::string depthQCGroup("DiagnosticFlags/DepthReject");
+  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));
+    writer.write_variable_level_qc("DEPTH_QC", depthQCData);
+  }
+  if (wholeReportPositionQCData.n_obs() != 0) {
+    writer.write_variable_surf_qc("POSITION_QC", wholeReportPositionQCData);
+  }
+  if (wholeReportTimeQCData.n_obs() != 0) {
+    writer.write_variable_surf_qc("JULD_QC", wholeReportTimeQCData);
   }
 }
 
@@ -306,9 +377,9 @@ feedback_io::MetaData NemoFeedback::setupMetaData(
   feedback_io::Data<double> lons(creator.create("MetaData", "longitude",
         static_cast<double>(0)));
   const std::string depthGroup = parameters_.depthGroup.value()
-    .value_or("MetaData");
+    .value_or(static_cast<std::string>(defaultDepthGroup));
   const std::string depthVariable = parameters_.depthVariable.value()
-    .value_or("depthBelowWaterSurface");
+    .value_or(static_cast<std::string>(defaultDepthVariable));
   feedback_io::Data<double> depths;
   if (obsdb_.has(depthGroup, depthVariable)) {
     depths = feedback_io::Data<double>(creator.create(depthGroup, depthVariable,

src/nemo-feedback/NemoFeedbackDataCreator.cc

@@ -62,11 +62,11 @@ feedback_io::Data<T> NemoFeedbackDataCreator::create_from_obsdb(
   oops::Log::trace() << NemoFeedbackDataCreator::className()
                      << ":create_from_obsdb " << obsGroup
                      << "/" << ufoName << std::endl;
-  std::vector<T> data;
+  std::vector<T> data(obsdb_.nlocs());
   obsdb_.get_db(obsGroup, ufoName, data);
   if (data.size() != obsdb_.nlocs())
     throw eckit::BadValue(NemoFeedbackDataCreator::className()
-        + ":create_from_obsdb no data ");
+        + ":create_from_obsdb no data ", Here());
   const T missingValue = util::missingValue<T>();
   // Convert oops missing values to NEMO missing value
   for_each(data.begin(), data.end(),
@@ -152,14 +152,14 @@ feedback_io::Data<int32_t> NemoFeedbackDataCreator::create_from_obsdb(const
                      << ":create_from_obsdb ufo::DiagnosticFlag "
                      << obsGroup << "/" << ufoName << std::endl;
   std::vector<int32_t> data;
-  std::vector<ufo::DiagnosticFlag> flagData;
+  std::vector<ufo::DiagnosticFlag> flagData(obsdb_.nlocs(), 0);
   obsdb_.get_db(obsGroup, ufoName, flagData);
   for (ufo::DiagnosticFlag flag : flagData) {
-    data.push_back(flag ? whenTrue : whenFalse);
+    data.emplace_back(flag ? whenTrue : whenFalse);
   }
   if (data.size() != obsdb_.nlocs())
     throw eckit::BadValue(NemoFeedbackDataCreator::className()
-        + ":create_from_obsdb ufo::DiagnosticFlag no data.");
+        + ":create_from_obsdb ufo::DiagnosticFlag no data.", Here());
 
   return feedback_io::Data<int32_t>(indexer_, std::move(data));
 }
@@ -234,7 +234,7 @@ std::tuple<std::string, feedback_io::Data<double>>
         feedback_io::Data<double>());
   }
 
-  std::vector<util::DateTime> datetimes;
+  std::vector<util::DateTime> datetimes(obsdb_.nlocs());
   obsdb_.get_db(obsGroup, ufoName, datetimes);
 
   std::vector<double> data;

src/tests/testinput/hofx_profiles_2var_writer.yaml

@@ -29,11 +29,6 @@ observations:
       observation vertical coordinate group: ObsValue
       interpolation method: linear
     obs filters:
-    - filter: Variable Assignment
-      assignments:
-      - name: DerivedObsValue/waterTemperature
-        type: float
-        value: 12.0
     - filter: NEMO Feedback Writer
       filename: testoutput/test_hofx3d_nc_prof_2vars_writer_out.nc
       reference date: 1950-01-01T00:00:00Z

src/tests/testinput/hofx_profiles_writer.yaml

@@ -39,6 +39,11 @@ observations:
         - name: waterPotentialTemperature
           feedback suffix: Hx
           ioda group: HofX
+      where:
+      - is_true: null
+        variable:
+          name: DiagnosticFlags/ObsToWrite/waterPotentialTemperature
+
     benchmarkFlag: 1000 # just to keep the ObsFilters test happy
     flaggedBenchmark: 0
     HofX: hofx

src/tests/testoutput/test_hofx3d_nc_prof_2vars_writer_out_ref.cdl

@@ -164,8 +164,8 @@ data:
   "    " ;
 
  DEPTH_QC =
-  _, _, _, _, _, _,
-  _, _, _, _, _, _ ;
+  1, 1, 1, 1, 4, 4,
+  1, 1, 1, 1, 4, 4 ;
 
  DEPTH_QC_FLAGS =
   _, _,
@@ -181,19 +181,19 @@ data:
   _, _,
   _, _ ;
 
- OBSERVATION_QC = 4, 4 ;
+ OBSERVATION_QC = 1, 1 ;
 
  OBSERVATION_QC_FLAGS =
   _, _,
   _, _ ;
 
- POSITION_QC = _, _ ;
+ POSITION_QC = 1, 1 ;
 
  POSITION_QC_FLAGS =
   _, _,
   _, _ ;
 
- JULD_QC = _, _ ;
+ JULD_QC = 4, 1 ;
 
  JULD_QC_FLAGS =
   _, _,
@@ -227,11 +227,11 @@ data:
   1, _,
   1, _ ;
 
- POTM_QC = _, _ ;
+ POTM_QC = 1, 1 ;
 
  POTM_LEVEL_QC =
-  _, _, _, _, _, _,
-  _, _, _, _, _, _ ;
+  1, 1, 1, 4, 1, 1,
+  1, 1, 1, 4, 1, 1 ;
 
  POTM_IOBSI = _, _ ;
 
@@ -269,11 +269,11 @@ data:
   1, _,
   1, _ ;
 
- PSAL_QC = _, _ ;
+ PSAL_QC = 1, 1 ;
 
  PSAL_LEVEL_QC =
-  _, _, _, _, _, _,
-  _, _, _, _, _, _ ;
+  1, 1, 1, 4, 1, 1,
+  1, 1, 1, 4, 1, 1 ;
 
  PSAL_IOBSI = _, _ ;