IODA Converter with In-Situ Observations Concatenation and ObsError Inflation

utils/obsproc/InsituAll2ioda.h

@@ -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
+

utils/obsproc/Rads2Ioda.h

@@ -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;

utils/obsproc/applications/gdas_obsprovider2ioda.h

@@ -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"
@@ -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;

utils/obsproc/util.h

@@ -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);

utils/test/CMakeLists.txt

@@ -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
 )
 
@@ -25,6 +26,7 @@ set( gdas_utils_test_ref
   testref/icecamsr2ioda.test
   testref/icecmirs2ioda.test
   testref/icecjpssrr2ioda.test
+  testref/insituall2ioda.test
   testref/viirsaod2ioda.test
 )
 
@@ -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)

utils/test/prepdata.sh

@@ -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

utils/test/testdata/insitu_profile_argo_1.cdl

@@ -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
+}