From ceedcce8eda13fc37b61f3c272a4c4f3995e986e Mon Sep 17 00:00:00 2001 From: Yaping Wang <49168260+ypwang19@users.noreply.github.com> Date: Mon, 18 Dec 2023 12:38:38 -0600 Subject: [PATCH] Code clean-up and a reference test for VIIRS AOD IODA converter (#808) ### Description: - Remove duplicate code (NetCDFToIodaConverter2d.h) and implement the channel dimension to the original constructors. - Update the VIIRS ioda converter ctest with reference test. --------- Co-authored-by: ypwang19 --- utils/obsproc/NetCDFToIodaConverter.h | 24 +- utils/obsproc/NetCDFToIodaConverter2D.h | 318 ------------------- utils/obsproc/Viirsaod2Ioda.h | 14 +- utils/obsproc/util.h | 22 +- utils/test/testinput/gdas_viirsaod2ioda.yaml | 6 + utils/test/testref/viirsaod2ioda.test | 26 ++ 6 files changed, 74 insertions(+), 336 deletions(-) delete mode 100644 utils/obsproc/NetCDFToIodaConverter2D.h create mode 100644 utils/test/testref/viirsaod2ioda.test diff --git a/utils/obsproc/NetCDFToIodaConverter.h b/utils/obsproc/NetCDFToIodaConverter.h index ecc6b6aa6..f3061472d 100644 --- a/utils/obsproc/NetCDFToIodaConverter.h +++ b/utils/obsproc/NetCDFToIodaConverter.h @@ -91,8 +91,26 @@ namespace gdasapp { // Update the group with the location dimension ioda::NewDimensionScales_t newDims {ioda::NewDimensionScale("Location", iodaVarsAll.location_)}; + + // Update the group with the location and channel dimension (if channel value is assigned) + if (iodaVars.channelValues_(0) > 0) { + newDims = { + ioda::NewDimensionScale("Location", iodaVarsAll.location_), + ioda::NewDimensionScale("Channel", iodaVarsAll.channel_) + }; + } ioda::ObsGroup ogrp = ioda::ObsGroup::generate(group, newDims); + // Create different dimensionScales for data w/wo channel info + std::vector dimensionScales { + ogrp.vars["Location"] + }; + if (iodaVars.channelValues_(0) > 0) { + dimensionScales = {ogrp.vars["Location"], ogrp.vars["Channel"]}; + ogrp.vars["Channel"].writeWithEigenRegular(iodaVars.channelValues_); + } + + // Set up the creation parameters ioda::VariableCreationParameters float_params = createVariableParams(); ioda::VariableCreationParameters int_params = createVariableParams(); @@ -112,14 +130,14 @@ namespace gdasapp { ioda::Variable iodaObsVal = ogrp.vars.createWithScales("ObsValue/"+variable_, - {ogrp.vars["Location"]}, float_params); + dimensionScales, float_params); ioda::Variable iodaObsErr = ogrp.vars.createWithScales("ObsError/"+variable_, - {ogrp.vars["Location"]}, float_params); + dimensionScales, float_params); ioda::Variable iodaPreQc = ogrp.vars.createWithScales("PreQC/"+variable_, - {ogrp.vars["Location"]}, int_params); + dimensionScales, int_params); // add input filenames to IODA file global attributes ogrp.atts.add("obs_source_files", inputFilenames_); diff --git a/utils/obsproc/NetCDFToIodaConverter2D.h b/utils/obsproc/NetCDFToIodaConverter2D.h deleted file mode 100644 index 0355fc4b9..000000000 --- a/utils/obsproc/NetCDFToIodaConverter2D.h +++ /dev/null @@ -1,318 +0,0 @@ -#pragma once - -#include -#include -#include -#include - -#include "eckit/config/LocalConfiguration.h" - -#include "ioda/Engines/HH.h" -#include "ioda/Group.h" -#include "ioda/ObsGroup.h" - -#include "oops/mpi/mpi.h" -#include "oops/util/DateTime.h" -#include "oops/util/Duration.h" -#include "oops/util/Logger.h" -#include "oops/util/missingValues.h" - -namespace gdasapp { - // A simple data structure to organize the info to provide to the ioda - // writter - struct IodaVars2d { - int location_; // Number of observation per variable - int nVars_; // number of obs variables, should be set to - // for now in the children classes - int channel_; // number of channels - 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 datetime_; // Epoch date in seconds - std::string referenceDate_; // Reference date for epoch time - Eigen::Array channelNumber_; - - // 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 floatMetadataName_; // String descriptor of the float metadata - Eigen::ArrayXXf intMetadata_; // Optional array of integer metadata - std::vector intMetadataName_; // String descriptor of the integer metadata - - // Optional global attributes - std::map strGlobalAttr_; - - // Constructor - explicit IodaVars2d(const int nobs = 0, - const int nchan = 1, - const std::vector fmnames = {}, - const std::vector imnames = {}) : - location_(nobs), nVars_(1), channel_(nchan), - nfMetadata_(fmnames.size()), niMetadata_(imnames.size()), - longitude_(location_), latitude_(location_), datetime_(location_), channelNumber_(channel_), - obsVal_(location_ * channel_), - obsError_(location_ * channel_), - preQc_(location_ * channel_), - floatMetadata_(location_, fmnames.size()), - floatMetadataName_(fmnames), - intMetadata_(location_, imnames.size()), - intMetadataName_(imnames) - { - oops::Log::trace() << "IodaVars2d::IodaVars2d created." << std::endl; - } - - // Append an other instance of IodaVars - void append(const IodaVars2d& 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_ * channel_ + other.location_ * other.channel_); - obsError_.conservativeResize(location_ * channel_ + other.location_ * other.channel_); - preQc_.conservativeResize(location_ * channel_ + other.location_ * other.channel_); - 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.channel_) = other.obsVal_; - obsError_.tail(other.location_ * other.channel_) = other.obsError_; - preQc_.tail(other.location_ * other.channel_) = other.preQc_; - floatMetadata_.bottomRows(other.location_) = other.floatMetadata_; - intMetadata_.bottomRows(other.location_) = other.intMetadata_; - - // Update obs count - location_ += other.location_; - oops::Log::trace() << "IodaVars2d::IodaVars2d done appending." << std::endl; - } - }; - - // Base class for the converters - class NetCDFToIodaConverter2d { - public: - // Constructor: Stores the configuration as a data members - explicit NetCDFToIodaConverter2d(const eckit::Configuration & fullConfig, - const eckit::mpi::Comm & comm): comm_(comm), - fullConfig_(fullConfig) { - // time window info - std::string winbegin; - std::string winend; - fullConfig.get("window begin", winbegin); - fullConfig.get("window end", winend); - windowBegin_ = util::DateTime(winbegin); - windowEnd_ = util::DateTime(winend); - variable_ = "None"; - oops::Log::info() << "--- Window begin: " << winbegin << std::endl; - oops::Log::info() << "--- Window end: " << winend << std::endl; - - // get input netcdf files - fullConfig.get("input files", inputFilenames_); - oops::Log::info() << "--- Input files: " << inputFilenames_ << std::endl; - - // ioda output file name - outputFilename_ = fullConfig.getString("output file"); - oops::Log::info() << "--- Output files: " << outputFilename_ << std::endl; - oops::Log::trace() << "NetCDFToIodaConverter::NetCDFToIodaConverter created." << std::endl; - } - - // Method to write out a IODA file (writter called in destructor) - void writeToIoda() { - // Extract ioda variables from the provider's files - int myrank = comm_.rank(); - int nobs(0); - - // Currently need the PE count to be less than the number of files - ASSERT(comm_.size() <= inputFilenames_.size()); - - // Read the provider's netcdf file - gdasapp::IodaVars2d 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; - oops::Log::info() << " obs count: " << iodaVars.location_ << std::endl; - } - nobs = iodaVars.location_; - - // Get the total number of obs across pe's - int nobsAll(0); - comm_.allReduce(nobs, nobsAll, eckit::mpi::sum()); - gdasapp::IodaVars2d iodaVarsAll(nobsAll, - iodaVars.channel_, - iodaVars.floatMetadataName_, - iodaVars.intMetadataName_); - - // Gather iodaVars arrays - gatherObs(iodaVars.longitude_, iodaVarsAll.longitude_); - gatherObs(iodaVars.latitude_, iodaVarsAll.latitude_); - gatherObs(iodaVars.datetime_, iodaVarsAll.datetime_); - gatherObs(iodaVars.obsVal_, iodaVarsAll.obsVal_); - gatherObs(iodaVars.obsError_, iodaVarsAll.obsError_); - gatherObs(iodaVars.preQc_, iodaVarsAll.preQc_); - - // Create empty group backed by HDF file - if (oops::mpi::world().rank() == 0) { - ioda::Group group = - ioda::Engines::HH::createFile(outputFilename_, - ioda::Engines::BackendCreateModes::Truncate_If_Exists); - - // Update the group with the location dimension - ioda::NewDimensionScales_t - newDims { - ioda::NewDimensionScale("Location", iodaVarsAll.location_), - ioda::NewDimensionScale("Channel", iodaVarsAll.channel_) - }; - - ioda::ObsGroup ogrp = ioda::ObsGroup::generate(group, newDims); - - ogrp.vars["Channel"].writeWithEigenRegular(iodaVars.channelNumber_); - - // Set up the creation parameters - ioda::VariableCreationParameters float_params = createVariableParams(); - ioda::VariableCreationParameters int_params = createVariableParams(); - ioda::VariableCreationParameters long_params = createVariableParams(); - - // Create the mendatory IODA variables - ioda::Variable iodaDatetime = - ogrp.vars.createWithScales("MetaData/dateTime", - {ogrp.vars["Location"]}, long_params); - iodaDatetime.atts.add("units", {iodaVars.referenceDate_}, {1}); - ioda::Variable iodaLat = - ogrp.vars.createWithScales("MetaData/latitude", - {ogrp.vars["Location"]}, float_params); - ioda::Variable iodaLon = - ogrp.vars.createWithScales("MetaData/longitude", - {ogrp.vars["Location"]}, float_params); - - ioda::Variable iodaObsVal = - ogrp.vars.createWithScales("ObsValue/"+variable_, - {ogrp.vars["Location"], - ogrp.vars["Channel"]}, float_params); - ioda::Variable iodaObsErr = - ogrp.vars.createWithScales("ObsError/"+variable_, - {ogrp.vars["Location"], - ogrp.vars["Channel"]}, float_params); - ioda::Variable iodaPreQc = - ogrp.vars.createWithScales("PreQC/"+variable_, - {ogrp.vars["Location"], - ogrp.vars["Channel"]}, int_params); - - // add input filenames to IODA file global attributes - ogrp.atts.add("obs_source_files", inputFilenames_); - - // add global attributes collected from the specific converter - for (const auto& globalAttr : iodaVars.strGlobalAttr_) { - ogrp.atts.add(globalAttr.first , globalAttr.second); - } - - // Create the optional IODA integer metadata - ioda::Variable tmpIntMeta; - int count = 0; - for (const std::string& strMeta : iodaVars.intMetadataName_) { - tmpIntMeta = ogrp.vars.createWithScales("MetaData/"+strMeta, - {ogrp.vars["Location"]}, int_params); - tmpIntMeta.writeWithEigenRegular(iodaVars.intMetadata_.col(count)); - count++; - } - - // Create the optional IODA float metadata - ioda::Variable tmpFloatMeta; - count = 0; - for (const std::string& strMeta : iodaVars.floatMetadataName_) { - tmpFloatMeta = ogrp.vars.createWithScales("MetaData/"+strMeta, - {ogrp.vars["Location"]}, float_params); - tmpFloatMeta.writeWithEigenRegular(iodaVars.floatMetadata_.col(count)); - count++; - } - - // Write obs info to group - oops::Log::info() << "Writing ioda file" << std::endl; - iodaLon.writeWithEigenRegular(iodaVarsAll.longitude_); - iodaLat.writeWithEigenRegular(iodaVarsAll.latitude_); - iodaDatetime.writeWithEigenRegular(iodaVarsAll.datetime_); - iodaObsVal.writeWithEigenRegular(iodaVarsAll.obsVal_); - iodaObsErr.writeWithEigenRegular(iodaVarsAll.obsError_); - iodaPreQc.writeWithEigenRegular(iodaVarsAll.preQc_); - } - } - - private: - // Virtual method that reads the provider's netcdf file and store the relevant - // info in a IodaVars struct - virtual gdasapp::IodaVars2d providerToIodaVars(const std::string fileName) = 0; - - // Gather for eigen array - template - void gatherObs(const Eigen::Array & obsPe, - Eigen::Array & obsAllPes) { - // define root pe - const size_t root = 0; - - // send pointer to the PE's data - std::vector send(obsPe.data(), obsPe.data() + obsPe.size()); - size_t ntasks = comm_.size(); - - // gather the sizes of the send buffers - int localnobs = send.size(); - std::vector sizes(ntasks); - comm_.allGather(localnobs, sizes.begin(), sizes.end()); - - // displacement for the received data - std::vector displs(ntasks); - size_t rcvsz = sizes[0]; - displs[0] = 0; - for (size_t jj = 1; jj < ntasks; ++jj) { - displs[jj] = displs[jj - 1] + sizes[jj - 1]; - rcvsz += sizes[jj]; - } - - // create receiving buffer - std::vector recv(0); - - // gather all send buffers - if (comm_.rank() == root) recv.resize(rcvsz); - comm_.barrier(); - comm_.gatherv(send, recv, sizes, displs, root); - - if (comm_.rank() == root) { - obsAllPes.segment(0, recv.size()) = - Eigen::Map>(recv.data(), recv.size()); - } - } - - // Short-cut to create type dependent VariableCreationParameters - template - ioda::VariableCreationParameters createVariableParams() { - ioda::VariableCreationParameters params; - params.chunk = true; // allow chunking - params.compressWithGZIP(); // compress using gzip - params.setFillValue(util::missingValue()); - - return params; - } - - protected: - util::DateTime windowBegin_; - util::DateTime windowEnd_; - std::vector inputFilenames_; - std::string outputFilename_; - std::string variable_; - const eckit::mpi::Comm & comm_; - const eckit::Configuration & fullConfig_; - }; -} // namespace gdasapp diff --git a/utils/obsproc/Viirsaod2Ioda.h b/utils/obsproc/Viirsaod2Ioda.h index 4bea0bd85..453a1d3cd 100644 --- a/utils/obsproc/Viirsaod2Ioda.h +++ b/utils/obsproc/Viirsaod2Ioda.h @@ -12,20 +12,20 @@ #include "ioda/Group.h" #include "ioda/ObsGroup.h" -#include "NetCDFToIodaConverter2D.h" +#include "NetCDFToIodaConverter.h" #include "superob.h" namespace gdasapp { - class Viirsaod2Ioda : public NetCDFToIodaConverter2d { + class Viirsaod2Ioda : public NetCDFToIodaConverter { public: explicit Viirsaod2Ioda(const eckit::Configuration & fullConfig, const eckit::mpi::Comm & comm) - : NetCDFToIodaConverter2d(fullConfig, comm) { + : NetCDFToIodaConverter(fullConfig, comm) { variable_ = "aerosolOpticalDepth"; } // Read netcdf file and populate iodaVars - gdasapp::IodaVars2d providerToIodaVars(const std::string fileName) final { + gdasapp::obsproc::iodavars::IodaVars providerToIodaVars(const std::string fileName) final { oops::Log::info() << "Processing files provided by VIIRSAOD" << std::endl; // Open the NetCDF file in read-only mode @@ -153,12 +153,12 @@ namespace gdasapp { int nchan(channelNumber.size()); oops::Log::info() << " number of channels " << nchan << std::endl; // Create instance of iodaVars object - gdasapp::IodaVars2d iodaVars(nobs, nchan, {}, {}); + gdasapp::obsproc::iodavars::IodaVars iodaVars(nobs, {}, {}); - oops::Log::info() << " eigen... " << std::endl; + oops::Log::info() << " eigen... " << dimRow << dimCol << std::endl; // Store into eigen arrays for (int k = 0; k < nchan; k++) { - iodaVars.channelNumber_(k) = channelNumber[k]; + iodaVars.channelValues_(k) = channelNumber[k]; int loc(0); for (int i = 0; i < dimRow; i++) { for (int j = 0; j < dimCol; j++) { diff --git a/utils/obsproc/util.h b/utils/obsproc/util.h index 13794a1ad..d44cc2356 100644 --- a/utils/obsproc/util.h +++ b/utils/obsproc/util.h @@ -92,6 +92,10 @@ namespace gdasapp { int nfMetadata_; // number of float metadata fields int niMetadata_; // number of int metadata fields + // Channels + int channel_; // Number of channels + Eigen::ArrayXi channelValues_; // Values specific to channels + // Non optional metadata Eigen::ArrayXf longitude_; // geo-location_ Eigen::ArrayXf latitude_; // " @@ -115,16 +119,18 @@ namespace gdasapp { // Constructor explicit IodaVars(const int nobs = 0, const std::vector fmnames = {}, - const std::vector imnames = {}) : + const std::vector imnames = {}): location_(nobs), nVars_(1), nfMetadata_(fmnames.size()), niMetadata_(imnames.size()), longitude_(location_), latitude_(location_), datetime_(location_), - obsVal_(location_), - obsError_(location_), - preQc_(location_), + obsVal_(location_*channel_), + obsError_(location_*channel_), + preQc_(location_*channel_), floatMetadata_(location_, fmnames.size()), floatMetadataName_(fmnames), intMetadata_(location_, imnames.size()), - intMetadataName_(imnames) + intMetadataName_(imnames), + channel_(1), + channelValues_(Eigen::ArrayXi::Constant(channel_, -1)) { oops::Log::trace() << "IodaVars::IodaVars created." << std::endl; } @@ -142,9 +148,9 @@ namespace gdasapp { 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_); + obsVal_.conservativeResize(location_ * channel_ + other.location_ * other.channel_); + obsError_.conservativeResize(location_ * channel_ + other.location_ * other.channel_); + preQc_.conservativeResize(location_ * channel_ + other.location_ * other.channel_); floatMetadata_.conservativeResize(location_ + other.location_, nfMetadata_); intMetadata_.conservativeResize(location_ + other.location_, niMetadata_); diff --git a/utils/test/testinput/gdas_viirsaod2ioda.yaml b/utils/test/testinput/gdas_viirsaod2ioda.yaml index ac99cf019..073b6f799 100644 --- a/utils/test/testinput/gdas_viirsaod2ioda.yaml +++ b/utils/test/testinput/gdas_viirsaod2ioda.yaml @@ -9,3 +9,9 @@ input files: thinning: threshold: 0 channel: 4 + +test: + reference filename: testref/viirsaod2ioda.test + test output filename: testoutput/viirsaod2ioda.test + float relative tolerance: 1e-6 + diff --git a/utils/test/testref/viirsaod2ioda.test b/utils/test/testref/viirsaod2ioda.test new file mode 100644 index 000000000..105d013d2 --- /dev/null +++ b/utils/test/testref/viirsaod2ioda.test @@ -0,0 +1,26 @@ +Reading: [viirs_aod_1.nc4,viirs_aod_2.nc4] +1970-01-01 00:00:00 +obsVal: + Min: 0.00219523 + Max: 0.163349 + Sum: 11.8086 +obsError: + Min: 0.00885619 + Max: 0.130309 + Sum: 8.75665 +preQc: + Min: 0 + Max: 2 + Sum: 105 +longitude: + Min: 141.747 + Max: 142.872 + Sum: 19953.1 +latitude: + Min: 50.2839 + Max: 51.3225 + Sum: 7140.63 +datetime: + Min: 1523765998 + Max: 1523765998 + Sum: 213327239720