Adding ocean basin in the Metadata (#732)

guillaumevernieres · web-flow · commit c607b20e573a · 2023-11-16T14:56:17.000-05:00
#### What was done - Reorganization of what's in `NetCDFToIodaConverter.h`: Moved structs and functions to `util.h` as an attempt to make the code slightly cleaner - Added a simple OceanMask struct that reads in regional ocean masks from [reccap2 ocean mask](https://reccap2-ocean.github.io/regions/). The struct has a method to return the integer mask at a specific location. - Fixed a bug in `trim`, related to the optional float metadata - implementation of the new feature in the RADS converter #### Notes It's currently not part of the `CI`, something to think about in the future. We could add `RECCAP2_region_masks_all_v20221025.nc` to the tarball or stage it on hpc somewhere ...
diff --git a/utils/obsproc/Ghrsst2Ioda.h b/utils/obsproc/Ghrsst2Ioda.h
@@ -26,7 +26,7 @@ namespace gdasapp {
     }
 
     // Read netcdf file and populate iodaVars
-    gdasapp::IodaVars providerToIodaVars(const std::string fileName) final {
+    gdasapp::obsproc::iodavars::IodaVars providerToIodaVars(const std::string fileName) final {
       oops::Log::info() << "Processing files provided by GHRSST" << std::endl;
 
       // Get the sst bounds from the configuration
@@ -182,7 +182,7 @@ namespace gdasapp {
       int nobs = sst_s.size() * sst_s[0].size();
 
       // Create instance of iodaVars object
-      gdasapp::IodaVars iodaVars(nobs, {}, {});
+      gdasapp::obsproc::iodavars::IodaVars iodaVars(nobs, {}, {});
 
       // Reference time is Jan 01 1981 00:00:00 GMT+0000
       iodaVars.referenceDate_ = refDate;
diff --git a/utils/obsproc/IcecAmsr2Ioda.h b/utils/obsproc/IcecAmsr2Ioda.h
@@ -26,7 +26,7 @@ namespace gdasapp {
     }
 
     // Read netcdf file and populate iodaVars
-    gdasapp::IodaVars providerToIodaVars(const std::string fileName) final {
+    gdasapp::obsproc::iodavars::IodaVars providerToIodaVars(const std::string fileName) final {
       oops::Log::info() << "Processing files provided by the AMSR2" << std::endl;
 
       // Open the NetCDF file in read-only mode
@@ -40,7 +40,7 @@ namespace gdasapp {
       int ntimes = dimxSize * dimySize * dimTimeSize;
 
       // Create instance of iodaVars object
-      gdasapp::IodaVars iodaVars(nobs, {}, {});
+      gdasapp::obsproc::iodavars::IodaVars iodaVars(nobs, {}, {});
 
       oops::Log::debug() << "--- iodaVars.location_: " << iodaVars.location_ << std::endl;
 
diff --git a/utils/obsproc/NetCDFToIodaConverter.h b/utils/obsproc/NetCDFToIodaConverter.h
@@ -1,7 +1,6 @@
 #pragma once
 
 #include <iostream>
-#include <map>
 #include <string>
 #include <vector>
 
@@ -17,162 +16,9 @@
 #include "oops/util/Logger.h"
 #include "oops/util/missingValues.h"
 
-namespace gdasapp {
-  namespace testutils {
-    template <typename Derived>
-    std::string checksum(const Eigen::ArrayBase<Derived>& arr, const std::string varname) {
-      std::stringstream result;
-      if (arr.size() == 0) {
-        result << varname << " is empty" << "\n";
-      } else {
-        auto minElement = arr.minCoeff();
-        auto maxElement = arr.maxCoeff();
-        auto sumElements = arr.sum();
-
-        result << varname << ":" << "\n";
-        result << "    Min: " << minElement << "\n";
-        result << "    Max: " << maxElement << "\n";
-        result << "    Sum: " << sumElements;
-      }
-      return result.str();
-    }
-  }  // namespace testutils
-
-  // A simple data structure to organize the info to provide to the ioda
-  // writter
-  struct IodaVars {
-    int location_;      // Number of observation per variable
-    int nVars_;         // number of obs variables, should be set to
-                       // for now in the children classes
-    int nfMetadata_;    // number of float metadata fields
-    int niMetadata_;    // number of int metadata fields
-
-    // Non optional metadata
-    Eigen::ArrayXf longitude_;  // geo-location_
-    Eigen::ArrayXf latitude_;   //     "
-    Eigen::Array<int64_t, Eigen::Dynamic, 1> datetime_;   // Epoch date in seconds
-    std::string referenceDate_;                           // Reference date for epoch time
-
-    // Obs info
-    Eigen::ArrayXf obsVal_;     // Observation value
-    Eigen::ArrayXf obsError_;   //      "      error
-    Eigen::ArrayXi preQc_;      // Quality control flag
-
-    // Optional metadata
-    Eigen::ArrayXXf floatMetadata_;                // Optional array of float metadata
-    std::vector<std::string> floatMetadataName_;  // String descriptor of the float metadata
-    Eigen::ArrayXXf intMetadata_;                  // Optional array of integer metadata
-    std::vector<std::string> intMetadataName_;    // String descriptor of the integer metadata
-
-    // Optional global attributes
-    std::map<std::string, std::string> strGlobalAttr_;
-
-    // Constructor
-    explicit IodaVars(const int nobs = 0,
-                      const std::vector<std::string> fmnames = {},
-                      const std::vector<std::string> imnames = {}) :
-      location_(nobs), nVars_(1), nfMetadata_(fmnames.size()), niMetadata_(imnames.size()),
-      longitude_(location_), latitude_(location_), datetime_(location_),
-      obsVal_(location_),
-      obsError_(location_),
-      preQc_(location_),
-      floatMetadata_(location_, fmnames.size()),
-      floatMetadataName_(fmnames),
-      intMetadata_(location_, imnames.size()),
-      intMetadataName_(imnames)
-    {
-      oops::Log::trace() << "IodaVars::IodaVars created." << std::endl;
-    }
-
-    // Append an other instance of IodaVars
-    void append(const IodaVars& other) {
-      // Check if the two instances can be concatenated
-      ASSERT(nVars_ == other.nVars_);
-      ASSERT(nfMetadata_ == other.nfMetadata_);
-      ASSERT(niMetadata_ == other.niMetadata_);
-      ASSERT(floatMetadataName_ == floatMetadataName_);
-      ASSERT(intMetadataName_ == intMetadataName_);
-
-      // Concatenate Eigen arrays and vectors
-      longitude_.conservativeResize(location_ + other.location_);
-      latitude_.conservativeResize(location_ + other.location_);
-      datetime_.conservativeResize(location_ + other.location_);
-      obsVal_.conservativeResize(location_ + other.location_);
-      obsError_.conservativeResize(location_ + other.location_);
-      preQc_.conservativeResize(location_ + other.location_);
-      floatMetadata_.conservativeResize(location_ + other.location_, nfMetadata_);
-      intMetadata_.conservativeResize(location_ + other.location_, niMetadata_);
-
-      // Copy data from the 'other' object to this object
-      longitude_.tail(other.location_) = other.longitude_;
-      latitude_.tail(other.location_) = other.latitude_;
-      datetime_.tail(other.location_) = other.datetime_;
-      obsVal_.tail(other.location_) = other.obsVal_;
-      obsError_.tail(other.location_) = other.obsError_;
-      preQc_.tail(other.location_) = other.preQc_;
-      floatMetadata_.bottomRows(other.location_) = other.floatMetadata_;
-      intMetadata_.bottomRows(other.location_) = other.intMetadata_;
-
-      // Update obs count
-      location_ += other.location_;
-      oops::Log::trace() << "IodaVars::IodaVars done appending." << std::endl;
-    }
-    void trim(const Eigen::Array<bool, Eigen::Dynamic, 1>& mask ) {
-      int newlocation = mask.count();
-
-      IodaVars iodaVarsMasked(newlocation,  floatMetadataName_, intMetadataName_);
-
-      // this feels downright crude, but apparently numpy-type masks are on the todo list for Eigen
-      int j = 0;
-      for (int i = 0; i < location_; i++) {
-        if (mask(i)) {
-          iodaVarsMasked.longitude_(j) = longitude_(i);
-          iodaVarsMasked.latitude_(j) = latitude_(i);
-          iodaVarsMasked.obsVal_(j) = obsVal_(i);
-          iodaVarsMasked.obsError_(j) = obsError_(i);
-          iodaVarsMasked.preQc_(j) = preQc_(i);
-          iodaVarsMasked.datetime_(j) = datetime_(i);
-          for (int k = 0; k < nfMetadata_; k++) {
-            iodaVarsMasked.intMetadata_(j, k) = intMetadata_(i, k);
-            }
-          for (int k = 0; k < niMetadata_; k++) {
-            iodaVarsMasked.intMetadata_(j, k) = intMetadata_(i, k);
-            }
-        j++;
-        }  // end if (mask(i))
-      }
-
-      longitude_ = iodaVarsMasked.longitude_;
-      latitude_ = iodaVarsMasked.latitude_;
-      datetime_ = iodaVarsMasked.datetime_;
-      obsVal_ = iodaVarsMasked.obsVal_;
-      obsError_ = iodaVarsMasked.obsError_;
-      preQc_ = iodaVarsMasked.preQc_;
-      floatMetadata_ = iodaVarsMasked.floatMetadata_;
-      intMetadata_ = iodaVarsMasked.intMetadata_;
-
-      // Update obs count
-      location_ = iodaVarsMasked.location_;
-      oops::Log::info() << "IodaVars::IodaVars done masking." << std::endl;
-    }
+# include "util.h"
 
-    // Testing
-    void testOutput() {
-      oops::Log::test() << referenceDate_ << std::endl;
-      oops::Log::test() <<
-        gdasapp::testutils::checksum(obsVal_, "obsVal") << std::endl;
-      oops::Log::test() <<
-        gdasapp::testutils::checksum(obsError_, "obsError") << std::endl;
-      oops::Log::test() <<
-        gdasapp::testutils::checksum(preQc_, "preQc") << std::endl;
-      oops::Log::test() <<
-        gdasapp::testutils::checksum(longitude_, "longitude") << std::endl;
-      oops::Log::test() <<
-        gdasapp::testutils::checksum(latitude_, "latitude") << std::endl;
-      oops::Log::test() <<
-        gdasapp::testutils::checksum(datetime_, "datetime") << std::endl;
-    }
-  };
+namespace gdasapp {
 
   // Base class for the converters
   class NetCDFToIodaConverter {
@@ -212,7 +58,7 @@ namespace gdasapp {
       ASSERT(comm_.size() <= inputFilenames_.size());
 
       // Read the provider's netcdf file
-      gdasapp::IodaVars iodaVars = providerToIodaVars(inputFilenames_[myrank]);
+      gdasapp::obsproc::iodavars::IodaVars iodaVars = providerToIodaVars(inputFilenames_[myrank]);
       for (int i = myrank + comm_.size(); i < inputFilenames_.size(); i += comm_.size()) {
         iodaVars.append(providerToIodaVars(inputFilenames_[i]));
         oops::Log::info() << " appending: " << inputFilenames_[i] << std::endl;
@@ -223,7 +69,7 @@ namespace gdasapp {
       // Get the total number of obs across pe's
       int nobsAll(0);
       comm_.allReduce(nobs, nobsAll, eckit::mpi::sum());
-      gdasapp::IodaVars iodaVarsAll(nobsAll,
+      gdasapp::obsproc::iodavars::IodaVars iodaVarsAll(nobsAll,
                                     iodaVars.floatMetadataName_,
                                     iodaVars.intMetadataName_);
 
@@ -316,7 +162,8 @@ namespace gdasapp {
    private:
     // Virtual method that reads the provider's netcdf file and store the relevant
     // info in a IodaVars struct
-    virtual gdasapp::IodaVars  providerToIodaVars(const std::string fileName) = 0;
+    virtual gdasapp::obsproc::iodavars::IodaVars providerToIodaVars(
+                                                        const std::string fileName) = 0;
 
     // Gather for eigen array
     template <typename T>
diff --git a/utils/obsproc/Rads2Ioda.h b/utils/obsproc/Rads2Ioda.h
@@ -25,7 +25,7 @@ namespace gdasapp {
     }
 
     // Read netcdf file and populate iodaVars
-    gdasapp::IodaVars providerToIodaVars(const std::string fileName) final {
+    gdasapp::obsproc::iodavars::IodaVars providerToIodaVars(const std::string fileName) final {
       oops::Log::info() << "Processing files provided by the RADS" << std::endl;
 
       // Open the NetCDF file in read-only mode
@@ -35,13 +35,13 @@ namespace gdasapp {
       int nobs = ncFile.getDim("time").getSize();
 
       // Set the int metadata names
-      std::vector<std::string> intMetadataNames = {"pass", "cycle", "mission"};
+      std::vector<std::string> intMetadataNames = {"pass", "cycle", "mission", "oceanBasin"};
 
       // Set the float metadata name
       std::vector<std::string> floatMetadataNames = {"mdt"};
 
       // Create instance of iodaVars object
-      gdasapp::IodaVars iodaVars(nobs, floatMetadataNames, intMetadataNames);
+      gdasapp::obsproc::iodavars::IodaVars iodaVars(nobs, floatMetadataNames, intMetadataNames);
 
       // Read non-optional metadata: datetime, longitude and latitude
       int lat[iodaVars.location_];  // NOLINT
@@ -89,8 +89,9 @@ namespace gdasapp {
       int cycle[iodaVars.location_];  // NOLINT
       ncFile.getVar("cycle").getVar(cycle);
 
+      // Store optional metadata, set ocean basins to -999 for now
       for (int i = 0; i < iodaVars.location_; i++) {
-        iodaVars.intMetadata_.row(i) << pass[i], cycle[i], mission_index;
+        iodaVars.intMetadata_.row(i) << pass[i], cycle[i], mission_index, -999;
       }
 
       // Get adt_egm2008 obs values and attributes
@@ -100,7 +101,7 @@ namespace gdasapp {
       ncFile.getVar("adt_egm2008").getAtt("scale_factor").getValues(&scaleFactor);
 
       // Read sla
-      int sla[iodaVars.location_];  // NOLINT
+      short sla[iodaVars.location_];  // NOLINT
       ncFile.getVar("sla").getVar(sla);
 
       // Update non-optional Eigen arrays
@@ -112,9 +113,28 @@ namespace gdasapp {
         iodaVars.obsError_(i) = 0.1;  // Do something for obs error
         iodaVars.preQc_(i) = 0;
         // Save MDT in optional floatMetadata
-        iodaVars.floatMetadata_(i, 0) = iodaVars.obsVal_(i) -
+        iodaVars.floatMetadata_.row(i) << iodaVars.obsVal_(i) -
                                         static_cast<float>(sla[i])*scaleFactor;
       }
+
+      // Basic QC
+      Eigen::Array<bool, Eigen::Dynamic, 1> boundsCheck =
+        (iodaVars.obsVal_ > -4.0 && iodaVars.obsVal_ < 4.0);
+      iodaVars.trim(boundsCheck);
+
+      // get ocean basin masks if asked in the config
+      obsproc::oceanmask::OceanMask* oceanMask = nullptr;
+      if (fullConfig_.has("ocean basin")) {
+          std::string fileName;
+          fullConfig_.get("ocean basin", fileName);
+          oceanMask = new obsproc::oceanmask::OceanMask(fileName);
+
+          for (int i = 0; i < iodaVars.location_; i++) {
+            iodaVars.intMetadata_.coeffRef(i, 3) =
+              oceanMask->getOceanMask(iodaVars.longitude_[i], iodaVars.latitude_[i]);
+          }
+      }
+
       return iodaVars;
     };
   };  // class Rads2Ioda
diff --git a/utils/obsproc/Smap2Ioda.h b/utils/obsproc/Smap2Ioda.h
@@ -25,7 +25,7 @@ namespace gdasapp {
     }
 
     // Read netcdf file and populate iodaVars
-    gdasapp::IodaVars providerToIodaVars(const std::string fileName) final {
+    gdasapp::obsproc::iodavars::IodaVars providerToIodaVars(const std::string fileName) final {
       oops::Log::info() << "Processing files provided by SMAP" << std::endl;
 
       // Open the NetCDF file in read-only mode
@@ -36,7 +36,7 @@ namespace gdasapp {
       int dim1  = ncFile.getDim("phony_dim_1").getSize();
       int nobs = dim0 * dim1;
 
-      gdasapp::IodaVars iodaVars(nobs, {}, {});
+      gdasapp::obsproc::iodavars::IodaVars iodaVars(nobs, {}, {});
 
       // TODO(AFE): these arrays can be done as 1D vectors, but those need proper ushorts in
       // the input files, at odd with the current ctests
diff --git a/utils/obsproc/Smos2Ioda.h b/utils/obsproc/Smos2Ioda.h
@@ -24,7 +24,7 @@ namespace gdasapp {
     }
 
     // Read netcdf file and populate iodaVars
-    gdasapp::IodaVars providerToIodaVars(const std::string fileName) final {
+    gdasapp::obsproc::iodavars::IodaVars providerToIodaVars(const std::string fileName) final {
       oops::Log::info() << "Processing files provided by SMOS" << std::endl;
 
       // Open the NetCDF file in read-only mode
@@ -43,7 +43,7 @@ namespace gdasapp {
       // std::vector<std::string> floatMetadataNames = {"mdt"};
       std::vector<std::string> floatMetadataNames = {};
       // Create instance of iodaVars object
-      gdasapp::IodaVars iodaVars(nobs, floatMetadataNames, intMetadataNames);
+      gdasapp::obsproc::iodavars::IodaVars iodaVars(nobs, floatMetadataNames, intMetadataNames);
 
       std::vector<float> lat(iodaVars.location_);
       ncFile.getVar("Latitude").getVar(lat.data());
diff --git a/utils/obsproc/util.h b/utils/obsproc/util.h
diff --git a/utils/test/testinput/gdas_rads2ioda.yaml b/utils/test/testinput/gdas_rads2ioda.yaml
