Skip to content

Commit

Permalink
IODA Converter with In-Situ Observations Concatenation and ObsError I…
Browse files Browse the repository at this point in the history
…nflation (#1472)

Concatenation and re-dating of ioda files.
  • Loading branch information
apchoiCMD authored Feb 13, 2025
1 parent 7b90ef7 commit d5083d5
Show file tree
Hide file tree
Showing 10 changed files with 507 additions and 4 deletions.
162 changes: 162 additions & 0 deletions utils/obsproc/InsituAll2ioda.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,162 @@
#pragma once

#include <cstdlib>
#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);

// TODO(Mindo): This is incomplete and needed to leverage ioda for the reading
// Check if the MetaData group is null
netCDF::NcGroup metaDataGroup = ncFile.getGroup("MetaData");
if (metaDataGroup.isNull()) {
oops::Log::error() << "Group 'MetaData' not found! Aborting execution..." << std::endl;
std::abort();
}

// 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::warning()
<< "WARNING: no depth found, assuming the observations are at the surface."
<< 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());

// Define and check obs groups
struct { const char* name; netCDF::NcGroup group; } groups[] = {
{"ObsValue", ncFile.getGroup("ObsValue")},
{"ObsError", ncFile.getGroup("ObsError")},
{"PreQC", ncFile.getGroup("PreQC")}
};

// Validate groups and abort if any is missing
for (const auto& g : groups) {
if (g.group.isNull()) {
oops::Log::error() << "Group '" << g.name
<< "' not found! Aborting execution..." << std::endl;
std::abort();
}
}

// Assign validated groups
netCDF::NcGroup& obsvalGroup = groups[0].group;
netCDF::NcGroup& obserrGroup = groups[1].group;
netCDF::NcGroup& preqcGroup = groups[2].group;

// 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 EpochTime String Format(1970-01-01T00:00:00Z)
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

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

// Extract EpochTime String Format(1858-11-17T00:00:00Z)
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
12 changes: 11 additions & 1 deletion utils/test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ list( APPEND utils_test_input
testinput/gdas_icecamsr2ioda.yaml
testinput/gdas_icecmirs2ioda.yaml
testinput/gdas_icecjpssrr2ioda.yaml
testinput/gdas_insituall2ioda.yaml
testinput/gdas_viirsaod2ioda.yaml
)

Expand All @@ -25,6 +26,7 @@ set( gdas_utils_test_ref
testref/icecamsr2ioda.test
testref/icecmirs2ioda.test
testref/icecjpssrr2ioda.test
testref/insituall2ioda.test
testref/viirsaod2ioda.test
)

Expand Down Expand Up @@ -174,10 +176,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)
2 changes: 2 additions & 0 deletions utils/test/prepdata.sh
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,8 @@ cdl2nc4 sss_smap_1.nc4 ${project_source_dir}/testdata/sss_smap_1.cdl
cdl2nc4 sss_smap_2.nc4 ${project_source_dir}/testdata/sss_smap_2.cdl
cdl2nc4 sss_smos_1.nc4 ${project_source_dir}/testdata/sss_smos_1.cdl
cdl2nc4 sss_smos_2.nc4 ${project_source_dir}/testdata/sss_smos_2.cdl
cdl2nc4 insitu_profile_argo_1.nc4 ${project_source_dir}/testdata/insitu_profile_argo_1.cdl
cdl2nc4 insitu_profile_argo_2.nc4 ${project_source_dir}/testdata/insitu_profile_argo_2.cdl
cdl2nc4 ghrsst_sst_ma_202103241540.nc4 ${project_source_dir}/testdata/ghrsst_sst_ma_202103241540.cdl
cdl2nc4 ghrsst_sst_ma_202103241550.nc4 ${project_source_dir}/testdata/ghrsst_sst_ma_202103241550.cdl
cdl2nc4 viirs_aod_1.nc4 ${project_source_dir}/testdata/viirs_aod_1.cdl
Expand Down
140 changes: 140 additions & 0 deletions utils/test/testdata/insitu_profile_argo_1.cdl
Original file line number Diff line number Diff line change
@@ -0,0 +1,140 @@
netcdf insitu_profile_argo_1 {
dimensions:
Location = UNLIMITED ; // (16 currently)
variables:
int64 Location(Location) ;
Location:suggested_chunk_dim = 10000LL ;

// global attributes:
string :datetimeRange = "1616598060", "1616622120" ;
string :dataProviderOrigin = "U.S. NOAA" ;
:_ioda_layout_version = 0 ;
string :sourceFiles = "2021032418-gdas.t18z.subpfl.tm00.bufr_d" ;
string :source = "NCEP data tank" ;
string :Converter = "BUFR to IODA Converter" ;
string :_ioda_layout = "ObsGroup" ;
string :platformLongDescription = "ARGO profiles from subpfl: temperature and salinity" ;
string :description = "6-hrly in situ ARGO profiles" ;
:history = "Thu Jan 30 19:31:00 2025: ncks -v MetaData/dateTime,MetaData/latitude,MetaData/longitude,MetaData/depth,MetaData/oceanBasin,MetaData/rcptdateTime,MetaData/sequenceNumber,MetaData/stationID,ObsValue/waterTemperature,ObsValue/salinity,ObsError/waterTemperature,ObsError/salinity,PreQC/waterTemperature,PreQC/salinity -d Location,3000,3015 gdas.t18z.insitu_profile_argo.2021032418.nc4 insitu_profile_argo_1.nc" ;
:NCO = "netCDF Operators version 5.0.6 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)" ;
data:

group: MetaData {
variables:
int64 dateTime(Location) ;
dateTime:_FillValue = 0LL ;
string dateTime:units = "seconds since 1970-01-01T00:00:00Z" ;
string dateTime:long_name = "Datetime" ;
float depth(Location) ;
depth:_FillValue = 2.147484e+09f ;
string depth:units = "m" ;
string depth:long_name = "Water depth" ;
float latitude(Location) ;
latitude:_FillValue = 3.402823e+38f ;
string latitude:units = "degrees_north" ;
latitude:valid_range = -90.f, 90.f ;
string latitude:long_name = "Latitude" ;
float longitude(Location) ;
longitude:_FillValue = 3.402823e+38f ;
string longitude:units = "degrees_east" ;
longitude:valid_range = -180.f, 180.f ;
string longitude:long_name = "Longitude" ;
int oceanBasin(Location) ;
oceanBasin:_FillValue = 999999 ;
string oceanBasin:long_name = "Ocean basin" ;
int64 rcptdateTime(Location) ;
rcptdateTime:_FillValue = 0LL ;
string rcptdateTime:units = "seconds since 1970-01-01T00:00:00Z" ;
string rcptdateTime:long_name = "receipt Datetime" ;
int sequenceNumber(Location) ;
sequenceNumber:_FillValue = 999999 ;
string sequenceNumber:long_name = "Sequence Number" ;
int stationID(Location) ;
stationID:_FillValue = 2147483647 ;
string stationID:long_name = "Station Identification" ;
data:

dateTime = 1616598720, 1616598720, 1616598720, 1616598720, 1616598720,
1616598720, 1616598720, 1616598720, 1616598720, 1616598720, 1616598720,
1616598720, 1616598720, 1616598720, 1616598720, 1616598720 ;

depth = 3, 4, 4.9, 6, 7, 8, 9, 10, 10.8, 12, 14, 16, 18.1, 20, 21.9, 24 ;

latitude = -8.42897, -8.42897, -8.42897, -8.42897, -8.42897, -8.42897,
-8.42897, -8.42897, -8.42897, -8.42897, -8.42897, -8.42897, -8.42897,
-8.42897, -8.42897, -8.42897 ;

longitude = -134.1183, -134.1183, -134.1183, -134.1183, -134.1183,
-134.1183, -134.1183, -134.1183, -134.1183, -134.1183, -134.1183,
-134.1183, -134.1183, -134.1183, -134.1183, -134.1183 ;

oceanBasin = 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 ;

rcptdateTime = 1616612820, 1616612820, 1616612820, 1616612820, 1616612820,
1616612820, 1616612820, 1616612820, 1616612820, 1616612820, 1616612820,
1616612820, 1616612820, 1616612820, 1616612820, 1616612820 ;

sequenceNumber = 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
13, 13 ;

stationID = 5905699, 5905699, 5905699, 5905699, 5905699, 5905699, 5905699,
5905699, 5905699, 5905699, 5905699, 5905699, 5905699, 5905699, 5905699,
5905699 ;
} // group MetaData

group: ObsError {
variables:
float salinity(Location) ;
salinity:_FillValue = 1.e+20f ;
string salinity:units = "psu" ;
string salinity:long_name = "ObsError" ;
float waterTemperature(Location) ;
waterTemperature:_FillValue = 1.e+20f ;
string waterTemperature:units = "degC" ;
string waterTemperature:long_name = "ObsError" ;
data:

salinity = 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01,
0.01, 0.01, 0.01, 0.01, 0.01, 0.01 ;

waterTemperature = 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02,
0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02 ;
} // group ObsError

group: ObsValue {
variables:
float salinity(Location) ;
salinity:_FillValue = 3.402823e+38f ;
string salinity:units = "psu" ;
salinity:valid_range = 0.f, 45.f ;
string salinity:long_name = "salinity" ;
float waterTemperature(Location) ;
waterTemperature:_FillValue = 3.402823e+38f ;
string waterTemperature:units = "degC" ;
waterTemperature:valid_range = -10.f, 50.f ;
string waterTemperature:long_name = "waterTemperature" ;
data:

salinity = 35.283, 35.284, 35.284, 35.284, 35.285, 35.283, 35.285, 35.283,
35.285, 35.283, 35.284, 35.285, 35.284, 35.284, 35.284, 35.284 ;

waterTemperature = 27.70999, 27.71401, 27.71401, 27.712, 27.71099,
27.71099, 27.712, 27.717, 27.71301, 27.716, 27.717, 27.716, 27.71801,
27.71801, 27.71801, 27.71801 ;
} // group ObsValue

group: PreQC {
variables:
int salinity(Location) ;
salinity:_FillValue = 999999 ;
string salinity:long_name = "PreQC" ;
int waterTemperature(Location) ;
waterTemperature:_FillValue = 999999 ;
string waterTemperature:long_name = "PreQC" ;
data:

salinity = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;

waterTemperature = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
} // group PreQC
}
Loading

0 comments on commit d5083d5

Please sign in to comment.