Unfinished example of a RADS to ioda converter

utils/CMakeLists.txt

@@ -18,3 +18,4 @@ find_package(soca REQUIRED)
 add_subdirectory(soca)
 add_subdirectory(ioda_example)
 add_subdirectory(test)
+add_subdirectory(obsproc)

utils/obsproc/CMakeLists.txt

@@ -0,0 +1 @@
+add_subdirectory(applications)

utils/obsproc/NetCDFToIodaConverter.h

@@ -0,0 +1,59 @@
+#pragma once
+
+#include <iostream>
+#include <string>
+#include <vector>
+
+#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"
+
+namespace gdasapp {
+  class NetCDFToIodaConverter {
+   public:
+    // Constructor that stores the configuration as a data member
+  explicit NetCDFToIodaConverter(const eckit::Configuration & fullConfig) {
+    // time window info
+    std::string winbegin;
+    std::string winend;
+    fullConfig.get("window begin", winbegin);
+    fullConfig.get("window end", winend);
+    this->windowBegin_ = util::DateTime(winbegin);
+    this->windowEnd_ = util::DateTime(winend);
+
+    // get input netcdf files
+    fullConfig.get("input files", this->inputFilenames_);
+
+    // ioda output file name
+    this->outputFilename_ = fullConfig.getString("output file");
+  }
+
+    // Method to write out a IODA file
+    void WriteToIoda() {
+      // Create empty group backed by HDF file
+      ioda::Group group =
+        ioda::Engines::HH::createFile(
+                                      this->outputFilename_,
+                                      ioda::Engines::BackendCreateModes::Truncate_If_Exists);
+      this->ReadFromNetCDF(group);
+    }
+
+   private:
+    // Virtual method to read a netcdf file and store the relevant info in a group
+    virtual void ReadFromNetCDF(ioda::Group &) = 0;
+
+
+   protected:
+    util::DateTime windowBegin_;
+    util::DateTime windowEnd_;
+    std::vector<std::string> inputFilenames_;
+    std::string outputFilename_;
+  };
+}  // namespace gdasapp

utils/obsproc/Rads2Ioda.h

@@ -0,0 +1,87 @@
+#pragma once
+
+#include <iostream>
+#include <netcdf>    // NOLINT (using C API)
+#include <string>
+
+#include "eckit/config/LocalConfiguration.h"
+
+#include <Eigen/Dense>    // NOLINT
+
+#include "ioda/Group.h"
+#include "ioda/ObsGroup.h"
+
+#include "NetCDFToIodaConverter.h"
+
+#include "oops/util/missingValues.h"
+
+namespace gdasapp {
+  class Rads2Ioda : public NetCDFToIodaConverter {
+   public:
+    explicit Rads2Ioda(const eckit::Configuration & fullConfig)
+    : NetCDFToIodaConverter(fullConfig) {
+    oops::Log::info() << "Window begin: " << this->windowBegin_ << std::endl;
+    oops::Log::info() << "Window end: " << this->windowEnd_ << std::endl;
+    oops::Log::info() << "IODA output file name: " << this->outputFilename_ << std::endl;
+    }
+
+    // Read netcdf file and populate group
+    void ReadFromNetCDF(ioda::Group & group) override {
+      std::string fileName = this->inputFilenames_[0];  // TODO(Guillaume): make it work for a list
+
+      // Open the NetCDF file in read-only mode
+      netCDF::NcFile ncFile(fileName, netCDF::NcFile::read);
+
+      // Get dimensions
+      int location = ncFile.getDim("time").getSize();
+      int nVars = 1;
+
+      // Get adt_egm2008 obs values and attributes
+      netCDF::NcVar adtNcVar = ncFile.getVar("adt_egm2008");
+      int adtObsVal[location];  // NOLINT (can't pass vector to getVar below)
+      adtNcVar.getVar(adtObsVal);
+      std::string units;
+      adtNcVar.getAtt("units").getValues(units);
+      float scaleFactor;
+      adtNcVar.getAtt("scale_factor").getValues(&scaleFactor);
+
+      // Update the group with 1 dimension: location
+      ioda::NewDimensionScales_t newDims {ioda::NewDimensionScale<int>("Location", location)};
+      ioda::ObsGroup ogrp = ioda::ObsGroup::generate(group, newDims);
+
+      // Set up the float creation parameters
+      ioda::VariableCreationParameters float_params;
+      float_params.chunk = true;               // allow chunking
+      float_params.compressWithGZIP();         // compress using gzip
+      float missing_value = util::missingValue(missing_value);
+      float_params.setFillValue<float>(missing_value);
+
+      // Create the IODA variables
+      ioda::Variable adtIodaDatetime =
+        ogrp.vars.createWithScales<float>("Metadata/dateTime",
+                                          {ogrp.vars["Location"]}, float_params);
+      // TODO(Mindo): Get the date info from the netcdf file
+      adtIodaDatetime.atts.add<std::string>("units", {"seconds since 9999-04-15T12:00:00Z"}, {1});
+
+      ioda::Variable adtIodaObsVal =
+        ogrp.vars.createWithScales<float>("ObsValue/absoluteDynamicTopography",
+                                          {ogrp.vars["Location"]}, float_params);
+      ioda::Variable adtIodaObsErr =
+        ogrp.vars.createWithScales<float>("ObsError/absoluteDynamicTopography",
+                                          {ogrp.vars["Location"]}, float_params);
+
+      // Write adt obs value to group
+      Eigen::ArrayXf tmpvar(location);
+      for (int i = 0; i <= location; i++) {
+        tmpvar[i] = static_cast<float>(adtObsVal[i])*scaleFactor;
+      }
+      adtIodaObsVal.writeWithEigenRegular(tmpvar);
+
+      // Write adt obs error to group
+      for (int i = 0; i <= location; i++) {
+        tmpvar[i] = 0.1;
+      }
+      adtIodaObsErr.writeWithEigenRegular(tmpvar);
+    };
+  };
+}  // namespace gdasapp

utils/obsproc/applications/CMakeLists.txt

@@ -0,0 +1,9 @@
+list( APPEND gdasapp_rads2ioda_src_files
+  gdas_rads2ioda.cc
+  gdas_rads2ioda.h
+  )
+ecbuild_add_executable( TARGET gdas_rads2ioda.x
+                        SOURCES ${gdasapp_rads2ioda_src_files} )
+
+target_compile_features( gdas_rads2ioda.x PUBLIC cxx_std_17)
+target_link_libraries( gdas_rads2ioda.x PUBLIC oops ioda NetCDF::NetCDF_CXX)

utils/obsproc/applications/gdas_rads2ioda.cc

@@ -0,0 +1,8 @@
+#include "gdas_rads2ioda.h"
+#include "oops/runs/Run.h"
+
+int main(int argc, char ** argv) {
+  oops::Run run(argc, argv);
+  gdasapp::Rads2IodaApp rads2iodaapp;
+  return run.execute(rads2iodaapp);
+}

utils/obsproc/applications/gdas_rads2ioda.h

@@ -0,0 +1,30 @@
+#pragma once
+
+#include <string>
+
+#include "eckit/config/LocalConfiguration.h"
+#include "oops/mpi/mpi.h"
+#include "oops/runs/Application.h"
+
+#include "../Rads2Ioda.h"
+
+namespace gdasapp {
+  class Rads2IodaApp : public oops::Application {
+   public:
+    explicit Rads2IodaApp(const eckit::mpi::Comm & comm = oops::mpi::world())
+      : Application(comm) {}
+    static const std::string classname() {return "gdasapp::Rads2IodaApp";}
+
+    int execute(const eckit::Configuration & fullConfig, bool /*validate*/) const {
+      Rads2Ioda rads2ioda(fullConfig);
+      rads2ioda.WriteToIoda();
+      return 0;
+    }
+    // -----------------------------------------------------------------------------
+   private:
+    std::string appname() const {
+      return "gdasapp::Rads2IodaApp";
+    }
+    // -----------------------------------------------------------------------------
+  };
+}  // namespace gdasapp

utils/test/CMakeLists.txt

@@ -1,6 +1,7 @@
 # Create Data directory for test input config and symlink all files
 list( APPEND utils_test_input
   testinput/gdas_meanioda.yaml
+  testinput/gdas_rads2ioda.yaml
 )
 
 file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/testinput)
@@ -22,3 +23,17 @@ ecbuild_add_test( TARGET  test_gdasapp_util_ioda_example
                   COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_meanioda.x
                   ARGS    "testinput/gdas_meanioda.yaml"
                   LIBS    gdas-utils)
+
+# Prepare data for the IODA converters
+file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc)
+ecbuild_add_test( TARGET  test_gdasapp_util_prepdata
+                  COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/../../utils/test/prepdata.sh
+                  ARGS    ${CMAKE_CURRENT_SOURCE_DIR}/../../utils/test
+                  WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc)
+
+# Test the RADS to IODA converter
+ecbuild_add_test( TARGET  test_gdasapp_util_rads2ioda
+                  COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_rads2ioda.x
+                  ARGS    "../testinput/gdas_rads2ioda.yaml"
+                  LIBS    gdas-utils
+                  WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc)

utils/test/prepdata.sh

@@ -0,0 +1,6 @@
+#!/bin/bash
+set -e
+
+project_source_dir=$1
+
+ncgen -o rads_adt_j3_2021182.nc4 ${project_source_dir}/testdata/rads_adt_j3_2021182.cdl

utils/test/testdata/rads_adt_j3_2021182.cdl

@@ -0,0 +1,110 @@
+netcdf output {
+dimensions:
+	time = UNLIMITED ; // (11 currently)
+variables:
+	int adt_egm2008(time) ;
+		adt_egm2008:_FillValue = 2147483647 ;
+		adt_egm2008:long_name = "absolute dynamic topography (EGM2008)" ;
+		adt_egm2008:standard_name = "absolute_dynamic_topography_egm2008" ;
+		adt_egm2008:units = "m" ;
+		adt_egm2008:scale_factor = 0.0001 ;
+		adt_egm2008:coordinates = "lon lat" ;
+	int adt_xgm2016(time) ;
+		adt_xgm2016:_FillValue = 2147483647 ;
+		adt_xgm2016:long_name = "absolute dynamic topography (XGM2016)" ;
+		adt_xgm2016:standard_name = "absolute_dynamic_topography_xgm2016" ;
+		adt_xgm2016:units = "m" ;
+		adt_xgm2016:scale_factor = 0.0001 ;
+		adt_xgm2016:coordinates = "lon lat" ;
+	int cycle(time) ;
+		cycle:_FillValue = 2147483647 ;
+		cycle:long_name = "cycle number" ;
+		cycle:field = 9905s ;
+	int lat(time) ;
+		lat:_FillValue = 2147483647 ;
+		lat:long_name = "latitude" ;
+		lat:standard_name = "latitude" ;
+		lat:units = "degrees_north" ;
+		lat:scale_factor = 1.e-06 ;
+		lat:field = 201s ;
+		lat:comment = "Positive latitude is North latitude, negative latitude is South latitude" ;
+	int lon(time) ;
+		lon:_FillValue = 2147483647 ;
+		lon:long_name = "longitude" ;
+		lon:standard_name = "longitude" ;
+		lon:units = "degrees_east" ;
+		lon:scale_factor = 1.e-06 ;
+		lon:field = 301s ;
+		lon:comment = "East longitude relative to Greenwich meridian" ;
+	int pass(time) ;
+		pass:_FillValue = 2147483647 ;
+		pass:long_name = "pass number" ;
+		pass:field = 9906s ;
+	short sla(time) ;
+		sla:_FillValue = 32767s ;
+		sla:long_name = "sea level anomaly" ;
+		sla:standard_name = "sea_surface_height_above_sea_level" ;
+		sla:units = "m" ;
+		sla:quality_flag = "swh sig0 range_rms range_numval flags swh_rms sig0_rms attitude" ;
+		sla:scale_factor = 0.0001 ;
+		sla:coordinates = "lon lat" ;
+		sla:field = 0s ;
+		sla:comment = "Sea level determined from satellite altitude - range - all altimetric corrections" ;
+	double time_dtg(time) ;
+		time_dtg:long_name = "time_dtg" ;
+		time_dtg:standard_name = "time_dtg" ;
+		time_dtg:units = "yyyymmddhhmmss" ;
+		time_dtg:coordinates = "lon lat" ;
+		time_dtg:comment = "UTC time formatted as yyyymmddhhmmss" ;
+	double time_mjd(time) ;
+		time_mjd:long_name = "Modified Julian Days" ;
+		time_mjd:standard_name = "time" ;
+		time_mjd:units = "days since 1858-11-17 00:00:00 UTC" ;
+		time_mjd:field = 105s ;
+		time_mjd:comment = "UTC time of measurement expressed in Modified Julian Days" ;
+
+// global attributes:
+		:Conventions = "CF-1.7" ;
+		:title = "RADS 4 pass file" ;
+		:institution = "EUMETSAT / NOAA / TU Delft" ;
+		:source = "radar altimeter" ;
+		:references = "RADS Data Manual, Version 4.2 or later" ;
+		:featureType = "trajectory" ;
+		:ellipsoid = "TOPEX" ;
+		:ellipsoid_axis = 6378136.3 ;
+		:ellipsoid_flattening = 0.00335281317789691 ;
+		:filename = "rads_adt_j3_2021182.nc" ;
+		:mission_name = "JASON-3" ;
+		:mission_phase = "a" ;
+		:log01 = "2021-07-03 | /Users/rads/bin/rads2nc --ymd=20210701000000,20210702000000 -C1,1000 -Sj3 -Vsla,adt_egm2008,adt_xgm2016,time_mjd,time_dtg,lon,lat,cycle,pass -X/Users/rads/cron/xgm2016 -X/Users/rads/cron/adt -X/Users/rads/cron/time_dtg -o/ftp/rads/adt//2021/rads_adt_j3_2021182.nc: RAW data from" ;
+		:history = "Thu Sep  7 14:43:07 2023: ncks -d time,0,10 rads_adt_j3_2021182.nc output.nc\n",
+			"2021-07-03 20:56:20 : /Users/rads/bin/rads2nc --ymd=20210701000000,20210702000000 -C1,1000 -Sj3 -Vsla,adt_egm2008,adt_xgm2016,time_mjd,time_dtg,lon,lat,cycle,pass -X/Users/rads/cron/xgm2016 -X/Users/rads/cron/adt -X/Users/rads/cron/time_dtg -o/ftp/rads/adt//2021/rads_adt_j3_2021182.nc" ;
+		:NCO = "netCDF Operators version 5.0.6 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)" ;
+data:
+
+ adt_egm2008 = -7884, -10580, -7180, -8899, -9341, -8404, -8400, -9468, 
+    -8810, -10000, -8592 ;
+
+ adt_xgm2016 = -8097, -10657, -7368, -9127, -9540, -8536, -8406, -9285, 
+    -8232, -9248, -7758 ;
+
+ cycle = 198, 198, 198, 198, 198, 198, 198, 198, 198, 198, 198 ;
+
+ lat = -65447761, -65484640, -65496721, -65508694, -65520561, -65532322, 
+    -65543975, -65555521, -65578291, -65589514, -65600630 ;
+
+ lon = -84965810, -84597998, -84475158, -84352200, -84229127, -84105938, 
+    -83982635, -83859219, -83612050, -83488298, -83364437 ;
+
+ pass = 184, 184, 184, 184, 184, 184, 184, 184, 184, 184, 184 ;
+
+ sla = 141, -2275, 1208, -498, -849, 25, 91, -893, -315, -1329, 126 ;
+
+ time_dtg = 20210701000000, 20210701000003, 20210701000004, 20210701000005, 
+    20210701000006, 20210701000007, 20210701000008, 20210701000009, 
+    20210701000011, 20210701000012, 20210701000013 ;
+
+ time_mjd = 59396.0000038319, 59396.0000392038, 59396.0000509944, 
+    59396.0000627851, 59396.0000745757, 59396.0000863663, 59396.0000981569, 
+    59396.0001099476, 59396.0001335288, 59396.0001453194, 59396.0001571101 ;
+}

utils/test/testinput/gdas_rads2ioda.yaml

@@ -0,0 +1,5 @@
+window begin: 2018-04-15T06:00:00Z
+window end: 2018-04-15T12:00:00Z
+output file: rads_adt_j3_2021182.ioda.nc
+input files:
+- rads_adt_j3_2021182.nc4