diff --git a/utils/CMakeLists.txt b/utils/CMakeLists.txt index beda28531..49290c0dc 100644 --- a/utils/CMakeLists.txt +++ b/utils/CMakeLists.txt @@ -18,3 +18,4 @@ find_package(soca REQUIRED) add_subdirectory(soca) add_subdirectory(ioda_example) add_subdirectory(test) +add_subdirectory(obsproc) diff --git a/utils/obsproc/CMakeLists.txt b/utils/obsproc/CMakeLists.txt new file mode 100644 index 000000000..4e4e3feb6 --- /dev/null +++ b/utils/obsproc/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(applications) diff --git a/utils/obsproc/NetCDFToIodaConverter.h b/utils/obsproc/NetCDFToIodaConverter.h new file mode 100644 index 000000000..face50faa --- /dev/null +++ b/utils/obsproc/NetCDFToIodaConverter.h @@ -0,0 +1,59 @@ +#pragma once + +#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" + +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 inputFilenames_; + std::string outputFilename_; + }; +} // namespace gdasapp diff --git a/utils/obsproc/Rads2Ioda.h b/utils/obsproc/Rads2Ioda.h new file mode 100644 index 000000000..6c1d683c3 --- /dev/null +++ b/utils/obsproc/Rads2Ioda.h @@ -0,0 +1,87 @@ +#pragma once + +#include +#include // NOLINT (using C API) +#include + +#include "eckit/config/LocalConfiguration.h" + +#include // 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("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(missing_value); + + // Create the IODA variables + ioda::Variable adtIodaDatetime = + ogrp.vars.createWithScales("Metadata/dateTime", + {ogrp.vars["Location"]}, float_params); + // TODO(Mindo): Get the date info from the netcdf file + adtIodaDatetime.atts.add("units", {"seconds since 9999-04-15T12:00:00Z"}, {1}); + + ioda::Variable adtIodaObsVal = + ogrp.vars.createWithScales("ObsValue/absoluteDynamicTopography", + {ogrp.vars["Location"]}, float_params); + ioda::Variable adtIodaObsErr = + ogrp.vars.createWithScales("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(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 diff --git a/utils/obsproc/applications/CMakeLists.txt b/utils/obsproc/applications/CMakeLists.txt new file mode 100644 index 000000000..f55107b6e --- /dev/null +++ b/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) diff --git a/utils/obsproc/applications/gdas_rads2ioda.cc b/utils/obsproc/applications/gdas_rads2ioda.cc new file mode 100644 index 000000000..39f8d2fd4 --- /dev/null +++ b/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); +} diff --git a/utils/obsproc/applications/gdas_rads2ioda.h b/utils/obsproc/applications/gdas_rads2ioda.h new file mode 100644 index 000000000..6e37f5f3d --- /dev/null +++ b/utils/obsproc/applications/gdas_rads2ioda.h @@ -0,0 +1,30 @@ +#pragma once + +#include + +#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 diff --git a/utils/test/CMakeLists.txt b/utils/test/CMakeLists.txt index 23a126909..d0de5c6de 100644 --- a/utils/test/CMakeLists.txt +++ b/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) diff --git a/utils/test/prepdata.sh b/utils/test/prepdata.sh new file mode 100755 index 000000000..e93f672ae --- /dev/null +++ b/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 diff --git a/utils/test/testdata/rads_adt_j3_2021182.cdl b/utils/test/testdata/rads_adt_j3_2021182.cdl new file mode 100644 index 000000000..742044652 --- /dev/null +++ b/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 ; +} diff --git a/utils/test/testinput/gdas_rads2ioda.yaml b/utils/test/testinput/gdas_rads2ioda.yaml new file mode 100644 index 000000000..bf64f8366 --- /dev/null +++ b/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