refactor: swap getValues for Interpolator

src/orca-jedi/CMakeLists.txt

@@ -13,10 +13,9 @@ errorcovariance/ErrorCovariance.h
 geometry/GeometryParameters.h
 geometry/Geometry.h
 geometry/Geometry.cc
-getvalues/GetValuesParameters.h
-getvalues/GetValues.h
-getvalues/GetValues.cc
-getvalues/LinearGetValues.h
+interpolator/InterpolatorParameters.h
+interpolator/Interpolator.h
+interpolator/Interpolator.cc
 increment/Increment.cc
 increment/Increment.h
 model/ModelBiasCovariance.h
@@ -31,9 +30,7 @@ nemo_io/NemoFieldReader.h
 nemo_io/NemoFieldReader.cc
 utilities/OrcaModelTraits.h
 variablechanges/VariableChangeParameters.h
-variablechanges/LinearVariableChangeParameters.h
 variablechanges/VariableChange.h
-variablechanges/LinearVariableChange.h
 )
 
 

src/orca-jedi/geometry/Geometry.cc

@@ -106,6 +106,18 @@ std::vector<size_t> Geometry::variableSizes(const oops::Variables & vars) const
   return varSizes;
 }
 
+void Geometry::latlon(std::vector<double> & lats, std::vector<double> & lons,
+    const bool halo) const {
+  const auto lonlat = atlas::array::make_view<double, 2>(funcSpace_.lonlat());
+  const size_t npts = funcSpace_.size();
+  lons.resize(npts);
+  lats.resize(npts);
+  for (size_t jj = 0; jj < npts; ++jj) {
+    lons[jj] = lonlat(jj, 0);
+    lats[jj] = lonlat(jj, 1);
+  }
+}
+
 const oops::Variables & Geometry::variables() const {
   return vars_;
 }

src/orca-jedi/geometry/Geometry.h

@@ -48,14 +48,14 @@ class Geometry : public util::Printable {
  public:
   typedef OrcaGeometryParameters Parameters__;
 
-  // Geometry(const Parameters_ &, const eckit::mpi::Comm &);
   Geometry(const eckit::Configuration &, const eckit::mpi::Comm &);
   ~Geometry();
 
   std::vector<size_t> variableSizes(const oops::Variables &) const;
   const eckit::mpi::Comm & getComm() const {return comm_;}
-
   const oops::Variables & variables() const;
+  void latlon(std::vector<double> & lats, std::vector<double> & lons,
+              const bool halo) const;
 
   const atlas::Mesh & mesh() const {return mesh_;}
   const atlas::functionspace::NodeColumns & funcSpace() const

src/orca-jedi/getvalues/GetValues.cc → src/orca-jedi/interpolator/Interpolator.cc

@@ -27,33 +27,29 @@
 #include "oops/util/Printable.h"
 #include "oops/util/missingValues.h"
 
-#include "ufo/GeoVaLs.h"
-#include "ufo/Locations.h"
-
 #include "orca-jedi/geometry/Geometry.h"
 #include "orca-jedi/state/State.h"
 
-#include "orca-jedi/getvalues/GetValues.h"
+#include "orca-jedi/interpolator/Interpolator.h"
 
 namespace eckit {
   class Configuration;
 }
 
 namespace ufo {
   class GeoVaLs;
-  class Locations;
 }
 
 namespace orcamodel {
   class State;
   class Geometry;
 
   atlas::functionspace::PointCloud atlasObsFuncSpaceFactory(
-      const ufo::Locations & locs) {
-      size_t nlocs = locs.size();
+      const std::vector<double>& locs) {
+      size_t nlocs = locs.size() / 2;
       if (nlocs == 0) {
         std::stringstream err_stream;
-        err_stream << "orcamodel::GetValues::"
+        err_stream << "orcamodel::Interpolator::"
                    << " Constructor called with no locations " << std::endl;
         err_stream << "         "
                    << "this might mean that there are no observations "
@@ -63,52 +59,49 @@ namespace orcamodel {
       }
 
       // Setup observation functionspace
-      oops::Log::trace() << "orcamodel::GetValues:: creating atlasObsFuncSpace "
+      oops::Log::trace() << "orcamodel::Interpolator:: creating atlasObsFuncSpace "
                          << "with nlocs = " << nlocs << std::endl;
-      auto lons = locs.lons();
-      auto lats = locs.lats();
       atlas::Field points("lonlat", atlas::array::make_datatype<double>(),
           atlas::array::make_shape(nlocs, 2));
       auto arrv_t = atlas::array::make_view<double, 2>(points);
-      for (atlas::idx_t i = 0; i < arrv_t.shape(0); ++i) {
-        arrv_t(i, 0) = lons[i];
-        arrv_t(i, 1) = lats[i];
+      for ( unsigned int j = 0, latIndex=0, lonIndex=1;
+        j < nlocs;
+        ++j, latIndex+=2, lonIndex+=2 ) {
+        arrv_t(j, 1) = locs[latIndex];
+        arrv_t(j, 0) = locs[lonIndex];
       }
-      oops::Log::trace() << "orcamodel::GetValues:: creating atlasObsFuncSpace "
+      oops::Log::trace() << "orcamodel::Interpolator:: creating atlasObsFuncSpace "
                          << "... done" << std::endl;
 
       return atlas::functionspace::PointCloud(std::move(points));
   }
 
-  GetValues::GetValues(const Geometry & geom, const ufo::Locations & locs,
-            const eckit::Configuration & conf) :
+  Interpolator::Interpolator(const eckit::Configuration & conf, const Geometry & geom,
+      const std::vector<double>& locs) :
+      nlocs_(locs.size() / 2),
       atlasObsFuncSpace_(std::move(atlasObsFuncSpaceFactory(locs))),
       interpolator_(eckit::LocalConfiguration(conf, "atlas-interpolator"),
                     geom.funcSpace(),
                     atlasObsFuncSpace_ ) {
     params_.validateAndDeserialize(conf);
-    oops::Log::trace() << "orcamodel::GetValues:: conf:" << conf
+    oops::Log::trace() << "orcamodel::Interpolator:: conf:" << conf
                        << std::endl;
-    oops::Log::debug() << "orcamodel::GetValues:: atlasObsFuncSpace_:"
+    oops::Log::debug() << "orcamodel::Interpolator:: atlasObsFuncSpace_:"
                        << atlasObsFuncSpace_ << std::endl;
   }
 
-  void GetValues::fillGeoVaLs(const State& state,
-      const util::DateTime& dt_begin, const util::DateTime& dt_end,
-      ufo::GeoVaLs& geovals) const {
-    std::size_t nlocs = geovals.nlocs();
+  void Interpolator::apply(const oops::Variables& vars, const State& state,
+      std::vector<double>& result) const {
 
-    oops::Log::trace() << "orcamodel::GetValues::fillGeoVaLs starting "
+    oops::Log::trace() << "orcamodel::Interpolator::apply starting "
                        << std::endl;
 
-    std::vector<double> vals(nlocs);
-    const oops::Variables geovalsVars = geovals.getVars();
-    const size_t nvars = geovalsVars.size();
+    const size_t nvars = vars.size();
     for (size_t j=0; j < nvars; ++j) {
-      if (!state.variables().has(geovalsVars[j])) {
+      if (!state.variables().has(vars[j])) {
         std::stringstream err_stream;
-        err_stream << "orcamodel::GetValues::fillGeoVals geovals varname \" "
-                   << "\" " << geovalsVars[j]
+        err_stream << "orcamodel::Interpolator::apply varname \" "
+                   << "\" " << vars[j]
                    << " not found in the model state." << std::endl;
         err_stream << "    add the variable to the state variables and "
                    << "add a mapping from the geometry to that variable."
@@ -117,12 +110,12 @@ namespace orcamodel {
       }
     }
     const std::vector<size_t> varSizes =
-      state.geometry()->variableSizes(geovalsVars);
+      state.geometry()->variableSizes(vars);
     for (size_t j=0; j < nvars; ++j) {
-      auto gv_varname = geovalsVars[j];
+      auto gv_varname = vars[j];
       if (varSizes[j] != 1) {
         std::stringstream err_stream;
-        err_stream << "orcamodel::GetValues::fillGeoVals interpolating "
+        err_stream << "orcamodel::Interpolator::apply interpolating "
                    << "data with levels > 1 not implemented." << std::endl;
         throw eckit::NotImplemented(err_stream.str(), Here());
       }
@@ -132,29 +125,22 @@ namespace orcamodel {
       auto field_view = atlas::array::make_view<double, 1>(tgt_field);
       atlas::field::MissingValue mv(state.stateFields()[gv_varname]);
       bool has_mv = static_cast<bool>(mv);
-      for (std::size_t i=0; i < nlocs; i++) {
+      for (std::size_t i=0; i < nlocs_; i++) {
         if (has_mv && mv(field_view(i))) {
-          vals[i] = util::missingValue(field_view(i));
+          result[i] = util::missingValue(field_view(i));
         } else {
-          vals[i] = field_view(i);
+          result[i] = field_view(i);
         }
       }
-      geovals.putAtLevel(vals, gv_varname, 0);
     }
-    oops::Log::debug() << "orcamodel::GetValues::fillGeoVaLs geovals print: "
-                       << geovals << std::endl;
-    oops::Log::trace() << "orcamodel::GetValues::fillGeoVaLs done "
+    oops::Log::trace() << "orcamodel::Interpolator::apply done "
                        << std::endl;
   }
 
-  void GetValues::print(std::ostream & os) const {
-    os << "orcamodel::GetValues: " << std::endl;
+  void Interpolator::print(std::ostream & os) const {
+    os << "orcamodel::Interpolator: " << std::endl;
     os << "  Obs function space " << atlasObsFuncSpace_ << std::endl;
     os << "  Interpolator " << interpolator_ << std::endl;
   }
 
 }  // namespace orcamodel
-
-
-
-

src/orca-jedi/getvalues/GetValues.h → src/orca-jedi/interpolator/Interpolator.h

@@ -25,10 +25,7 @@
 #include "oops/util/ObjectCounter.h"
 #include "oops/util/Printable.h"
 
-#include "ufo/GeoVaLs.h"
-#include "ufo/Locations.h"
-
-#include "orca-jedi/getvalues/GetValuesParameters.h"
+#include "orca-jedi/interpolator/InterpolatorParameters.h"
 #include "orca-jedi/geometry/Geometry.h"
 #include "orca-jedi/state/State.h"
 
@@ -37,42 +34,40 @@ namespace eckit {
   class Configuration;
 }
 
-namespace ufo {
-  class GeoVaLs;
-  class Locations;
-}
-
 namespace orcamodel {
   class State;
   class Geometry;
+  class Increment;
 
 atlas::functionspace::PointCloud atlasObsFuncSpaceFactory(
-    const ufo::Locations & locs);
+    const std::vector<double> & locs);
 
-class GetValues : public util::Printable,
-  private util::ObjectCounter<GetValues> {
+class Interpolator : public util::Printable,
+  private util::ObjectCounter<Interpolator> {
  public:
-  static const std::string classname() {return "orcamodel::GetValues";}
+  static const std::string classname() {return "orcamodel::Interpolator";}
 
-  typedef OrcaGetValuesParameters Parameters_;
+  typedef OrcaInterpolatorParameters Parameters_;
 
-  GetValues(const Geometry & geom, const ufo::Locations & locs,
-    const eckit::Configuration & conf);
+  Interpolator(const eckit::Configuration & conf, const Geometry & geom,
+      const std::vector<double>& locs);
 
-  virtual ~GetValues() {}
+  virtual ~Interpolator() {}
 
-  void fillGeoVaLs(const State& state, const util::DateTime& dt_begin,
-    const util::DateTime& dt_end, ufo::GeoVaLs& geovals) const;
+  void apply(const oops::Variables& vars, const State& state, std::vector<double>& result) const;
+  void apply(const oops::Variables& vars, const Increment& inc, std::vector<double>& result) const {
+    throw eckit::NotImplemented("Increment interpolation not implemented", Here());
+  }
+  void applyAD(const oops::Variables& vars, const Increment& inc, const std::vector<double> &) const {
+    throw eckit::NotImplemented("Adjoint interpolation not implemented", Here());
+  }
 
  private:
   void print(std::ostream &) const override;
+  int64_t nlocs_;
   atlas::functionspace::PointCloud atlasObsFuncSpace_;
   atlas::Interpolation interpolator_;
   Parameters_ params_;
 };
 
 }  // namespace orcamodel
-
-
-
-

src/orca-jedi/getvalues/GetValuesParameters.h → src/orca-jedi/interpolator/InterpolatorParameters.h

@@ -26,8 +26,8 @@ class OrcaAtlasInterpolatorParameters : public oops::Parameters {
     max_fraction_elems_to_try{"max_fraction_elems_to_try", this};
 };
 
-class OrcaGetValuesParameters : public oops::Parameters {
-  OOPS_CONCRETE_PARAMETERS(OrcaGetValuesParameters, oops::Parameters)
+class OrcaInterpolatorParameters : public oops::Parameters {
+  OOPS_CONCRETE_PARAMETERS(OrcaInterpolatorParameters, oops::Parameters)
 
  public:
   oops::RequiredParameter<OrcaAtlasInterpolatorParameters> atlasInterpolator{