From 4335e71d9c0282018879c4fc48c98be48d7b2623 Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Tue, 23 Apr 2024 13:16:28 -0500 Subject: [PATCH 01/18] added an empty test rtofsinsitu --- utils/obsproc/RTOFSInSitu.h | 227 ++++++++++++++++++ .../applications/gdas_obsprovider2ioda.h | 4 + utils/test/CMakeLists.txt | 9 + utils/test/testinput/gdas_rtofsinsitu.yaml | 12 + 4 files changed, 252 insertions(+) create mode 100644 utils/obsproc/RTOFSInSitu.h create mode 100644 utils/test/testinput/gdas_rtofsinsitu.yaml diff --git a/utils/obsproc/RTOFSInSitu.h b/utils/obsproc/RTOFSInSitu.h new file mode 100644 index 000000000..5c8aca9d8 --- /dev/null +++ b/utils/obsproc/RTOFSInSitu.h @@ -0,0 +1,227 @@ +#pragma once + +#include +#include // NOLINT (using C API) +#include +#include +#include + +#include "eckit/config/LocalConfiguration.h" + +#include // NOLINT + +#include "ioda/Group.h" +#include "ioda/ObsGroup.h" + +#include "NetCDFToIodaConverter.h" +#include "superob.h" + + +namespace gdasapp { + + class RTOFSInSitu : public NetCDFToIodaConverter { + public: + explicit RTOFSInSitu(const eckit::Configuration & fullConfig, const eckit::mpi::Comm & comm) + : NetCDFToIodaConverter(fullConfig, comm) { + variable_ = "seaSurfaceTemperature"; + } + + // Read netcdf file and populate iodaVars + gdasapp::obsproc::iodavars::IodaVars providerToIodaVars(const std::string fileName) final { + oops::Log::info() << "Processing files provided by RTOFS" << std::endl; + + // Get the sst bounds from the configuration + float sstMin; + fullConfig_.get("bounds.min", sstMin); + float sstMax; + fullConfig_.get("bounds.max", sstMax); + + /* + // Open the NetCDF file in read-only mode + netCDF::NcFile ncFile(fileName, netCDF::NcFile::read); + oops::Log::info() << "Reading... " << fileName << std::endl; + // Get number of obs + int dimLon = ncFile.getDim("lon").getSize(); + int dimLat = ncFile.getDim("lat").getSize(); + int dimTime = ncFile.getDim("time").getSize(); + + // Read non-optional metadata: datetime, longitude and latitude + // latitude + std::vector lat(dimLat); + ncFile.getVar("lat").getVar(lat.data()); + + // longitude + std::vector lon(dimLon); + ncFile.getVar("lon").getVar(lon.data()); + + // Generate the lat-lon grid + std::vector> lon2d(dimLat, std::vector(dimLon)); + std::vector> lat2d(dimLat, std::vector(dimLon)); + for (int i = 0; i < dimLat; ++i) { + for (int j = 0; j < dimLon; ++j) { + lon2d[i][j] = lon[j]; + lat2d[i][j] = lat[i]; + } + } + + // Read the reference time + std::vector refTime(dimTime); + ncFile.getVar("time").getVar(refTime.data()); + std::string refDate; + ncFile.getVar("time").getAtt("units").getValues(refDate); + + // Reformat the reference time + std::regex dateRegex(R"(\b\d{4}-\d{2}-\d{2} \d{2}:\d{2}:\d{2}\b)"); + std::smatch match; + std::regex_search(refDate, match, dateRegex); + refDate = match.str(); + std::tm tmStruct = {}; + std::istringstream ss(refDate); + ss >> std::get_time(&tmStruct, "%Y-%m-%d %H:%M:%S"); + std::ostringstream isoFormatted; + isoFormatted << std::put_time(&tmStruct, "seconds since %Y-%m-%dT%H:%M:%SZ"); + refDate = isoFormatted.str(); + + // Read sst_dtime to add to the reference time + // TODO(AMG): What's below does not read the field the same way python does + std::vector sstdTime(dimTime*dimLat*dimLon); + ncFile.getVar("sst_dtime").getVar(sstdTime.data()); + float dtimeScaleFactor; + ncFile.getVar("sst_dtime").getAtt("scale_factor").getValues(&dtimeScaleFactor); + int32_t dtimeFillValue; + ncFile.getVar("sst_dtime").getAtt("_FillValue").getValues(&dtimeFillValue); + + // Read SST ObsValue + std::vector sstObsVal(dimTime*dimLat*dimLon); + ncFile.getVar("sea_surface_temperature").getVar(sstObsVal.data()); + float sstOffSet; + ncFile.getVar("sea_surface_temperature").getAtt("add_offset").getValues(&sstOffSet); + float sstScaleFactor; + ncFile.getVar("sea_surface_temperature").getAtt("scale_factor").getValues(&sstScaleFactor); + + // Read SST Bias + std::vector sstObsBias(dimTime*dimLat*dimLon); + ncFile.getVar("sses_bias").getVar(sstObsBias.data()); + float biasScaleFactor; + ncFile.getVar("sses_bias").getAtt("scale_factor").getValues(&biasScaleFactor); + + // Read Error + std::vector sstObsErr(dimTime*dimLat*dimLon); + ncFile.getVar("sses_standard_deviation").getVar(sstObsErr.data()); + float errOffSet; + ncFile.getVar("sses_standard_deviation").getAtt("add_offset").getValues(&errOffSet); + float errScaleFactor; + ncFile.getVar("sses_standard_deviation").getAtt("scale_factor").getValues(&errScaleFactor); + + // Read preQc + signed char sstPreQC[dimTime][dimLat][dimLon]; + ncFile.getVar("quality_level").getVar(sstPreQC); + + // Apply scaling/unit change and compute the necessary fields + std::vector> mask(dimLat, std::vector(dimLon)); + std::vector> sst(dimLat, std::vector(dimLon)); + std::vector> obserror(dimLat, std::vector(dimLon)); + std::vector> preqc(dimLat, std::vector(dimLon)); + std::vector> seconds(dimLat, std::vector(dimLon)); + size_t index = 0; + for (int i = 0; i < dimLat; i++) { + for (int j = 0; j < dimLon; j++) { + // provider's QC flag + // Note: the qc flags in GDS2.0 run from 0 to 5, with higher numbers being better. + // IODA typically expects 0 to be good, and larger numbers to be worse so the + // provider's QC is flipped + preqc[i][j] = 5 - static_cast(sstPreQC[0][i][j]); + + // bias corrected sst, regressed to the drifter depth + // Remove added sstOffSet for Celsius + sst[i][j] = static_cast(sstObsVal[index]) * sstScaleFactor + - static_cast(sstObsBias[index]) * biasScaleFactor; + // mask + if (sst[i][j] >= sstMin && sst[i][j] <= sstMax && preqc[i][j] ==0) { + mask[i][j] = 1; + } else { + mask[i][j] = 0; + } + + // obs error + // TODO(Somebody): add sampled std. dev. of sst to the total obs error + obserror[i][j] = static_cast(sstObsErr[index]) * errScaleFactor + errOffSet; + + // epoch time in seconds and avoid to use FillValue in calculation + if (sstdTime[index] == dtimeFillValue) { + seconds[i][j] = 0; + } else { + seconds[i][j] = static_cast(sstdTime[index]) * dtimeScaleFactor + + static_cast(refTime[0]); + } + index++; + } + } + + // Superobing + // TODO(Guillaume): Save the sampling std dev of sst so it can be used + // as a proxi for obs error + std::vector> sst_s; + std::vector> lon2d_s; + std::vector> lat2d_s; + std::vector> obserror_s; + std::vector> seconds_s; + if ( fullConfig_.has("binning") ) { + sst_s = gdasapp::superobutils::subsample2D(sst, mask, fullConfig_); + lon2d_s = gdasapp::superobutils::subsample2D(lon2d, mask, fullConfig_); + lat2d_s = gdasapp::superobutils::subsample2D(lat2d, mask, fullConfig_); + obserror_s = gdasapp::superobutils::subsample2D(obserror, mask, fullConfig_); + seconds_s = gdasapp::superobutils::subsample2D(seconds, mask, fullConfig_); + } else { + sst_s = sst; + lon2d_s = lon2d; + lat2d_s = lat2d; + obserror_s = obserror; + seconds_s = seconds; + } +*/ + + // number of obs after subsampling + // int nobs = sst_s.size() * sst_s[0].size(); + int nobs = 0; + + // Set the int metadata names + // std::vector intMetadataNames = {"oceanBasin"}; + std::vector intMetadataNames = {"dummyBasin"}; + + // Set the float metadata name + std::vector floatMetadataNames = {}; + + // Create instance of iodaVars object + gdasapp::obsproc::iodavars::IodaVars iodaVars(nobs, floatMetadataNames, intMetadataNames); + +/* + // Reference time is Jan 01 1981 00:00:00 GMT+0000 + iodaVars.referenceDate_ = refDate; + + // Store into eigen arrays + int loc(0); + for (int i = 0; i < sst_s.size(); i++) { + for (int j = 0; j < sst_s[0].size(); j++) { + iodaVars.longitude_(loc) = lon2d_s[i][j]; + iodaVars.latitude_(loc) = lat2d_s[i][j]; + iodaVars.obsVal_(loc) = sst_s[i][j]; + iodaVars.obsError_(loc) = obserror_s[i][j]; + iodaVars.preQc_(loc) = 0; + iodaVars.datetime_(loc) = seconds_s[i][j]; + // Store optional metadata, set ocean basins to -999 for now + iodaVars.intMetadata_.row(loc) << -999; + loc += 1; + } + } + + // Basic QC + Eigen::Array boundsCheck = + (iodaVars.obsVal_ > sstMin && iodaVars.obsVal_ < sstMax && iodaVars.datetime_ > 0.0); + iodaVars.trim(boundsCheck); + */ + + return iodaVars; + }; + }; // class RTOFSInSitu +} // namespace gdasapp diff --git a/utils/obsproc/applications/gdas_obsprovider2ioda.h b/utils/obsproc/applications/gdas_obsprovider2ioda.h index 2fef42183..c14c9a1a2 100644 --- a/utils/obsproc/applications/gdas_obsprovider2ioda.h +++ b/utils/obsproc/applications/gdas_obsprovider2ioda.h @@ -8,6 +8,7 @@ #include "../Ghrsst2Ioda.h" #include "../IcecAmsr2Ioda.h" +#include "../RTOFSInSitu.h" #include "../Rads2Ioda.h" #include "../Smap2Ioda.h" #include "../Smos2Ioda.h" @@ -31,6 +32,9 @@ namespace gdasapp { } else if (provider == "GHRSST") { Ghrsst2Ioda conv2ioda(fullConfig, this->getComm()); conv2ioda.writeToIoda(); + } else if (provider == "RTOFS") { + RTOFSInSitu conv2ioda(fullConfig, this->getComm()); + conv2ioda.writeToIoda(); } else if (provider == "SMAP") { Smap2Ioda conv2ioda(fullConfig, this->getComm()); conv2ioda.writeToIoda(); diff --git a/utils/test/CMakeLists.txt b/utils/test/CMakeLists.txt index fcbb26e37..4861e8f45 100644 --- a/utils/test/CMakeLists.txt +++ b/utils/test/CMakeLists.txt @@ -3,6 +3,7 @@ list( APPEND utils_test_input testinput/gdas_meanioda.yaml testinput/gdas_rads2ioda.yaml testinput/gdas_ghrsst2ioda.yaml + testinput/gdas_rtofsinsitu.yaml testinput/gdas_smap2ioda.yaml testinput/gdas_smos2ioda.yaml testinput/gdas_icecamsr2ioda.yaml @@ -11,6 +12,7 @@ list( APPEND utils_test_input set( gdas_utils_test_ref testref/ghrsst2ioda.test + # testref/rtofsinsitu.test testref/rads2ioda.test testref/smap2ioda.test testref/smos2ioda.test @@ -63,6 +65,13 @@ ecbuild_add_test( TARGET test_gdasapp_util_ghrsst2ioda LIBS gdas-utils WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc) +# Test the RTOFS in Situ converter +ecbuild_add_test( TARGET test_gdasapp_util_rtofsinsitu + COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x + ARGS "../testinput/gdas_rtofsinsitu.yaml" + LIBS gdas-utils + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc) + # Test the SMAP to IODA converter ecbuild_add_test( TARGET test_gdasapp_util_smap2ioda COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x diff --git a/utils/test/testinput/gdas_rtofsinsitu.yaml b/utils/test/testinput/gdas_rtofsinsitu.yaml new file mode 100644 index 000000000..8b6a129e6 --- /dev/null +++ b/utils/test/testinput/gdas_rtofsinsitu.yaml @@ -0,0 +1,12 @@ +provider: RTOFS +window begin: 2021-03-24T15:00:00Z +window end: 2021-03-24T21:00:00Z +output file: rtofsinsitu_2024032600.profile +input files: +- 2024032600.profile +- 2024032600.sfc + +# test: + # reference filename: testref/rtofsinsitu.test + # test output filename: testoutput/rtofsinsitu.test + # float relative tolerance: 1e-6 From 547f381242d135ac839fda2aff3349d7642817eb Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Wed, 15 May 2024 10:47:43 -0500 Subject: [PATCH 02/18] basic working code --- utils/obsproc/rtofs/CMakeLists.txt | 38 +++++ utils/obsproc/rtofs/README | 11 ++ utils/obsproc/rtofs/RTOFSDataFile.cc | 168 ++++++++++++++++++++++ utils/obsproc/rtofs/RTOFSDataFile.h | 33 +++++ utils/obsproc/rtofs/RTOFSOb.cc | 161 +++++++++++++++++++++ utils/obsproc/rtofs/RTOFSOb.h | 63 +++++++++ utils/obsproc/rtofs/util.cc | 203 +++++++++++++++++++++++++++ utils/obsproc/rtofs/util.h | 24 ++++ 8 files changed, 701 insertions(+) create mode 100644 utils/obsproc/rtofs/CMakeLists.txt create mode 100644 utils/obsproc/rtofs/README create mode 100644 utils/obsproc/rtofs/RTOFSDataFile.cc create mode 100644 utils/obsproc/rtofs/RTOFSDataFile.h create mode 100644 utils/obsproc/rtofs/RTOFSOb.cc create mode 100644 utils/obsproc/rtofs/RTOFSOb.h create mode 100644 utils/obsproc/rtofs/util.cc create mode 100644 utils/obsproc/rtofs/util.h diff --git a/utils/obsproc/rtofs/CMakeLists.txt b/utils/obsproc/rtofs/CMakeLists.txt new file mode 100644 index 000000000..63488dc49 --- /dev/null +++ b/utils/obsproc/rtofs/CMakeLists.txt @@ -0,0 +1,38 @@ +set( + rtofs_src_files + RTOFSDataFile.h + RTOFSDataFile.cc + RTOFSOb.h + RTOFSOb.cc + util.cc + util.h +) + +list( + APPEND gdasapp_provider2ioda_src_files + ${rtofs_src_files} +) + + +add_library(rtofs STATIC + ${rtofs_src_files} +) + +target_compile_features( + rtofs + PUBLIC cxx_std_17 +) + +# target_link_libraries( + # rtofs + # PUBLIC oops ioda NetCDF::NetCDF_CXX +# ) + + +# ecbuild_add_executable( + # TARGET gdas_obsprovider2ioda.x + # SOURCES ${gdasapp_provider2ioda_src_files} +# ) + +# target_compile_features( gdas_obsprovider2ioda.x PUBLIC cxx_std_17) +# target_link_libraries( gdas_obsprovider2ioda.x PUBLIC oops ioda NetCDF::NetCDF_CXX) diff --git a/utils/obsproc/rtofs/README b/utils/obsproc/rtofs/README new file mode 100644 index 000000000..c420242f2 --- /dev/null +++ b/utils/obsproc/rtofs/README @@ -0,0 +1,11 @@ +Author: E. Givelberg +May 15, 2024 + +The code in this rtofs directory reads a binary Fortran +file and generates a data object that can be ingested into ioda. + +The fortran binary file consists of records, including additional +information that needs to be skipped. +There are typically 8 bytes to be skipped between 2 data arrays. +The read data is converted from the big endian to the host +format, which is little endian in linux. diff --git a/utils/obsproc/rtofs/RTOFSDataFile.cc b/utils/obsproc/rtofs/RTOFSDataFile.cc new file mode 100644 index 000000000..e6620e1e6 --- /dev/null +++ b/utils/obsproc/rtofs/RTOFSDataFile.cc @@ -0,0 +1,168 @@ +#include +#include +#include + + +#include +using std::cerr; +using std::endl; +using std::cout; +#include // std::get_time + +#include "RTOFSDataFile.h" +#include "RTOFSOb.h" +#include "util.h" + + +namespace rtofs +{ + +typedef char char12[12]; +typedef char char7[7]; + + +std::time_t + SecondsSinceReferenceTime(char12 time) +{ + std::tm referenceTime = {}; + referenceTime.tm_year = 70; // 1970 + referenceTime.tm_mon = 0; // January (months are 0-based) + referenceTime.tm_mday = 1; // 1st day of the month + referenceTime.tm_hour = 0; + referenceTime.tm_min = 0; + referenceTime.tm_sec = 0; + std::time_t referenceTimestamp = std::mktime(&referenceTime); + + std::tm t = {}; + std::istringstream ss(time); + ss >> std::get_time(&t, "%Y%m%d%H%M"); + + std::time_t timestamp = std::mktime(&t); + return std::difftime(timestamp, referenceTimestamp); + +} // SecondsSinceReferenceTime + + + +RTOFSDataFile:: +RTOFSDataFile(std::string filename): + filename(filename) +{ + + if (! file_exists(filename)) + { + cerr << "File not found" << endl; + exit(1); + } + + + const char * fn = filename.c_str(); + f = fopen(fn, "rb"); + if (!f) + { + cerr << "Error opening file " << fn << endl; + exit(1); + } + + read_file(); + +} // RTOFSDataFile::RTOFSDataFile + + +void +RTOFSDataFile:: + read_file() +{ + fseek(f, 4, SEEK_CUR); + + int n_read = read_int(f); + int mx_lvl = read_int(f); + int vrsn = read_int(f); + + ob = new RTOFSOb(n_read, mx_lvl, vrsn); + nobs = n_read; + + ob->read(f); + + skip8bytes(f); + + char12 * ob_dtg = new char12[n_read]; + fread(ob_dtg, sizeof(char[12]), n_read, f); + for (int i = 0; i < n_read; i ++) + ob->dtg[i] = SecondsSinceReferenceTime(ob_dtg[i]); + +/* + { + for (int i = 0; i < n_read; i ++) + { + for (int j = 0; j < 12; j ++) + cout << ob_dtg[i][j]; + + cout + << " = " + << SecondsSinceReferenceTime(ob_dtg[i]) + << " sec" + << endl + ; + cout << endl; + } + } +*/ + + skip8bytes(f); + + char12 * ob_rct = new char12[n_read]; + fread(ob_rct, sizeof(char[12]), n_read, f); + +/* + { + for (int i = 0; i < n_read; i ++) + { + for (int j = 0; j < 12; j ++) + cout << ob_rct[i][j]; + + cout + << " = " + << SecondsSinceReferenceTime(ob_rct[i]) + << " sec" + << endl + ; + cout << endl; + } + } +*/ + + + skip8bytes(f); + + char7 * ob_sgn = new char7[n_read]; + fread(ob_sgn, sizeof(char[7]), n_read, f); + +/* + for (int i = 0; i < n_read; i ++) + { + // for (int j = 0; j < 12; j ++) + // cout << ob_dtg[i][j]; + // cout << (unsigned int) ob_dtg[i][j] << " "; + // cout << ob_rct[i][j]; + for (int j = 0; j < 7; j ++) + cout << ob_sgn[i][j]; + cout << endl; + } +*/ + + if (vrsn == 2) + { + float ** glb_sal = new float * [n_read]; + for (int i = 0; i < n_read; i ++) + { + int k = ob->lt[i]; + glb_sal[i] = alloc_read_float_array(f, k); + } + } + + +} // read_file + + +} // namespace rtofs diff --git a/utils/obsproc/rtofs/RTOFSDataFile.h b/utils/obsproc/rtofs/RTOFSDataFile.h new file mode 100644 index 000000000..ea68e4852 --- /dev/null +++ b/utils/obsproc/rtofs/RTOFSDataFile.h @@ -0,0 +1,33 @@ +#ifndef RTOFSDataFile__H +#define RTOFSDataFile__H + +#include + + +namespace rtofs +{ + +class RTOFSOb; + +class RTOFSDataFile +{ +public: + RTOFSDataFile(std::string filename); + int NumberOfObservations() { return nobs; } + RTOFSOb & observations() { return * ob; } + + +private: + std::string filename; + int nobs; + FILE * f; + RTOFSOb * ob; + + void read_file(); + +}; // class RTOFSDataFile + +} // namespace rtofs + + +#endif // RTOFSDataFile__H diff --git a/utils/obsproc/rtofs/RTOFSOb.cc b/utils/obsproc/rtofs/RTOFSOb.cc new file mode 100644 index 000000000..3e97fe5cc --- /dev/null +++ b/utils/obsproc/rtofs/RTOFSOb.cc @@ -0,0 +1,161 @@ +#include + +#include +using std::string; + +#include +using std::endl; +using std::cerr; + +#include +using std::ofstream; + + +#include "RTOFSOb.h" +#include "util.h" + + +namespace rtofs +{ + + +RTOFSOb:: + RTOFSOb(int n, int mx_lvl, int vrsn): + n(n), + mx_lvl(mx_lvl), + vrsn(vrsn) +{ + dtg = new std::time_t[n]; + allocate(); +} + + + +void +RTOFSOb:: + allocate() +{ + lat = new float[n]; + lon = new float[n]; + btm = new float[n]; + ls = new int[n]; + lt = new int[n]; + sal_typ = new int[n]; + sqc = new float[n]; + tmp_typ = new int[n]; + tqc = new float[n]; + + lvl = new float * [n]; + sal = new float * [n]; + sal_err = new float * [n]; + sprb = new float * [n]; + tmp = new float * [n]; + tmp_err = new float * [n]; + clm_sal = new float * [n]; // skipped? + tprb = new float * [n]; + cssd = new float * [n]; + clm_tmp = new float * [n]; // skipped? + ctsd = new float * [n]; + flg = new int * [n]; +} + + +void +RTOFSOb:: + allocate2d() +{ + for (int i = 0; i < n; i ++) + { + int k = lt[i]; + + lvl[i] = new float[k]; + sal[i] = new float[k]; + sal_err[i] = new float[k]; + sprb[i] = new float[k]; + tmp[i] = new float[k]; + tmp_err[i] = new float[k]; + tprb[i] = new float[k]; + clm_sal[i] = new float[k]; + cssd[i] = new float[k]; + clm_tmp[i] = new float[k]; + ctsd[i] = new float[k]; + flg[i] = new int[k]; + } +} + + + +void +RTOFSOb:: + read(FILE * f) +{ + read_float_array(f, btm, n); + read_float_array(f, lat, n); + read_float_array(f, lon, n); + read_int_array(f, ls, n); + read_int_array(f, lt, n); + read_int_array(f, sal_typ, n); + read_float_array(f, sqc, n); + read_int_array(f, tmp_typ, n); + read_float_array(f, tqc, n); + + allocate2d(); + + for (int i = 0; i < n; i ++) + { + int k = lt[i]; + + read_float_array(f, lvl[i], k); + read_float_array(f, sal[i], k); + read_float_array(f, sal_err[i], k); + read_float_array(f, sprb[i], k); + read_float_array(f, tmp[i], k); + read_float_array(f, tmp_err[i], k); + read_float_array(f, tprb[i], k); + read_float_array(f, clm_sal[i], k); + read_float_array(f, cssd[i], k); + read_float_array(f, clm_tmp[i], k); + read_float_array(f, ctsd[i], k); + read_int_array(f, flg[i], k); + } +} + + +void +RTOFSOb:: + print(std::string DirectoryName) +{ + if (! file_exists(DirectoryName)) + { + cerr << "Directory " << DirectoryName << "doesn't exist" << endl; + exit(1); + } + print_float_array(DirectoryName + "/latitude", lat, n); + print_float_array(DirectoryName + "/longitude", lon, n); + print_float_array(DirectoryName + "/btm", btm, n); + print_float_array(DirectoryName + "/tqc", tqc, n); + print_float_array(DirectoryName + "/sqc", sqc, n); + print_int_array(DirectoryName + "/lt", lt, n); + print_int_array(DirectoryName + "/ls", ls, n); + print_int_array(DirectoryName + "/sal_typ", sal_typ, n); + print_int_array(DirectoryName + "/tmp_typ", sal_typ, n); + + print_2d_float_array(DirectoryName + "/tmp", tmp, n, lt); + print_2d_float_array(DirectoryName + "/sal", sal, n, lt); + + // print lvl2d array + { + ofstream o; + o.open(DirectoryName + "/lvl2d"); + for (int i = 0; i < n; i ++) + for (int j = 0; j < lt[i]; j ++) + o + << j + << endl + ; + o.close(); + } +} + + +} // namespace rtofs diff --git a/utils/obsproc/rtofs/RTOFSOb.h b/utils/obsproc/rtofs/RTOFSOb.h new file mode 100644 index 000000000..1e8429ac5 --- /dev/null +++ b/utils/obsproc/rtofs/RTOFSOb.h @@ -0,0 +1,63 @@ +#ifndef RTOFSOb__H +#define RTOFSOb__H + + +#include +#include + +#include +using std::string; + +namespace rtofs +{ + + +class RTOFSOb +{ +public: + RTOFSOb(int n, int mx_lvl, int vrsn); + + void read(FILE * f); + int NumberOfObservations() { return n; } + void print(std::string DirectoryName); + + + float * btm; + float * lat; + float * lon; + int * ls; + int * lt; + int * sal_typ; + float * sqc; + int * tmp_typ; + float * tqc; + + float ** lvl; + float ** sal; + float ** sal_err; + float ** sprb; + float ** tmp; + float ** tmp_err; + float ** tprb; + float ** clm_sal; + float ** cssd; + float ** clm_tmp; + float ** ctsd; + int ** flg; + + std::time_t * dtg; + +private: + int n; + int mx_lvl; + int vrsn; + + void allocate(); + void allocate2d(); + +}; // class RTOFSOb + + +} // namespace rtofs + +#endif // RTOFSOb__H diff --git a/utils/obsproc/rtofs/util.cc b/utils/obsproc/rtofs/util.cc new file mode 100644 index 000000000..7a49531de --- /dev/null +++ b/utils/obsproc/rtofs/util.cc @@ -0,0 +1,203 @@ +// #define _BSD_SOURCE +#include +#include // used in file_exists + +#include +using std::cout; +using std::endl; + +#include +using std::ofstream; + +#include +using std::setprecision; + +#include "util.h" + + +namespace rtofs +{ + + +void skip8bytes(FILE * f) +{ + fseek(f, 8, SEEK_CUR); +} + + +void + read_float_array(FILE * f, float * a, int n) +{ + skip8bytes(f); + + fread(a, sizeof(float), n, f); + + uint32_t * data = (uint32_t *) a; + + for (int i = 0; i < n; i ++) + data[i] = be32toh(data[i]); + +} // read_float_array + + +void + read_int_array(FILE * f, int * a, int n) +{ + skip8bytes(f); + + fread(a, sizeof(int), n, f); + + uint32_t * data = (uint32_t *) a; + + for (int i = 0; i < n; i ++) + data[i] = be32toh(data[i]); + +} // read_int_array + + + +void + print_level(std::string filename, float ** a, int n, int level) +{ + ofstream o(filename); + o + << std::fixed + << std::setprecision(9) + ; + for (int i = 0; i < n; i ++) + o + << a[i] [level] + << endl + ; + o.close(); +} + + +void + print_int_array(std::string filename, int * a, int n) +{ + ofstream o(filename); + for (int i = 0; i < n; i ++) + o + << a[i] + << endl + ; + o.close(); +} + + +void + print_float_array(std::string filename, float * a, int n) +{ + ofstream o(filename); + o + << std::fixed + << std::setprecision(9) + ; + for (int i = 0; i < n; i ++) + o + << a[i] + << endl + ; +} + +void + print_2d_float_array(std::string filename, float ** a, int n, int * dim2) +{ + ofstream o(filename); + o + << std::fixed + << std::setprecision(9) + ; + for (int i = 0; i < n; i ++) + for (int j = 0; j < dim2[i]; j ++) + o + << a[i] [j] + << endl + ; +} + + +bool + file_exists (const std::string& name) +{ + struct stat buffer; + return (stat (name.c_str(), &buffer) == 0); +} + + +int + read_int(FILE * f) +{ + int dummy; + fread(&dummy, sizeof(int), 1, f); + dummy = (int) be32toh(dummy); + return dummy; + +} // read_int + + + +void + read_real_array(FILE * f, float * a, int n) +{ + fread(a, sizeof(float), n, f); + + uint32_t * data = (uint32_t *) a; + + for (int i = 0; i < n; i ++) + data[i] = be32toh(data[i]); + +} // read_real_array + + +float * + alloc_read_float_array(FILE * f, int n) +{ + float * a = new float[n]; + skip8bytes(f); + read_real_array(f, a, n); + return a; +} + + +void + print_int_array(int * a, int n) +{ + for (int i = 0; i < n; i ++) + cout + // << setw(12) + // << i + 1 + // << setw(12) + << a[i] + << endl + ; +} + + +int * + alloc_read_int_array(FILE * f, int n) +{ + int * a = new int[n]; + skip8bytes(f); + read_int_array(f, a, n); + return a; +} + + +void + print_real_array(float * a, int n) +{ + cout + << std::fixed + << std::setprecision(9) + ; + for (int i = 0; i < n; i ++) + cout + << a[i] + << endl + ; +} + + +} // namespace rtofs diff --git a/utils/obsproc/rtofs/util.h b/utils/obsproc/rtofs/util.h new file mode 100644 index 000000000..a4b1ebbb4 --- /dev/null +++ b/utils/obsproc/rtofs/util.h @@ -0,0 +1,24 @@ +#ifndef RTOFS_util__H +#define RTOFS_util__H + +#include + +namespace rtofs +{ + +bool file_exists (const std::string& name); + +void skip8bytes(FILE * f); +int read_int(FILE * f); +void read_float_array(FILE * f, float * a, int n); +void read_int_array(FILE * f, int * a, int n); + +void print_float_array(std::string filename, float * a, int n); +void print_int_array(std::string filename, int * a, int n); +void print_2d_float_array(std::string filename, float ** a, int n, int * dim2); +float * alloc_read_float_array(FILE * f, int n); + +} // namespace rtofs + + +#endif // RTOFS_util__H From f7fbda3df291c1a80f3a31b281d62e89fb0feb2c Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Wed, 15 May 2024 10:48:34 -0500 Subject: [PATCH 03/18] basic working code --- utils/obsproc/CMakeLists.txt | 24 ++ utils/obsproc/RTOFSInSitu.h | 292 +++++++--------------- utils/obsproc/applications/CMakeLists.txt | 13 +- 3 files changed, 122 insertions(+), 207 deletions(-) diff --git a/utils/obsproc/CMakeLists.txt b/utils/obsproc/CMakeLists.txt index 4e4e3feb6..5ea5aeea3 100644 --- a/utils/obsproc/CMakeLists.txt +++ b/utils/obsproc/CMakeLists.txt @@ -1 +1,25 @@ + +add_subdirectory(rtofs) + +# add_library(obsproc STATIC + # Ghrsst2Ioda.h + # IcecAmsr2Ioda.h + # NetCDFToIodaConverter.h + # Rads2Ioda.h + # RTOFSInSitu.h + # Smap2Ioda.h + # Smos2Ioda.h + # superob.h + # util.h + # Viirsaod2Ioda.h +# ) +# target_compile_features( + # onsproc + # PUBLIC cxx_std_17 +# ) +# target_link_libraries( + # obsproc + # PUBLIC oops ioda NetCDF::NetCDF_CXX +# ) + add_subdirectory(applications) diff --git a/utils/obsproc/RTOFSInSitu.h b/utils/obsproc/RTOFSInSitu.h index 5c8aca9d8..f8cff5432 100644 --- a/utils/obsproc/RTOFSInSitu.h +++ b/utils/obsproc/RTOFSInSitu.h @@ -14,214 +14,98 @@ #include "ioda/ObsGroup.h" #include "NetCDFToIodaConverter.h" -#include "superob.h" - - -namespace gdasapp { - - class RTOFSInSitu : public NetCDFToIodaConverter { - public: - explicit RTOFSInSitu(const eckit::Configuration & fullConfig, const eckit::mpi::Comm & comm) - : NetCDFToIodaConverter(fullConfig, comm) { - variable_ = "seaSurfaceTemperature"; - } - - // Read netcdf file and populate iodaVars - gdasapp::obsproc::iodavars::IodaVars providerToIodaVars(const std::string fileName) final { - oops::Log::info() << "Processing files provided by RTOFS" << std::endl; - - // Get the sst bounds from the configuration - float sstMin; - fullConfig_.get("bounds.min", sstMin); - float sstMax; - fullConfig_.get("bounds.max", sstMax); - - /* - // Open the NetCDF file in read-only mode - netCDF::NcFile ncFile(fileName, netCDF::NcFile::read); - oops::Log::info() << "Reading... " << fileName << std::endl; - // Get number of obs - int dimLon = ncFile.getDim("lon").getSize(); - int dimLat = ncFile.getDim("lat").getSize(); - int dimTime = ncFile.getDim("time").getSize(); - - // Read non-optional metadata: datetime, longitude and latitude - // latitude - std::vector lat(dimLat); - ncFile.getVar("lat").getVar(lat.data()); - - // longitude - std::vector lon(dimLon); - ncFile.getVar("lon").getVar(lon.data()); - - // Generate the lat-lon grid - std::vector> lon2d(dimLat, std::vector(dimLon)); - std::vector> lat2d(dimLat, std::vector(dimLon)); - for (int i = 0; i < dimLat; ++i) { - for (int j = 0; j < dimLon; ++j) { - lon2d[i][j] = lon[j]; - lat2d[i][j] = lat[i]; - } - } - - // Read the reference time - std::vector refTime(dimTime); - ncFile.getVar("time").getVar(refTime.data()); - std::string refDate; - ncFile.getVar("time").getAtt("units").getValues(refDate); - - // Reformat the reference time - std::regex dateRegex(R"(\b\d{4}-\d{2}-\d{2} \d{2}:\d{2}:\d{2}\b)"); - std::smatch match; - std::regex_search(refDate, match, dateRegex); - refDate = match.str(); - std::tm tmStruct = {}; - std::istringstream ss(refDate); - ss >> std::get_time(&tmStruct, "%Y-%m-%d %H:%M:%S"); - std::ostringstream isoFormatted; - isoFormatted << std::put_time(&tmStruct, "seconds since %Y-%m-%dT%H:%M:%SZ"); - refDate = isoFormatted.str(); - - // Read sst_dtime to add to the reference time - // TODO(AMG): What's below does not read the field the same way python does - std::vector sstdTime(dimTime*dimLat*dimLon); - ncFile.getVar("sst_dtime").getVar(sstdTime.data()); - float dtimeScaleFactor; - ncFile.getVar("sst_dtime").getAtt("scale_factor").getValues(&dtimeScaleFactor); - int32_t dtimeFillValue; - ncFile.getVar("sst_dtime").getAtt("_FillValue").getValues(&dtimeFillValue); - - // Read SST ObsValue - std::vector sstObsVal(dimTime*dimLat*dimLon); - ncFile.getVar("sea_surface_temperature").getVar(sstObsVal.data()); - float sstOffSet; - ncFile.getVar("sea_surface_temperature").getAtt("add_offset").getValues(&sstOffSet); - float sstScaleFactor; - ncFile.getVar("sea_surface_temperature").getAtt("scale_factor").getValues(&sstScaleFactor); - - // Read SST Bias - std::vector sstObsBias(dimTime*dimLat*dimLon); - ncFile.getVar("sses_bias").getVar(sstObsBias.data()); - float biasScaleFactor; - ncFile.getVar("sses_bias").getAtt("scale_factor").getValues(&biasScaleFactor); - - // Read Error - std::vector sstObsErr(dimTime*dimLat*dimLon); - ncFile.getVar("sses_standard_deviation").getVar(sstObsErr.data()); - float errOffSet; - ncFile.getVar("sses_standard_deviation").getAtt("add_offset").getValues(&errOffSet); - float errScaleFactor; - ncFile.getVar("sses_standard_deviation").getAtt("scale_factor").getValues(&errScaleFactor); - - // Read preQc - signed char sstPreQC[dimTime][dimLat][dimLon]; - ncFile.getVar("quality_level").getVar(sstPreQC); - - // Apply scaling/unit change and compute the necessary fields - std::vector> mask(dimLat, std::vector(dimLon)); - std::vector> sst(dimLat, std::vector(dimLon)); - std::vector> obserror(dimLat, std::vector(dimLon)); - std::vector> preqc(dimLat, std::vector(dimLon)); - std::vector> seconds(dimLat, std::vector(dimLon)); - size_t index = 0; - for (int i = 0; i < dimLat; i++) { - for (int j = 0; j < dimLon; j++) { - // provider's QC flag - // Note: the qc flags in GDS2.0 run from 0 to 5, with higher numbers being better. - // IODA typically expects 0 to be good, and larger numbers to be worse so the - // provider's QC is flipped - preqc[i][j] = 5 - static_cast(sstPreQC[0][i][j]); - - // bias corrected sst, regressed to the drifter depth - // Remove added sstOffSet for Celsius - sst[i][j] = static_cast(sstObsVal[index]) * sstScaleFactor - - static_cast(sstObsBias[index]) * biasScaleFactor; - // mask - if (sst[i][j] >= sstMin && sst[i][j] <= sstMax && preqc[i][j] ==0) { - mask[i][j] = 1; - } else { - mask[i][j] = 0; - } - - // obs error - // TODO(Somebody): add sampled std. dev. of sst to the total obs error - obserror[i][j] = static_cast(sstObsErr[index]) * errScaleFactor + errOffSet; - - // epoch time in seconds and avoid to use FillValue in calculation - if (sstdTime[index] == dtimeFillValue) { - seconds[i][j] = 0; - } else { - seconds[i][j] = static_cast(sstdTime[index]) * dtimeScaleFactor - + static_cast(refTime[0]); - } - index++; - } - } - - // Superobing - // TODO(Guillaume): Save the sampling std dev of sst so it can be used - // as a proxi for obs error - std::vector> sst_s; - std::vector> lon2d_s; - std::vector> lat2d_s; - std::vector> obserror_s; - std::vector> seconds_s; - if ( fullConfig_.has("binning") ) { - sst_s = gdasapp::superobutils::subsample2D(sst, mask, fullConfig_); - lon2d_s = gdasapp::superobutils::subsample2D(lon2d, mask, fullConfig_); - lat2d_s = gdasapp::superobutils::subsample2D(lat2d, mask, fullConfig_); - obserror_s = gdasapp::superobutils::subsample2D(obserror, mask, fullConfig_); - seconds_s = gdasapp::superobutils::subsample2D(seconds, mask, fullConfig_); - } else { - sst_s = sst; - lon2d_s = lon2d; - lat2d_s = lat2d; - obserror_s = obserror; - seconds_s = seconds; - } + +#include "rtofs/RTOFSDataFile.h" +#include "rtofs/RTOFSOb.h" + + +namespace gdasapp +{ + + +class RTOFSInSitu: + public NetCDFToIodaConverter +{ +public: + explicit RTOFSInSitu( + const eckit::Configuration & fullConfig, + const eckit::mpi::Comm & comm + ): + NetCDFToIodaConverter(fullConfig, comm) + { + variable_ = "temperature"; + } + + // Read binary file and populate iodaVars + gdasapp::obsproc::iodavars::IodaVars + providerToIodaVars(const std::string fileName) final + { + oops::Log::info() + << "Processing files provided by RTOFS" + << std::endl + ; + + // Set the int metadata names + std::vector intMetadataNames = {"temperature"}; + + // Set the float metadata name + std::vector floatMetadataNames = {}; + + // read the file + oops::Log::info() + << "Processing file " + << fileName + << std::endl + ; + + rtofs::RTOFSDataFile RTOFSFile(fileName); + int n = RTOFSFile.NumberOfObservations(); + rtofs::RTOFSOb & ob = RTOFSFile.observations(); + + int NumberOfTemperatureValues = 0; + for (int i = 0; i < n; i ++) + NumberOfTemperatureValues += ob.lt[i]; + + gdasapp::obsproc::iodavars::IodaVars iodaVars( + NumberOfTemperatureValues, + floatMetadataNames, + intMetadataNames + ); + iodaVars.referenceDate_ = "seconds since 1970-01-01T00:00:00Z"; + + // debugging output: + // std::string DataDirectory = "/home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/build/gdas-utils/test/obsproc/RTOFSdata"; + // ob.print(DataDirectory); + +// oops::Log::debug() << "--- iodaVars.location_: " << iodaVars.location_ << std::endl; + + int k = 0; + for (int i = 0; i < n; i ++) + for (int j = 0; j < ob.lt[i]; j ++) + { + iodaVars.longitude_(k) = ob.lon[i]; + iodaVars.latitude_(k) = ob.lat[i]; + iodaVars.obsVal_(k) = ob.tmp[i][j]; + iodaVars.obsError_(k) = ob.tmp_err[i][j]; + iodaVars.datetime_(k) = ob.dtg[i]; + // iodaVars.preQc_(k) = oneDimFlagsVal[i]; + // iodaVars.intMetadata_.row(k) << -999; + + k ++; + } + + // basic test for iodaVars.trim +/* + Eigen::Array mask + = (iodaVars.obsVal_ > 0.0); + iodaVars.trim(mask); */ - // number of obs after subsampling - // int nobs = sst_s.size() * sst_s[0].size(); - int nobs = 0; - // Set the int metadata names - // std::vector intMetadataNames = {"oceanBasin"}; - std::vector intMetadataNames = {"dummyBasin"}; + return iodaVars; - // Set the float metadata name - std::vector floatMetadataNames = {}; + }; // providerToIodaVars - // Create instance of iodaVars object - gdasapp::obsproc::iodavars::IodaVars iodaVars(nobs, floatMetadataNames, intMetadataNames); +}; // class RTOFSInSitu -/* - // Reference time is Jan 01 1981 00:00:00 GMT+0000 - iodaVars.referenceDate_ = refDate; - - // Store into eigen arrays - int loc(0); - for (int i = 0; i < sst_s.size(); i++) { - for (int j = 0; j < sst_s[0].size(); j++) { - iodaVars.longitude_(loc) = lon2d_s[i][j]; - iodaVars.latitude_(loc) = lat2d_s[i][j]; - iodaVars.obsVal_(loc) = sst_s[i][j]; - iodaVars.obsError_(loc) = obserror_s[i][j]; - iodaVars.preQc_(loc) = 0; - iodaVars.datetime_(loc) = seconds_s[i][j]; - // Store optional metadata, set ocean basins to -999 for now - iodaVars.intMetadata_.row(loc) << -999; - loc += 1; - } - } - - // Basic QC - Eigen::Array boundsCheck = - (iodaVars.obsVal_ > sstMin && iodaVars.obsVal_ < sstMax && iodaVars.datetime_ > 0.0); - iodaVars.trim(boundsCheck); - */ - - return iodaVars; - }; - }; // class RTOFSInSitu } // namespace gdasapp diff --git a/utils/obsproc/applications/CMakeLists.txt b/utils/obsproc/applications/CMakeLists.txt index 04c7ad3b7..6fd956b89 100644 --- a/utils/obsproc/applications/CMakeLists.txt +++ b/utils/obsproc/applications/CMakeLists.txt @@ -1,9 +1,16 @@ list( APPEND gdasapp_provider2ioda_src_files - gdas_obsprovider2ioda.cc - gdas_obsprovider2ioda.h - ) + gdas_obsprovider2ioda.cc + gdas_obsprovider2ioda.h +) + ecbuild_add_executable( TARGET gdas_obsprovider2ioda.x SOURCES ${gdasapp_provider2ioda_src_files} ) target_compile_features( gdas_obsprovider2ioda.x PUBLIC cxx_std_17) + target_link_libraries( gdas_obsprovider2ioda.x PUBLIC oops ioda NetCDF::NetCDF_CXX) + +# message("CMAKE_SOURCE_DIR = ${CMAKE_SOURCE_DIR}") +link_directories(${CMAKE_SOURCE_DIR}/rtofs) +target_link_libraries( gdas_obsprovider2ioda.x PRIVATE rtofs) + From 483430b3a83dc85489aff3251ce2255d5e1e8321 Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Wed, 15 May 2024 10:53:02 -0500 Subject: [PATCH 04/18] basic working code --- utils/test/testinput/gdas_rtofsinsitu.yaml | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/utils/test/testinput/gdas_rtofsinsitu.yaml b/utils/test/testinput/gdas_rtofsinsitu.yaml index 8b6a129e6..4c8684d9d 100644 --- a/utils/test/testinput/gdas_rtofsinsitu.yaml +++ b/utils/test/testinput/gdas_rtofsinsitu.yaml @@ -1,10 +1,11 @@ provider: RTOFS window begin: 2021-03-24T15:00:00Z window end: 2021-03-24T21:00:00Z -output file: rtofsinsitu_2024032600.profile +output file: rtofsinsitu_2024032600.profile.nc4 input files: -- 2024032600.profile -- 2024032600.sfc +- rtofsinsitu_2024032600.profile +# - 2024032600.profile +# - 2024032600.sfc # test: # reference filename: testref/rtofsinsitu.test From 56b71c61fb45b73ebfd29a3e4edd9a5eb9ebd06a Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Wed, 15 May 2024 17:15:41 -0500 Subject: [PATCH 05/18] fixed almost all style issues --- utils/obsproc/RTOFSInSitu.h | 144 ++++++----- .../applications/gdas_obsprovider2ioda.h | 2 +- utils/obsproc/rtofs/RTOFSDataFile.cc | 212 ++++++---------- utils/obsproc/rtofs/RTOFSDataFile.h | 31 ++- utils/obsproc/rtofs/RTOFSOb.cc | 229 +++++++++--------- utils/obsproc/rtofs/RTOFSOb.h | 97 ++++---- utils/obsproc/rtofs/util.cc | 225 ++++++++--------- utils/obsproc/rtofs/util.h | 11 +- 8 files changed, 440 insertions(+), 511 deletions(-) diff --git a/utils/obsproc/RTOFSInSitu.h b/utils/obsproc/RTOFSInSitu.h index f8cff5432..c5f3a92fa 100644 --- a/utils/obsproc/RTOFSInSitu.h +++ b/utils/obsproc/RTOFSInSitu.h @@ -24,88 +24,84 @@ namespace gdasapp class RTOFSInSitu: - public NetCDFToIodaConverter + public NetCDFToIodaConverter { -public: - explicit RTOFSInSitu( - const eckit::Configuration & fullConfig, - const eckit::mpi::Comm & comm - ): - NetCDFToIodaConverter(fullConfig, comm) - { - variable_ = "temperature"; - } - - // Read binary file and populate iodaVars - gdasapp::obsproc::iodavars::IodaVars - providerToIodaVars(const std::string fileName) final - { - oops::Log::info() - << "Processing files provided by RTOFS" - << std::endl - ; - - // Set the int metadata names - std::vector intMetadataNames = {"temperature"}; - - // Set the float metadata name - std::vector floatMetadataNames = {}; - - // read the file - oops::Log::info() - << "Processing file " - << fileName - << std::endl - ; - - rtofs::RTOFSDataFile RTOFSFile(fileName); - int n = RTOFSFile.NumberOfObservations(); - rtofs::RTOFSOb & ob = RTOFSFile.observations(); - - int NumberOfTemperatureValues = 0; - for (int i = 0; i < n; i ++) - NumberOfTemperatureValues += ob.lt[i]; - - gdasapp::obsproc::iodavars::IodaVars iodaVars( - NumberOfTemperatureValues, - floatMetadataNames, - intMetadataNames - ); - iodaVars.referenceDate_ = "seconds since 1970-01-01T00:00:00Z"; - - // debugging output: - // std::string DataDirectory = "/home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/build/gdas-utils/test/obsproc/RTOFSdata"; - // ob.print(DataDirectory); + public: + explicit RTOFSInSitu( + const eckit::Configuration & fullConfig, + const eckit::mpi::Comm & comm): + NetCDFToIodaConverter(fullConfig, comm) + { + variable_ = "temperature"; + } + + // Read binary file and populate iodaVars + gdasapp::obsproc::iodavars::IodaVars + providerToIodaVars(const std::string fileName) final + { + oops::Log::info() + << "Processing files provided by RTOFS" + << std::endl; + + // Set the int metadata names + std::vector intMetadataNames = {"temperature"}; + + // Set the float metadata name + std::vector floatMetadataNames = {}; + + // read the file + oops::Log::info() + << "Processing file " + << fileName + << std::endl; + + rtofs::RTOFSDataFile RTOFSFile(fileName); + int n = RTOFSFile.NumberOfObservations(); + rtofs::RTOFSOb & ob = RTOFSFile.observations(); + + int NumberOfTemperatureValues = 0; + for (int i = 0; i < n; i ++) + NumberOfTemperatureValues += ob.lt[i]; + + gdasapp::obsproc::iodavars::IodaVars iodaVars( + NumberOfTemperatureValues, + floatMetadataNames, + intMetadataNames); + iodaVars.referenceDate_ = "seconds since 1970-01-01T00:00:00Z"; + + // debugging output: + // std::string DataDirectory = + // "/home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/build/gdas-utils/test/obsproc/RTOFSdata"; + // ob.print(DataDirectory); // oops::Log::debug() << "--- iodaVars.location_: " << iodaVars.location_ << std::endl; - int k = 0; - for (int i = 0; i < n; i ++) - for (int j = 0; j < ob.lt[i]; j ++) - { - iodaVars.longitude_(k) = ob.lon[i]; - iodaVars.latitude_(k) = ob.lat[i]; - iodaVars.obsVal_(k) = ob.tmp[i][j]; - iodaVars.obsError_(k) = ob.tmp_err[i][j]; - iodaVars.datetime_(k) = ob.dtg[i]; - // iodaVars.preQc_(k) = oneDimFlagsVal[i]; - // iodaVars.intMetadata_.row(k) << -999; - - k ++; - } - - // basic test for iodaVars.trim + int k = 0; + for (int i = 0; i < n; i ++) + for (int j = 0; j < ob.lt[i]; j ++) + { + iodaVars.longitude_(k) = ob.lon[i]; + iodaVars.latitude_(k) = ob.lat[i]; + iodaVars.obsVal_(k) = ob.tmp[i][j]; + iodaVars.obsError_(k) = ob.tmp_err[i][j]; + iodaVars.datetime_(k) = ob.dtg[i]; + // iodaVars.preQc_(k) = oneDimFlagsVal[i]; + // iodaVars.intMetadata_.row(k) << -999; + + k++; + } + + // basic test for iodaVars.trim /* - Eigen::Array mask - = (iodaVars.obsVal_ > 0.0); - iodaVars.trim(mask); + Eigen::Array mask + = (iodaVars.obsVal_ > 0.0); + iodaVars.trim(mask); */ - return iodaVars; - - }; // providerToIodaVars - + return iodaVars; + }; // providerToIodaVars }; // class RTOFSInSitu + } // namespace gdasapp diff --git a/utils/obsproc/applications/gdas_obsprovider2ioda.h b/utils/obsproc/applications/gdas_obsprovider2ioda.h index c14c9a1a2..b3bcd5c89 100644 --- a/utils/obsproc/applications/gdas_obsprovider2ioda.h +++ b/utils/obsproc/applications/gdas_obsprovider2ioda.h @@ -8,8 +8,8 @@ #include "../Ghrsst2Ioda.h" #include "../IcecAmsr2Ioda.h" -#include "../RTOFSInSitu.h" #include "../Rads2Ioda.h" +#include "../RTOFSInSitu.h" #include "../Smap2Ioda.h" #include "../Smos2Ioda.h" #include "../Viirsaod2Ioda.h" diff --git a/utils/obsproc/rtofs/RTOFSDataFile.cc b/utils/obsproc/rtofs/RTOFSDataFile.cc index e6620e1e6..79398d854 100644 --- a/utils/obsproc/rtofs/RTOFSDataFile.cc +++ b/utils/obsproc/rtofs/RTOFSDataFile.cc @@ -1,15 +1,16 @@ -#include -#include -#include +#include "RTOFSDataFile.h" +#include #include using std::cerr; using std::endl; using std::cout; -#include // std::get_time +#include // std::get_time + +#include +#include -#include "RTOFSDataFile.h" #include "RTOFSOb.h" #include "util.h" @@ -21,148 +22,93 @@ typedef char char12[12]; typedef char char7[7]; -std::time_t - SecondsSinceReferenceTime(char12 time) +std::time_t + SecondsSinceReferenceTime(char12 time) { - std::tm referenceTime = {}; - referenceTime.tm_year = 70; // 1970 - referenceTime.tm_mon = 0; // January (months are 0-based) - referenceTime.tm_mday = 1; // 1st day of the month - referenceTime.tm_hour = 0; - referenceTime.tm_min = 0; - referenceTime.tm_sec = 0; - std::time_t referenceTimestamp = std::mktime(&referenceTime); + std::tm referenceTime = {}; + referenceTime.tm_year = 70; // 1970 + referenceTime.tm_mon = 0; // January (months are 0-based) + referenceTime.tm_mday = 1; // 1st day of the month + referenceTime.tm_hour = 0; + referenceTime.tm_min = 0; + referenceTime.tm_sec = 0; + std::time_t referenceTimestamp = std::mktime(&referenceTime); - std::tm t = {}; - std::istringstream ss(time); - ss >> std::get_time(&t, "%Y%m%d%H%M"); + std::tm t = {}; + std::istringstream ss(time); + ss >> std::get_time(&t, "%Y%m%d%H%M"); - std::time_t timestamp = std::mktime(&t); - return std::difftime(timestamp, referenceTimestamp); - -} // SecondsSinceReferenceTime + std::time_t timestamp = std::mktime(&t); + return std::difftime(timestamp, referenceTimestamp); +} // SecondsSinceReferenceTime RTOFSDataFile:: RTOFSDataFile(std::string filename): - filename(filename) + filename(filename) { + if (!file_exists(filename)) + { + cerr << "File not found" << endl; + exit(1); + } - if (! file_exists(filename)) - { - cerr << "File not found" << endl; - exit(1); - } - - - const char * fn = filename.c_str(); - f = fopen(fn, "rb"); - if (!f) - { - cerr << "Error opening file " << fn << endl; - exit(1); - } - read_file(); + const char * fn = filename.c_str(); + f = fopen(fn, "rb"); + if (!f) + { + cerr << "Error opening file " << fn << endl; + exit(1); + } -} // RTOFSDataFile::RTOFSDataFile + read_file(); +} // RTOFSDataFile::RTOFSDataFile void RTOFSDataFile:: - read_file() + read_file() { - fseek(f, 4, SEEK_CUR); - - int n_read = read_int(f); - int mx_lvl = read_int(f); - int vrsn = read_int(f); - - ob = new RTOFSOb(n_read, mx_lvl, vrsn); - nobs = n_read; - - ob->read(f); - - skip8bytes(f); - - char12 * ob_dtg = new char12[n_read]; - fread(ob_dtg, sizeof(char[12]), n_read, f); - for (int i = 0; i < n_read; i ++) - ob->dtg[i] = SecondsSinceReferenceTime(ob_dtg[i]); - -/* - { - for (int i = 0; i < n_read; i ++) - { - for (int j = 0; j < 12; j ++) - cout << ob_dtg[i][j]; - - cout - << " = " - << SecondsSinceReferenceTime(ob_dtg[i]) - << " sec" - << endl - ; - cout << endl; - } - } -*/ - - skip8bytes(f); - - char12 * ob_rct = new char12[n_read]; - fread(ob_rct, sizeof(char[12]), n_read, f); - -/* - { - for (int i = 0; i < n_read; i ++) - { - for (int j = 0; j < 12; j ++) - cout << ob_rct[i][j]; - - cout - << " = " - << SecondsSinceReferenceTime(ob_rct[i]) - << " sec" - << endl - ; - cout << endl; - } - } -*/ - - - skip8bytes(f); - - char7 * ob_sgn = new char7[n_read]; - fread(ob_sgn, sizeof(char[7]), n_read, f); - -/* - for (int i = 0; i < n_read; i ++) - { - // for (int j = 0; j < 12; j ++) - // cout << ob_dtg[i][j]; - // cout << (unsigned int) ob_dtg[i][j] << " "; - // cout << ob_rct[i][j]; - for (int j = 0; j < 7; j ++) - cout << ob_sgn[i][j]; - cout << endl; - } -*/ - - if (vrsn == 2) - { - float ** glb_sal = new float * [n_read]; - for (int i = 0; i < n_read; i ++) - { - int k = ob->lt[i]; - glb_sal[i] = alloc_read_float_array(f, k); - } - } - - -} // read_file - - -} // namespace rtofs + fseek(f, 4, SEEK_CUR); + + int n_read = read_int(f); + int mx_lvl = read_int(f); + int vrsn = read_int(f); + + ob = new RTOFSOb(n_read, mx_lvl, vrsn); + nobs = n_read; + + ob->read(f); + + skip8bytes(f); + + char12 * ob_dtg = new char12[n_read]; + fread(ob_dtg, sizeof(char[12]), n_read, f); + for (int i = 0; i < n_read; i ++) + ob->dtg[i] = SecondsSinceReferenceTime(ob_dtg[i]); + + skip8bytes(f); + + char12 * ob_rct = new char12[n_read]; + fread(ob_rct, sizeof(char[12]), n_read, f); + + skip8bytes(f); + + char7 * ob_sgn = new char7[n_read]; + fread(ob_sgn, sizeof(char[7]), n_read, f); + + if (vrsn == 2) + { + float ** glb_sal = new float * [n_read]; + for (int i = 0; i < n_read; i ++) + { + int k = ob->lt[i]; + glb_sal[i] = alloc_read_float_array(f, k); + } + } +} // read_file + + +} // namespace rtofs diff --git a/utils/obsproc/rtofs/RTOFSDataFile.h b/utils/obsproc/rtofs/RTOFSDataFile.h index ea68e4852..40a5b45bf 100644 --- a/utils/obsproc/rtofs/RTOFSDataFile.h +++ b/utils/obsproc/rtofs/RTOFSDataFile.h @@ -1,5 +1,5 @@ -#ifndef RTOFSDataFile__H -#define RTOFSDataFile__H +#ifndef BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_RTOFSDATAFILE_H_ +#define BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_RTOFSDATAFILE_H_ #include @@ -11,23 +11,22 @@ class RTOFSOb; class RTOFSDataFile { -public: - RTOFSDataFile(std::string filename); - int NumberOfObservations() { return nobs; } - RTOFSOb & observations() { return * ob; } + public: + explicit RTOFSDataFile(std::string filename); + int NumberOfObservations() { return nobs; } + RTOFSOb & observations() { return * ob; } -private: - std::string filename; - int nobs; - FILE * f; - RTOFSOb * ob; + private: + std::string filename; + int nobs; + FILE * f; + RTOFSOb * ob; - void read_file(); + void read_file(); +}; // class RTOFSDataFile -}; // class RTOFSDataFile +} // namespace rtofs -} // namespace rtofs - -#endif // RTOFSDataFile__H +#endif // BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_RTOFSDATAFILE_H_ diff --git a/utils/obsproc/rtofs/RTOFSOb.cc b/utils/obsproc/rtofs/RTOFSOb.cc index 3e97fe5cc..66c057acf 100644 --- a/utils/obsproc/rtofs/RTOFSOb.cc +++ b/utils/obsproc/rtofs/RTOFSOb.cc @@ -1,3 +1,5 @@ +#include "RTOFSOb.h" + #include #include @@ -10,8 +12,6 @@ using std::cerr; #include using std::ofstream; - -#include "RTOFSOb.h" #include "util.h" @@ -20,142 +20,141 @@ namespace rtofs RTOFSOb:: - RTOFSOb(int n, int mx_lvl, int vrsn): - n(n), - mx_lvl(mx_lvl), - vrsn(vrsn) + RTOFSOb(int n, int mx_lvl, int vrsn): + n(n), + mx_lvl(mx_lvl), + vrsn(vrsn) { - dtg = new std::time_t[n]; - allocate(); + dtg = new std::time_t[n]; + allocate(); } void RTOFSOb:: - allocate() + allocate() { - lat = new float[n]; - lon = new float[n]; - btm = new float[n]; - ls = new int[n]; - lt = new int[n]; - sal_typ = new int[n]; - sqc = new float[n]; - tmp_typ = new int[n]; - tqc = new float[n]; - - lvl = new float * [n]; - sal = new float * [n]; - sal_err = new float * [n]; - sprb = new float * [n]; - tmp = new float * [n]; - tmp_err = new float * [n]; - clm_sal = new float * [n]; // skipped? - tprb = new float * [n]; - cssd = new float * [n]; - clm_tmp = new float * [n]; // skipped? - ctsd = new float * [n]; - flg = new int * [n]; + lat = new float[n]; + lon = new float[n]; + btm = new float[n]; + ls = new int[n]; + lt = new int[n]; + sal_typ = new int[n]; + sqc = new float[n]; + tmp_typ = new int[n]; + tqc = new float[n]; + + lvl = new float * [n]; + sal = new float * [n]; + sal_err = new float * [n]; + sprb = new float * [n]; + tmp = new float * [n]; + tmp_err = new float * [n]; + clm_sal = new float * [n]; // skipped? + tprb = new float * [n]; + cssd = new float * [n]; + clm_tmp = new float * [n]; // skipped? + ctsd = new float * [n]; + flg = new int * [n]; } void RTOFSOb:: - allocate2d() + allocate2d() { - for (int i = 0; i < n; i ++) - { - int k = lt[i]; - - lvl[i] = new float[k]; - sal[i] = new float[k]; - sal_err[i] = new float[k]; - sprb[i] = new float[k]; - tmp[i] = new float[k]; - tmp_err[i] = new float[k]; - tprb[i] = new float[k]; - clm_sal[i] = new float[k]; - cssd[i] = new float[k]; - clm_tmp[i] = new float[k]; - ctsd[i] = new float[k]; - flg[i] = new int[k]; - } + for (int i = 0; i < n; i ++) + { + int k = lt[i]; + + lvl[i] = new float[k]; + sal[i] = new float[k]; + sal_err[i] = new float[k]; + sprb[i] = new float[k]; + tmp[i] = new float[k]; + tmp_err[i] = new float[k]; + tprb[i] = new float[k]; + clm_sal[i] = new float[k]; + cssd[i] = new float[k]; + clm_tmp[i] = new float[k]; + ctsd[i] = new float[k]; + flg[i] = new int[k]; + } } void RTOFSOb:: - read(FILE * f) + read(FILE * f) { - read_float_array(f, btm, n); - read_float_array(f, lat, n); - read_float_array(f, lon, n); - read_int_array(f, ls, n); - read_int_array(f, lt, n); - read_int_array(f, sal_typ, n); - read_float_array(f, sqc, n); - read_int_array(f, tmp_typ, n); - read_float_array(f, tqc, n); - - allocate2d(); - - for (int i = 0; i < n; i ++) - { - int k = lt[i]; - - read_float_array(f, lvl[i], k); - read_float_array(f, sal[i], k); - read_float_array(f, sal_err[i], k); - read_float_array(f, sprb[i], k); - read_float_array(f, tmp[i], k); - read_float_array(f, tmp_err[i], k); - read_float_array(f, tprb[i], k); - read_float_array(f, clm_sal[i], k); - read_float_array(f, cssd[i], k); - read_float_array(f, clm_tmp[i], k); - read_float_array(f, ctsd[i], k); - read_int_array(f, flg[i], k); - } + read_float_array(f, btm, n); + read_float_array(f, lat, n); + read_float_array(f, lon, n); + read_int_array(f, ls, n); + read_int_array(f, lt, n); + read_int_array(f, sal_typ, n); + read_float_array(f, sqc, n); + read_int_array(f, tmp_typ, n); + read_float_array(f, tqc, n); + + allocate2d(); + + for (int i = 0; i < n; i ++) + { + int k = lt[i]; + + read_float_array(f, lvl[i], k); + read_float_array(f, sal[i], k); + read_float_array(f, sal_err[i], k); + read_float_array(f, sprb[i], k); + read_float_array(f, tmp[i], k); + read_float_array(f, tmp_err[i], k); + read_float_array(f, tprb[i], k); + read_float_array(f, clm_sal[i], k); + read_float_array(f, cssd[i], k); + read_float_array(f, clm_tmp[i], k); + read_float_array(f, ctsd[i], k); + read_int_array(f, flg[i], k); + } } -void +void RTOFSOb:: - print(std::string DirectoryName) + print(std::string DirectoryName) { - if (! file_exists(DirectoryName)) - { - cerr << "Directory " << DirectoryName << "doesn't exist" << endl; - exit(1); - } - print_float_array(DirectoryName + "/latitude", lat, n); - print_float_array(DirectoryName + "/longitude", lon, n); - print_float_array(DirectoryName + "/btm", btm, n); - print_float_array(DirectoryName + "/tqc", tqc, n); - print_float_array(DirectoryName + "/sqc", sqc, n); - print_int_array(DirectoryName + "/lt", lt, n); - print_int_array(DirectoryName + "/ls", ls, n); - print_int_array(DirectoryName + "/sal_typ", sal_typ, n); - print_int_array(DirectoryName + "/tmp_typ", sal_typ, n); - - print_2d_float_array(DirectoryName + "/tmp", tmp, n, lt); - print_2d_float_array(DirectoryName + "/sal", sal, n, lt); - - // print lvl2d array - { - ofstream o; - o.open(DirectoryName + "/lvl2d"); - for (int i = 0; i < n; i ++) - for (int j = 0; j < lt[i]; j ++) - o - << j - << endl - ; - o.close(); - } -} - - -} // namespace rtofs + if (!file_exists(DirectoryName)) + { + cerr << "Directory " << DirectoryName << "doesn't exist" << endl; + exit(1); + } + print_float_array(DirectoryName + "/latitude", lat, n); + print_float_array(DirectoryName + "/longitude", lon, n); + print_float_array(DirectoryName + "/btm", btm, n); + print_float_array(DirectoryName + "/tqc", tqc, n); + print_float_array(DirectoryName + "/sqc", sqc, n); + print_int_array(DirectoryName + "/lt", lt, n); + print_int_array(DirectoryName + "/ls", ls, n); + print_int_array(DirectoryName + "/sal_typ", sal_typ, n); + print_int_array(DirectoryName + "/tmp_typ", sal_typ, n); + + print_2d_float_array(DirectoryName + "/tmp", tmp, n, lt); + print_2d_float_array(DirectoryName + "/sal", sal, n, lt); + + // print lvl2d array + { + ofstream o; + o.open(DirectoryName + "/lvl2d"); + for (int i = 0; i < n; i ++) + for (int j = 0; j < lt[i]; j ++) + o + << j + << endl; + o.close(); + } +} // RTOFSOb::print + + +} // namespace rtofs diff --git a/utils/obsproc/rtofs/RTOFSOb.h b/utils/obsproc/rtofs/RTOFSOb.h index 1e8429ac5..fdd14df0c 100644 --- a/utils/obsproc/rtofs/RTOFSOb.h +++ b/utils/obsproc/rtofs/RTOFSOb.h @@ -1,5 +1,5 @@ -#ifndef RTOFSOb__H -#define RTOFSOb__H +#ifndef BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_RTOFSOB_H_ +#define BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_RTOFSOB_H_ #include @@ -14,50 +14,49 @@ namespace rtofs class RTOFSOb { -public: - RTOFSOb(int n, int mx_lvl, int vrsn); - - void read(FILE * f); - int NumberOfObservations() { return n; } - void print(std::string DirectoryName); - - - float * btm; - float * lat; - float * lon; - int * ls; - int * lt; - int * sal_typ; - float * sqc; - int * tmp_typ; - float * tqc; - - float ** lvl; - float ** sal; - float ** sal_err; - float ** sprb; - float ** tmp; - float ** tmp_err; - float ** tprb; - float ** clm_sal; - float ** cssd; - float ** clm_tmp; - float ** ctsd; - int ** flg; - - std::time_t * dtg; - -private: - int n; - int mx_lvl; - int vrsn; - - void allocate(); - void allocate2d(); - -}; // class RTOFSOb - - -} // namespace rtofs - -#endif // RTOFSOb__H + public: + RTOFSOb(int n, int mx_lvl, int vrsn); + + void read(FILE * f); + int NumberOfObservations() { return n; } + void print(std::string DirectoryName); + + + float * btm; + float * lat; + float * lon; + int * ls; + int * lt; + int * sal_typ; + float * sqc; + int * tmp_typ; + float * tqc; + + float ** lvl; + float ** sal; + float ** sal_err; + float ** sprb; + float ** tmp; + float ** tmp_err; + float ** tprb; + float ** clm_sal; + float ** cssd; + float ** clm_tmp; + float ** ctsd; + int ** flg; + + std::time_t * dtg; + + private: + int n; + int mx_lvl; + int vrsn; + + void allocate(); + void allocate2d(); +}; // class RTOFSOb + + +} // namespace rtofs + +#endif // BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_RTOFSOB_H_ diff --git a/utils/obsproc/rtofs/util.cc b/utils/obsproc/rtofs/util.cc index 7a49531de..f30de85e0 100644 --- a/utils/obsproc/rtofs/util.cc +++ b/utils/obsproc/rtofs/util.cc @@ -1,3 +1,5 @@ +#include "util.h" + // #define _BSD_SOURCE #include #include // used in file_exists @@ -12,8 +14,6 @@ using std::ofstream; #include using std::setprecision; -#include "util.h" - namespace rtofs { @@ -21,183 +21,172 @@ namespace rtofs void skip8bytes(FILE * f) { - fseek(f, 8, SEEK_CUR); + fseek(f, 8, SEEK_CUR); } void - read_float_array(FILE * f, float * a, int n) + read_float_array(FILE * f, float * a, int n) { - skip8bytes(f); + skip8bytes(f); - fread(a, sizeof(float), n, f); + fread(a, sizeof(float), n, f); - uint32_t * data = (uint32_t *) a; + // uint32_t * data = (uint32_t *) a; + uint32_t * data = reinterpret_cast(a); - for (int i = 0; i < n; i ++) - data[i] = be32toh(data[i]); - -} // read_float_array + for (int i = 0; i < n; i ++) + data[i] = be32toh(data[i]); +} // read_float_array void - read_int_array(FILE * f, int * a, int n) + read_int_array(FILE * f, int * a, int n) { - skip8bytes(f); - - fread(a, sizeof(int), n, f); + skip8bytes(f); - uint32_t * data = (uint32_t *) a; + fread(a, sizeof(int), n, f); - for (int i = 0; i < n; i ++) - data[i] = be32toh(data[i]); + // uint32_t * data = (uint32_t *) a; + uint32_t * data = reinterpret_cast(a); -} // read_int_array + for (int i = 0; i < n; i ++) + data[i] = be32toh(data[i]); +} // read_int_array -void - print_level(std::string filename, float ** a, int n, int level) +void + print_level(std::string filename, float ** a, int n, int level) { - ofstream o(filename); - o - << std::fixed - << std::setprecision(9) - ; - for (int i = 0; i < n; i ++) - o - << a[i] [level] - << endl - ; - o.close(); + ofstream o(filename); + o + << std::fixed + << std::setprecision(9); + + for (int i = 0; i < n; i ++) + o + << a[i] [level] + << endl; + + o.close(); } -void - print_int_array(std::string filename, int * a, int n) +void + print_int_array(std::string filename, int * a, int n) { - ofstream o(filename); - for (int i = 0; i < n; i ++) - o - << a[i] - << endl - ; - o.close(); + ofstream o(filename); + for (int i = 0; i < n; i ++) + o + << a[i] + << endl; + o.close(); } -void - print_float_array(std::string filename, float * a, int n) +void + print_float_array(std::string filename, float * a, int n) { - ofstream o(filename); - o - << std::fixed - << std::setprecision(9) - ; - for (int i = 0; i < n; i ++) - o - << a[i] - << endl - ; + ofstream o(filename); + o + << std::fixed + << std::setprecision(9); + for (int i = 0; i < n; i ++) + o + << a[i] + << endl; } -void - print_2d_float_array(std::string filename, float ** a, int n, int * dim2) +void + print_2d_float_array(std::string filename, float ** a, int n, int * dim2) { - ofstream o(filename); - o - << std::fixed - << std::setprecision(9) - ; - for (int i = 0; i < n; i ++) - for (int j = 0; j < dim2[i]; j ++) - o - << a[i] [j] - << endl - ; + ofstream o(filename); + o + << std::fixed + << std::setprecision(9); + for (int i = 0; i < n; i ++) + for (int j = 0; j < dim2[i]; j ++) + o + << a[i] [j] + << endl; } -bool - file_exists (const std::string& name) +bool + file_exists(const std::string& name) { - struct stat buffer; - return (stat (name.c_str(), &buffer) == 0); + struct stat buffer; + return (stat (name.c_str(), &buffer) == 0); } int - read_int(FILE * f) + read_int(FILE * f) { - int dummy; - fread(&dummy, sizeof(int), 1, f); - dummy = (int) be32toh(dummy); - return dummy; - -} // read_int + int dummy; + fread(&dummy, sizeof(int), 1, f); + // dummy = (int) be32toh(dummy); + dummy = static_cast(be32toh(dummy)); + return dummy; +} // read_int void - read_real_array(FILE * f, float * a, int n) + read_real_array(FILE * f, float * a, int n) { - fread(a, sizeof(float), n, f); - - uint32_t * data = (uint32_t *) a; + fread(a, sizeof(float), n, f); - for (int i = 0; i < n; i ++) - data[i] = be32toh(data[i]); + // uint32_t * data = (uint32_t *) a; + uint32_t * data = reinterpret_cast(a); -} // read_real_array + for (int i = 0; i < n; i ++) + data[i] = be32toh(data[i]); +} // read_real_array -float * - alloc_read_float_array(FILE * f, int n) +float * + alloc_read_float_array(FILE * f, int n) { - float * a = new float[n]; - skip8bytes(f); - read_real_array(f, a, n); - return a; + float * a = new float[n]; + skip8bytes(f); + read_real_array(f, a, n); + return a; } -void - print_int_array(int * a, int n) +void + print_int_array(int * a, int n) { - for (int i = 0; i < n; i ++) - cout - // << setw(12) - // << i + 1 - // << setw(12) - << a[i] - << endl - ; + for (int i = 0; i < n; i ++) + cout + << a[i] + << endl; } -int * - alloc_read_int_array(FILE * f, int n) +int * + alloc_read_int_array(FILE * f, int n) { - int * a = new int[n]; - skip8bytes(f); - read_int_array(f, a, n); - return a; + int * a = new int[n]; + skip8bytes(f); + read_int_array(f, a, n); + return a; } -void - print_real_array(float * a, int n) +void + print_real_array(float * a, int n) { - cout - << std::fixed - << std::setprecision(9) - ; - for (int i = 0; i < n; i ++) - cout - << a[i] - << endl - ; + cout + << std::fixed + << std::setprecision(9); + for (int i = 0; i < n; i ++) + cout + << a[i] + << endl; } -} // namespace rtofs +} // namespace rtofs diff --git a/utils/obsproc/rtofs/util.h b/utils/obsproc/rtofs/util.h index a4b1ebbb4..6dfae46f7 100644 --- a/utils/obsproc/rtofs/util.h +++ b/utils/obsproc/rtofs/util.h @@ -1,12 +1,13 @@ -#ifndef RTOFS_util__H -#define RTOFS_util__H +#ifndef BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_UTIL_H_ +#define BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_UTIL_H_ #include +#include namespace rtofs { -bool file_exists (const std::string& name); +bool file_exists(const std::string& name); void skip8bytes(FILE * f); int read_int(FILE * f); @@ -18,7 +19,7 @@ void print_int_array(std::string filename, int * a, int n); void print_2d_float_array(std::string filename, float ** a, int n, int * dim2); float * alloc_read_float_array(FILE * f, int n); -} // namespace rtofs +} // namespace rtofs -#endif // RTOFS_util__H +#endif // BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_UTIL_H_ From bc8cfa3a5998a5764baaffd6269f83c1a862f287 Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Tue, 21 May 2024 12:55:27 -0500 Subject: [PATCH 06/18] removed the rtofs dir/lib; everhthing is in one header file; passes tests --- utils/obsproc/CMakeLists.txt | 25 +- utils/obsproc/RTOFSInSitu.h | 563 ++++++++++++++++++++-- utils/obsproc/applications/CMakeLists.txt | 7 +- 3 files changed, 533 insertions(+), 62 deletions(-) diff --git a/utils/obsproc/CMakeLists.txt b/utils/obsproc/CMakeLists.txt index 5ea5aeea3..9cce213e1 100644 --- a/utils/obsproc/CMakeLists.txt +++ b/utils/obsproc/CMakeLists.txt @@ -1,25 +1,2 @@ - -add_subdirectory(rtofs) - -# add_library(obsproc STATIC - # Ghrsst2Ioda.h - # IcecAmsr2Ioda.h - # NetCDFToIodaConverter.h - # Rads2Ioda.h - # RTOFSInSitu.h - # Smap2Ioda.h - # Smos2Ioda.h - # superob.h - # util.h - # Viirsaod2Ioda.h -# ) -# target_compile_features( - # onsproc - # PUBLIC cxx_std_17 -# ) -# target_link_libraries( - # obsproc - # PUBLIC oops ioda NetCDF::NetCDF_CXX -# ) - +# add_subdirectory(rtofs) add_subdirectory(applications) diff --git a/utils/obsproc/RTOFSInSitu.h b/utils/obsproc/RTOFSInSitu.h index c5f3a92fa..32b9b967f 100644 --- a/utils/obsproc/RTOFSInSitu.h +++ b/utils/obsproc/RTOFSInSitu.h @@ -1,28 +1,152 @@ #pragma once +#include +// #define _BSD_SOURCE +#include +#include // used in file_exists + +#include +#include +using std::ofstream; +#include // std::get_time +using std::setprecision; #include +using std::cerr; +using std::endl; +using std::cout; #include // NOLINT (using C API) -#include #include +using std::string; +#include #include -#include "eckit/config/LocalConfiguration.h" +#include "NetCDFToIodaConverter.h" +#include "oops/util/dateFunctions.h" +#include "oops/util/DateTime.h" -#include // NOLINT -#include "ioda/Group.h" -#include "ioda/ObsGroup.h" -#include "NetCDFToIodaConverter.h" -#include "rtofs/RTOFSDataFile.h" -#include "rtofs/RTOFSOb.h" + + + +namespace rtofs +{ + +bool file_exists(const std::string& name); + +void skip8bytes(FILE * f); +int read_int(FILE * f); +void read_float_array(FILE * f, float * a, int n); +void read_int_array(FILE * f, int * a, int n); + +void print_float_array(std::string filename, float * a, int n); +void print_int_array(std::string filename, int * a, int n); +void print_2d_float_array(std::string filename, float ** a, int n, int * dim2); +float * alloc_read_float_array(FILE * f, int n); + + + +class RTOFSOb +{ + public: + RTOFSOb(int n, int mx_lvl, int vrsn); + + typedef char char12[12]; + typedef char char7[7]; + + void read(FILE * f); + int NumberOfObservations() { return n; } + void print(std::string DirectoryName); + + + float * btm; // bottom depth + float * lat; // latitude + float * lon; // longitude + int * ls; // ?? + int * lt; // array of dimensions, the number of levels + int * sal_typ; // ?? salinity type ?? + float * sqc; // salinity qc + int * tmp_typ; // ?? temperature type ?? + float * tqc; // temperature qc + + // observations per level: + float ** lvl; // depth + float ** sal; // salinity (PSU) + float ** sal_err; // salinity error (std deviation, PSU) + float ** sprb; // ?? salinity ... ?? + float ** tmp; // in situ temperature (C) + float ** tmp_err; // in situ temperature error (C) + float ** tprb; // ?? temperature ... ?? + float ** clm_sal; // climatology of salinity + float ** cssd; // climatological std deviation for salinity + float ** clm_tmp; // climatology of temperature + float ** ctsd; // climatological std deviation for temperature + int ** flg; // ?? qc flags ?? + + char12 * dtg; // date (Julian, seconds) + // std::time_t * dtg; // date (Julian, seconds) + + private: + int n; + int mx_lvl; + int vrsn; + + void allocate(); + void allocate2d(); +}; // class RTOFSOb + + + +class RTOFSDataFile +{ + public: + explicit RTOFSDataFile(std::string filename); + int NumberOfObservations() { return nobs; } + RTOFSOb & observations() { return * ob; } + + + private: + std::string filename; + int nobs; + FILE * f; + RTOFSOb * ob; + + void read_file(); +}; // class RTOFSDataFile + +} // namespace rtofs + + namespace gdasapp { +int64_t + SecondsSinceReferenceTime(rtofs::RTOFSOb::char12 time) +{ + // parse the date + std::string s(time); + int year = std::stoi(s.substr(0, 4)); + int month = std::stoi(s.substr(4, 2)); + int day = std::stoi(s.substr(6, 2)); + int hour = std::stoi(s.substr(8, 2)); + int minute = std::stoi(s.substr(10, 2)); + int second = 0; + + uint64_t julianDate = + util::datefunctions::dateToJulian(year, month, day); + + // 2440588 = Julian date for January 1, 1970. + int daysSinceEpoch = julianDate - 2440588; + int secondsOffset = util::datefunctions::hmsToSeconds(hour, minute, second); + return static_cast(daysSinceEpoch * 86400.0f) + secondsOffset; +} // SecondsSinceReferenceTime + + + class RTOFSInSitu: public NetCDFToIodaConverter { @@ -32,7 +156,7 @@ class RTOFSInSitu: const eckit::mpi::Comm & comm): NetCDFToIodaConverter(fullConfig, comm) { - variable_ = "temperature"; + variable_ = "waterTemperature"; } // Read binary file and populate iodaVars @@ -49,12 +173,12 @@ class RTOFSInSitu: // Set the float metadata name std::vector floatMetadataNames = {}; - // read the file oops::Log::info() << "Processing file " << fileName << std::endl; + // read the file rtofs::RTOFSDataFile RTOFSFile(fileName); int n = RTOFSFile.NumberOfObservations(); rtofs::RTOFSOb & ob = RTOFSFile.observations(); @@ -69,35 +193,24 @@ class RTOFSInSitu: intMetadataNames); iodaVars.referenceDate_ = "seconds since 1970-01-01T00:00:00Z"; - // debugging output: - // std::string DataDirectory = - // "/home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/build/gdas-utils/test/obsproc/RTOFSdata"; - // ob.print(DataDirectory); - -// oops::Log::debug() << "--- iodaVars.location_: " << iodaVars.location_ << std::endl; - int k = 0; for (int i = 0; i < n; i ++) - for (int j = 0; j < ob.lt[i]; j ++) { - iodaVars.longitude_(k) = ob.lon[i]; - iodaVars.latitude_(k) = ob.lat[i]; - iodaVars.obsVal_(k) = ob.tmp[i][j]; - iodaVars.obsError_(k) = ob.tmp_err[i][j]; - iodaVars.datetime_(k) = ob.dtg[i]; - // iodaVars.preQc_(k) = oneDimFlagsVal[i]; - // iodaVars.intMetadata_.row(k) << -999; - - k++; - } + int64_t date = SecondsSinceReferenceTime(ob.dtg[i]); - // basic test for iodaVars.trim -/* - Eigen::Array mask - = (iodaVars.obsVal_ > 0.0); - iodaVars.trim(mask); -*/ + for (int j = 0; j < ob.lt[i]; j ++) + { + iodaVars.longitude_(k) = ob.lon[i]; + iodaVars.latitude_(k) = ob.lat[i]; + iodaVars.obsVal_(k) = ob.tmp[i][j]; + iodaVars.obsError_(k) = ob.tmp_err[i][j]; + iodaVars.datetime_(k) = date; + // iodaVars.preQc_(k) = oneDimFlagsVal[i]; + // iodaVars.intMetadata_.row(k) << -999; + k++; + } + } return iodaVars; }; // providerToIodaVars @@ -105,3 +218,385 @@ class RTOFSInSitu: } // namespace gdasapp + + + +namespace rtofs +{ + +// the variables marked "skipped" are arrays which are present +// in the binary file, but were skipped by the fortran code +// that was reading the files. + +RTOFSOb:: + RTOFSOb(int n, int mx_lvl, int vrsn): + n(n), + mx_lvl(mx_lvl), + vrsn(vrsn) +{ + allocate(); +} + + + +void +RTOFSOb:: + allocate() +{ + dtg = new char12[n]; + + lat = new float[n]; + lon = new float[n]; + btm = new float[n]; + ls = new int[n]; + lt = new int[n]; + sal_typ = new int[n]; + sqc = new float[n]; + tmp_typ = new int[n]; + tqc = new float[n]; + + lvl = new float * [n]; + sal = new float * [n]; + sal_err = new float * [n]; + sprb = new float * [n]; + tmp = new float * [n]; + tmp_err = new float * [n]; + clm_sal = new float * [n]; // skipped? + tprb = new float * [n]; + cssd = new float * [n]; + clm_tmp = new float * [n]; // skipped? + ctsd = new float * [n]; + flg = new int * [n]; +} + + +void +RTOFSOb:: + allocate2d() +{ + for (int i = 0; i < n; i ++) + { + int k = lt[i]; + + lvl[i] = new float[k]; + sal[i] = new float[k]; + sal_err[i] = new float[k]; + sprb[i] = new float[k]; + tmp[i] = new float[k]; + tmp_err[i] = new float[k]; + tprb[i] = new float[k]; + clm_sal[i] = new float[k]; + cssd[i] = new float[k]; + clm_tmp[i] = new float[k]; + ctsd[i] = new float[k]; + flg[i] = new int[k]; + } +} + + + +void +RTOFSOb:: + read(FILE * f) +{ + read_float_array(f, btm, n); + read_float_array(f, lat, n); + read_float_array(f, lon, n); + read_int_array(f, ls, n); + read_int_array(f, lt, n); + read_int_array(f, sal_typ, n); + read_float_array(f, sqc, n); + read_int_array(f, tmp_typ, n); + read_float_array(f, tqc, n); + + allocate2d(); + + for (int i = 0; i < n; i ++) + { + int k = lt[i]; + + read_float_array(f, lvl[i], k); + read_float_array(f, sal[i], k); + read_float_array(f, sal_err[i], k); + read_float_array(f, sprb[i], k); + read_float_array(f, tmp[i], k); + read_float_array(f, tmp_err[i], k); + read_float_array(f, tprb[i], k); + read_float_array(f, clm_sal[i], k); + read_float_array(f, cssd[i], k); + read_float_array(f, clm_tmp[i], k); + read_float_array(f, ctsd[i], k); + read_int_array(f, flg[i], k); + } +} + + +void +RTOFSOb:: + print(std::string DirectoryName) +{ + if (!file_exists(DirectoryName)) + { + cerr << "Directory " << DirectoryName << "doesn't exist" << endl; + exit(1); + } + print_float_array(DirectoryName + "/latitude", lat, n); + print_float_array(DirectoryName + "/longitude", lon, n); + print_float_array(DirectoryName + "/btm", btm, n); + print_float_array(DirectoryName + "/tqc", tqc, n); + print_float_array(DirectoryName + "/sqc", sqc, n); + print_int_array(DirectoryName + "/lt", lt, n); + print_int_array(DirectoryName + "/ls", ls, n); + print_int_array(DirectoryName + "/sal_typ", sal_typ, n); + print_int_array(DirectoryName + "/tmp_typ", sal_typ, n); + + print_2d_float_array(DirectoryName + "/tmp", tmp, n, lt); + print_2d_float_array(DirectoryName + "/sal", sal, n, lt); + + // print lvl2d array + { + ofstream o; + o.open(DirectoryName + "/lvl2d"); + for (int i = 0; i < n; i ++) + for (int j = 0; j < lt[i]; j ++) + o + << j + << endl; + o.close(); + } +} // RTOFSOb::print + + + +RTOFSDataFile:: +RTOFSDataFile(std::string filename): + filename(filename) +{ + if (!file_exists(filename)) + { + cerr << "File not found" << endl; + exit(1); + } + + + const char * fn = filename.c_str(); + f = fopen(fn, "rb"); + if (!f) + { + cerr << "Error opening file " << fn << endl; + exit(1); + } + + read_file(); +} // RTOFSDataFile::RTOFSDataFile + + +void +RTOFSDataFile:: + read_file() +{ + fseek(f, 4, SEEK_CUR); + + int n_read = read_int(f); + int mx_lvl = read_int(f); + int vrsn = read_int(f); + + ob = new RTOFSOb(n_read, mx_lvl, vrsn); + nobs = n_read; + + ob->read(f); + + skip8bytes(f); + + fread(ob->dtg, sizeof(RTOFSOb::char12), n_read, f); + + skip8bytes(f); + + RTOFSOb::char12 * ob_rct = new RTOFSOb::char12[n_read]; + fread(ob_rct, sizeof(RTOFSOb::char12), n_read, f); + + skip8bytes(f); + + RTOFSOb::char7 * ob_sgn = new RTOFSOb::char7[n_read]; + fread(ob_sgn, sizeof(RTOFSOb::char7), n_read, f); + + if (vrsn == 2) + { + float ** glb_sal = new float * [n_read]; + for (int i = 0; i < n_read; i ++) + { + int k = ob->lt[i]; + glb_sal[i] = alloc_read_float_array(f, k); + } + } +} // read_file + + + +void skip8bytes(FILE * f) +{ + fseek(f, 8, SEEK_CUR); +} + + +void + read_float_array(FILE * f, float * a, int n) +{ + skip8bytes(f); + + fread(a, sizeof(float), n, f); + + uint32_t * data = reinterpret_cast(a); + + for (int i = 0; i < n; i ++) + data[i] = be32toh(data[i]); +} // read_float_array + + +void + read_int_array(FILE * f, int * a, int n) +{ + skip8bytes(f); + + fread(a, sizeof(int), n, f); + + uint32_t * data = reinterpret_cast(a); + + for (int i = 0; i < n; i ++) + data[i] = be32toh(data[i]); +} // read_int_array + + + +void + print_level(std::string filename, float ** a, int n, int level) +{ + ofstream o(filename); + o + << std::fixed + << std::setprecision(9); + + for (int i = 0; i < n; i ++) + o + << a[i] [level] + << endl; + + o.close(); +} + + +void + print_int_array(std::string filename, int * a, int n) +{ + ofstream o(filename); + for (int i = 0; i < n; i ++) + o + << a[i] + << endl; + o.close(); +} + + +void + print_float_array(std::string filename, float * a, int n) +{ + ofstream o(filename); + o + << std::fixed + << std::setprecision(9); + for (int i = 0; i < n; i ++) + o + << a[i] + << endl; +} + +void + print_2d_float_array(std::string filename, float ** a, int n, int * dim2) +{ + ofstream o(filename); + o + << std::fixed + << std::setprecision(9); + for (int i = 0; i < n; i ++) + for (int j = 0; j < dim2[i]; j ++) + o + << a[i] [j] + << endl; +} + + +bool + file_exists(const std::string& name) +{ + struct stat buffer; + return (stat (name.c_str(), &buffer) == 0); +} + + +int + read_int(FILE * f) +{ + int dummy; + fread(&dummy, sizeof(int), 1, f); + dummy = static_cast(be32toh(dummy)); + return dummy; +} // read_int + + + +void + read_real_array(FILE * f, float * a, int n) +{ + fread(a, sizeof(float), n, f); + + uint32_t * data = reinterpret_cast(a); + + for (int i = 0; i < n; i ++) + data[i] = be32toh(data[i]); +} // read_real_array + + +float * + alloc_read_float_array(FILE * f, int n) +{ + float * a = new float[n]; + skip8bytes(f); + read_real_array(f, a, n); + return a; +} + + +void + print_int_array(int * a, int n) +{ + for (int i = 0; i < n; i ++) + cout + << a[i] + << endl; +} + + +int * + alloc_read_int_array(FILE * f, int n) +{ + int * a = new int[n]; + skip8bytes(f); + read_int_array(f, a, n); + return a; +} + + +void + print_real_array(float * a, int n) +{ + cout + << std::fixed + << std::setprecision(9); + for (int i = 0; i < n; i ++) + cout + << a[i] + << endl; +} + + +} // namespace rtofs diff --git a/utils/obsproc/applications/CMakeLists.txt b/utils/obsproc/applications/CMakeLists.txt index 6fd956b89..f03631855 100644 --- a/utils/obsproc/applications/CMakeLists.txt +++ b/utils/obsproc/applications/CMakeLists.txt @@ -10,7 +10,6 @@ target_compile_features( gdas_obsprovider2ioda.x PUBLIC cxx_std_17) target_link_libraries( gdas_obsprovider2ioda.x PUBLIC oops ioda NetCDF::NetCDF_CXX) -# message("CMAKE_SOURCE_DIR = ${CMAKE_SOURCE_DIR}") -link_directories(${CMAKE_SOURCE_DIR}/rtofs) -target_link_libraries( gdas_obsprovider2ioda.x PRIVATE rtofs) - +# to be used when rtofs is a static library: +# link_directories(${CMAKE_SOURCE_DIR}/rtofs) +# target_link_libraries( gdas_obsprovider2ioda.x PRIVATE rtofs) From bc4474d3b9c042cdf261f6594c6631960f1d4b70 Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Tue, 21 May 2024 13:01:34 -0500 Subject: [PATCH 07/18] removed the rtofs dir/lib --- utils/obsproc/rtofs/CMakeLists.txt | 38 ------ utils/obsproc/rtofs/README | 11 -- utils/obsproc/rtofs/RTOFSDataFile.cc | 114 ---------------- utils/obsproc/rtofs/RTOFSDataFile.h | 32 ----- utils/obsproc/rtofs/RTOFSOb.cc | 160 ---------------------- utils/obsproc/rtofs/RTOFSOb.h | 62 --------- utils/obsproc/rtofs/util.cc | 192 --------------------------- utils/obsproc/rtofs/util.h | 25 ---- 8 files changed, 634 deletions(-) delete mode 100644 utils/obsproc/rtofs/CMakeLists.txt delete mode 100644 utils/obsproc/rtofs/README delete mode 100644 utils/obsproc/rtofs/RTOFSDataFile.cc delete mode 100644 utils/obsproc/rtofs/RTOFSDataFile.h delete mode 100644 utils/obsproc/rtofs/RTOFSOb.cc delete mode 100644 utils/obsproc/rtofs/RTOFSOb.h delete mode 100644 utils/obsproc/rtofs/util.cc delete mode 100644 utils/obsproc/rtofs/util.h diff --git a/utils/obsproc/rtofs/CMakeLists.txt b/utils/obsproc/rtofs/CMakeLists.txt deleted file mode 100644 index 63488dc49..000000000 --- a/utils/obsproc/rtofs/CMakeLists.txt +++ /dev/null @@ -1,38 +0,0 @@ -set( - rtofs_src_files - RTOFSDataFile.h - RTOFSDataFile.cc - RTOFSOb.h - RTOFSOb.cc - util.cc - util.h -) - -list( - APPEND gdasapp_provider2ioda_src_files - ${rtofs_src_files} -) - - -add_library(rtofs STATIC - ${rtofs_src_files} -) - -target_compile_features( - rtofs - PUBLIC cxx_std_17 -) - -# target_link_libraries( - # rtofs - # PUBLIC oops ioda NetCDF::NetCDF_CXX -# ) - - -# ecbuild_add_executable( - # TARGET gdas_obsprovider2ioda.x - # SOURCES ${gdasapp_provider2ioda_src_files} -# ) - -# target_compile_features( gdas_obsprovider2ioda.x PUBLIC cxx_std_17) -# target_link_libraries( gdas_obsprovider2ioda.x PUBLIC oops ioda NetCDF::NetCDF_CXX) diff --git a/utils/obsproc/rtofs/README b/utils/obsproc/rtofs/README deleted file mode 100644 index c420242f2..000000000 --- a/utils/obsproc/rtofs/README +++ /dev/null @@ -1,11 +0,0 @@ -Author: E. Givelberg -May 15, 2024 - -The code in this rtofs directory reads a binary Fortran -file and generates a data object that can be ingested into ioda. - -The fortran binary file consists of records, including additional -information that needs to be skipped. -There are typically 8 bytes to be skipped between 2 data arrays. -The read data is converted from the big endian to the host -format, which is little endian in linux. diff --git a/utils/obsproc/rtofs/RTOFSDataFile.cc b/utils/obsproc/rtofs/RTOFSDataFile.cc deleted file mode 100644 index 79398d854..000000000 --- a/utils/obsproc/rtofs/RTOFSDataFile.cc +++ /dev/null @@ -1,114 +0,0 @@ -#include "RTOFSDataFile.h" - -#include - -#include -using std::cerr; -using std::endl; -using std::cout; -#include // std::get_time - -#include -#include - -#include "RTOFSOb.h" -#include "util.h" - - -namespace rtofs -{ - -typedef char char12[12]; -typedef char char7[7]; - - -std::time_t - SecondsSinceReferenceTime(char12 time) -{ - std::tm referenceTime = {}; - referenceTime.tm_year = 70; // 1970 - referenceTime.tm_mon = 0; // January (months are 0-based) - referenceTime.tm_mday = 1; // 1st day of the month - referenceTime.tm_hour = 0; - referenceTime.tm_min = 0; - referenceTime.tm_sec = 0; - std::time_t referenceTimestamp = std::mktime(&referenceTime); - - std::tm t = {}; - std::istringstream ss(time); - ss >> std::get_time(&t, "%Y%m%d%H%M"); - - std::time_t timestamp = std::mktime(&t); - return std::difftime(timestamp, referenceTimestamp); -} // SecondsSinceReferenceTime - - - -RTOFSDataFile:: -RTOFSDataFile(std::string filename): - filename(filename) -{ - if (!file_exists(filename)) - { - cerr << "File not found" << endl; - exit(1); - } - - - const char * fn = filename.c_str(); - f = fopen(fn, "rb"); - if (!f) - { - cerr << "Error opening file " << fn << endl; - exit(1); - } - - read_file(); -} // RTOFSDataFile::RTOFSDataFile - - -void -RTOFSDataFile:: - read_file() -{ - fseek(f, 4, SEEK_CUR); - - int n_read = read_int(f); - int mx_lvl = read_int(f); - int vrsn = read_int(f); - - ob = new RTOFSOb(n_read, mx_lvl, vrsn); - nobs = n_read; - - ob->read(f); - - skip8bytes(f); - - char12 * ob_dtg = new char12[n_read]; - fread(ob_dtg, sizeof(char[12]), n_read, f); - for (int i = 0; i < n_read; i ++) - ob->dtg[i] = SecondsSinceReferenceTime(ob_dtg[i]); - - skip8bytes(f); - - char12 * ob_rct = new char12[n_read]; - fread(ob_rct, sizeof(char[12]), n_read, f); - - skip8bytes(f); - - char7 * ob_sgn = new char7[n_read]; - fread(ob_sgn, sizeof(char[7]), n_read, f); - - if (vrsn == 2) - { - float ** glb_sal = new float * [n_read]; - for (int i = 0; i < n_read; i ++) - { - int k = ob->lt[i]; - glb_sal[i] = alloc_read_float_array(f, k); - } - } -} // read_file - - -} // namespace rtofs diff --git a/utils/obsproc/rtofs/RTOFSDataFile.h b/utils/obsproc/rtofs/RTOFSDataFile.h deleted file mode 100644 index 40a5b45bf..000000000 --- a/utils/obsproc/rtofs/RTOFSDataFile.h +++ /dev/null @@ -1,32 +0,0 @@ -#ifndef BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_RTOFSDATAFILE_H_ -#define BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_RTOFSDATAFILE_H_ - -#include - - -namespace rtofs -{ - -class RTOFSOb; - -class RTOFSDataFile -{ - public: - explicit RTOFSDataFile(std::string filename); - int NumberOfObservations() { return nobs; } - RTOFSOb & observations() { return * ob; } - - - private: - std::string filename; - int nobs; - FILE * f; - RTOFSOb * ob; - - void read_file(); -}; // class RTOFSDataFile - -} // namespace rtofs - - -#endif // BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_RTOFSDATAFILE_H_ diff --git a/utils/obsproc/rtofs/RTOFSOb.cc b/utils/obsproc/rtofs/RTOFSOb.cc deleted file mode 100644 index 66c057acf..000000000 --- a/utils/obsproc/rtofs/RTOFSOb.cc +++ /dev/null @@ -1,160 +0,0 @@ -#include "RTOFSOb.h" - -#include - -#include -using std::string; - -#include -using std::endl; -using std::cerr; - -#include -using std::ofstream; - -#include "util.h" - - -namespace rtofs -{ - - -RTOFSOb:: - RTOFSOb(int n, int mx_lvl, int vrsn): - n(n), - mx_lvl(mx_lvl), - vrsn(vrsn) -{ - dtg = new std::time_t[n]; - allocate(); -} - - - -void -RTOFSOb:: - allocate() -{ - lat = new float[n]; - lon = new float[n]; - btm = new float[n]; - ls = new int[n]; - lt = new int[n]; - sal_typ = new int[n]; - sqc = new float[n]; - tmp_typ = new int[n]; - tqc = new float[n]; - - lvl = new float * [n]; - sal = new float * [n]; - sal_err = new float * [n]; - sprb = new float * [n]; - tmp = new float * [n]; - tmp_err = new float * [n]; - clm_sal = new float * [n]; // skipped? - tprb = new float * [n]; - cssd = new float * [n]; - clm_tmp = new float * [n]; // skipped? - ctsd = new float * [n]; - flg = new int * [n]; -} - - -void -RTOFSOb:: - allocate2d() -{ - for (int i = 0; i < n; i ++) - { - int k = lt[i]; - - lvl[i] = new float[k]; - sal[i] = new float[k]; - sal_err[i] = new float[k]; - sprb[i] = new float[k]; - tmp[i] = new float[k]; - tmp_err[i] = new float[k]; - tprb[i] = new float[k]; - clm_sal[i] = new float[k]; - cssd[i] = new float[k]; - clm_tmp[i] = new float[k]; - ctsd[i] = new float[k]; - flg[i] = new int[k]; - } -} - - - -void -RTOFSOb:: - read(FILE * f) -{ - read_float_array(f, btm, n); - read_float_array(f, lat, n); - read_float_array(f, lon, n); - read_int_array(f, ls, n); - read_int_array(f, lt, n); - read_int_array(f, sal_typ, n); - read_float_array(f, sqc, n); - read_int_array(f, tmp_typ, n); - read_float_array(f, tqc, n); - - allocate2d(); - - for (int i = 0; i < n; i ++) - { - int k = lt[i]; - - read_float_array(f, lvl[i], k); - read_float_array(f, sal[i], k); - read_float_array(f, sal_err[i], k); - read_float_array(f, sprb[i], k); - read_float_array(f, tmp[i], k); - read_float_array(f, tmp_err[i], k); - read_float_array(f, tprb[i], k); - read_float_array(f, clm_sal[i], k); - read_float_array(f, cssd[i], k); - read_float_array(f, clm_tmp[i], k); - read_float_array(f, ctsd[i], k); - read_int_array(f, flg[i], k); - } -} - - -void -RTOFSOb:: - print(std::string DirectoryName) -{ - if (!file_exists(DirectoryName)) - { - cerr << "Directory " << DirectoryName << "doesn't exist" << endl; - exit(1); - } - print_float_array(DirectoryName + "/latitude", lat, n); - print_float_array(DirectoryName + "/longitude", lon, n); - print_float_array(DirectoryName + "/btm", btm, n); - print_float_array(DirectoryName + "/tqc", tqc, n); - print_float_array(DirectoryName + "/sqc", sqc, n); - print_int_array(DirectoryName + "/lt", lt, n); - print_int_array(DirectoryName + "/ls", ls, n); - print_int_array(DirectoryName + "/sal_typ", sal_typ, n); - print_int_array(DirectoryName + "/tmp_typ", sal_typ, n); - - print_2d_float_array(DirectoryName + "/tmp", tmp, n, lt); - print_2d_float_array(DirectoryName + "/sal", sal, n, lt); - - // print lvl2d array - { - ofstream o; - o.open(DirectoryName + "/lvl2d"); - for (int i = 0; i < n; i ++) - for (int j = 0; j < lt[i]; j ++) - o - << j - << endl; - o.close(); - } -} // RTOFSOb::print - - -} // namespace rtofs diff --git a/utils/obsproc/rtofs/RTOFSOb.h b/utils/obsproc/rtofs/RTOFSOb.h deleted file mode 100644 index fdd14df0c..000000000 --- a/utils/obsproc/rtofs/RTOFSOb.h +++ /dev/null @@ -1,62 +0,0 @@ -#ifndef BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_RTOFSOB_H_ -#define BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_RTOFSOB_H_ - - -#include -#include - -#include -using std::string; - -namespace rtofs -{ - - -class RTOFSOb -{ - public: - RTOFSOb(int n, int mx_lvl, int vrsn); - - void read(FILE * f); - int NumberOfObservations() { return n; } - void print(std::string DirectoryName); - - - float * btm; - float * lat; - float * lon; - int * ls; - int * lt; - int * sal_typ; - float * sqc; - int * tmp_typ; - float * tqc; - - float ** lvl; - float ** sal; - float ** sal_err; - float ** sprb; - float ** tmp; - float ** tmp_err; - float ** tprb; - float ** clm_sal; - float ** cssd; - float ** clm_tmp; - float ** ctsd; - int ** flg; - - std::time_t * dtg; - - private: - int n; - int mx_lvl; - int vrsn; - - void allocate(); - void allocate2d(); -}; // class RTOFSOb - - -} // namespace rtofs - -#endif // BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_RTOFSOB_H_ diff --git a/utils/obsproc/rtofs/util.cc b/utils/obsproc/rtofs/util.cc deleted file mode 100644 index f30de85e0..000000000 --- a/utils/obsproc/rtofs/util.cc +++ /dev/null @@ -1,192 +0,0 @@ -#include "util.h" - -// #define _BSD_SOURCE -#include -#include // used in file_exists - -#include -using std::cout; -using std::endl; - -#include -using std::ofstream; - -#include -using std::setprecision; - - -namespace rtofs -{ - - -void skip8bytes(FILE * f) -{ - fseek(f, 8, SEEK_CUR); -} - - -void - read_float_array(FILE * f, float * a, int n) -{ - skip8bytes(f); - - fread(a, sizeof(float), n, f); - - // uint32_t * data = (uint32_t *) a; - uint32_t * data = reinterpret_cast(a); - - for (int i = 0; i < n; i ++) - data[i] = be32toh(data[i]); -} // read_float_array - - -void - read_int_array(FILE * f, int * a, int n) -{ - skip8bytes(f); - - fread(a, sizeof(int), n, f); - - // uint32_t * data = (uint32_t *) a; - uint32_t * data = reinterpret_cast(a); - - for (int i = 0; i < n; i ++) - data[i] = be32toh(data[i]); -} // read_int_array - - - -void - print_level(std::string filename, float ** a, int n, int level) -{ - ofstream o(filename); - o - << std::fixed - << std::setprecision(9); - - for (int i = 0; i < n; i ++) - o - << a[i] [level] - << endl; - - o.close(); -} - - -void - print_int_array(std::string filename, int * a, int n) -{ - ofstream o(filename); - for (int i = 0; i < n; i ++) - o - << a[i] - << endl; - o.close(); -} - - -void - print_float_array(std::string filename, float * a, int n) -{ - ofstream o(filename); - o - << std::fixed - << std::setprecision(9); - for (int i = 0; i < n; i ++) - o - << a[i] - << endl; -} - -void - print_2d_float_array(std::string filename, float ** a, int n, int * dim2) -{ - ofstream o(filename); - o - << std::fixed - << std::setprecision(9); - for (int i = 0; i < n; i ++) - for (int j = 0; j < dim2[i]; j ++) - o - << a[i] [j] - << endl; -} - - -bool - file_exists(const std::string& name) -{ - struct stat buffer; - return (stat (name.c_str(), &buffer) == 0); -} - - -int - read_int(FILE * f) -{ - int dummy; - fread(&dummy, sizeof(int), 1, f); - // dummy = (int) be32toh(dummy); - dummy = static_cast(be32toh(dummy)); - return dummy; -} // read_int - - - -void - read_real_array(FILE * f, float * a, int n) -{ - fread(a, sizeof(float), n, f); - - // uint32_t * data = (uint32_t *) a; - uint32_t * data = reinterpret_cast(a); - - for (int i = 0; i < n; i ++) - data[i] = be32toh(data[i]); -} // read_real_array - - -float * - alloc_read_float_array(FILE * f, int n) -{ - float * a = new float[n]; - skip8bytes(f); - read_real_array(f, a, n); - return a; -} - - -void - print_int_array(int * a, int n) -{ - for (int i = 0; i < n; i ++) - cout - << a[i] - << endl; -} - - -int * - alloc_read_int_array(FILE * f, int n) -{ - int * a = new int[n]; - skip8bytes(f); - read_int_array(f, a, n); - return a; -} - - -void - print_real_array(float * a, int n) -{ - cout - << std::fixed - << std::setprecision(9); - for (int i = 0; i < n; i ++) - cout - << a[i] - << endl; -} - - -} // namespace rtofs diff --git a/utils/obsproc/rtofs/util.h b/utils/obsproc/rtofs/util.h deleted file mode 100644 index 6dfae46f7..000000000 --- a/utils/obsproc/rtofs/util.h +++ /dev/null @@ -1,25 +0,0 @@ -#ifndef BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_UTIL_H_ -#define BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_UTIL_H_ - -#include -#include - -namespace rtofs -{ - -bool file_exists(const std::string& name); - -void skip8bytes(FILE * f); -int read_int(FILE * f); -void read_float_array(FILE * f, float * a, int n); -void read_int_array(FILE * f, int * a, int n); - -void print_float_array(std::string filename, float * a, int n); -void print_int_array(std::string filename, int * a, int n); -void print_2d_float_array(std::string filename, float ** a, int n, int * dim2); -float * alloc_read_float_array(FILE * f, int n); - -} // namespace rtofs - - -#endif // BUNDLE_GDAS_UTILS_OBSPROC_RTOFS_UTIL_H_ From 24a188285c4ed6bd3a85d741895eea54312e0dd8 Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Tue, 21 May 2024 14:44:23 -0500 Subject: [PATCH 08/18] removed the test (cause it uses binary file) --- utils/test/CMakeLists.txt | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/utils/test/CMakeLists.txt b/utils/test/CMakeLists.txt index 4861e8f45..2d912c2d2 100644 --- a/utils/test/CMakeLists.txt +++ b/utils/test/CMakeLists.txt @@ -65,12 +65,14 @@ ecbuild_add_test( TARGET test_gdasapp_util_ghrsst2ioda LIBS gdas-utils WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc) +#TODO: +# add binary files for the test # Test the RTOFS in Situ converter -ecbuild_add_test( TARGET test_gdasapp_util_rtofsinsitu - COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x - ARGS "../testinput/gdas_rtofsinsitu.yaml" - LIBS gdas-utils - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc) +# ecbuild_add_test( TARGET test_gdasapp_util_rtofsinsitu + # COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x + # ARGS "../testinput/gdas_rtofsinsitu.yaml" + # LIBS gdas-utils + # WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc) # Test the SMAP to IODA converter ecbuild_add_test( TARGET test_gdasapp_util_smap2ioda From 18f4e816ec1018a9c8b78af6185daa6fd928b117 Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Fri, 31 May 2024 16:17:40 -0500 Subject: [PATCH 09/18] two converters for tmp, sal --- utils/obsproc/CMakeLists.txt | 1 - utils/obsproc/RTOFSInSitu.h | 276 ++++++++++++++---- utils/obsproc/RTOFSSalinity.h | 35 +++ utils/obsproc/RTOFSTemperature.h | 34 +++ .../applications/gdas_obsprovider2ioda.h | 13 +- 5 files changed, 294 insertions(+), 65 deletions(-) create mode 100644 utils/obsproc/RTOFSSalinity.h create mode 100644 utils/obsproc/RTOFSTemperature.h diff --git a/utils/obsproc/CMakeLists.txt b/utils/obsproc/CMakeLists.txt index 9cce213e1..4e4e3feb6 100644 --- a/utils/obsproc/CMakeLists.txt +++ b/utils/obsproc/CMakeLists.txt @@ -1,2 +1 @@ -# add_subdirectory(rtofs) add_subdirectory(applications) diff --git a/utils/obsproc/RTOFSInSitu.h b/utils/obsproc/RTOFSInSitu.h index 32b9b967f..afd09deaf 100644 --- a/utils/obsproc/RTOFSInSitu.h +++ b/utils/obsproc/RTOFSInSitu.h @@ -19,6 +19,7 @@ using std::cout; using std::string; #include #include +using std::vector; #include "NetCDFToIodaConverter.h" #include "oops/util/dateFunctions.h" @@ -29,6 +30,9 @@ using std::string; +namespace gdasapp +{ + namespace rtofs { @@ -43,8 +47,12 @@ void read_int_array(FILE * f, int * a, int n); void print_float_array(std::string filename, float * a, int n); void print_int_array(std::string filename, int * a, int n); void print_2d_float_array(std::string filename, float ** a, int n, int * dim2); +template +void print_2d_array(std::string filename, NUMBER ** a, int n, int * dim2); +void print_level(std::string filename, float ** a, int n, int level, int * dim2); float * alloc_read_float_array(FILE * f, int n); +std::vector find_unique_values(int * instrument_list, int n); class RTOFSOb @@ -57,6 +65,7 @@ class RTOFSOb void read(FILE * f); int NumberOfObservations() { return n; } + int TotalNumberOfValues(); void print(std::string DirectoryName); @@ -85,7 +94,6 @@ class RTOFSOb int ** flg; // ?? qc flags ?? char12 * dtg; // date (Julian, seconds) - // std::time_t * dtg; // date (Julian, seconds) private: int n; @@ -97,12 +105,22 @@ class RTOFSOb }; // class RTOFSOb +int +RTOFSOb:: + TotalNumberOfValues() +{ + int NumberOfValues = 0; + for (int i = 0; i < n; i ++) + NumberOfValues += lt[i]; + return NumberOfValues; +} + + class RTOFSDataFile { public: explicit RTOFSDataFile(std::string filename); - int NumberOfObservations() { return nobs; } RTOFSOb & observations() { return * ob; } @@ -115,13 +133,6 @@ class RTOFSDataFile void read_file(); }; // class RTOFSDataFile -} // namespace rtofs - - - - -namespace gdasapp -{ int64_t @@ -146,7 +157,6 @@ int64_t } // SecondsSinceReferenceTime - class RTOFSInSitu: public NetCDFToIodaConverter { @@ -155,74 +165,173 @@ class RTOFSInSitu: const eckit::Configuration & fullConfig, const eckit::mpi::Comm & comm): NetCDFToIodaConverter(fullConfig, comm) - { - variable_ = "waterTemperature"; - } + {} + + + protected: + void ProcessFile(std::string filename); + gdasapp::obsproc::iodavars::IodaVars TemperatureIodaVars(); + gdasapp::obsproc::iodavars::IodaVars SalinityIodaVars(); + + private: // Read binary file and populate iodaVars - gdasapp::obsproc::iodavars::IodaVars - providerToIodaVars(const std::string fileName) final + virtual gdasapp::obsproc::iodavars::IodaVars + providerToIodaVars(const std::string filename) = 0; + + rtofs::RTOFSOb * pob; +}; // class RTOFSInSitu + + +void +RTOFSInSitu:: + ProcessFile(std::string filename) +{ + oops::Log::info() + << "Processing RTOFS file " + << filename + << std::endl; + + // read the file + rtofs::RTOFSDataFile RTOFSFile(filename); + pob = & RTOFSFile.observations(); + rtofs::RTOFSOb & ob = * pob; + + + int n = ob.NumberOfObservations(); + + std::vector tmp_instrument = + rtofs::find_unique_values(ob.tmp_typ, n); + + std::vector sal_instrument = + rtofs::find_unique_values(ob.sal_typ, n); +} // RTOFSInSitu::ProcessFile + + + +gdasapp::obsproc::iodavars::IodaVars +RTOFSInSitu:: + TemperatureIodaVars() +{ + // Set the int metadata names + std::vector intMetadataNames = {"tmp_typ"}; + + // Set the float metadata name + std::vector floatMetadataNames = {"level"}; + rtofs::RTOFSOb & ob = * pob; + + gdasapp::obsproc::iodavars::IodaVars iodaVars( + ob.TotalNumberOfValues(), + floatMetadataNames, + intMetadataNames); + + iodaVars.referenceDate_ = "seconds since 1970-01-01T00:00:00Z"; + iodaVars.niMetadata_ = 1; + iodaVars.nfMetadata_ = 1; + + int n = ob.NumberOfObservations(); + + int k = 0; + for (int i = 0; i < n; i ++) { - oops::Log::info() - << "Processing files provided by RTOFS" - << std::endl; + int64_t date = SecondsSinceReferenceTime(ob.dtg[i]); - // Set the int metadata names - std::vector intMetadataNames = {"temperature"}; + for (int j = 0; j < ob.lt[i]; j ++) + { + iodaVars.longitude_(k) = ob.lon[i]; + iodaVars.latitude_(k) = ob.lat[i]; + iodaVars.obsVal_(k) = ob.tmp[i][j]; + iodaVars.obsError_(k) = ob.tmp_err[i][j]; + iodaVars.datetime_(k) = date; + iodaVars.preQc_(k) = ob.flg[i][j]; + iodaVars.intMetadata_.row(k) << ob.tmp_typ[i]; + iodaVars.floatMetadata_.row(k) << ob.lvl[i][j]; + + k++; + } + } - // Set the float metadata name - std::vector floatMetadataNames = {}; + return iodaVars; +} // RTOFSInSitu::TemperatureIodaVars - oops::Log::info() - << "Processing file " - << fileName - << std::endl; - // read the file - rtofs::RTOFSDataFile RTOFSFile(fileName); - int n = RTOFSFile.NumberOfObservations(); - rtofs::RTOFSOb & ob = RTOFSFile.observations(); - int NumberOfTemperatureValues = 0; - for (int i = 0; i < n; i ++) - NumberOfTemperatureValues += ob.lt[i]; - gdasapp::obsproc::iodavars::IodaVars iodaVars( - NumberOfTemperatureValues, - floatMetadataNames, - intMetadataNames); - iodaVars.referenceDate_ = "seconds since 1970-01-01T00:00:00Z"; - int k = 0; - for (int i = 0; i < n; i ++) + +gdasapp::obsproc::iodavars::IodaVars +RTOFSInSitu:: + SalinityIodaVars() +{ + // Set the int metadata names + std::vector intMetadataNames = {"sal_typ"}; + + // Set the float metadata name + std::vector floatMetadataNames = {"level"}; + + rtofs::RTOFSOb & ob = * pob; + + gdasapp::obsproc::iodavars::IodaVars iodaVars( + ob.TotalNumberOfValues(), + floatMetadataNames, + intMetadataNames); + + iodaVars.referenceDate_ = "seconds since 1970-01-01T00:00:00Z"; + iodaVars.niMetadata_ = 1; + iodaVars.nfMetadata_ = 1; + + int n = ob.NumberOfObservations(); + + int k = 0; + for (int i = 0; i < n; i ++) + { + int64_t date = SecondsSinceReferenceTime(ob.dtg[i]); + + for (int j = 0; j < ob.lt[i]; j ++) { - int64_t date = SecondsSinceReferenceTime(ob.dtg[i]); - - for (int j = 0; j < ob.lt[i]; j ++) - { - iodaVars.longitude_(k) = ob.lon[i]; - iodaVars.latitude_(k) = ob.lat[i]; - iodaVars.obsVal_(k) = ob.tmp[i][j]; - iodaVars.obsError_(k) = ob.tmp_err[i][j]; - iodaVars.datetime_(k) = date; - // iodaVars.preQc_(k) = oneDimFlagsVal[i]; - // iodaVars.intMetadata_.row(k) << -999; - - k++; - } + iodaVars.longitude_(k) = ob.lon[i]; + iodaVars.latitude_(k) = ob.lat[i]; + iodaVars.obsVal_(k) = ob.sal[i][j]; + iodaVars.obsError_(k) = ob.sal_err[i][j]; + iodaVars.datetime_(k) = date; + iodaVars.preQc_(k) = ob.flg[i][j]; + iodaVars.intMetadata_.row(k) << ob.sal_typ[i]; + iodaVars.floatMetadata_.row(k) << ob.lvl[i][j]; + + k++; } + } - return iodaVars; - }; // providerToIodaVars -}; // class RTOFSInSitu + return iodaVars; +} // RTOFSInSitu::SalinityIodaVars -} // namespace gdasapp -namespace rtofs + + +std::vector + find_unique_values(int * instrument_list, int n) { + std::vector instrument; + + // Function to insert a value while keeping + // the vector sorted and without duplicates + auto insertUniqueSorted = [&](int value) + { + auto it = std::lower_bound(instrument.begin(), instrument.end(), value); + if (it == instrument.end() || *it != value) + instrument.insert(it, value); + }; + + for (int i = 0; i < n; i ++) + insertUniqueSorted(instrument_list[i]); + + return instrument; +} + + // the variables marked "skipped" are arrays which are present // in the binary file, but were skipped by the fortran code @@ -348,10 +457,24 @@ RTOFSOb:: print_int_array(DirectoryName + "/lt", lt, n); print_int_array(DirectoryName + "/ls", ls, n); print_int_array(DirectoryName + "/sal_typ", sal_typ, n); - print_int_array(DirectoryName + "/tmp_typ", sal_typ, n); + print_int_array(DirectoryName + "/tmp_typ", tmp_typ, n); print_2d_float_array(DirectoryName + "/tmp", tmp, n, lt); print_2d_float_array(DirectoryName + "/sal", sal, n, lt); + print_2d_float_array(DirectoryName + "/lvl", lvl, n, lt); + print_2d_float_array(DirectoryName + "/cssd", cssd, n, lt); + print_2d_float_array(DirectoryName + "/ctsd", ctsd, n, lt); + print_2d_float_array(DirectoryName + "/tprb", tprb, n, lt); + print_2d_float_array(DirectoryName + "/sprb", sprb, n, lt); + print_2d_float_array(DirectoryName + "/clm_tmp", clm_tmp, n, lt); + + print_level(DirectoryName + "/tmp0", tmp, n, 0, lt); + print_level(DirectoryName + "/clm_tmp0", clm_tmp, n, 0, lt); + print_level(DirectoryName + "/tmp10", tmp, n, 10, lt); + print_level(DirectoryName + "/clm_tmp10", clm_tmp, n, 10, lt); + print_level(DirectoryName + "/tprb0", tprb, n, 0, lt); + + print_2d_array(DirectoryName + "/flg", flg, n, lt); // print lvl2d array { @@ -420,6 +543,18 @@ RTOFSDataFile:: RTOFSOb::char7 * ob_sgn = new RTOFSOb::char7[n_read]; fread(ob_sgn, sizeof(RTOFSOb::char7), n_read, f); + // we don't know what sgn is. + { + ofstream o("sgn"); + for (int i = 0; i < n_read; i ++) + { + for (int k = 0; k < 7; k ++) + o << ob_sgn[i][k]; + o << std::endl; + } + } + + // we never had an input file with vrsn == 2 if (vrsn == 2) { float ** glb_sal = new float * [n_read]; @@ -469,7 +604,7 @@ void void - print_level(std::string filename, float ** a, int n, int level) + print_level(std::string filename, float ** a, int n, int level, int * dim2) { ofstream o(filename); o @@ -477,6 +612,7 @@ void << std::setprecision(9); for (int i = 0; i < n; i ++) + if (level < dim2[i]) o << a[i] [level] << endl; @@ -525,6 +661,22 @@ void } +template +void + print_2d_array(std::string filename, NUMBER ** a, int n, int * dim2) +{ + ofstream o(filename); + o + << std::fixed + << std::setprecision(9); + for (int i = 0; i < n; i ++) + for (int j = 0; j < dim2[i]; j ++) + o + << a[i] [j] + << endl; +} + + bool file_exists(const std::string& name) { @@ -600,3 +752,5 @@ void } // namespace rtofs + +} // namespace gdasapp diff --git a/utils/obsproc/RTOFSSalinity.h b/utils/obsproc/RTOFSSalinity.h new file mode 100644 index 000000000..4789dfa2f --- /dev/null +++ b/utils/obsproc/RTOFSSalinity.h @@ -0,0 +1,35 @@ +#pragma once + + +#include + +#include "RTOFSInSitu.h" + + +namespace gdasapp +{ + + +class RTOFSSalinity: + public rtofs::RTOFSInSitu +{ + public: + explicit RTOFSSalinity( + const eckit::Configuration & fullConfig, + const eckit::mpi::Comm & comm): + RTOFSInSitu(fullConfig, comm) + { + variable_ = "waterSalinity"; + } + + // Read binary file and populate iodaVars + gdasapp::obsproc::iodavars::IodaVars + providerToIodaVars(const std::string filename) final + { + ProcessFile(filename); + return SalinityIodaVars(); + } +}; // class RTOFSSalinity + + +} // namespace gdasapp diff --git a/utils/obsproc/RTOFSTemperature.h b/utils/obsproc/RTOFSTemperature.h new file mode 100644 index 000000000..add6a7a48 --- /dev/null +++ b/utils/obsproc/RTOFSTemperature.h @@ -0,0 +1,34 @@ +#pragma once + +#include + +#include "RTOFSInSitu.h" + + +namespace gdasapp +{ + + +class RTOFSTemperature: + public rtofs::RTOFSInSitu +{ + public: + explicit RTOFSTemperature( + const eckit::Configuration & fullConfig, + const eckit::mpi::Comm & comm): + RTOFSInSitu(fullConfig, comm) + { + variable_ = "waterTemperature"; + } + + // Read binary file and populate iodaVars + gdasapp::obsproc::iodavars::IodaVars + providerToIodaVars(const std::string filename) final + { + ProcessFile(filename); + return TemperatureIodaVars(); + } +}; // class RTOFSTemperature + + +} // namespace gdasapp diff --git a/utils/obsproc/applications/gdas_obsprovider2ioda.h b/utils/obsproc/applications/gdas_obsprovider2ioda.h index b3bcd5c89..eb53fce0b 100644 --- a/utils/obsproc/applications/gdas_obsprovider2ioda.h +++ b/utils/obsproc/applications/gdas_obsprovider2ioda.h @@ -9,7 +9,8 @@ #include "../Ghrsst2Ioda.h" #include "../IcecAmsr2Ioda.h" #include "../Rads2Ioda.h" -#include "../RTOFSInSitu.h" +#include "../RTOFSSalinity.h" +#include "../RTOFSTemperature.h" #include "../Smap2Ioda.h" #include "../Smos2Ioda.h" #include "../Viirsaod2Ioda.h" @@ -19,8 +20,10 @@ namespace gdasapp { public: explicit ObsProvider2IodaApp(const eckit::mpi::Comm & comm = oops::mpi::world()) : Application(comm) {} + static const std::string classname() {return "gdasapp::ObsProvider2IodaApp";} + int execute(const eckit::Configuration & fullConfig, bool /*validate*/) const { // Get the file provider string identifier from the config std::string provider; @@ -32,8 +35,12 @@ namespace gdasapp { } else if (provider == "GHRSST") { Ghrsst2Ioda conv2ioda(fullConfig, this->getComm()); conv2ioda.writeToIoda(); - } else if (provider == "RTOFS") { - RTOFSInSitu conv2ioda(fullConfig, this->getComm()); + } else if (provider == "RTOFStmp") { + // RTOFSInSitu conv2ioda(fullConfig, this->getComm()); + RTOFSTemperature conv2ioda(fullConfig, this->getComm()); + conv2ioda.writeToIoda(); + } else if (provider == "RTOFSsal") { + RTOFSSalinity conv2ioda(fullConfig, this->getComm()); conv2ioda.writeToIoda(); } else if (provider == "SMAP") { Smap2Ioda conv2ioda(fullConfig, this->getComm()); From 4f75398db42cf9c7f1f71419a150a046f77d135e Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Wed, 5 Jun 2024 19:21:31 -0500 Subject: [PATCH 10/18] rtofs directory; lint exclude --- utils/obsproc/CMakeLists.txt | 1 + utils/obsproc/RTOFSInSitu.h | 698 +-------------------------- utils/obsproc/RTOFSSalinity.h | 56 ++- utils/obsproc/RTOFSTemperature.h | 53 +- utils/obsproc/rtofs/CMakeLists.txt | 44 ++ utils/obsproc/rtofs/README | 14 + utils/obsproc/rtofs/RTOFSDataFile.cc | 95 ++++ utils/obsproc/rtofs/RTOFSDataFile.h | 29 ++ utils/obsproc/rtofs/RTOFSOb.cc | 218 +++++++++ utils/obsproc/rtofs/RTOFSOb.h | 66 +++ utils/obsproc/rtofs/j | 33 ++ utils/obsproc/rtofs/util.cc | 216 +++++++++ utils/obsproc/rtofs/util.h | 61 +++ utils/test/CMakeLists.txt | 23 +- 14 files changed, 906 insertions(+), 701 deletions(-) create mode 100644 utils/obsproc/rtofs/CMakeLists.txt create mode 100644 utils/obsproc/rtofs/README create mode 100644 utils/obsproc/rtofs/RTOFSDataFile.cc create mode 100644 utils/obsproc/rtofs/RTOFSDataFile.h create mode 100644 utils/obsproc/rtofs/RTOFSOb.cc create mode 100644 utils/obsproc/rtofs/RTOFSOb.h create mode 100644 utils/obsproc/rtofs/j create mode 100644 utils/obsproc/rtofs/util.cc create mode 100644 utils/obsproc/rtofs/util.h diff --git a/utils/obsproc/CMakeLists.txt b/utils/obsproc/CMakeLists.txt index 4e4e3feb6..c092a8758 100644 --- a/utils/obsproc/CMakeLists.txt +++ b/utils/obsproc/CMakeLists.txt @@ -1 +1,2 @@ +add_subdirectory(rtofs) add_subdirectory(applications) diff --git a/utils/obsproc/RTOFSInSitu.h b/utils/obsproc/RTOFSInSitu.h index afd09deaf..ccbd44eb7 100644 --- a/utils/obsproc/RTOFSInSitu.h +++ b/utils/obsproc/RTOFSInSitu.h @@ -1,32 +1,12 @@ #pragma once -#include -// #define _BSD_SOURCE -#include -#include // used in file_exists - -#include -#include -using std::ofstream; -#include // std::get_time -using std::setprecision; -#include -using std::cerr; -using std::endl; -using std::cout; -#include // NOLINT (using C API) #include -using std::string; -#include #include -using std::vector; #include "NetCDFToIodaConverter.h" -#include "oops/util/dateFunctions.h" -#include "oops/util/DateTime.h" - - - +#include "rtofs/RTOFSDataFile.h" +#include "rtofs/RTOFSOb.h" +#include "rtofs/util.h" @@ -37,126 +17,6 @@ namespace gdasapp namespace rtofs { -bool file_exists(const std::string& name); - -void skip8bytes(FILE * f); -int read_int(FILE * f); -void read_float_array(FILE * f, float * a, int n); -void read_int_array(FILE * f, int * a, int n); - -void print_float_array(std::string filename, float * a, int n); -void print_int_array(std::string filename, int * a, int n); -void print_2d_float_array(std::string filename, float ** a, int n, int * dim2); -template -void print_2d_array(std::string filename, NUMBER ** a, int n, int * dim2); -void print_level(std::string filename, float ** a, int n, int level, int * dim2); -float * alloc_read_float_array(FILE * f, int n); - -std::vector find_unique_values(int * instrument_list, int n); - - -class RTOFSOb -{ - public: - RTOFSOb(int n, int mx_lvl, int vrsn); - - typedef char char12[12]; - typedef char char7[7]; - - void read(FILE * f); - int NumberOfObservations() { return n; } - int TotalNumberOfValues(); - void print(std::string DirectoryName); - - - float * btm; // bottom depth - float * lat; // latitude - float * lon; // longitude - int * ls; // ?? - int * lt; // array of dimensions, the number of levels - int * sal_typ; // ?? salinity type ?? - float * sqc; // salinity qc - int * tmp_typ; // ?? temperature type ?? - float * tqc; // temperature qc - - // observations per level: - float ** lvl; // depth - float ** sal; // salinity (PSU) - float ** sal_err; // salinity error (std deviation, PSU) - float ** sprb; // ?? salinity ... ?? - float ** tmp; // in situ temperature (C) - float ** tmp_err; // in situ temperature error (C) - float ** tprb; // ?? temperature ... ?? - float ** clm_sal; // climatology of salinity - float ** cssd; // climatological std deviation for salinity - float ** clm_tmp; // climatology of temperature - float ** ctsd; // climatological std deviation for temperature - int ** flg; // ?? qc flags ?? - - char12 * dtg; // date (Julian, seconds) - - private: - int n; - int mx_lvl; - int vrsn; - - void allocate(); - void allocate2d(); -}; // class RTOFSOb - - -int -RTOFSOb:: - TotalNumberOfValues() -{ - int NumberOfValues = 0; - for (int i = 0; i < n; i ++) - NumberOfValues += lt[i]; - return NumberOfValues; -} - - - -class RTOFSDataFile -{ - public: - explicit RTOFSDataFile(std::string filename); - RTOFSOb & observations() { return * ob; } - - - private: - std::string filename; - int nobs; - FILE * f; - RTOFSOb * ob; - - void read_file(); -}; // class RTOFSDataFile - - - -int64_t - SecondsSinceReferenceTime(rtofs::RTOFSOb::char12 time) -{ - // parse the date - std::string s(time); - int year = std::stoi(s.substr(0, 4)); - int month = std::stoi(s.substr(4, 2)); - int day = std::stoi(s.substr(6, 2)); - int hour = std::stoi(s.substr(8, 2)); - int minute = std::stoi(s.substr(10, 2)); - int second = 0; - - uint64_t julianDate = - util::datefunctions::dateToJulian(year, month, day); - - // 2440588 = Julian date for January 1, 1970. - int daysSinceEpoch = julianDate - 2440588; - int secondsOffset = util::datefunctions::hmsToSeconds(hour, minute, second); - return static_cast(daysSinceEpoch * 86400.0f) + secondsOffset; -} // SecondsSinceReferenceTime - - class RTOFSInSitu: public NetCDFToIodaConverter { @@ -171,18 +31,17 @@ class RTOFSInSitu: protected: void ProcessFile(std::string filename); - gdasapp::obsproc::iodavars::IodaVars TemperatureIodaVars(); - gdasapp::obsproc::iodavars::IodaVars SalinityIodaVars(); + rtofs::RTOFSOb * pob; private: // Read binary file and populate iodaVars virtual gdasapp::obsproc::iodavars::IodaVars providerToIodaVars(const std::string filename) = 0; - - rtofs::RTOFSOb * pob; }; // class RTOFSInSitu +// read the binary file and possibly do additional filtering +// to prepare an rtofs::RTOFSOb to be converted to iodavars object void RTOFSInSitu:: ProcessFile(std::string filename) @@ -197,9 +56,9 @@ RTOFSInSitu:: pob = & RTOFSFile.observations(); rtofs::RTOFSOb & ob = * pob; - int n = ob.NumberOfObservations(); + // to be used with selection by instrument: std::vector tmp_instrument = rtofs::find_unique_values(ob.tmp_typ, n); @@ -208,549 +67,6 @@ RTOFSInSitu:: } // RTOFSInSitu::ProcessFile - -gdasapp::obsproc::iodavars::IodaVars -RTOFSInSitu:: - TemperatureIodaVars() -{ - // Set the int metadata names - std::vector intMetadataNames = {"tmp_typ"}; - - // Set the float metadata name - std::vector floatMetadataNames = {"level"}; - rtofs::RTOFSOb & ob = * pob; - - gdasapp::obsproc::iodavars::IodaVars iodaVars( - ob.TotalNumberOfValues(), - floatMetadataNames, - intMetadataNames); - - iodaVars.referenceDate_ = "seconds since 1970-01-01T00:00:00Z"; - iodaVars.niMetadata_ = 1; - iodaVars.nfMetadata_ = 1; - - int n = ob.NumberOfObservations(); - - int k = 0; - for (int i = 0; i < n; i ++) - { - int64_t date = SecondsSinceReferenceTime(ob.dtg[i]); - - for (int j = 0; j < ob.lt[i]; j ++) - { - iodaVars.longitude_(k) = ob.lon[i]; - iodaVars.latitude_(k) = ob.lat[i]; - iodaVars.obsVal_(k) = ob.tmp[i][j]; - iodaVars.obsError_(k) = ob.tmp_err[i][j]; - iodaVars.datetime_(k) = date; - iodaVars.preQc_(k) = ob.flg[i][j]; - iodaVars.intMetadata_.row(k) << ob.tmp_typ[i]; - iodaVars.floatMetadata_.row(k) << ob.lvl[i][j]; - - k++; - } - } - - return iodaVars; -} // RTOFSInSitu::TemperatureIodaVars - - - - - - -gdasapp::obsproc::iodavars::IodaVars -RTOFSInSitu:: - SalinityIodaVars() -{ - // Set the int metadata names - std::vector intMetadataNames = {"sal_typ"}; - - // Set the float metadata name - std::vector floatMetadataNames = {"level"}; - - rtofs::RTOFSOb & ob = * pob; - - gdasapp::obsproc::iodavars::IodaVars iodaVars( - ob.TotalNumberOfValues(), - floatMetadataNames, - intMetadataNames); - - iodaVars.referenceDate_ = "seconds since 1970-01-01T00:00:00Z"; - iodaVars.niMetadata_ = 1; - iodaVars.nfMetadata_ = 1; - - int n = ob.NumberOfObservations(); - - int k = 0; - for (int i = 0; i < n; i ++) - { - int64_t date = SecondsSinceReferenceTime(ob.dtg[i]); - - for (int j = 0; j < ob.lt[i]; j ++) - { - iodaVars.longitude_(k) = ob.lon[i]; - iodaVars.latitude_(k) = ob.lat[i]; - iodaVars.obsVal_(k) = ob.sal[i][j]; - iodaVars.obsError_(k) = ob.sal_err[i][j]; - iodaVars.datetime_(k) = date; - iodaVars.preQc_(k) = ob.flg[i][j]; - iodaVars.intMetadata_.row(k) << ob.sal_typ[i]; - iodaVars.floatMetadata_.row(k) << ob.lvl[i][j]; - - k++; - } - } - - return iodaVars; -} // RTOFSInSitu::SalinityIodaVars - - - - - - - -std::vector - find_unique_values(int * instrument_list, int n) -{ - std::vector instrument; - - // Function to insert a value while keeping - // the vector sorted and without duplicates - auto insertUniqueSorted = [&](int value) - { - auto it = std::lower_bound(instrument.begin(), instrument.end(), value); - if (it == instrument.end() || *it != value) - instrument.insert(it, value); - }; - - for (int i = 0; i < n; i ++) - insertUniqueSorted(instrument_list[i]); - - return instrument; -} - - - -// the variables marked "skipped" are arrays which are present -// in the binary file, but were skipped by the fortran code -// that was reading the files. - -RTOFSOb:: - RTOFSOb(int n, int mx_lvl, int vrsn): - n(n), - mx_lvl(mx_lvl), - vrsn(vrsn) -{ - allocate(); -} - - - -void -RTOFSOb:: - allocate() -{ - dtg = new char12[n]; - - lat = new float[n]; - lon = new float[n]; - btm = new float[n]; - ls = new int[n]; - lt = new int[n]; - sal_typ = new int[n]; - sqc = new float[n]; - tmp_typ = new int[n]; - tqc = new float[n]; - - lvl = new float * [n]; - sal = new float * [n]; - sal_err = new float * [n]; - sprb = new float * [n]; - tmp = new float * [n]; - tmp_err = new float * [n]; - clm_sal = new float * [n]; // skipped? - tprb = new float * [n]; - cssd = new float * [n]; - clm_tmp = new float * [n]; // skipped? - ctsd = new float * [n]; - flg = new int * [n]; -} - - -void -RTOFSOb:: - allocate2d() -{ - for (int i = 0; i < n; i ++) - { - int k = lt[i]; - - lvl[i] = new float[k]; - sal[i] = new float[k]; - sal_err[i] = new float[k]; - sprb[i] = new float[k]; - tmp[i] = new float[k]; - tmp_err[i] = new float[k]; - tprb[i] = new float[k]; - clm_sal[i] = new float[k]; - cssd[i] = new float[k]; - clm_tmp[i] = new float[k]; - ctsd[i] = new float[k]; - flg[i] = new int[k]; - } -} - - - -void -RTOFSOb:: - read(FILE * f) -{ - read_float_array(f, btm, n); - read_float_array(f, lat, n); - read_float_array(f, lon, n); - read_int_array(f, ls, n); - read_int_array(f, lt, n); - read_int_array(f, sal_typ, n); - read_float_array(f, sqc, n); - read_int_array(f, tmp_typ, n); - read_float_array(f, tqc, n); - - allocate2d(); - - for (int i = 0; i < n; i ++) - { - int k = lt[i]; - - read_float_array(f, lvl[i], k); - read_float_array(f, sal[i], k); - read_float_array(f, sal_err[i], k); - read_float_array(f, sprb[i], k); - read_float_array(f, tmp[i], k); - read_float_array(f, tmp_err[i], k); - read_float_array(f, tprb[i], k); - read_float_array(f, clm_sal[i], k); - read_float_array(f, cssd[i], k); - read_float_array(f, clm_tmp[i], k); - read_float_array(f, ctsd[i], k); - read_int_array(f, flg[i], k); - } -} - - -void -RTOFSOb:: - print(std::string DirectoryName) -{ - if (!file_exists(DirectoryName)) - { - cerr << "Directory " << DirectoryName << "doesn't exist" << endl; - exit(1); - } - print_float_array(DirectoryName + "/latitude", lat, n); - print_float_array(DirectoryName + "/longitude", lon, n); - print_float_array(DirectoryName + "/btm", btm, n); - print_float_array(DirectoryName + "/tqc", tqc, n); - print_float_array(DirectoryName + "/sqc", sqc, n); - print_int_array(DirectoryName + "/lt", lt, n); - print_int_array(DirectoryName + "/ls", ls, n); - print_int_array(DirectoryName + "/sal_typ", sal_typ, n); - print_int_array(DirectoryName + "/tmp_typ", tmp_typ, n); - - print_2d_float_array(DirectoryName + "/tmp", tmp, n, lt); - print_2d_float_array(DirectoryName + "/sal", sal, n, lt); - print_2d_float_array(DirectoryName + "/lvl", lvl, n, lt); - print_2d_float_array(DirectoryName + "/cssd", cssd, n, lt); - print_2d_float_array(DirectoryName + "/ctsd", ctsd, n, lt); - print_2d_float_array(DirectoryName + "/tprb", tprb, n, lt); - print_2d_float_array(DirectoryName + "/sprb", sprb, n, lt); - print_2d_float_array(DirectoryName + "/clm_tmp", clm_tmp, n, lt); - - print_level(DirectoryName + "/tmp0", tmp, n, 0, lt); - print_level(DirectoryName + "/clm_tmp0", clm_tmp, n, 0, lt); - print_level(DirectoryName + "/tmp10", tmp, n, 10, lt); - print_level(DirectoryName + "/clm_tmp10", clm_tmp, n, 10, lt); - print_level(DirectoryName + "/tprb0", tprb, n, 0, lt); - - print_2d_array(DirectoryName + "/flg", flg, n, lt); - - // print lvl2d array - { - ofstream o; - o.open(DirectoryName + "/lvl2d"); - for (int i = 0; i < n; i ++) - for (int j = 0; j < lt[i]; j ++) - o - << j - << endl; - o.close(); - } -} // RTOFSOb::print - - - -RTOFSDataFile:: -RTOFSDataFile(std::string filename): - filename(filename) -{ - if (!file_exists(filename)) - { - cerr << "File not found" << endl; - exit(1); - } - - - const char * fn = filename.c_str(); - f = fopen(fn, "rb"); - if (!f) - { - cerr << "Error opening file " << fn << endl; - exit(1); - } - - read_file(); -} // RTOFSDataFile::RTOFSDataFile - - -void -RTOFSDataFile:: - read_file() -{ - fseek(f, 4, SEEK_CUR); - - int n_read = read_int(f); - int mx_lvl = read_int(f); - int vrsn = read_int(f); - - ob = new RTOFSOb(n_read, mx_lvl, vrsn); - nobs = n_read; - - ob->read(f); - - skip8bytes(f); - - fread(ob->dtg, sizeof(RTOFSOb::char12), n_read, f); - - skip8bytes(f); - - RTOFSOb::char12 * ob_rct = new RTOFSOb::char12[n_read]; - fread(ob_rct, sizeof(RTOFSOb::char12), n_read, f); - - skip8bytes(f); - - RTOFSOb::char7 * ob_sgn = new RTOFSOb::char7[n_read]; - fread(ob_sgn, sizeof(RTOFSOb::char7), n_read, f); - - // we don't know what sgn is. - { - ofstream o("sgn"); - for (int i = 0; i < n_read; i ++) - { - for (int k = 0; k < 7; k ++) - o << ob_sgn[i][k]; - o << std::endl; - } - } - - // we never had an input file with vrsn == 2 - if (vrsn == 2) - { - float ** glb_sal = new float * [n_read]; - for (int i = 0; i < n_read; i ++) - { - int k = ob->lt[i]; - glb_sal[i] = alloc_read_float_array(f, k); - } - } -} // read_file - - - -void skip8bytes(FILE * f) -{ - fseek(f, 8, SEEK_CUR); -} - - -void - read_float_array(FILE * f, float * a, int n) -{ - skip8bytes(f); - - fread(a, sizeof(float), n, f); - - uint32_t * data = reinterpret_cast(a); - - for (int i = 0; i < n; i ++) - data[i] = be32toh(data[i]); -} // read_float_array - - -void - read_int_array(FILE * f, int * a, int n) -{ - skip8bytes(f); - - fread(a, sizeof(int), n, f); - - uint32_t * data = reinterpret_cast(a); - - for (int i = 0; i < n; i ++) - data[i] = be32toh(data[i]); -} // read_int_array - - - -void - print_level(std::string filename, float ** a, int n, int level, int * dim2) -{ - ofstream o(filename); - o - << std::fixed - << std::setprecision(9); - - for (int i = 0; i < n; i ++) - if (level < dim2[i]) - o - << a[i] [level] - << endl; - - o.close(); -} - - -void - print_int_array(std::string filename, int * a, int n) -{ - ofstream o(filename); - for (int i = 0; i < n; i ++) - o - << a[i] - << endl; - o.close(); -} - - -void - print_float_array(std::string filename, float * a, int n) -{ - ofstream o(filename); - o - << std::fixed - << std::setprecision(9); - for (int i = 0; i < n; i ++) - o - << a[i] - << endl; -} - -void - print_2d_float_array(std::string filename, float ** a, int n, int * dim2) -{ - ofstream o(filename); - o - << std::fixed - << std::setprecision(9); - for (int i = 0; i < n; i ++) - for (int j = 0; j < dim2[i]; j ++) - o - << a[i] [j] - << endl; -} - - -template -void - print_2d_array(std::string filename, NUMBER ** a, int n, int * dim2) -{ - ofstream o(filename); - o - << std::fixed - << std::setprecision(9); - for (int i = 0; i < n; i ++) - for (int j = 0; j < dim2[i]; j ++) - o - << a[i] [j] - << endl; -} - - -bool - file_exists(const std::string& name) -{ - struct stat buffer; - return (stat (name.c_str(), &buffer) == 0); -} - - -int - read_int(FILE * f) -{ - int dummy; - fread(&dummy, sizeof(int), 1, f); - dummy = static_cast(be32toh(dummy)); - return dummy; -} // read_int - - - -void - read_real_array(FILE * f, float * a, int n) -{ - fread(a, sizeof(float), n, f); - - uint32_t * data = reinterpret_cast(a); - - for (int i = 0; i < n; i ++) - data[i] = be32toh(data[i]); -} // read_real_array - - -float * - alloc_read_float_array(FILE * f, int n) -{ - float * a = new float[n]; - skip8bytes(f); - read_real_array(f, a, n); - return a; -} - - -void - print_int_array(int * a, int n) -{ - for (int i = 0; i < n; i ++) - cout - << a[i] - << endl; -} - - -int * - alloc_read_int_array(FILE * f, int n) -{ - int * a = new int[n]; - skip8bytes(f); - read_int_array(f, a, n); - return a; -} - - -void - print_real_array(float * a, int n) -{ - cout - << std::fixed - << std::setprecision(9); - for (int i = 0; i < n; i ++) - cout - << a[i] - << endl; -} - - } // namespace rtofs } // namespace gdasapp diff --git a/utils/obsproc/RTOFSSalinity.h b/utils/obsproc/RTOFSSalinity.h index 4789dfa2f..d13073c31 100644 --- a/utils/obsproc/RTOFSSalinity.h +++ b/utils/obsproc/RTOFSSalinity.h @@ -2,7 +2,9 @@ #include +#include +#include "rtofs/RTOFSOb.h" #include "RTOFSInSitu.h" @@ -24,12 +26,58 @@ class RTOFSSalinity: // Read binary file and populate iodaVars gdasapp::obsproc::iodavars::IodaVars - providerToIodaVars(const std::string filename) final + providerToIodaVars(const std::string filename) final; +}; // class RTOFSSalinity + + +// Read binary file and populate iodaVars +gdasapp::obsproc::iodavars::IodaVars +RTOFSSalinity:: + providerToIodaVars(const std::string filename) +{ + ProcessFile(filename); + + // Set the int metadata names + std::vector intMetadataNames = {"sal_typ"}; + + // Set the float metadata name + std::vector floatMetadataNames = {"level"}; + + rtofs::RTOFSOb & ob = * pob; + + gdasapp::obsproc::iodavars::IodaVars iodaVars( + ob.TotalNumberOfValues(), + floatMetadataNames, + intMetadataNames); + + iodaVars.referenceDate_ = "seconds since 1970-01-01T00:00:00Z"; + iodaVars.niMetadata_ = 1; + iodaVars.nfMetadata_ = 1; + + int n = ob.NumberOfObservations(); + + int k = 0; + for (int i = 0; i < n; i ++) { - ProcessFile(filename); - return SalinityIodaVars(); + int64_t date = rtofs::RTOFSOb::SecondsSinceReferenceTime(ob.dtg[i]); + + for (int j = 0; j < ob.lt[i]; j ++) + { + iodaVars.longitude_(k) = ob.lon[i]; + iodaVars.latitude_(k) = ob.lat[i]; + iodaVars.obsVal_(k) = ob.sal[i][j]; + iodaVars.obsError_(k) = ob.sal_err[i][j]; + iodaVars.datetime_(k) = date; + iodaVars.preQc_(k) = ob.flg[i][j]; + iodaVars.intMetadata_.row(k) << ob.sal_typ[i]; + iodaVars.floatMetadata_.row(k) << ob.lvl[i][j]; + + k++; + } } -}; // class RTOFSSalinity + + return iodaVars; +} // RTOFSInSitu::providerToIodaVars } // namespace gdasapp diff --git a/utils/obsproc/RTOFSTemperature.h b/utils/obsproc/RTOFSTemperature.h index add6a7a48..c66b39dc2 100644 --- a/utils/obsproc/RTOFSTemperature.h +++ b/utils/obsproc/RTOFSTemperature.h @@ -1,7 +1,9 @@ #pragma once #include +#include +#include "rtofs/RTOFSOb.h" #include "RTOFSInSitu.h" @@ -23,11 +25,56 @@ class RTOFSTemperature: // Read binary file and populate iodaVars gdasapp::obsproc::iodavars::IodaVars - providerToIodaVars(const std::string filename) final + providerToIodaVars(const std::string filename) final; +}; // class RTOFSTemperature + + + +gdasapp::obsproc::iodavars::IodaVars +RTOFSTemperature:: + providerToIodaVars(const std::string filename) +{ + ProcessFile(filename); + + // Set the int metadata names + std::vector intMetadataNames = {"tmp_typ"}; + + // Set the float metadata name + std::vector floatMetadataNames = {"level"}; + rtofs::RTOFSOb & ob = * pob; + + gdasapp::obsproc::iodavars::IodaVars iodaVars( + ob.TotalNumberOfValues(), + floatMetadataNames, + intMetadataNames); + + iodaVars.referenceDate_ = "seconds since 1970-01-01T00:00:00Z"; + iodaVars.niMetadata_ = 1; + iodaVars.nfMetadata_ = 1; + + int n = ob.NumberOfObservations(); + + int k = 0; + for (int i = 0; i < n; i ++) { - ProcessFile(filename); - return TemperatureIodaVars(); + int64_t date = rtofs::RTOFSOb::SecondsSinceReferenceTime(ob.dtg[i]); + + for (int j = 0; j < ob.lt[i]; j ++) + { + iodaVars.longitude_(k) = ob.lon[i]; + iodaVars.latitude_(k) = ob.lat[i]; + iodaVars.obsVal_(k) = ob.tmp[i][j]; + iodaVars.obsError_(k) = ob.tmp_err[i][j]; + iodaVars.datetime_(k) = date; + iodaVars.preQc_(k) = ob.flg[i][j]; + iodaVars.intMetadata_.row(k) << ob.tmp_typ[i]; + iodaVars.floatMetadata_.row(k) << ob.lvl[i][j]; + + k++; + } } + + return iodaVars; }; // class RTOFSTemperature diff --git a/utils/obsproc/rtofs/CMakeLists.txt b/utils/obsproc/rtofs/CMakeLists.txt new file mode 100644 index 000000000..25f22474f --- /dev/null +++ b/utils/obsproc/rtofs/CMakeLists.txt @@ -0,0 +1,44 @@ +# message("RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR") +# message("Current source directory: ${CMAKE_CURRENT_SOURCE_DIR}") +# message("Current binary directory: ${CMAKE_CURRENT_BINARY_DIR}") +# message("CMAKE_INCLUDE_PATH=${CMAKE_INCLUDE_PATH}") +# message("C_INCLUDE_PATH=${C_INCLUDE_PATH}") +# message("CPLUS_INCLUDE_PATH=${CPLUS_INCLUDE_PATH}") +# message("RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR") + +# include_directories(${CMAKE_SOURCE_DIR}) + +set( + rtofs_src_files + RTOFSDataFile.h + RTOFSDataFile.cc + RTOFSOb.h + RTOFSOb.cc + util.cc + util.h +) + +list( + APPEND gdasapp_provider2ioda_src_files + ${rtofs_src_files} +) + + +add_library(rtofs STATIC + ${rtofs_src_files} +) +# ecbuild_add_library( + # TARGET rtofs + # TYPE STATIC + # SOURCES ${rtofs_src_files} +# ) + +target_compile_features( + rtofs + PUBLIC cxx_std_17 +) + +target_link_libraries( rtofs PUBLIC oops ioda NetCDF::NetCDF_CXX ) + +# target_link_libraries(rtofs PUBLIC oops) +target_include_directories( rtofs PUBLIC ${oops_INCLUDE_DIRS} ) diff --git a/utils/obsproc/rtofs/README b/utils/obsproc/rtofs/README new file mode 100644 index 000000000..32cfbb5a4 --- /dev/null +++ b/utils/obsproc/rtofs/README @@ -0,0 +1,14 @@ +Author: E. Givelberg +Tue Jun 4 13:33:48 CDT 2024 + +The code in this rtofs directory reads a binary Fortran +file and generates a data object that can be ingested into ioda. + +The fortran binary file consists of records, including additional +information that needs to be skipped. +There are typically 8 bytes to be skipped between 2 data arrays. +The read data is converted from the big endian to the host +format, which is little endian in linux. + +There are currently several data objects that are read in RTOFSDataFile +and are discarded. diff --git a/utils/obsproc/rtofs/RTOFSDataFile.cc b/utils/obsproc/rtofs/RTOFSDataFile.cc new file mode 100644 index 000000000..cd7c97d01 --- /dev/null +++ b/utils/obsproc/rtofs/RTOFSDataFile.cc @@ -0,0 +1,95 @@ +#include +using std::cerr; +using std::endl; + +#include +using std::ofstream; + +#include "RTOFSDataFile.h" +#include "util.h" + + +namespace gdasapp +{ +namespace rtofs +{ + + +RTOFSDataFile:: +RTOFSDataFile(std::string filename) +{ + if (!file_exists(filename)) + { + cerr << "File not found" << endl; + exit(1); + } + + + const char * fn = filename.c_str(); + f = fopen(fn, "rb"); + if (!f) + { + cerr << "Error opening file " << fn << endl; + exit(1); + } + + read_file(); +} // RTOFSDataFile::RTOFSDataFile + + +void +RTOFSDataFile:: + read_file() +{ + fseek(f, 4, SEEK_CUR); + + int n_read = read_int(f); + int mx_lvl = read_int(f); + int vrsn = read_int(f); + + ob = new RTOFSOb(n_read, mx_lvl, vrsn); + + ob->read(f); + + skip8bytes(f); + + fread(ob->dtg, sizeof(RTOFSOb::char12), n_read, f); + + skip8bytes(f); + + RTOFSOb::char12 * ob_rct = new RTOFSOb::char12[n_read]; + fread(ob_rct, sizeof(RTOFSOb::char12), n_read, f); + + skip8bytes(f); + + RTOFSOb::char7 * ob_sgn = new RTOFSOb::char7[n_read]; + fread(ob_sgn, sizeof(RTOFSOb::char7), n_read, f); + +#ifdef RTOFS_OUTPUT_SGN + // we don't know what sgn is. + // by default, it is not output + { + ofstream o("sgn"); + for (int i = 0; i < n_read; i ++) + { + for (int k = 0; k < 7; k ++) + o << ob_sgn[i][k]; + o << std::endl; + } + } +#endif // RTOFS_OUTPUT_SGN + + // we never had an input file with vrsn == 2 + if (vrsn == 2) + { + float ** glb_sal = new float * [n_read]; + for (int i = 0; i < n_read; i ++) + { + int k = ob->lt[i]; + glb_sal[i] = alloc_read_float_array(f, k); + } + } +} // RTOFSDataFile::read_file + +} // namespace rtofs +} // namespace gdasapp diff --git a/utils/obsproc/rtofs/RTOFSDataFile.h b/utils/obsproc/rtofs/RTOFSDataFile.h new file mode 100644 index 000000000..39d344ac3 --- /dev/null +++ b/utils/obsproc/rtofs/RTOFSDataFile.h @@ -0,0 +1,29 @@ +#pragma once + +#include +#include + +#include "RTOFSOb.h" + + +namespace gdasapp +{ +namespace rtofs +{ + +class RTOFSDataFile +{ + public: + explicit RTOFSDataFile(std::string filename); + RTOFSOb & observations() { return * ob; } + + private: + FILE * f; + RTOFSOb * ob; + + void read_file(); +}; // class RTOFSDataFile + +} // namespace rtofs + +} // namespace gdasapp diff --git a/utils/obsproc/rtofs/RTOFSOb.cc b/utils/obsproc/rtofs/RTOFSOb.cc new file mode 100644 index 000000000..43bdabea8 --- /dev/null +++ b/utils/obsproc/rtofs/RTOFSOb.cc @@ -0,0 +1,218 @@ +// NOLINT +#include +using std::cerr; +using std::endl; +using std::cout; + +#include +using std::ofstream; + +#include "RTOFSOb.h" +#include "util.h" + + +#include "oops/util/dateFunctions.h" + + +namespace gdasapp +{ +namespace rtofs +{ + + + +int64_t +RTOFSOb:: + SecondsSinceReferenceTime(rtofs::RTOFSOb::char12 time) +{ + // parse the date + std::string s(time); + int year = std::stoi(s.substr(0, 4)); + int month = std::stoi(s.substr(4, 2)); + int day = std::stoi(s.substr(6, 2)); + int hour = std::stoi(s.substr(8, 2)); + int minute = std::stoi(s.substr(10, 2)); + int second = 0; + + uint64_t julianDate = + util::datefunctions::dateToJulian(year, month, day); + + // 2440588 = Julian date for January 1, 1970. + int daysSinceEpoch = julianDate - 2440588; + int secondsOffset = util::datefunctions::hmsToSeconds(hour, minute, second); + return static_cast(daysSinceEpoch * 86400.0f) + secondsOffset; +} // SecondsSinceReferenceTime + + + +// the variables marked "skipped" are arrays which are present +// in the binary file, but were skipped by the fortran code +// that was reading the files. + +RTOFSOb:: + RTOFSOb(int n, int mx_lvl, int vrsn): + n(n), + mx_lvl(mx_lvl), + vrsn(vrsn) +{ + allocate(); +} + + +int +RTOFSOb:: + TotalNumberOfValues() +{ + int NumberOfValues = 0; + for (int i = 0; i < n; i ++) + NumberOfValues += lt[i]; + return NumberOfValues; +} + + + +void +RTOFSOb:: + allocate() +{ + dtg = new char12[n]; + + lat = new float[n]; + lon = new float[n]; + btm = new float[n]; + ls = new int[n]; + lt = new int[n]; + sal_typ = new int[n]; + sqc = new float[n]; + tmp_typ = new int[n]; + tqc = new float[n]; + + lvl = new float * [n]; + sal = new float * [n]; + sal_err = new float * [n]; + sprb = new float * [n]; + tmp = new float * [n]; + tmp_err = new float * [n]; + clm_sal = new float * [n]; // skipped + tprb = new float * [n]; + cssd = new float * [n]; + clm_tmp = new float * [n]; // skipped + ctsd = new float * [n]; + flg = new int * [n]; +} + + +void +RTOFSOb:: + allocate2d() +{ + for (int i = 0; i < n; i ++) + { + int k = lt[i]; + + lvl[i] = new float[k]; + sal[i] = new float[k]; + sal_err[i] = new float[k]; + sprb[i] = new float[k]; + tmp[i] = new float[k]; + tmp_err[i] = new float[k]; + tprb[i] = new float[k]; + clm_sal[i] = new float[k]; + cssd[i] = new float[k]; + clm_tmp[i] = new float[k]; + ctsd[i] = new float[k]; + flg[i] = new int[k]; + } +} + + + +void +RTOFSOb:: + read(FILE * f) +{ + read_float_array(f, btm, n); + read_float_array(f, lat, n); + read_float_array(f, lon, n); + read_int_array(f, ls, n); + read_int_array(f, lt, n); + read_int_array(f, sal_typ, n); + read_float_array(f, sqc, n); + read_int_array(f, tmp_typ, n); + read_float_array(f, tqc, n); + + allocate2d(); + + for (int i = 0; i < n; i ++) + { + int k = lt[i]; + + read_float_array(f, lvl[i], k); + read_float_array(f, sal[i], k); + read_float_array(f, sal_err[i], k); + read_float_array(f, sprb[i], k); + read_float_array(f, tmp[i], k); + read_float_array(f, tmp_err[i], k); + read_float_array(f, tprb[i], k); + read_float_array(f, clm_sal[i], k); + read_float_array(f, cssd[i], k); + read_float_array(f, clm_tmp[i], k); + read_float_array(f, ctsd[i], k); + read_int_array(f, flg[i], k); + } +} + + +// helpful function to dump ascii files that can be visualized +// with python scripts +void +RTOFSOb:: + print(std::string DirectoryName) +{ + if (!file_exists(DirectoryName)) + { + cerr << "Directory " << DirectoryName << "doesn't exist" << endl; + exit(1); + } + print_float_array(DirectoryName + "/latitude", lat, n); + print_float_array(DirectoryName + "/longitude", lon, n); + print_float_array(DirectoryName + "/btm", btm, n); + print_float_array(DirectoryName + "/tqc", tqc, n); + print_float_array(DirectoryName + "/sqc", sqc, n); + print_int_array(DirectoryName + "/lt", lt, n); + print_int_array(DirectoryName + "/ls", ls, n); + print_int_array(DirectoryName + "/sal_typ", sal_typ, n); + print_int_array(DirectoryName + "/tmp_typ", tmp_typ, n); + + print_2d_array(DirectoryName + "/tmp", tmp, n, lt); + print_2d_array(DirectoryName + "/sal", sal, n, lt); + print_2d_array(DirectoryName + "/lvl", lvl, n, lt); + print_2d_array(DirectoryName + "/cssd", cssd, n, lt); + print_2d_array(DirectoryName + "/ctsd", ctsd, n, lt); + print_2d_array(DirectoryName + "/tprb", tprb, n, lt); + print_2d_array(DirectoryName + "/sprb", sprb, n, lt); + print_2d_array(DirectoryName + "/clm_tmp", clm_tmp, n, lt); + + print_level(DirectoryName + "/tmp0", tmp, n, 0, lt); + print_level(DirectoryName + "/clm_tmp0", clm_tmp, n, 0, lt); + print_level(DirectoryName + "/tmp10", tmp, n, 10, lt); + print_level(DirectoryName + "/clm_tmp10", clm_tmp, n, 10, lt); + print_level(DirectoryName + "/tprb0", tprb, n, 0, lt); + + print_2d_array(DirectoryName + "/flg", flg, n, lt); + + // print lvl2d array + { + ofstream o; + o.open(DirectoryName + "/lvl2d"); + for (int i = 0; i < n; i ++) + for (int j = 0; j < lt[i]; j ++) + o + << j + << endl; + o.close(); + } +} // RTOFSOb::print + +} // namespace rtofs +} // namespace gdasapp diff --git a/utils/obsproc/rtofs/RTOFSOb.h b/utils/obsproc/rtofs/RTOFSOb.h new file mode 100644 index 000000000..fc7cc6c3b --- /dev/null +++ b/utils/obsproc/rtofs/RTOFSOb.h @@ -0,0 +1,66 @@ +#pragma once + +#include +#include + + +namespace gdasapp +{ + +namespace rtofs +{ + +class RTOFSOb +{ + public: + RTOFSOb(int n, int mx_lvl, int vrsn); + + typedef char char12[12]; + typedef char char7[7]; + + static int64_t SecondsSinceReferenceTime(char12 time); + + void read(FILE * f); + int NumberOfObservations() { return n; } + int TotalNumberOfValues(); + void print(std::string DirectoryName); + + + float * btm; // bottom depth + float * lat; // latitude + float * lon; // longitude + int * ls; // ?? + int * lt; // array of dimensions, the number of levels + int * sal_typ; // ?? salinity type ?? + float * sqc; // salinity qc + int * tmp_typ; // ?? temperature type ?? + float * tqc; // temperature qc + + // observations per level: + float ** lvl; // depth + float ** sal; // salinity (PSU) + float ** sal_err; // salinity error (std deviation, PSU) + float ** sprb; // ?? salinity ... ?? + float ** tmp; // in situ temperature (C) + float ** tmp_err; // in situ temperature error (C) + float ** tprb; // ?? temperature ... ?? + float ** clm_sal; // climatology of salinity + float ** cssd; // climatological std deviation for salinity + float ** clm_tmp; // climatology of temperature + float ** ctsd; // climatological std deviation for temperature + int ** flg; // ?? qc flags ?? + + char12 * dtg; // date (Julian, seconds) + + private: + int n; + int mx_lvl; + int vrsn; + + void allocate(); + void allocate2d(); +}; // class RTOFSOb + +} // namespace rtofs + +} // namespace gdasapp diff --git a/utils/obsproc/rtofs/j b/utils/obsproc/rtofs/j new file mode 100644 index 000000000..4eb98b594 --- /dev/null +++ b/utils/obsproc/rtofs/j @@ -0,0 +1,33 @@ +UpdateCTestConfiguration from :/home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/build/gdas-utils/DartConfiguration.tcl +UpdateCTestConfiguration from :/home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/build/gdas-utils/DartConfiguration.tcl +Test project /home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/build/gdas-utils +Constructing a list of tests +Done constructing a list of tests +Updating test list for fixtures +Added 0 tests to meet fixture requirements +Checking test dependency graph... +Checking test dependency graph end +test 1 + Start 1: test_gdasapp_util_coding_norms + +1: Test command: /home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/build/bin/gdas-utils_cpplint.py "--quiet" "--recursive" "/home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/bundle/gdas-utils/test/../" +1: Environment variables: +1: OMP_NUM_THREADS=1 +1: Test timeout computed to be: 10000000 +1: /home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/bundle/gdas-utils/test/../obsproc/rtofs/RTOFSDataFile.cc:1: bundle/gdas-utils/obsproc/rtofs/RTOFSDataFile.cc should include its header file bundle/gdas-utils/obsproc/rtofs/RTOFSDataFile.h [build/include] [5] +1: /home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/bundle/gdas-utils/test/../obsproc/rtofs/RTOFSOb.cc:2: bundle/gdas-utils/obsproc/rtofs/RTOFSOb.cc should include its header file bundle/gdas-utils/obsproc/rtofs/RTOFSOb.h [build/include] [5] +1/1 Test #1: test_gdasapp_util_coding_norms ...***Failed 2.31 sec + +0% tests passed, 1 tests failed out of 1 + +Label Time Summary: +gdas-utils = 2.31 sec*proc (1 test) +script = 2.31 sec*proc (1 test) + +Total Test time (real) = 2.32 sec + +The following tests FAILED: + 1 - test_gdasapp_util_coding_norms (Failed) +Errors while running CTest +Output from these tests are in: /home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/build/gdas-utils/Testing/Temporary/LastTest.log +Use "--rerun-failed --output-on-failure" to re-run the failed cases verbosely. diff --git a/utils/obsproc/rtofs/util.cc b/utils/obsproc/rtofs/util.cc new file mode 100644 index 000000000..ae96fc655 --- /dev/null +++ b/utils/obsproc/rtofs/util.cc @@ -0,0 +1,216 @@ +#include "util.h" + +#include // used in file_exists + +#include +using std::ofstream; + +#include +using std::setprecision; + +#include +using std::cerr; +using std::endl; +using std::cout; + + + +namespace gdasapp +{ + +namespace rtofs +{ + + +void skip8bytes(FILE * f) +{ + fseek(f, 8, SEEK_CUR); +} + + +void + read_float_array(FILE * f, float * a, int n) +{ + skip8bytes(f); + + fread(a, sizeof(float), n, f); + + uint32_t * data = reinterpret_cast(a); + + for (int i = 0; i < n; i ++) + data[i] = be32toh(data[i]); +} // read_float_array + + +void + read_int_array(FILE * f, int * a, int n) +{ + skip8bytes(f); + + fread(a, sizeof(int), n, f); + + uint32_t * data = reinterpret_cast(a); + + for (int i = 0; i < n; i ++) + data[i] = be32toh(data[i]); +} // read_int_array + + + +void + print_level(std::string filename, float ** a, int n, int level, int * dim2) +{ + ofstream o(filename); + o + << std::fixed + << std::setprecision(9); + + for (int i = 0; i < n; i ++) + if (level < dim2[i]) + o + << a[i] [level] + << endl; + + o.close(); +} + + +void + print_int_array(std::string filename, int * a, int n) +{ + ofstream o(filename); + for (int i = 0; i < n; i ++) + o + << a[i] + << endl; + o.close(); +} + + +void + print_float_array(std::string filename, float * a, int n) +{ + ofstream o(filename); + o + << std::fixed + << std::setprecision(9); + for (int i = 0; i < n; i ++) + o + << a[i] + << endl; +} + +void + print_2d_float_array(std::string filename, float ** a, int n, int * dim2) +{ + ofstream o(filename); + o + << std::fixed + << std::setprecision(9); + for (int i = 0; i < n; i ++) + for (int j = 0; j < dim2[i]; j ++) + o + << a[i] [j] + << endl; +} + + + +bool + file_exists(const std::string& name) +{ + struct stat buffer; + return (stat (name.c_str(), &buffer) == 0); +} + + +int + read_int(FILE * f) +{ + int dummy; + fread(&dummy, sizeof(int), 1, f); + dummy = static_cast(be32toh(dummy)); + return dummy; +} // read_int + + + +void + read_real_array(FILE * f, float * a, int n) +{ + fread(a, sizeof(float), n, f); + + uint32_t * data = reinterpret_cast(a); + + for (int i = 0; i < n; i ++) + data[i] = be32toh(data[i]); +} // read_real_array + + +float * + alloc_read_float_array(FILE * f, int n) +{ + float * a = new float[n]; + skip8bytes(f); + read_real_array(f, a, n); + return a; +} + + +void + print_int_array(int * a, int n) +{ + for (int i = 0; i < n; i ++) + cout + << a[i] + << endl; +} + + +int * + alloc_read_int_array(FILE * f, int n) +{ + int * a = new int[n]; + skip8bytes(f); + read_int_array(f, a, n); + return a; +} + + +void + print_real_array(float * a, int n) +{ + cout + << std::fixed + << std::setprecision(9); + for (int i = 0; i < n; i ++) + cout + << a[i] + << endl; +} + + +std::vector + find_unique_values(int * instrument_list, int n) +{ + std::vector instrument; + + // Function to insert a value while keeping + // the vector sorted and without duplicates + auto insertUniqueSorted = [&](int value) + { + auto it = std::lower_bound(instrument.begin(), instrument.end(), value); + if (it == instrument.end() || *it != value) + instrument.insert(it, value); + }; + + for (int i = 0; i < n; i ++) + insertUniqueSorted(instrument_list[i]); + + return instrument; +} + + +} // namespace rtofs + +} // namespace gdasapp diff --git a/utils/obsproc/rtofs/util.h b/utils/obsproc/rtofs/util.h new file mode 100644 index 000000000..5dee3d6ea --- /dev/null +++ b/utils/obsproc/rtofs/util.h @@ -0,0 +1,61 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + + +namespace gdasapp +{ + + +namespace rtofs +{ + +bool file_exists(const std::string& name); + +void skip8bytes(FILE * f); +int read_int(FILE * f); +void read_float_array(FILE * f, float * a, int n); +void read_int_array(FILE * f, int * a, int n); + +void print_float_array(std::string filename, float * a, int n); +void print_int_array(std::string filename, int * a, int n); + + +// void print_2d_float_array(std::string filename, float ** a, int n, int * dim2); + +// template +// void print_2d_array(std::string filename, NUMBER ** a, int n, int * dim2); + +void print_level(std::string filename, float ** a, int n, int level, int * dim2); +float * alloc_read_float_array(FILE * f, int n); + + +std::vector find_unique_values(int * instrument_list, int n); + + + + +template +void + print_2d_array(std::string filename, NUMBER ** a, int n, int * dim2) +{ + std::ofstream o(filename); + o + << std::fixed + << std::setprecision(9); + for (int i = 0; i < n; i ++) + for (int j = 0; j < dim2[i]; j ++) + o + << a[i] [j] + << std::endl; +} + + +} // namespace rtofs + +} // namespace gdasapp diff --git a/utils/test/CMakeLists.txt b/utils/test/CMakeLists.txt index 2d912c2d2..3eee8010d 100644 --- a/utils/test/CMakeLists.txt +++ b/utils/test/CMakeLists.txt @@ -3,7 +3,8 @@ list( APPEND utils_test_input testinput/gdas_meanioda.yaml testinput/gdas_rads2ioda.yaml testinput/gdas_ghrsst2ioda.yaml - testinput/gdas_rtofsinsitu.yaml + testinput/gdas_rtofstmp.yaml + testinput/gdas_rtofssal.yaml testinput/gdas_smap2ioda.yaml testinput/gdas_smos2ioda.yaml testinput/gdas_icecamsr2ioda.yaml @@ -12,7 +13,8 @@ list( APPEND utils_test_input set( gdas_utils_test_ref testref/ghrsst2ioda.test - # testref/rtofsinsitu.test + testref/rtofstmp.test + testref/rtofssal.test testref/rads2ioda.test testref/smap2ioda.test testref/smos2ioda.test @@ -32,10 +34,13 @@ CREATE_SYMLINK( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ${utils_ execute_process( COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/cpplint.py ${CMAKE_BINARY_DIR}/bin/${PROJECT_NAME}_cpplint.py) # add linter for the utils +# TODO: the linter complains about local .h files included in the .cc +# in the rtofs directory +# the exclude argument needs to be removed ecbuild_add_test( TARGET test_gdasapp_util_coding_norms TYPE SCRIPT COMMAND ${CMAKE_BINARY_DIR}/bin/${PROJECT_NAME}_cpplint.py - ARGS --quiet --recursive ${CMAKE_CURRENT_SOURCE_DIR}/../ + ARGS --quiet --recursive --exclude=${CMAKE_CURRENT_SOURCE_DIR}/../obsproc/rtofs/*.cc ${CMAKE_CURRENT_SOURCE_DIR}/../ WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/bin ) # Test example IODA utility that computes the mean of a variable @@ -74,6 +79,18 @@ ecbuild_add_test( TARGET test_gdasapp_util_ghrsst2ioda # LIBS gdas-utils # WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc) +ecbuild_add_test( TARGET test_gdasapp_util_rtofstmp + COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x + ARGS "../testinput/gdas_rtofstmp.yaml" + LIBS gdas-utils + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc) + +ecbuild_add_test( TARGET test_gdasapp_util_rtofssal + COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x + ARGS "../testinput/gdas_rtofssal.yaml" + LIBS gdas-utils + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc) + # Test the SMAP to IODA converter ecbuild_add_test( TARGET test_gdasapp_util_smap2ioda COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x From 821dd36737def5bec0bd12505de8d682b4f4ac33 Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Wed, 5 Jun 2024 19:26:06 -0500 Subject: [PATCH 11/18] junk file j deleted --- utils/obsproc/rtofs/j | 33 --------------------------------- 1 file changed, 33 deletions(-) delete mode 100644 utils/obsproc/rtofs/j diff --git a/utils/obsproc/rtofs/j b/utils/obsproc/rtofs/j deleted file mode 100644 index 4eb98b594..000000000 --- a/utils/obsproc/rtofs/j +++ /dev/null @@ -1,33 +0,0 @@ -UpdateCTestConfiguration from :/home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/build/gdas-utils/DartConfiguration.tcl -UpdateCTestConfiguration from :/home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/build/gdas-utils/DartConfiguration.tcl -Test project /home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/build/gdas-utils -Constructing a list of tests -Done constructing a list of tests -Updating test list for fixtures -Added 0 tests to meet fixture requirements -Checking test dependency graph... -Checking test dependency graph end -test 1 - Start 1: test_gdasapp_util_coding_norms - -1: Test command: /home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/build/bin/gdas-utils_cpplint.py "--quiet" "--recursive" "/home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/bundle/gdas-utils/test/../" -1: Environment variables: -1: OMP_NUM_THREADS=1 -1: Test timeout computed to be: 10000000 -1: /home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/bundle/gdas-utils/test/../obsproc/rtofs/RTOFSDataFile.cc:1: bundle/gdas-utils/obsproc/rtofs/RTOFSDataFile.cc should include its header file bundle/gdas-utils/obsproc/rtofs/RTOFSDataFile.h [build/include] [5] -1: /home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/bundle/gdas-utils/test/../obsproc/rtofs/RTOFSOb.cc:2: bundle/gdas-utils/obsproc/rtofs/RTOFSOb.cc should include its header file bundle/gdas-utils/obsproc/rtofs/RTOFSOb.h [build/include] [5] -1/1 Test #1: test_gdasapp_util_coding_norms ...***Failed 2.31 sec - -0% tests passed, 1 tests failed out of 1 - -Label Time Summary: -gdas-utils = 2.31 sec*proc (1 test) -script = 2.31 sec*proc (1 test) - -Total Test time (real) = 2.32 sec - -The following tests FAILED: - 1 - test_gdasapp_util_coding_norms (Failed) -Errors while running CTest -Output from these tests are in: /home/edwardg/da/edwardg/global-workflow/sorc/gdas.cd/build/gdas-utils/Testing/Temporary/LastTest.log -Use "--rerun-failed --output-on-failure" to re-run the failed cases verbosely. From e34f55231d9c904d93863112348fc9513e761958 Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Thu, 6 Jun 2024 13:22:54 -0500 Subject: [PATCH 12/18] added cmake code for staging rtofs input files --- utils/obsproc/applications/CMakeLists.txt | 5 +- .../applications/gdas_obsprovider2ioda.h | 1 - utils/obsproc/rtofs/CMakeLists.txt | 20 +------ utils/test/CMakeLists.txt | 53 +++++++++++++++---- 4 files changed, 45 insertions(+), 34 deletions(-) diff --git a/utils/obsproc/applications/CMakeLists.txt b/utils/obsproc/applications/CMakeLists.txt index f03631855..99b517aa2 100644 --- a/utils/obsproc/applications/CMakeLists.txt +++ b/utils/obsproc/applications/CMakeLists.txt @@ -10,6 +10,5 @@ target_compile_features( gdas_obsprovider2ioda.x PUBLIC cxx_std_17) target_link_libraries( gdas_obsprovider2ioda.x PUBLIC oops ioda NetCDF::NetCDF_CXX) -# to be used when rtofs is a static library: -# link_directories(${CMAKE_SOURCE_DIR}/rtofs) -# target_link_libraries( gdas_obsprovider2ioda.x PRIVATE rtofs) +link_directories(${CMAKE_SOURCE_DIR}/rtofs) +target_link_libraries(gdas_obsprovider2ioda.x PRIVATE rtofs) diff --git a/utils/obsproc/applications/gdas_obsprovider2ioda.h b/utils/obsproc/applications/gdas_obsprovider2ioda.h index eb53fce0b..ac08b19c6 100644 --- a/utils/obsproc/applications/gdas_obsprovider2ioda.h +++ b/utils/obsproc/applications/gdas_obsprovider2ioda.h @@ -36,7 +36,6 @@ namespace gdasapp { Ghrsst2Ioda conv2ioda(fullConfig, this->getComm()); conv2ioda.writeToIoda(); } else if (provider == "RTOFStmp") { - // RTOFSInSitu conv2ioda(fullConfig, this->getComm()); RTOFSTemperature conv2ioda(fullConfig, this->getComm()); conv2ioda.writeToIoda(); } else if (provider == "RTOFSsal") { diff --git a/utils/obsproc/rtofs/CMakeLists.txt b/utils/obsproc/rtofs/CMakeLists.txt index 25f22474f..d524b22a3 100644 --- a/utils/obsproc/rtofs/CMakeLists.txt +++ b/utils/obsproc/rtofs/CMakeLists.txt @@ -1,13 +1,3 @@ -# message("RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR") -# message("Current source directory: ${CMAKE_CURRENT_SOURCE_DIR}") -# message("Current binary directory: ${CMAKE_CURRENT_BINARY_DIR}") -# message("CMAKE_INCLUDE_PATH=${CMAKE_INCLUDE_PATH}") -# message("C_INCLUDE_PATH=${C_INCLUDE_PATH}") -# message("CPLUS_INCLUDE_PATH=${CPLUS_INCLUDE_PATH}") -# message("RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR") - -# include_directories(${CMAKE_SOURCE_DIR}) - set( rtofs_src_files RTOFSDataFile.h @@ -27,18 +17,10 @@ list( add_library(rtofs STATIC ${rtofs_src_files} ) -# ecbuild_add_library( - # TARGET rtofs - # TYPE STATIC - # SOURCES ${rtofs_src_files} -# ) target_compile_features( rtofs PUBLIC cxx_std_17 ) -target_link_libraries( rtofs PUBLIC oops ioda NetCDF::NetCDF_CXX ) - -# target_link_libraries(rtofs PUBLIC oops) -target_include_directories( rtofs PUBLIC ${oops_INCLUDE_DIRS} ) +target_link_libraries(rtofs PUBLIC oops ioda NetCDF::NetCDF_CXX) diff --git a/utils/test/CMakeLists.txt b/utils/test/CMakeLists.txt index 3eee8010d..87c8ffd00 100644 --- a/utils/test/CMakeLists.txt +++ b/utils/test/CMakeLists.txt @@ -70,26 +70,57 @@ ecbuild_add_test( TARGET test_gdasapp_util_ghrsst2ioda LIBS gdas-utils WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc) -#TODO: -# add binary files for the test -# Test the RTOFS in Situ converter -# ecbuild_add_test( TARGET test_gdasapp_util_rtofsinsitu - # COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x - # ARGS "../testinput/gdas_rtofsinsitu.yaml" - # LIBS gdas-utils - # WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc) - -ecbuild_add_test( TARGET test_gdasapp_util_rtofstmp + + +# copy rtofs binary input files to the testing area +# and generate the tests +execute_process( + COMMAND hostname + OUTPUT_VARIABLE HOSTNAME OUTPUT_STRIP_TRAILING_WHITESPACE +) + +set(RTOFS_INPUT_FILE1 "rtofsinsitu_2024032600.profile") +set(RTOFS_INPUT_FILE2 "rtofsinsitu_2024032700.profile") + +if (${HOSTNAME} MATCHES "^([Oo][rR]ion)") + set(RTOFS_FILES_PATH + "/work/noaa/da/marineda/gfs-marine/data/obs/ci/obs" + ) +endif() +if (${HOSTNAME} MATCHES "^(Hera)") + set(RTOFS_FILES_PATH + "/scratch1/NCEPDEV/da/common/ci/obs" + ) +endif() + +set(RTOFS_FILE1 "${RTOFS_FILES_PATH}/${RTOFS_INPUT_FILE1}") +set(RTOFS_FILE2 "${RTOFS_FILES_PATH}/${RTOFS_INPUT_FILE2}") +set(DESTINATION_DIR "${CMAKE_CURRENT_BINARY_DIR}/obsproc") + +if (EXISTS ${RTOFS_FILE1} AND EXISTS ${RTOFS_FILE2}) + message("Found RTOFS files; generating tests") + + file(COPY ${RTOFS_FILE1} DESTINATION ${DESTINATION_DIR}) + file(COPY ${RTOFS_FILE2} DESTINATION ${DESTINATION_DIR}) + + # Test the RTOFStmp to IODA converter + ecbuild_add_test( TARGET test_gdasapp_util_rtofstmp COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x ARGS "../testinput/gdas_rtofstmp.yaml" LIBS gdas-utils WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc) -ecbuild_add_test( TARGET test_gdasapp_util_rtofssal + # Test the RTOFSsal to IODA converter + ecbuild_add_test( TARGET test_gdasapp_util_rtofssal COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x ARGS "../testinput/gdas_rtofssal.yaml" LIBS gdas-utils WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc) +else() + message("Error: RTOFS input files not found; no test generated.") +endif() + + # Test the SMAP to IODA converter ecbuild_add_test( TARGET test_gdasapp_util_smap2ioda From c625e187b8042af71d2579ccb79a1943721d2d09 Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Fri, 7 Jun 2024 01:13:38 -0500 Subject: [PATCH 13/18] cmake test checks rtofs dir location on orion/hera --- utils/test/CMakeLists.txt | 47 +++++++++++++++++++-------------------- 1 file changed, 23 insertions(+), 24 deletions(-) diff --git a/utils/test/CMakeLists.txt b/utils/test/CMakeLists.txt index 87c8ffd00..21347131d 100644 --- a/utils/test/CMakeLists.txt +++ b/utils/test/CMakeLists.txt @@ -74,50 +74,49 @@ ecbuild_add_test( TARGET test_gdasapp_util_ghrsst2ioda # copy rtofs binary input files to the testing area # and generate the tests -execute_process( - COMMAND hostname - OUTPUT_VARIABLE HOSTNAME OUTPUT_STRIP_TRAILING_WHITESPACE -) set(RTOFS_INPUT_FILE1 "rtofsinsitu_2024032600.profile") set(RTOFS_INPUT_FILE2 "rtofsinsitu_2024032700.profile") -if (${HOSTNAME} MATCHES "^([Oo][rR]ion)") - set(RTOFS_FILES_PATH - "/work/noaa/da/marineda/gfs-marine/data/obs/ci/obs" - ) -endif() -if (${HOSTNAME} MATCHES "^(Hera)") +# rtofs directory on orion: +set(RTOFS_FILES_PATH + "/work/noaa/da/marineda/gfs-marine/data/obs/ci/obs" +) +if (NOT EXISTS ${RTOFS_FILES_PATH}) + # rtofs directory on hera: set(RTOFS_FILES_PATH "/scratch1/NCEPDEV/da/common/ci/obs" ) endif() +if (NOT EXISTS ${RTOFS_FILES_PATH}) + message("Error: staging directory for RTOFS test files not found") +else() + set(RTOFS_FILE1 "${RTOFS_FILES_PATH}/${RTOFS_INPUT_FILE1}") + set(RTOFS_FILE2 "${RTOFS_FILES_PATH}/${RTOFS_INPUT_FILE2}") + set(DESTINATION_DIR "${CMAKE_CURRENT_BINARY_DIR}/obsproc") -set(RTOFS_FILE1 "${RTOFS_FILES_PATH}/${RTOFS_INPUT_FILE1}") -set(RTOFS_FILE2 "${RTOFS_FILES_PATH}/${RTOFS_INPUT_FILE2}") -set(DESTINATION_DIR "${CMAKE_CURRENT_BINARY_DIR}/obsproc") - -if (EXISTS ${RTOFS_FILE1} AND EXISTS ${RTOFS_FILE2}) - message("Found RTOFS files; generating tests") + if (EXISTS ${RTOFS_FILE1} AND EXISTS ${RTOFS_FILE2}) + message("Found RTOFS files; generating tests") - file(COPY ${RTOFS_FILE1} DESTINATION ${DESTINATION_DIR}) - file(COPY ${RTOFS_FILE2} DESTINATION ${DESTINATION_DIR}) + file(COPY ${RTOFS_FILE1} DESTINATION ${DESTINATION_DIR}) + file(COPY ${RTOFS_FILE2} DESTINATION ${DESTINATION_DIR}) - # Test the RTOFStmp to IODA converter - ecbuild_add_test( TARGET test_gdasapp_util_rtofstmp + # Test the RTOFStmp to IODA converter + ecbuild_add_test( TARGET test_gdasapp_util_rtofstmp COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x ARGS "../testinput/gdas_rtofstmp.yaml" LIBS gdas-utils WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc) - # Test the RTOFSsal to IODA converter - ecbuild_add_test( TARGET test_gdasapp_util_rtofssal + # Test the RTOFSsal to IODA converter + ecbuild_add_test( TARGET test_gdasapp_util_rtofssal COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x ARGS "../testinput/gdas_rtofssal.yaml" LIBS gdas-utils WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc) -else() - message("Error: RTOFS input files not found; no test generated.") + else() + message("Error: RTOFS input files not found; no test generated.") + endif() endif() From 542b72d1b30f34b7674312fa5ae19ac6074afc0f Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Fri, 7 Jun 2024 10:57:47 -0500 Subject: [PATCH 14/18] trying to add changeds from the develop branch --- .github/workflows/norms.yaml | 2 +- CMakeLists.txt | 2 ++ sorc/fv3 | 2 +- sorc/fv3-jedi | 2 +- sorc/gsibec | 2 +- sorc/gsw | 2 +- sorc/ioda | 2 +- sorc/iodaconv | 2 +- sorc/land-imsproc | 2 +- sorc/land-jediincr | 2 +- sorc/oops | 2 +- sorc/saber | 2 +- sorc/soca | 2 +- sorc/ufo | 2 +- sorc/vader | 2 +- utils/obsproc/NetCDFToIodaConverter.h | 8 +++++++ utils/test/testinput/gdas_rtofsinsitu.yaml | 13 ----------- utils/test/testinput/gdas_rtofssal.yaml | 12 ++++++++++ utils/test/testinput/gdas_rtofstmp.yaml | 12 ++++++++++ utils/test/testref/rtofssal.test | 26 ++++++++++++++++++++++ utils/test/testref/rtofstmp.test | 26 ++++++++++++++++++++++ 21 files changed, 100 insertions(+), 27 deletions(-) delete mode 100644 utils/test/testinput/gdas_rtofsinsitu.yaml create mode 100644 utils/test/testinput/gdas_rtofssal.yaml create mode 100644 utils/test/testinput/gdas_rtofstmp.yaml create mode 100644 utils/test/testref/rtofssal.test create mode 100644 utils/test/testref/rtofstmp.test diff --git a/.github/workflows/norms.yaml b/.github/workflows/norms.yaml index 7d180f88f..aca4a4596 100644 --- a/.github/workflows/norms.yaml +++ b/.github/workflows/norms.yaml @@ -25,4 +25,4 @@ jobs: - name: Run C++ linter on utils run: | cd $GITHUB_WORKSPACE/GDASApp/utils/test/ - ./cpplint.py --quiet --recursive $GITHUB_WORKSPACE/GDASApp/utils + ./cpplint.py --quiet --recursive $GITHUB_WORKSPACE/GDASApp/utils --exclude $GITHUB_WORKSPACE/GDASApp/utils/obsproc/rtofs diff --git a/CMakeLists.txt b/CMakeLists.txt index f4f2fd98f..57c149629 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,6 +8,8 @@ cmake_minimum_required( VERSION 3.20 FATAL_ERROR ) find_package(ecbuild 3.5 REQUIRED HINTS ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/../ecbuild) project(GDASApp VERSION 1.0.0 LANGUAGES C CXX Fortran ) +# include_directories(${CMAKE_SOURCE_DIR}) +include_directories(${CMAKE_CURRENT_SOURCE_DIR}) include(GNUInstallDirs) enable_testing() diff --git a/sorc/fv3 b/sorc/fv3 index 61450b4e3..47b766295 160000 --- a/sorc/fv3 +++ b/sorc/fv3 @@ -1 +1 @@ -Subproject commit 61450b4e3e80bb96b26c5f3808ce60b5e5cb4207 +Subproject commit 47b766295682df2e53b59d0dd07b5c511b84dd67 diff --git a/sorc/fv3-jedi b/sorc/fv3-jedi index 35a2d4a84..e960612c0 160000 --- a/sorc/fv3-jedi +++ b/sorc/fv3-jedi @@ -1 +1 @@ -Subproject commit 35a2d4a845227de2a70781474e39617727c1f4e6 +Subproject commit e960612c08e81c02c075286d1833094e99b357ea diff --git a/sorc/gsibec b/sorc/gsibec index c8ac58d9b..e6e1944e1 160000 --- a/sorc/gsibec +++ b/sorc/gsibec @@ -1 +1 @@ -Subproject commit c8ac58d9b43eb8f890e565d12c88f1b0579c9ccd +Subproject commit e6e1944e1cb8f339b4eb79acffb5470bb7723b42 diff --git a/sorc/gsw b/sorc/gsw index 1a02ebaf6..697cbeb76 160000 --- a/sorc/gsw +++ b/sorc/gsw @@ -1 +1 @@ -Subproject commit 1a02ebaf6f7a4e9f2c2d2dd973fb050e697bcc74 +Subproject commit 697cbeb7605d70ed3857664c5f54a5c05346e31f diff --git a/sorc/ioda b/sorc/ioda index 91eb2c764..2967b99a5 160000 --- a/sorc/ioda +++ b/sorc/ioda @@ -1 +1 @@ -Subproject commit 91eb2c7643b33c1af2470d289eefb7aec2667387 +Subproject commit 2967b99a56f60652dba232128f1d6c739d3f4cb7 diff --git a/sorc/iodaconv b/sorc/iodaconv index c0330c8a7..cd8b6be4a 160000 --- a/sorc/iodaconv +++ b/sorc/iodaconv @@ -1 +1 @@ -Subproject commit c0330c8a7d8ba51f74b98c27ce1ddb2db2a276ab +Subproject commit cd8b6be4abff376c632a8270b03f5decc9e2e869 diff --git a/sorc/land-imsproc b/sorc/land-imsproc index 2ac7c0d90..a3b233446 160000 --- a/sorc/land-imsproc +++ b/sorc/land-imsproc @@ -1 +1 @@ -Subproject commit 2ac7c0d90d3ebd449ad11bce155a71d381b4534a +Subproject commit a3b233446002af7643956bc31a7e68b89bf1a9d9 diff --git a/sorc/land-jediincr b/sorc/land-jediincr index 79576b79d..7c1f6a3f8 160000 --- a/sorc/land-jediincr +++ b/sorc/land-jediincr @@ -1 +1 @@ -Subproject commit 79576b79df6a3f1ab8e7dde4425821f2915006b3 +Subproject commit 7c1f6a3f8f949e376786eef7dba55a9e10e9778d diff --git a/sorc/oops b/sorc/oops index a12ec3c30..f69f4b038 160000 --- a/sorc/oops +++ b/sorc/oops @@ -1 +1 @@ -Subproject commit a12ec3c305ebe77bb5862f475c7c523988e18b14 +Subproject commit f69f4b0389602c2c68f46bf935c2aac03427c70c diff --git a/sorc/saber b/sorc/saber index 12b5a2602..188fa533a 160000 --- a/sorc/saber +++ b/sorc/saber @@ -1 +1 @@ -Subproject commit 12b5a2602d7fa74c32c4992f805ce724fe76c89b +Subproject commit 188fa533ac56a83c8a2f05195676ce64c644fd67 diff --git a/sorc/soca b/sorc/soca index 18f7813dc..82ec943fa 160000 --- a/sorc/soca +++ b/sorc/soca @@ -1 +1 @@ -Subproject commit 18f7813dc4eb5033dbba7bc0da65e3c9166fb5a4 +Subproject commit 82ec943fa332840b426fd97370e13a7697adf3c4 diff --git a/sorc/ufo b/sorc/ufo index 1f5a4ce81..fc956f44b 160000 --- a/sorc/ufo +++ b/sorc/ufo @@ -1 +1 @@ -Subproject commit 1f5a4ce81e10295cfc911ebb93c6c431d18e0314 +Subproject commit fc956f44b6bf044ba835a9c320e97d1a5a373060 diff --git a/sorc/vader b/sorc/vader index 7688c0700..f5125e96c 160000 --- a/sorc/vader +++ b/sorc/vader @@ -1 +1 @@ -Subproject commit 7688c07004912b964c479b4387bd3084cea31024 +Subproject commit f5125e96caa1920efe79f91fa71b623369d252c7 diff --git a/utils/obsproc/NetCDFToIodaConverter.h b/utils/obsproc/NetCDFToIodaConverter.h index 89b55b227..214fa25bf 100644 --- a/utils/obsproc/NetCDFToIodaConverter.h +++ b/utils/obsproc/NetCDFToIodaConverter.h @@ -49,6 +49,9 @@ namespace gdasapp { oops::Log::trace() << "NetCDFToIodaConverter::NetCDFToIodaConverter created." << std::endl; } + + + // Method to write out a IODA file (writter called in destructor) void writeToIoda() { // Extract ioda variables from the provider's files @@ -60,6 +63,8 @@ namespace gdasapp { // Read the provider's netcdf file gdasapp::obsproc::iodavars::IodaVars iodaVars = providerToIodaVars(inputFilenames_[myrank]); + + for (int i = myrank + comm_.size(); i < inputFilenames_.size(); i += comm_.size()) { iodaVars.append(providerToIodaVars(inputFilenames_[i])); oops::Log::info() << " appending: " << inputFilenames_[i] << std::endl; @@ -200,6 +205,9 @@ namespace gdasapp { } } + + + private: // Virtual method that reads the provider's netcdf file and store the relevant // info in a IodaVars struct diff --git a/utils/test/testinput/gdas_rtofsinsitu.yaml b/utils/test/testinput/gdas_rtofsinsitu.yaml deleted file mode 100644 index 4c8684d9d..000000000 --- a/utils/test/testinput/gdas_rtofsinsitu.yaml +++ /dev/null @@ -1,13 +0,0 @@ -provider: RTOFS -window begin: 2021-03-24T15:00:00Z -window end: 2021-03-24T21:00:00Z -output file: rtofsinsitu_2024032600.profile.nc4 -input files: -- rtofsinsitu_2024032600.profile -# - 2024032600.profile -# - 2024032600.sfc - -# test: - # reference filename: testref/rtofsinsitu.test - # test output filename: testoutput/rtofsinsitu.test - # float relative tolerance: 1e-6 diff --git a/utils/test/testinput/gdas_rtofssal.yaml b/utils/test/testinput/gdas_rtofssal.yaml new file mode 100644 index 000000000..8bf0165eb --- /dev/null +++ b/utils/test/testinput/gdas_rtofssal.yaml @@ -0,0 +1,12 @@ +provider: RTOFSsal +window begin: 2021-03-24T15:00:00Z +window end: 2021-03-24T21:00:00Z +output file: rtofssal_2024032600.profile.nc4 +input files: +- rtofsinsitu_2024032600.profile +- rtofsinsitu_2024032700.profile + +test: + reference filename: testref/rtofssal.test + test output filename: testoutput/rtofssal.test + float relative tolerance: 1e-6 diff --git a/utils/test/testinput/gdas_rtofstmp.yaml b/utils/test/testinput/gdas_rtofstmp.yaml new file mode 100644 index 000000000..2b92319bb --- /dev/null +++ b/utils/test/testinput/gdas_rtofstmp.yaml @@ -0,0 +1,12 @@ +provider: RTOFStmp +window begin: 2021-03-24T15:00:00Z +window end: 2021-03-24T21:00:00Z +output file: rtofstmp_2024032600.profile.nc4 +input files: +- rtofsinsitu_2024032600.profile +- rtofsinsitu_2024032700.profile + +test: + reference filename: testref/rtofstmp.test + test output filename: testoutput/rtofstmp.test + float relative tolerance: 1e-6 diff --git a/utils/test/testref/rtofssal.test b/utils/test/testref/rtofssal.test new file mode 100644 index 000000000..a63dacdb9 --- /dev/null +++ b/utils/test/testref/rtofssal.test @@ -0,0 +1,26 @@ +Reading: [rtofsinsitu_2024032600.profile,rtofsinsitu_2024032700.profile] +seconds since 1970-01-01T00:00:00Z +obsVal: + Min: -999 + Max: 43.2123 + Sum: 2.51883e+07 +obsError: + Min: 0.0101875 + Max: 1.74092 + Sum: 32320.7 +preQc: + Min: 0 + Max: 1042 + Sum: 3100218 +longitude: + Min: -179.959 + Max: 179.732 + Sum: -2.48938e+07 +latitude: + Min: -70.1871 + Max: 82.324 + Sum: 1.25108e+06 +datetime: + Min: 1711411200 + Max: 1711583940 + Sum: 1247390658115740 diff --git a/utils/test/testref/rtofstmp.test b/utils/test/testref/rtofstmp.test new file mode 100644 index 000000000..4c65720fb --- /dev/null +++ b/utils/test/testref/rtofstmp.test @@ -0,0 +1,26 @@ +Reading: [rtofsinsitu_2024032600.profile,rtofsinsitu_2024032700.profile] +seconds since 1970-01-01T00:00:00Z +obsVal: + Min: -1.98999 + Max: 34.2 + Sum: 5.97847e+06 +obsError: + Min: 0.0146336 + Max: 1.35821 + Sum: 86632.5 +preQc: + Min: 0 + Max: 1042 + Sum: 3100218 +longitude: + Min: -179.959 + Max: 179.732 + Sum: -2.48938e+07 +latitude: + Min: -70.1871 + Max: 82.324 + Sum: 1.25108e+06 +datetime: + Min: 1711411200 + Max: 1711583940 + Sum: 1247390658115740 From dc240641736a8454bdb2962e16904d7ed00dc12a Mon Sep 17 00:00:00 2001 From: Ed Givelberg Date: Sun, 9 Jun 2024 17:41:59 -0400 Subject: [PATCH 15/18] Update norms.yaml excluding *.cc in rtofs from the linter --- .github/workflows/norms.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/norms.yaml b/.github/workflows/norms.yaml index 7d180f88f..2e721e7af 100644 --- a/.github/workflows/norms.yaml +++ b/.github/workflows/norms.yaml @@ -25,4 +25,4 @@ jobs: - name: Run C++ linter on utils run: | cd $GITHUB_WORKSPACE/GDASApp/utils/test/ - ./cpplint.py --quiet --recursive $GITHUB_WORKSPACE/GDASApp/utils + ./cpplint.py --quiet --recursive --exclude=$GITHUB_WORKSPACE/GDASApp/utils/obsproc/rtofs/*.cc $GITHUB_WORKSPACE/GDASApp/utils From 7d5c4033d59b77b05b89d361bda6e96049b784a2 Mon Sep 17 00:00:00 2001 From: Ed Givelberg Date: Sun, 9 Jun 2024 19:50:04 -0400 Subject: [PATCH 16/18] Update RTOFSInSitu.h --- utils/obsproc/RTOFSInSitu.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/obsproc/RTOFSInSitu.h b/utils/obsproc/RTOFSInSitu.h index a387064ed..ccbd44eb7 100644 --- a/utils/obsproc/RTOFSInSitu.h +++ b/utils/obsproc/RTOFSInSitu.h @@ -69,4 +69,4 @@ RTOFSInSitu:: } // namespace rtofs -} // namespace gdasapp \ No newline at end of file +} // namespace gdasapp From 2949b9a0ed8a244b0e0aa7f3cb47175d214b7b77 Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Mon, 10 Jun 2024 08:36:23 -0500 Subject: [PATCH 17/18] fixed yaml --- utils/test/testinput/gdas_rtofsinsitu.yaml | 13 ------------- 1 file changed, 13 deletions(-) delete mode 100644 utils/test/testinput/gdas_rtofsinsitu.yaml diff --git a/utils/test/testinput/gdas_rtofsinsitu.yaml b/utils/test/testinput/gdas_rtofsinsitu.yaml deleted file mode 100644 index 4c8684d9d..000000000 --- a/utils/test/testinput/gdas_rtofsinsitu.yaml +++ /dev/null @@ -1,13 +0,0 @@ -provider: RTOFS -window begin: 2021-03-24T15:00:00Z -window end: 2021-03-24T21:00:00Z -output file: rtofsinsitu_2024032600.profile.nc4 -input files: -- rtofsinsitu_2024032600.profile -# - 2024032600.profile -# - 2024032600.sfc - -# test: - # reference filename: testref/rtofsinsitu.test - # test output filename: testoutput/rtofsinsitu.test - # float relative tolerance: 1e-6 From 6c9afd10abb5966461c6ebba43631452db78b6f5 Mon Sep 17 00:00:00 2001 From: Edward Givelberg Date: Mon, 10 Jun 2024 11:48:37 -0500 Subject: [PATCH 18/18] conflict syntax --- utils/test/CMakeLists.txt | 5 ----- 1 file changed, 5 deletions(-) diff --git a/utils/test/CMakeLists.txt b/utils/test/CMakeLists.txt index abff93352..ce9eae1ba 100644 --- a/utils/test/CMakeLists.txt +++ b/utils/test/CMakeLists.txt @@ -118,11 +118,6 @@ else() endif() -<<<<<<< HEAD -======= - - ->>>>>>> 7d5c4033d59b77b05b89d361bda6e96049b784a2 # Test the SMAP to IODA converter ecbuild_add_test( TARGET test_gdasapp_util_smap2ioda COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x