Add snow depth error inflation based on elevation difference from model

src/ufo/filters/CMakeLists.txt

@@ -323,6 +323,8 @@ set ( filters_files
       obsfunctions/ObsErrorFactorSfcPressure.h
       obsfunctions/ObsErrorFactorConventional.cc
       obsfunctions/ObsErrorFactorConventional.h
+      obsfunctions/ObsErrorFactorSnowElevationDiff.cc
+      obsfunctions/ObsErrorFactorSnowElevationDiff.h
       obsfunctions/PotentialTemperatureFromTemperature.cc
       obsfunctions/PotentialTemperatureFromTemperature.h
       obsfunctions/NearSSTRetCheckIR.cc

src/ufo/filters/obsfunctions/ObsErrorFactorSnowElevationDiff.cc

@@ -0,0 +1,113 @@
+/*
+ * (C) Copyright 2020 UCAR
+ *
+ * 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.
+ */
+
+#include "ufo/filters/obsfunctions/ObsErrorFactorSnowElevationDiff.h"
+
+#include <cmath>
+
+#include "eckit/exception/Exceptions.h"
+#include "ioda/ObsDataVector.h"
+#include "oops/util/Logger.h"
+#include "oops/util/missingValues.h"
+#include "ufo/filters/ObsFilterData.h"
+
+namespace ufo {
+
+static ObsFunctionMaker<ObsErrorFactorSnowElevationDiff> 
+    makerObsErrorFactorSnowElevationDiff_("ObsErrorFactorSnowElevationDiff");
+
+// -----------------------------------------------------------------------------
+
+ObsErrorFactorSnowElevationDiff::ObsErrorFactorSnowElevationDiff(
+    const eckit::Configuration &config)
+  : invars_() {
+  oops::Log::trace() << "ObsErrorFactorSnowElevationDiff constructor" << std::endl;
+  oops::Log::debug() << "ObsErrorFactorSnowElevationDiff: config = " << config << std::endl;
+
+  // Initialize options
+  options_.reset(new ObsErrorFactorSnowElevationDiffParameters());
+  options_->deserialize(config);
+
+  // Include observation elevation variable
+  const std::string obs_elev_var = options_->obs_elevation_var.value();
+  invars_ += Variable(obs_elev_var);
+
+  // Include model elevation variable
+  const std::string model_elev_var = options_->model_elevation_var.value();
+  invars_ += Variable(model_elev_var);
+}
+
+// -----------------------------------------------------------------------------
+
+ObsErrorFactorSnowElevationDiff::~ObsErrorFactorSnowElevationDiff() {
+  oops::Log::trace() << "ObsErrorFactorSnowElevationDiff destructor" << std::endl;
+}
+
+// -----------------------------------------------------------------------------
+
+void ObsErrorFactorSnowElevationDiff::compute(
+    const ObsFilterData & data,
+    ioda::ObsDataVector<float> & obserr) const {
+  oops::Log::trace() << "ObsErrorFactorSnowElevationDiff compute start" << std::endl;
+
+  const float missing = util::missingValue<float>();
+  const float h_scale = options_->elevation_scale_m.value();
+
+  // If no observations on this processor then nothing to do
+  if (data.nlocs() == 0) return;
+
+  // Ensure that only one output variable is expected
+  ASSERT(obserr.nvars() == 1);
+
+  // Get dimensions
+  size_t nlocs = data.nlocs();
+
+  // Get observation elevation variable
+  std::vector<float> ob_elevation(nlocs);
+  const std::string obs_elev_var = options_->obs_elevation_var.value();
+  data.get(Variable(obs_elev_var), ob_elevation);
+
+  // Get model surface elevation variable
+  std::vector<float> model_elevation(nlocs);
+  const std::string model_elev_var = options_->model_elevation_var.value();
+  data.get(Variable(model_elev_var), model_elevation);
+
+  // Compute inflation factor for each observation
+  float dz, inflation_factor_1, inflation_factor;
+  int iv = 0;
+
+  for (size_t iloc = 0; iloc < nlocs; ++iloc) {
+    // If missing observation or model elevation, set obserror to missing
+    if (ob_elevation[iloc] == missing || model_elevation[iloc] == missing) {
+      obserr[iv][iloc] = missing;
+    } else {
+      // Compute elevation difference
+      dz = std::abs(model_elevation[iloc] - ob_elevation[iloc]);
+
+      // Compute inflation_factor_1 = exp(-1 * dz^2 / (h^2))
+      inflation_factor_1 = std::exp(-1.0f * dz * dz / (h_scale * h_scale));
+
+      // Compute output inflation_factor = 1 / inflation_factor_1
+      inflation_factor = 1.0f / inflation_factor_1;
+
+      // Inflate error by the inflation factor
+      obserr[iv][iloc] *= inflation_factor;
+    }
+  }
+
+  oops::Log::trace() << "ObsErrorFactorSnowElevationDiff compute complete" << std::endl;
+}
+
+// -----------------------------------------------------------------------------
+
+const ufo::Variables & ObsErrorFactorSnowElevationDiff::requiredVariables() const {
+  return invars_;
+}
+
+// -----------------------------------------------------------------------------
+
+}  // namespace ufo

src/ufo/filters/obsfunctions/ObsErrorFactorSnowElevationDiff.h

@@ -0,0 +1,78 @@
+/*
+ * (C) Copyright 2020 UCAR
+ *
+ * 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.
+ */
+
+#ifndef UFO_FILTERS_OBSFUNCTIONS_OBSERRORFACTORSNOWELEVATIONDIFF_H_
+#define UFO_FILTERS_OBSFUNCTIONS_OBSERRORFACTORSNOWELEVATIONDIFF_H_
+
+#include <memory>
+#include <string>
+#include <vector>
+
+#include "oops/util/parameters/Parameter.h"
+#include "oops/util/parameters/Parameters.h"
+#include "oops/util/parameters/RequiredParameter.h"
+
+#include "ufo/filters/obsfunctions/ObsFunctionBase.h"
+#include "ufo/filters/Variables.h"
+
+namespace ufo {
+
+/// \brief Options controlling ObsErrorFactorSnowElevationDiff ObsFunction
+class ObsErrorFactorSnowElevationDiffParameters : public oops::Parameters {
+  OOPS_CONCRETE_PARAMETERS(ObsErrorFactorSnowElevationDiffParameters, Parameters)
+
+ public:
+  oops::Parameter<std::string> obs_elevation_var{"observation_elevation", "MetaData/stationElevation", this};
+  oops::Parameter<std::string> model_elevation_var{"model_elevation", "GeoVaLs/filtered_orography", this};
+  oops::Parameter<float> elevation_scale_m{"elevation_scale_m", 800., this};
+};
+
+// -----------------------------------------------------------------------------
+
+/// \brief Inflate observation error based on elevation difference between model and observation.
+///
+/// This routine computes an observation error inflation factor based on the elevation difference
+/// between the model surface elevation and the observed station elevation.
+/// The inflation factor is computed as: 1 / exp(-1 * dz^2 / (h^2))
+/// where dz = |model_elevation - obs_elevation| 
+/// and h is the elevation_scale parameter in m (Elevation difference at which the obs-error will be inflated by a factor of exp(1)).
+///
+/// Authors: Tseganeh Z Gichamo and Gihub Copilot. All the ML generated code sections have been revised and tested by Tseganeh Z. Gichamo
+/// 
+/// ### example configurations for application of this filter: ###
+///
+///     - filter: Perform Action
+///       filter variables:
+///       - name: totalSnowDepth
+///       action:
+///         name: inflate error
+///         inflation variable:
+///           name: ObsFunction/ObsErrorFactorSnowElevationDiff
+///           options:
+///             observation_elevation: MetaData/stationElevation
+///             model_elevation: GeoVaLs/filtered_orography
+///             elevation_scale_m: 800.0   // (Unit m) 
+///
+class ObsErrorFactorSnowElevationDiff : public ObsFunctionBase<float> {
+ public:
+  static const std::string classname() {return "ObsErrorFactorSnowElevationDiff";}
+
+  explicit ObsErrorFactorSnowElevationDiff(const eckit::Configuration &config);
+  ~ObsErrorFactorSnowElevationDiff();
+
+  void compute(const ObsFilterData &, ioda::ObsDataVector<float> &) const;
+  const ufo::Variables & requiredVariables() const;
+ private:
+  ufo::Variables invars_;
+  std::unique_ptr<ObsErrorFactorSnowElevationDiffParameters> options_;
+};
+
+// -----------------------------------------------------------------------------
+
+}  // namespace ufo
+
+#endif  // UFO_FILTERS_OBSFUNCTIONS_OBSERRORFACTORSNOWELEVATIONDIFF_H_