Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

IODA Converter with In-Situ Observations Concatenation and ObsError Inflation #1472

Merged
merged 15 commits into from
Feb 13, 2025
Merged
148 changes: 148 additions & 0 deletions utils/obsproc/InsituAll2ioda.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,148 @@
#pragma once

#include <iostream>
#include <map>
#include <netcdf> // NOLINT (using C API)
#include <string>
#include <vector>

#include "eckit/config/LocalConfiguration.h"

#include <Eigen/Dense> // NOLINT

#include "ioda/Group.h"
#include "ioda/ObsGroup.h"

#include "NetCDFToIodaConverter.h"

namespace gdasapp {

class InsituAll2Ioda : public NetCDFToIodaConverter {
public:
explicit InsituAll2Ioda(const eckit::Configuration &fullConfig, const eckit::mpi::Comm &comm)
: NetCDFToIodaConverter(fullConfig, comm) {
ASSERT(fullConfig_.has("variable"));
fullConfig_.get("variable", variable_);
}

// Read NetCDF file and populate data based on YAML configuration
gdasapp::obsproc::iodavars::IodaVars providerToIodaVars(const std::string fileName) final {
oops::Log::info() << "Processing files provided from ALL in-situ files" << std::endl;

// Abort the case where the 'error ratio' key is not found
ASSERT(fullConfig_.has("error ratio"));

// Get the obs. error ratio from the configuration (unit per day)
float errRatio;
fullConfig_.get("error ratio", errRatio);
// Convert errRatio from meters per day to its unit per second
errRatio /= 86400.0;

// Open the NetCDF file in read-only mode
netCDF::NcFile ncFile(fileName, netCDF::NcFile::read);
oops::Log::info() << "Reading... " << fileName << std::endl;

// Get the number of obs in the file
int nobs = ncFile.getDim("Location").getSize();

// Set the int metadata names
std::vector<std::string> intMetadataNames = {"oceanBasin"};

// Set the float metadata name
std::vector<std::string> floatMetadataNames = {"depth"};

// Create instance of iodaVars object
gdasapp::obsproc::iodavars::IodaVars iodaVars(nobs, floatMetadataNames, intMetadataNames);

// Check if the MetaData group is null
netCDF::NcGroup metaDataGroup = ncFile.getGroup("MetaData");
if (metaDataGroup.isNull()) {
oops::Log::debug() << "Group 'MetaData' not found!" << std::endl;
}

// Read non-optional metadata: datetime, longitude, latitude and optional: others
netCDF::NcVar latitudeVar = metaDataGroup.getVar("latitude");
std::vector<float> latitudeData(iodaVars.location_);
latitudeVar.getVar(latitudeData.data());

netCDF::NcVar longitudeVar = metaDataGroup.getVar("longitude");
std::vector<float> longitudeData(iodaVars.location_);
longitudeVar.getVar(longitudeData.data());

netCDF::NcVar datetimeVar = metaDataGroup.getVar("dateTime");
std::vector<int64_t> datetimeData(iodaVars.location_);
datetimeVar.getVar(datetimeData.data());
iodaVars.referenceDate_ = "seconds since 1970-01-01T00:00:00Z"; // Applied to All in-situ obs

netCDF::NcVar depthVar = metaDataGroup.getVar("depth");
std::vector<float> depthData(iodaVars.location_, 0); // Initialize with surface value 0

if (!depthVar.isNull()) { // Checking from surface in-situ obs
oops::Log::info() << "Variable 'depth' found and Reading!" << std::endl;
depthVar.getVar(depthData.data());
} else {
oops::Log::info() << "Variable 'depth' NOT found in 'MetaData'!" << std::endl;
}

// Save in optional floatMetadata
for (int i = 0; i < iodaVars.location_; i++) {
iodaVars.floatMetadata_.row(i) << depthData[i];
}

netCDF::NcVar oceanbasinVar = metaDataGroup.getVar("oceanBasin");
std::vector<int> oceanbasinData(iodaVars.location_);
oceanbasinVar.getVar(oceanbasinData.data());

// Check if the groups are null
netCDF::NcGroup obsvalGroup = ncFile.getGroup("ObsValue");
netCDF::NcGroup obserrGroup = ncFile.getGroup("ObsError");
netCDF::NcGroup preqcGroup = ncFile.getGroup("PreQC");
if (obsvalGroup.isNull()) {
oops::Log::debug() << "Group 'ObsValue' not found!" << std::endl;
} else if (obserrGroup.isNull()) {
oops::Log::debug() << "Group 'ObsError' not found!" << std::endl;
} else if (preqcGroup.isNull()) {
oops::Log::debug() << "Group 'PreQC' not found!" << std::endl;
}

// Get obs values, errors and preqc
netCDF::NcVar obsvalVar = obsvalGroup.getVar(variable_);
std::vector<float> obsvalData(iodaVars.location_);
obsvalVar.getVar(obsvalData.data());

netCDF::NcVar obserrVar = obserrGroup.getVar(variable_);
std::vector<float> obserrData(iodaVars.location_);
obserrVar.getVar(obserrData.data());

netCDF::NcVar preqcVar = preqcGroup.getVar(variable_);
std::vector<int> preqcData(iodaVars.location_);
preqcVar.getVar(preqcData.data());

// Update non-optional Eigen arrays
for (int i = 0; i < iodaVars.location_; i++) {
iodaVars.longitude_(i) = longitudeData[i];
iodaVars.latitude_(i) = latitudeData[i];
iodaVars.datetime_(i) = datetimeData[i];
iodaVars.obsVal_(i) = obsvalData[i];
iodaVars.obsError_(i) = obserrData[i];
iodaVars.preQc_(i) = preqcData[i];
// Save in optional intMetadata
iodaVars.intMetadata_.row(i) << oceanbasinData[i];
}

// Extract epochdate format from referenceDate
std::string extractedDate = iodaVars.referenceDate_.substr(14);

// Redating and adjusting Errors
if (iodaVars.datetime_.size() == 0) {
oops::Log::info() << "datetime_ is empty" << std::endl;
} else {
// Redating and Adjusting Error
iodaVars.reDate(windowBegin_, windowEnd_, extractedDate, errRatio);
}

return iodaVars;
};
}; // class InsituAll2Ioda
} // namespace gdasapp

4 changes: 3 additions & 1 deletion utils/obsproc/Rads2Ioda.h
Original file line number Diff line number Diff line change
Expand Up @@ -132,12 +132,14 @@ namespace gdasapp {
(iodaVars.obsVal_ > -4.0 && iodaVars.obsVal_ < 4.0);
iodaVars.trim(boundsCheck);

std::string extractedDate = iodaVars.referenceDate_.substr(14);

// Redating and adjusting Errors
if (iodaVars.datetime_.size() == 0) {
oops::Log::info() << "datetime_ is empty" << std::endl;
} else {
// Redating and Adjusting Error
iodaVars.reDate(windowBegin_, windowEnd_, errRatio);
iodaVars.reDate(windowBegin_, windowEnd_, extractedDate, errRatio);
}

return iodaVars;
Expand Down
4 changes: 4 additions & 0 deletions utils/obsproc/applications/gdas_obsprovider2ioda.h
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
#include "../IcecAmsr2Ioda.h"
#include "../IcecJpssrr2Ioda.h"
#include "../IcecMirs2Ioda.h"
#include "../InsituAll2ioda.h"
#include "../Rads2Ioda.h"
#include "../RTOFSSalinity.h"
#include "../RTOFSTemperature.h"
Expand Down Expand Up @@ -65,6 +66,9 @@ namespace gdasapp {
} else if (provider == "VIIRSAOD") {
Viirsaod2Ioda conv2ioda(fullConfig, this->getComm());
conv2ioda.writeToIoda();
} else if (provider == "INSITUOBS") {
InsituAll2Ioda conv2ioda(fullConfig, this->getComm());
conv2ioda.writeToIoda();
} else {
oops::Log::info() << "Provider not implemented" << std::endl;
return 1;
Expand Down
4 changes: 2 additions & 2 deletions utils/obsproc/util.h
Original file line number Diff line number Diff line change
Expand Up @@ -225,11 +225,11 @@ namespace gdasapp {

// Changing the date and Adjusting Errors
void reDate(const util::DateTime & windowBegin, const util::DateTime & windowEnd,
float errRatio) {
const std::string &epochDate, float errRatio) {
// windowBegin and End into DAwindowTimes
std::vector<util::DateTime> DAwindowTimes = {windowBegin, windowEnd};
// Epoch DateTime from Provider
util::DateTime epochDtime("1858-11-17T00:00:00Z");
util::DateTime epochDtime(epochDate);
// Convert DA Window DateTime objects to epoch time offsets in seconds
std::vector<int64_t> timeOffsets
= ioda::convertDtimeToTimeOffsets(epochDtime, DAwindowTimes);
Expand Down
10 changes: 9 additions & 1 deletion utils/test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -174,10 +174,18 @@ ecbuild_add_test( TARGET test_gdasapp_util_icecmirs2ioda
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc
TEST_DEPENDS test_gdasapp_util_prepdata)

# Test the JPSSRR to IODA converter
# Test the JPSSRR to IODA converter
ecbuild_add_test( TARGET test_gdasapp_util_icecjpssrr2ioda
COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x
ARGS "../testinput/gdas_icecjpssrr2ioda.yaml"
LIBS gdas-utils
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc
TEST_DEPENDS test_gdasapp_util_prepdata)

# Test the INSITU to IODA converter
ecbuild_add_test( TARGET test_gdasapp_util_insituall2ioda
COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x
ARGS "../testinput/gdas_insituall2ioda.yaml"
LIBS gdas-utils
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc
TEST_DEPENDS test_gdasapp_util_prepdata)
19 changes: 19 additions & 0 deletions utils/test/testinput/gdas_insituall2ioda.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
provider: INSITUOBS
window begin: 2021-03-24T21:00:00Z
window end: 2021-03-25T03:00:00Z
variable: waterTemperature
#variable: salinity
error ratio: 0.4
output file: insitu_profile_argo_waterTemperature.ioda.nc
#output file: insitu_profile_argo_salinity.ioda.nc
input files:
- gdas.t00z.insitu_profile_argo.2021032400.nc4
- gdas.t06z.insitu_profile_argo.2021032406.nc4
- gdas.t12z.insitu_profile_argo.2021032412.nc4
- gdas.t18z.insitu_profile_argo.2021032418.nc4
- gdas.t00z.insitu_profile_argo.2021032500.nc4
- gdas.t06z.insitu_profile_argo.2021032506.nc4
- gdas.t12z.insitu_profile_argo.2021032512.nc4
- gdas.t18z.insitu_profile_argo.2021032518.nc4
- gdas.t00z.insitu_profile_argo.2021032600.nc4