Skip to content

Commit

Permalink
Adding ocean basin in the Metadata (#732)
Browse files Browse the repository at this point in the history 
#### 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 ...
  • Loading branch information
@guillaumevernieres
guillaumevernieres authored Nov 16, 2023
1 parent edbe844 commit c607b20
Show file tree
Hide file tree
Showing 8 changed files with 259 additions and 173 deletions.
4 changes: 2 additions & 2 deletions utils/obsproc/Ghrsst2Ioda.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down
4 changes: 2 additions & 2 deletions utils/obsproc/IcecAmsr2Ioda.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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;

Expand Down
165 changes: 6 additions & 159 deletions utils/obsproc/NetCDFToIodaConverter.h
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
#pragma once

#include <iostream>
#include <map>
#include <string>
#include <vector>

Expand All @@ -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 {
Expand Down Expand Up @@ -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;
Expand All @@ -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_);

Expand Down Expand Up @@ -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>
Expand Down
32 changes: 26 additions & 6 deletions utils/obsproc/Rads2Ioda.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
4 changes: 2 additions & 2 deletions utils/obsproc/Smap2Ioda.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
4 changes: 2 additions & 2 deletions utils/obsproc/Smos2Ioda.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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());
Expand Down
Loading

0 comments on commit c607b20

Please sign in to comment.