CalcHofX.h

Go to the documentation of this file.

/*
 * (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 OOPS_ASSIMILATION_CALCHOFX_H_
#define OOPS_ASSIMILATION_CALCHOFX_H_
 
#include <memory>
#include <string>
#include <vector>
 
#include "eckit/config/LocalConfiguration.h"
 
#include "oops/base/ObsAuxControls.h"
#include "oops/base/Observations.h"
#include "oops/base/ObsFilters.h"
#include "oops/base/ObsSpaces.h"
#include "oops/base/ObsVector.h"
#include "oops/base/Variables.h"
#include "oops/interface/GeoVaLs.h"
#include "oops/interface/Locations.h"
#include "oops/interface/ObsDataVector.h"
#include "oops/interface/ObsDiagnostics.h"
#include "oops/interface/ObsOperator.h"
#include "oops/util/Logger.h"
#include "oops/util/parameters/Parameter.h"
#include "oops/util/parameters/Parameters.h"
#include "oops/util/parameters/RequiredParameter.h"
 
namespace oops {
 
template <typename OBS>
class CalcHofXParameters : public Parameters {
  OOPS_CONCRETE_PARAMETERS(CalcHofXParameters, Parameters)
 
  typedef typename OBS::ObsOperator::Parameters_ ObsOperatorParameters_;
 
 public:
  oops::RequiredParameter<ObsOperatorParameters_> obsOperator{"obs operator", this};
  oops::Parameter<std::vector<ObsFilterParametersWrapper<OBS>>> obsFilters{"obs filters", {}, this};
};
 
// -----------------------------------------------------------------------------
 
/// \brief Computes observation operator (from GeoVaLs), applies bias correction
///        and runs QC filters
template <typename OBS>
class CalcHofX {
  typedef GeoVaLs<OBS>               GeoVaLs_;
  typedef Locations<OBS>             Locations_;
  typedef ObsAuxControls<OBS>        ObsAuxCtrls_;
  typedef ObsDiagnostics<OBS>        ObsDiags_;
  typedef Observations<OBS>          Observations_;
  typedef ObsFilters<OBS>            ObsFilters_;
  typedef ObsOperator<OBS>           ObsOperator_;
  typedef ObsSpaces<OBS>             ObsSpaces_;
  typedef ObsVector<OBS>             ObsVector_;
  template <typename DATA> using ObsData_ = ObsDataVector<OBS, DATA>;
  template <typename DATA> using ObsDataVec_ = std::vector<std::shared_ptr<ObsData_<DATA>>>;
 
  typedef std::vector<std::unique_ptr<GeoVaLs_>>       GeoVaLsVec_;
  typedef std::vector<std::unique_ptr<Locations_>>     LocationsVec_;
  typedef std::vector<std::unique_ptr<ObsFilters_>>    ObsFiltersVec_;
  typedef std::vector<std::unique_ptr<ObsOperator_>>   ObsOperatorVec_;
  typedef std::vector<std::shared_ptr<ObsVector_>>     ObsVectorVec_;
  typedef std::vector<Variables>                       VariablesVec_;
 
 public:
/// \brief Initializes ObsOperators, Locations, and QC data
  CalcHofX(const ObsSpaces_ &, const eckit::Configuration &);
 
/// \brief Initializes variables, obs bias, obs filters (could be different for
/// different iterations
  void initialize(const ObsAuxCtrls_ &, const int iteration = 0);
 
/// \brief Computes H(x) from the filled in GeoVaLs
  Observations_ compute(const GeoVaLsVec_ &);
 
/// \brief accessor to the locations
  const LocationsVec_ & locations() const { return locations_; }
/// \brief accessor to variables required from the model
  const VariablesVec_ & requiredVars() const { return geovars_; }
/// \brief accessor to QC flags
  const ObsDataVec_<int> & qcflags() const { return qcflags_; }
 
/// \brief read QC flags from \p qcname variable from file
  void readQcFlags(const std::string & qcname);
/// \brief reset QC flags and ObsErrors
  void resetQc();
/// \brief save QC flags to ObsSpaces
  void saveQcFlags(const std::string &) const;
/// \brief save obs error variances (modified in QC) to ObsSpaces
  void saveObsErrors(const std::string &) const;
/// \brief mask obs errors with QC flags
  void maskObsErrors();
 
 private:
  eckit::LocalConfiguration obsconfig_;
  const ObsSpaces_ &   obspaces_;   // ObsSpaces used in H(x)
  ObsOperatorVec_      obsops_;     // Obs operators
  LocationsVec_        locations_;  // locations
  const ObsAuxCtrls_ * biascoeff_;  // bias coefficients
  ObsDataVec_<int>     qcflags_;    // QC flags
  ObsVectorVec_        obserrs_;    // Obs error variances (used in QC filters)
  ObsFiltersVec_       filters_;    // QC filters
  VariablesVec_        geovars_;    // variables required from the model
};
 
// -----------------------------------------------------------------------------
 
template <typename OBS>
CalcHofX<OBS>::CalcHofX(const ObsSpaces_ & obspaces,
                        const eckit::Configuration & config) :
  obsconfig_(config), obspaces_(obspaces), obsops_(), locations_(),
  biascoeff_(nullptr), qcflags_(), obserrs_(), filters_(),
  geovars_(obspaces_.size())
{
  std::vector<eckit::LocalConfiguration> obsconfs = config.getSubConfigurations();
  for (size_t jj = 0; jj < obspaces_.size(); ++jj) {
    // obsconfs[jj] contains not only options controlling the obs operator and filters (known to
    // CalcHofX) but also those controlling the obs space (unknown to it). So we can't call
    // validateAndDeserialize() here, since "obs space" would be treated as an unrecognized
    // keyword. In the long term the code constructing the CalcHofX will probably need to split
    // the contents of the "observations" vector into two vectors, one containing the "obs space"
    // sections and the other the "obs operator" and "obs filters" sections, and pass the former to
    // the constructor of ObsSpaces and the latter to the constructor of CalcHofX.
    CalcHofXParameters<OBS> observerParams;
    observerParams.deserialize(obsconfs[jj]);
    /// Set up observation operators
    obsops_.emplace_back(new ObsOperator_(obspaces_[jj], observerParams.obsOperator));
    locations_.emplace_back(new Locations_(obsops_[jj]->locations()));
 
    /// Allocate QC flags
    qcflags_.emplace_back(std::make_shared<ObsData_<int>>(obspaces[jj],
                            obspaces[jj].obsvariables()));
    /// Allocate and read initial obs error
    obserrs_.emplace_back(std::make_shared<ObsVector_>(obspaces[jj], "ObsError"));
  }
  Log::trace() << "CalcHofX<OBS> constructed" << std::endl;
}
 
// -----------------------------------------------------------------------------
 
template <typename OBS>
void CalcHofX<OBS>::initialize(const ObsAuxCtrls_ & obsaux, const int iteration) {
  std::vector<eckit::LocalConfiguration> obsconfs = obsconfig_.getSubConfigurations();
  biascoeff_ = &obsaux;
  filters_.clear();
  for (size_t jj = 0; jj < obspaces_.size(); ++jj) {
    CalcHofXParameters<OBS> observerParams;
    observerParams.deserialize(obsconfs[jj]);
    /// Set up QC filters and run preprocess
    filters_.emplace_back(new ObsFilters_(obspaces_[jj], observerParams.obsFilters,
                                          qcflags_[jj], *obserrs_[jj], iteration));
    filters_[jj]->preProcess();
 
    /// Set up variables requested from the model
    geovars_[jj] += obsops_[jj]->requiredVars();
    geovars_[jj] += (*biascoeff_)[jj].requiredVars();
    geovars_[jj] += filters_[jj]->requiredVars();
  }
  Log::trace() << "CalcHofX<OBS>::initialize done" << std::endl;
}
 
// -----------------------------------------------------------------------------
 
template <typename OBS>
Observations<OBS> CalcHofX<OBS>::compute(const GeoVaLsVec_ & geovals) {
  oops::Log::trace() << "CalcHofX<OBS>::compute start" << std::endl;
 
  Observations<OBS> yobs(obspaces_);
  Observations<OBS> ybias(obspaces_);
  ybias.zero();
  for (size_t jj = 0; jj < obspaces_.size(); ++jj) {
    /// call prior filters
    filters_[jj]->priorFilter(*geovals[jj]);
    /// compute H(x)
    oops::Variables vars;
    vars += filters_[jj]->requiredHdiagnostics();
    vars += (*biascoeff_)[jj].requiredHdiagnostics();
    ObsDiags_ ydiags(obspaces_[jj], *locations_[jj], vars);
    obsops_[jj]->simulateObs(*geovals[jj], yobs[jj], (*biascoeff_)[jj], ybias[jj], ydiags);
    /// call posterior filters
    filters_[jj]->postFilter(yobs[jj], ybias[jj], ydiags);
  }
  oops::Log::trace() << "CalcHofX<OBS>::compute done" << std::endl;
  return yobs;
}
 
// -----------------------------------------------------------------------------
 
template <typename OBS>
void CalcHofX<OBS>::readQcFlags(const std::string & qcname) {
  for (const auto & qcflag : qcflags_) {
    qcflag->read(qcname);
  }
}
 
// -----------------------------------------------------------------------------
 
template <typename OBS>
void CalcHofX<OBS>::resetQc() {
  for (const auto & qcflag : qcflags_) {
    qcflag->zero();
  }
  for (const auto & obserr : obserrs_) {
    obserr->read("ObsError");
  }
}
 
// -----------------------------------------------------------------------------
 
template <typename OBS>
void CalcHofX<OBS>::saveQcFlags(const std::string & name) const {
  for (const auto & qcflag : qcflags_) {
    qcflag->save(name);
  }
}
 
// -----------------------------------------------------------------------------
 
template <typename OBS>
void CalcHofX<OBS>::saveObsErrors(const std::string & name) const {
  for (const auto & obserr : obserrs_) {
    obserr->save(name);
  }
}
 
// -----------------------------------------------------------------------------
 
template <typename OBS>
void CalcHofX<OBS>::maskObsErrors() {
  for (size_t jj = 0; jj < obserrs_.size(); ++jj) {
    obserrs_[jj]->mask(*qcflags_[jj]);
  }
}
 
// -----------------------------------------------------------------------------
 
}  // namespace oops
 
#endif  // OOPS_ASSIMILATION_CALCHOFX_H_