diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index edb94ee4..b214792c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -13,6 +13,7 @@ add_subdirectory(ssec) add_subdirectory(wrfda_ncdiag) add_subdirectory(single_observation) add_subdirectory(mrms) +add_subdirectory(tomorrow_io) # Optional components if(iodaconv_gnssro_ENABLED) diff --git a/src/tomorrow_io/CMakeLists.txt b/src/tomorrow_io/CMakeLists.txt new file mode 100644 index 00000000..4a5128c1 --- /dev/null +++ b/src/tomorrow_io/CMakeLists.txt @@ -0,0 +1,23 @@ +# (C) Copyright 2024 The Tomorrow Companies, Inc. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + +ecbuild_add_executable( TARGET convert_tio.x + SOURCES + creation.cc + creation_tio_sat_inst_specs.cc + copy.h + copy_function.cc + copy_helpers.cc + copy_tms.cc + copy_tms_datetime.cc + detect_tms_type.cc + product.cc + product.h + main_tio.cc + LIBS + ioda + ) +set_target_properties(convert_tio.x PROPERTIES CXX_STANDARD 17) + diff --git a/src/tomorrow_io/README.md b/src/tomorrow_io/README.md new file mode 100644 index 00000000..0e71b535 --- /dev/null +++ b/src/tomorrow_io/README.md @@ -0,0 +1,12 @@ +# The Tomorrow Microwave Sounder (TMS) converter + +The source codes here provide the ability to convert the Tomorrow.io TMS L1B-TC product +and the NASA TROPICS L1B product into IODA format. + +Usage: `convert_tio.x input_file_1 [input_file_2 ...] output_file` + +This converter is written in C++ and uses IODA to read and write the observation data files. +It provides an example of how to convert data to ioda using the C++ interface. + +For details, contact [Ryan Honeyager](mailto:ryan.honeyager@tomorrow.io) (@rhoneyager-tomorrow). + diff --git a/src/tomorrow_io/copy.h b/src/tomorrow_io/copy.h new file mode 100644 index 00000000..c302db27 --- /dev/null +++ b/src/tomorrow_io/copy.h @@ -0,0 +1,63 @@ +#pragma once +/* + * (C) Copyright 2024 The Tomorrow Companiec, Inc. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ +#include +#include +#include +#include + +#include "ioda/ObsGroup.h" + +namespace tio_converter { + +using DimensionRanges_t = std::vector>>; + +struct VariableInfo { + mutable ioda::Variable var; ///< The variable. + /// Optionally defines a subset of indices along each axis. Bounds are inclusive on both sides. + DimensionRanges_t range; +}; + +struct VariableDerivedInfo { + /// Supplementary information about the dimensions of the variable. + DimensionRanges_t range; + /// The dimensions of the variable. + ioda::Dimensions dims; + /// The type of the variable's data. Ex: unsigned little-endian 32-bit integer. + ioda::Type type; + /// The starting indices for a hyperslab selection. + std::vector selection_start; + /// The span along each axis for a hyperslab selection. + std::vector selection_count; + /// The number of data elements in this selection. + size_t selection_num_elements; + /// The size, in bytes, needed to store the variable's data, accounting for the selection. + size_t size_bytes; + /// Selection from the file + ioda::Selection selection_ioda; + /// Selection within memory (starts at 0,0,0,...) + ioda::Selection selection_membuf; + + VariableDerivedInfo(); + VariableDerivedInfo(const VariableInfo &); + VariableDerivedInfo(const VariableInfo &, const DimensionRanges_t &); +}; + +/// @brief Get the fill value assigned to a variable as a vector of bytes, and optionally +/// convert to a different type representation. +/// @param var is the variable to be queried. +/// @param as_type is the desired return value's data type. Normally this is the source +/// variable's data type, but optionally you can convert to a different representation. +/// Useful when converting between differing source and destination types. +std::vector get_fill_value(const ioda::Variable &var, std::optional as_type = {}); + +/// @brief Generic function to copy a hyperslab of data from one variable to another. +/// @param from is the specification of the source data. This includes the variable and the hyperslab. +/// @param to is the specification of the destination location. +void copy(const VariableInfo &from, const VariableInfo &to); + +} // namespace tio_converter diff --git a/src/tomorrow_io/copy_function.cc b/src/tomorrow_io/copy_function.cc new file mode 100644 index 00000000..2c48bdc3 --- /dev/null +++ b/src/tomorrow_io/copy_function.cc @@ -0,0 +1,56 @@ +/* + * (C) Copyright 2024 The Tomorrow Companiec, Inc. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ +#include +#include +#include +#include +#include +#include + +#include "copy.h" +#include "hdf5.h" + +namespace tio_converter { + +void copy(const VariableInfo &from, const VariableInfo &to) { + using namespace ioda; + using std::byte; + using std::max; + using std::memcmp; + using std::memcpy; + using std::vector; + + VariableDerivedInfo from_info(from); + VariableDerivedInfo to_info(to); + + vector buffer(to_info.size_bytes); + from.var.read( + gsl::make_span(buffer.data(), buffer.size()), // Read into the buffer + to_info.type, // Convert data into the destination data type (e.g. int, float, ...) + from_info.selection_membuf, // Needed to tell ioda how the data should be mapped into memory + from_info.selection_ioda // The hyperslab being read + ); + + vector from_fill = get_fill_value(from.var, to_info.type); + vector to_fill = get_fill_value(to.var); + const size_t buffer_size_of_element_bytes = to_info.type.getSize(); + // Iterate over buffer in buffer_size_of_element_bytes increments. + // If we match from_fill_as_bytes_in_dest_representation, replace with the + // contents of to_fill_as_bytes_in_dest_representation. + for (size_t i = 0; i < buffer.size(); i += buffer_size_of_element_bytes) { + if (!memcmp(buffer.data() + i, from_fill.data(), buffer_size_of_element_bytes)) + memcpy(buffer.data() + i, to_fill.data(), buffer_size_of_element_bytes); + } + + to.var.write(gsl::make_span(buffer.data(), buffer.size()), // Write from this buffer + to_info.type, // Output variable type + to_info.selection_membuf, // Data mapping in memory + to_info.selection_ioda // The hyperslab being written + ); +} + +} // namespace tio_converter diff --git a/src/tomorrow_io/copy_helpers.cc b/src/tomorrow_io/copy_helpers.cc new file mode 100644 index 00000000..4be43c70 --- /dev/null +++ b/src/tomorrow_io/copy_helpers.cc @@ -0,0 +1,119 @@ +/* + * (C) Copyright 2024 The Tomorrow Companiec, Inc. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ +#include + +#include "copy.h" +#include "hdf5.h" + +namespace tio_converter { + +namespace { +// This only returns predefined HDF5 types, and those hid_ts are static objects. +// NOTE (RH): ioda really should be extended to return the endianness of data. +// We assume little endian until this is fixed. +hid_t get_hdf5_type(const ioda::Type &typ) { + using namespace ioda; + using std::logic_error; + const size_t len = typ.getSize(); + const TypeClass cls = typ.getClass(); + if (cls == TypeClass::Integer) { + bool sgn = typ.isTypeSigned(); + if (sgn && (len == 1)) return H5T_STD_I8LE; + if (!sgn && (len == 1)) return H5T_STD_U8LE; + if (sgn && (len == 2)) return H5T_STD_I16LE; + if (!sgn && (len == 2)) return H5T_STD_U16LE; + if (sgn && (len == 4)) return H5T_STD_I32LE; + if (!sgn && (len == 4)) return H5T_STD_U32LE; + if (sgn && (len == 8)) return H5T_STD_I64LE; + if (!sgn && (len == 8)) return H5T_STD_U64LE; + } else if (cls == TypeClass::Float) { +#ifdef H5T_NATIVE_FLOAT16 // Introduced in recent HDF5 versions + if (len == 2) return H5T_IEEE_F16LE; +#endif + if (len == 4) return H5T_IEEE_F32LE; + if (len == 8) return H5T_IEEE_F64LE; + } + throw logic_error("Unsupported object type."); +} +} // namespace + +VariableDerivedInfo::VariableDerivedInfo() = default; +VariableDerivedInfo::VariableDerivedInfo(const VariableInfo &vi) + : VariableDerivedInfo(vi, vi.range) {} +VariableDerivedInfo::VariableDerivedInfo(const VariableInfo &vi, const DimensionRanges_t &di) { + using ioda::Dimensions_t; + using ioda::SelectionOperator; + using std::accumulate; + using std::multiplies; + using std::vector; + + range = di; + dims = vi.var.getDimensions(); + type = vi.var.getType(); + + vector zero_starts(dims.dimensionality); + selection_start.resize(dims.dimensionality); + selection_count.resize(dims.dimensionality); + + for (size_t i = 0; i < dims.dimensionality; ++i) { + if (range.size() > i && range[i]) { + selection_start[i] = range[i]->first; + selection_count[i] = range[i]->second - range[i]->first + 1; + } else { + selection_start[i] = 0; + selection_count[i] = dims.dimsCur[i]; + } + } + selection_ioda.extent(dims.dimsCur) + .select({SelectionOperator::SET, selection_start, selection_count}); + selection_membuf.extent(selection_count) + .select({SelectionOperator::SET, zero_starts, selection_count}); + + // Determine the size of a buffer needed to read this data in its entirety. + // This is just selection_count. + selection_num_elements + = accumulate(selection_count.begin(), selection_count.end(), 1, multiplies()); + const size_t size_of_element_bytes = type.getSize(); + size_bytes = selection_num_elements * size_of_element_bytes; +} + +std::vector get_fill_value(const ioda::Variable &var, std::optional as_type) { + using std::logic_error; + using std::max; + using std::memcpy; + using std::vector; + const size_t len_bytes_src = var.getType().getSize(); + vector fill_bytes(len_bytes_src); + // BUG (RH): ioda's getFillValue is very slightly buggy in that it reports a + // spurious warning when reading the fill value of the TMS L1B MultiMask variable, + // which has an unsigned char data type. + // "ioda::Variable: hdf and netcdf fill value specifications do not match" + // In this case, we can just read the _FillValue attribute directly. + if (var.atts.exists("_FillValue")) { + ioda::Attribute fvAttr = var.atts.open("_FillValue"); + fvAttr.read(gsl::make_span(fill_bytes.data(), fill_bytes.size()), fvAttr.getType()); + } else { + const auto src_fill = var.getFillValue(); + memcpy(fill_bytes.data(), &(src_fill.fillValue_.ui64), len_bytes_src); + } + + if (!as_type) return fill_bytes; + + const size_t len_bytes_to = as_type->getSize(); + const hid_t h5type_from = get_hdf5_type(var.getType()); + const hid_t h5type_to = get_hdf5_type(*as_type); + fill_bytes.resize(max(len_bytes_src, len_bytes_to)); + + herr_t cvt_res = H5Tconvert(h5type_from, h5type_to, + 1, // Only one 'element' to be converted + fill_bytes.data(), nullptr, H5P_DEFAULT); + if (cvt_res < 0) throw logic_error("Fill value type conversion failed."); + fill_bytes.resize(len_bytes_to); + return fill_bytes; +} + +} // namespace tio_converter diff --git a/src/tomorrow_io/copy_tms.cc b/src/tomorrow_io/copy_tms.cc new file mode 100644 index 00000000..94e0c087 --- /dev/null +++ b/src/tomorrow_io/copy_tms.cc @@ -0,0 +1,279 @@ +/* + * (C) Copyright 2024 The Tomorrow Companiec, Inc. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ +#include +#include +#include +#include +#include +#include + +#include "copy.h" +#include "product.h" + +namespace tio_converter { + +void expandAndCopyCalQualityFlag(const ioda::Variable &from, + std::vector, ioda::Variable>> &bits_to) { + // Source var has dims of channels x scans x spots + // Destination vars have dims of Location x Channel. + // A bit of remapping magic must occur to keep the values in line. + // Let's read and write one channel at a time. + using namespace ioda; + using std::valarray; + using std::vector; + Dimensions from_dims = from.getDimensions(); + Type mem_type = from.getType(); + size_t num_channels = from_dims.dimsCur[0]; + size_t num_locations_per_channel = from_dims.dimsCur[1] * from_dims.dimsCur[2]; + size_t buf_typesz = mem_type.getSize(); + size_t buf_sz = num_locations_per_channel * buf_typesz; + + vector buf(buf_sz); + + for (size_t ich = 0; ich < num_channels; ++ich) { + Selection from_file_selection; + vector from_starts{(Dimensions_t)ich, 0, 0}, + from_counts{1, from_dims.dimsCur[1], from_dims.dimsCur[2]}; + from_file_selection.extent(from_dims.dimsCur) + .select({SelectionOperator::SET, from_starts, from_counts}); + + Selection buf_selection; + vector buf_starts{0, 0}, buf_counts{from_dims.dimsCur[1], from_dims.dimsCur[2]}; + buf_selection.extent(buf_counts).select({SelectionOperator::SET, buf_starts, buf_counts}); + + from.read(gsl::make_span(buf.data(), buf.size()), mem_type, buf_selection, from_file_selection); + valarray buf_from(buf.data(), buf_sz); + + // Break up the input data according to the bit mapping in bits_to. + for (auto &bit_to : bits_to) { + const auto &bit_mask = bit_to.first; + Variable to = bit_to.second; + + valarray buf_to(buf_sz); + buf_to = buf_from & (char)bit_mask.to_ulong(); + + vector buf_to2(std::begin(buf_to), std::end(buf_to)); + replace_if( + buf_to2.begin(), buf_to2.end(), [](char val) { return val > 0; }, 1); + + Dimensions to_dims = to.getDimensions(); + + Selection to_file_selection; + vector to_starts{0, (Dimensions_t)ich}, + to_counts{(Dimensions_t)num_locations_per_channel, 1}; + to_file_selection.extent(to_dims.dimsCur) + .select({SelectionOperator::SET, to_starts, to_counts}); + + to.write(gsl::make_span(buf_to2.data(), buf_to2.size()), mem_type, buf_selection, + to_file_selection); + } + } +} + +void convert_tropics(Converter_Params &p) { + using namespace std; + using namespace ioda; + + const DimensionRanges_t offset_start{{{p.locationOffset, p.locationOffset + p.inNumLocs - 1}}}; + + // Diagnostic flags + + auto makeBits = [](std::initializer_list bits) -> std::bitset<8> { + std::bitset<8> res; + for (auto bit : bits) res.set(bit); + return res; + }; + vector, Variable>> calmap{ + {makeBits({0}), p.out.vars["MetaData/nonOceanFlag_legacy"]}, + {makeBits({1}), p.out.vars["MetaData/intrusionFlag"]}, + {makeBits({2}), p.out.vars["MetaData/spacecraftManeuverFlag"]}, + {makeBits({3}), p.out.vars["MetaData/coldCalibrationFlag"]}, + {makeBits({4}), p.out.vars["MetaData/hotCalibrationFlag"]}, + {makeBits({5}), p.out.vars["MetaData/satelliteAscendingFlag"]}, + {makeBits({6}), p.out.vars["MetaData/dayOrNightQualifier"]}, + {makeBits({7}), p.out.vars["MetaData/payloadOrientationFlag"]}, + {makeBits({1, 2, 3, 4}), p.out.vars["MetaData/compositePreQCFlag"]}}; + expandAndCopyCalQualityFlag(p.in.vars["calQualityFlag"], calmap); + + copy(VariableInfo{.var = p.in.vars["NonOceanFlag"]}, + VariableInfo{.var = p.out.vars["MetaData/nonOceanFlag"], .range = offset_start}); + + // TROPICS TB has channel first, and we need to change to a var that is location first. + for (size_t ich = 0; ich < 12; ++ich) { + const DimensionRanges_t single_ch{{{ich, ich}}, {}, {}}; + const DimensionRanges_t dest_range{{{p.locationOffset, p.locationOffset + p.inNumLocs - 1}}, + {{ich, ich}}}; + + copy(VariableInfo{.var = p.in.vars["brightness_temperature"], .range = single_ch}, + VariableInfo{.var = p.out.vars["ObsValue/brightnessTemperature"], .range = dest_range}); + } +} + +void convert_tio(Converter_Params &p) { + using namespace std; + using namespace ioda; + + const DimensionRanges_t only_ch1{{}, {}, {{0, 0}}}; + const DimensionRanges_t offset_start{{{p.locationOffset, p.locationOffset + p.inNumLocs - 1}}}; + + // Diagnostic flags + + // NOTE (RH): many of the same variables are common to TIO and TROPICS TMS satellites, + // but they come from different source variables. Several of these are distinct variables + // in TIO, but in TROPICS they are gathered into a single "calQualityFlag" variable. + copy(VariableInfo{.var = p.in.vars["MultiMask"]}, + VariableInfo{.var = p.out.vars["MetaData/multiMask"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["clear_fraction"]}, + VariableInfo{.var = p.out.vars["MetaData/cloudClearFraction"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["flagICTCal"]}, + VariableInfo{.var = p.out.vars["MetaData/ICTCalibrationFlag"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["flagLunarIntrusion"]}, + VariableInfo{.var = p.out.vars["MetaData/lunarIntrusionFlag"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["flagNDCal"]}, + VariableInfo{.var = p.out.vars["MetaData/NDCalibrationFlag"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["flagSolarIntrusion"]}, + VariableInfo{.var = p.out.vars["MetaData/solarIntrusionFlag"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["land_fraction"]}, + VariableInfo{.var = p.out.vars["MetaData/landAreaFraction"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["flagAscDesc"], .range = only_ch1}, + VariableInfo{.var = p.out.vars["MetaData/satelliteAscendingFlag"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["flagColdCal"]}, + VariableInfo{.var = p.out.vars["MetaData/coldCalibrationFlag"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["flagDayNight"], .range = only_ch1}, + VariableInfo{.var = p.out.vars["MetaData/dayOrNightQualifier"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["flagManeuver"], .range = only_ch1}, + VariableInfo{.var = p.out.vars["MetaData/spacecraftManeuverFlag"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["flagNonOcean"], .range = only_ch1}, + VariableInfo{.var = p.out.vars["MetaData/nonOceanFlag"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["flagPLOrientation"], .range = only_ch1}, + VariableInfo{.var = p.out.vars["MetaData/payloadOrientationFlag"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["brightness_temperature"]}, + VariableInfo{.var = p.out.vars["ObsValue/brightnessTemperature"], .range = offset_start}); + + // Note (RH): + // For user convenience, we could introduce a "compositePreQCFlag" that merges the flags in + // ICTCalibrationFlag, NDCalibrationFlag, coldCalibrationFlag, spacecraftManeuverFlag, + // lunarIntrusionFlag, and solarIntrusionFlag. + // This flag is constructed by reading in the relevant variables and running a bitwise OR on the contents. + /* + { + VariableInfo vartest{.var = p.out.vars["MetaData/spacecraftManeuverFlag"], + .range = offset_start}; + VariableDerivedInfo vartest_info(vartest); + const size_t sz = vartest_info.selection_num_elements; + valarray compositePreQCFlag_data(sz); + vector buffer(sz); + + const vector qcvars{ + {"MetaData/ICTCalibrationFlag"}, {"MetaData/lunarIntrusionFlag"}, + {"MetaData/NDCalibrationFlag"}, {"MetaData/solarIntrusionFlag"}, + {"MetaData/coldCalibrationFlag"}, {"MetaData/spacecraftManeuverFlag"}}; + for (const auto &qcvar : qcvars) { + std::cout << qcvar << std::endl; + p.out.vars[qcvar].read( + gsl::make_span(buffer.data(), sz), // Read into this buffer + // No need for a type conversion. These QC flags are just a stream of single-byte fields. + vartest_info.type, + vartest_info.selection_membuf, // Where should this be put in "buffer"? + vartest_info.selection_ioda // The hyperslab to be read + ); + compositePreQCFlag_data = compositePreQCFlag_data & valarray(buffer.data(), sz); + } + vector outvar(begin(compositePreQCFlag_data), end(compositePreQCFlag_data)); + p.out.vars["MetaData/compositePreQCFlag"].write( + gsl::make_span(outvar.data(), sz), vartest_info.type, vartest_info.selection_membuf, + vartest_info.selection_ioda); + } + */ +} + +void convert(TMS_Type typ, Converter_Params &p) { + using namespace std; + using namespace ioda; + + // Note that the ordering of the variable dimensions is inconsistent between + // TROPICS and TMS products. + // TROPICS: channels x scans x spots + // TMS v6: spots x scans x channels + // Outgoing: Location x Channel + // So long as all of the variables are treated in + // the same way within a product, we are okay. TMS will be written with locations in + // spot x scan ordering, while tropics will be written in scan x spot. + // The subsidiary functions will handle any discrepancies when relevant. This + // only applies for manually reconstructed fields like dateTime. + + // This specifies that only data for the first channel or band should be copied. + const DimensionRanges_t only_ch1_TMS{{}, {}, {{0, 0}}}; + const DimensionRanges_t only_band1_TROPICS{{{0, 0}}, {}, {}}; + const DimensionRanges_t single_ch{(typ == TMS_Type::TIO) ? only_ch1_TMS : only_band1_TROPICS}; + + // This specifies that we are copying to a certain range in the output variable. + const DimensionRanges_t offset_start{{{p.locationOffset, p.locationOffset + p.inNumLocs - 1}}}; + + // MetaData + + // datetime + convert_datetime(typ, p); + + // Satellite instrument specs are populated elsewhere. No copying needed. + + // Orbit geometry + + copy(VariableInfo{.var = p.in.vars["latitude"], .range = single_ch}, + VariableInfo{.var = p.out.vars["MetaData/latitude"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["longitude"], .range = single_ch}, + VariableInfo{.var = p.out.vars["MetaData/longitude"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["sensor_azimuth_angle"], .range = single_ch}, + VariableInfo{.var = p.out.vars["MetaData/sensorAzimuthAngle"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["sensor_view_angle"], .range = single_ch}, + VariableInfo{.var = p.out.vars["MetaData/sensorViewAngle"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["sensor_zenith_angle"], .range = single_ch}, + VariableInfo{.var = p.out.vars["MetaData/sensorZenithAngle"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["solar_azimuth_angle"], .range = single_ch}, + VariableInfo{.var = p.out.vars["MetaData/solarAzimuthAngle"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["solar_zenith_angle"], .range = single_ch}, + VariableInfo{.var = p.out.vars["MetaData/solarZenithAngle"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["lunar_azimuth_angle"], .range = single_ch}, + VariableInfo{.var = p.out.vars["MetaData/lunarAzimuthAngle"], .range = offset_start}); + + copy(VariableInfo{.var = p.in.vars["lunar_zenith_angle"], .range = single_ch}, + VariableInfo{.var = p.out.vars["MetaData/lunarZenithAngle"], .range = offset_start}); + + // Diagnostic flags + + // ObsValue is handled in the TIO and TOPICS-specific conversion section. + + // No ObsError assignments + + // TIO and TROPICS-specific conversions + if (typ == TMS_Type::TIO) + convert_tio(p); + else + convert_tropics(p); +} + +} // namespace tio_converter diff --git a/src/tomorrow_io/copy_tms_datetime.cc b/src/tomorrow_io/copy_tms_datetime.cc new file mode 100644 index 00000000..89cb6a68 --- /dev/null +++ b/src/tomorrow_io/copy_tms_datetime.cc @@ -0,0 +1,122 @@ +/* + * (C) Copyright 2024 The Tomorrow Companiec, Inc. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ +#include +#include // setenv +#include // std::put_time +#include +#include +#include + +#include "product.h" + +#if __cplusplus < 202002L +/// \brief Stringify a time point. +/// \note This code will be superseded with the switch to C++20, as +/// C++20 adds ostream operators for all time types. +template +std::ostream &operator<<(std::ostream &os, const std::chrono::time_point &tp) { + const std::time_t t = Clock_T::to_time_t(tp); + os << std::put_time(std::gmtime(&t), ioda::Types::Chrono_Time_Format); + return os; +} +#endif + +typedef int64_t Chrono_Time_Rep_t; +typedef std::ratio<1, 1> Chrono_Time_Period_t; +typedef std::chrono::system_clock Chrono_Clock_t; +typedef std::chrono::duration Chrono_Duration_t; +typedef std::chrono::time_point Chrono_Time_Point_t; + +namespace tio_converter { + +namespace { +Chrono_Time_Point_t makeTimePoint(uint16_t year, uint16_t month, uint16_t dayofmonth, uint16_t hour, + uint16_t min, uint16_t sec) { + struct std::tm timeTemp {}; + timeTemp.tm_year = year - 1900; + timeTemp.tm_mon = month - 1; + timeTemp.tm_mday = dayofmonth; + timeTemp.tm_hour = hour; + timeTemp.tm_min = min; + timeTemp.tm_sec = sec; + timeTemp.tm_isdst = 0; + std::time_t tt = std::mktime(&timeTemp); + return std::chrono::time_point_cast( + std::chrono::system_clock::from_time_t(tt)); +} + +std::string stringifyEpoch(const Chrono_Time_Point_t &epoch) { + auto tt = std::chrono::system_clock::to_time_t(epoch); + struct std::tm *epochTime; + epochTime = std::gmtime(&tt); + char epochString[21]; + size_t ret = std::strftime(epochString, 21, "%Y-%m-%dT%H:%M:%SZ", epochTime); + if (!ret) throw std::logic_error("Invalid time"); + return std::string("seconds since ") + std::string(epochString); +} + +} // namespace + +void convert_datetime(TMS_Type typ, Converter_Params &p) { + using namespace ioda; + using namespace std; + + // Force UTC time zone. Needed for time calculations before C++20. + setenv("TZ", "UTC", 1); + + Variable datetime = p.out.vars["MetaData/dateTime"]; + auto epoch = makeTimePoint(1970, 1, 1, 0, 0, 0); + + // dateTime + vector scantimestamps(p.inNumScans); + { + vector hour, minute, second, month, day; + vector year; + // IODA's type system does not support 8-bit unsigned integers because these are interpreted as + // 'unsigned char' in C++, and not all compilers have support. + p.in.vars["Hour"].read(hour); + p.in.vars["Minute"].read(minute); + p.in.vars["Second"].read(second); + p.in.vars["Month"].read(month); + p.in.vars["Day"].read(day); + p.in.vars["Year"].read(year); + for (size_t i = 0; i < p.inNumScans; ++i) { + scantimestamps[i] = makeTimePoint(year[i], month[i], day[i], hour[i], minute[i], second[i]); + + // For debugging: + // if (!i) { + // std::cout << "First time point is " << year[0] << "-" << month[0] << "-" << day[0] << " " + // << hour[0] << ":" << minute[0] << ":" << second[0] << std::endl; + // } + } + + vector locationtimestamps_as_seconds(p.inNumLocs); + for (size_t scan = 0; scan < p.inNumScans; ++scan) + for (size_t spot = 0; spot < p.inNumSpots; ++spot) { + // The locations need to align to typical variable coordinate ordering in the file. + // We use the brightness_temperature variable as the basis for this ordering. + // TMS dims are spot x scan x ch. + // TROPICS dims are ch x scan x spot. + // The relative ordering of spot and scan are key. + size_t i = (typ == TMS_Type::TIO) ? (spot * p.inNumScans) + scan // T.io TMS + : (scan * p.inNumSpots) + spot; // TROPICS TMS + locationtimestamps_as_seconds[i] + = std::chrono::duration_cast(scantimestamps[scan] - epoch).count(); + } + + // time stamps are for each scan line, but need to map to each location + vector sel_zero{0}, sel_start{p.locationOffset}, sel_count{p.inNumLocs}; + + ioda::Selection selection_ioda, selection_membuf; + selection_ioda.extent(p.out.vars["Location"].getDimensions().dimsCur) + .select({SelectionOperator::SET, sel_start, sel_count}); + selection_membuf.extent({p.inNumLocs}).select({SelectionOperator::SET, sel_zero, sel_count}); + + datetime.write(locationtimestamps_as_seconds, selection_membuf, selection_ioda); + } +} +} // namespace tio_converter diff --git a/src/tomorrow_io/creation.cc b/src/tomorrow_io/creation.cc new file mode 100644 index 00000000..dea23727 --- /dev/null +++ b/src/tomorrow_io/creation.cc @@ -0,0 +1,230 @@ +/* + * (C) Copyright 2024 The Tomorrow Companiec, Inc. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ +#include +#include +#include + +#include "product.h" + +namespace tio_converter { + +void create_vars_tropics(ioda::ObsGroup &out, ioda::Variable Location, ioda::Variable Channel, + ioda::VariableCreationParameters uint8_params, + ioda::VariableCreationParameters float_params) { + using namespace ioda; + using namespace std; + // The TROPICS calQualityFlag is channels * scans * spots. + // We need to break up the variable's contents into several sub-variables. + // Several of these are covered in the general TIO/TROPICS function. + // | IODA Name | TROPICS Bit Pos | Meaning | + // | ------------------- | --------------- | ------- | + // | flagNonOcean_legacy | 1 | Is this ocean or not? | + // | flagIntrusion | 2 | Solar or Lunar intrusion (TROPICS only) | + // | flagManeuver | 3 | Is spacecraft repositioning? | + // | flagColdCal | 4 | Cold caibration inconsistency | + // | flagHotCal | 5 | Hot calibration inconsistency | + // | flagAscDesc | 6 | Ascending / Descending flag. 0: asc, 1: desc. | + // | flagDayNight | 7 | Day / night flag. 0: day, 1: night. | + // | flagPLOrientation | 8 | Payload orientation (0: fore / 1: aft) | + + Variable flagNonOcean_legacy = out.vars.createWithScales( + "MetaData/nonOceanFlag_legacy", {Location, Channel}, uint8_params); + + Variable flagIntrusion = out.vars.createWithScales("MetaData/intrusionFlag", + {Location, Channel}, uint8_params); + flagIntrusion.atts.add("description", string("Set to 1 if potential solar/lunar intrusion.")); + + Variable flagHotCal = out.vars.createWithScales("MetaData/hotCalibrationFlag", + {Location, Channel}, uint8_params); + flagHotCal.atts.add("description", string("Set to 1 if hot calibration is inconsistent.")); + + Variable compositePreQCFlag = out.vars.createWithScales( + "MetaData/compositePreQCFlag", {Location, Channel}, uint8_params); + compositePreQCFlag.atts.add( + "description", + std::string("Set to 1 if this observation has failed one or more upstream checks.")); +} + +void create_vars_tio(ioda::ObsGroup &out, ioda::Variable Location, ioda::Variable Channel, + ioda::VariableCreationParameters uint8_params, + ioda::VariableCreationParameters float_params) { + using namespace ioda; + using namespace std; + Variable MultiMask + = out.vars.createWithScales("MetaData/multiMask", {Location}, uint8_params); + MultiMask.atts.add("description", + std::string("1=clear land; 2=cloudy ocean; 3=clear ocean, 4=cloudy land")); + MultiMask.atts.add("long_name", + std::string("Combined land/ocean/cloud mask derived from channel 1.")); + + Variable clear_fraction = out.vars.createWithScales("MetaData/cloudClearFraction", + {Location, Channel}, float_params); + clear_fraction.atts.add( + "description", std::string("Cloud-clear fraction of each spot from time-matched geostationary " + "satellite cloud mask, weighted by each channel's antenna pattern")); + clear_fraction.atts.add("long_name", std::string("Cloud-clear fraction")); + clear_fraction.atts.add("units", std::string("1")); + + Variable flagICTCal = out.vars.createWithScales("MetaData/ICTCalibrationFlag", + {Location, Channel}, uint8_params); + flagICTCal.atts.add("description", + std::string("Outlier detection flag for internal calibration target spots.")); + flagICTCal.atts.add("long_name", std::string("Internal Calibration Target Spot Flag")); + + Variable flagLunarIntrusion = out.vars.createWithScales( + "MetaData/lunarIntrusionFlag", {Location, Channel}, uint8_params); + flagLunarIntrusion.atts.add("description", + std::string("True if there is a lunar intrusion into the cold space " + "or noise diode calibration sectors.")); + flagLunarIntrusion.atts.add("long_name", std::string("Lunar intrusion flag")); + + Variable flagNDCal = out.vars.createWithScales("MetaData/NDCalibrationFlag", + {Location, Channel}, uint8_params); + flagNDCal.atts.add("description", + std::string("Outlier detection flag for noise diode calibration spots.")); + flagNDCal.atts.add("long_name", std::string("Noise Diode Calibration Spot Flag")); + + Variable flagSolarIntrusion = out.vars.createWithScales( + "MetaData/solarIntrusionFlag", {Location, Channel}, uint8_params); + flagSolarIntrusion.atts.add("description", + std::string("True if there is a solar intrusion into the cold space " + "or noise diode calibration sectors.")); + flagSolarIntrusion.atts.add("long_name", std::string("Solar intrusion flag")); + + Variable land_fraction = out.vars.createWithScales("MetaData/landAreaFraction", + {Location, Channel}, float_params); + land_fraction.atts.add("description", + std::string("Land fraction weighted by each channel's antenna pattern")); + land_fraction.atts.add("long_name", std::string("Land fraction")); + land_fraction.atts.add("units", std::string("1")); +} + +void create_vars(TMS_Type typ, int sat_bufr_id, int sat_sub_bufr_id, int inst_bufr_id, + ioda::ObsGroup &out, ioda::Variable Location, ioda::Variable Channel) { + using namespace ioda; + VariableCreationParameters datetime_params; + VariableCreationParameters float_params; + VariableCreationParameters uint8_params; + + datetime_params.chunk = true; + datetime_params.compressWithGZIP(); + datetime_params.setFillValue(0); // Must be set for ioda to read this variable. + + float_params.chunk = true; + float_params.compressWithGZIP(); + float_params.setFillValue(-999); + + uint8_params.chunk = true; + uint8_params.compressWithGZIP(); + uint8_params.setFillValue(255); + + // MetaData + + // datetime + + Variable datetime + = out.vars.createWithScales("MetaData/dateTime", {Location}, datetime_params); + datetime.atts.add("units", std::string("seconds since 1970-01-01T00:00:00Z")); + + // Orbit geometry + + Variable latitude + = out.vars.createWithScales("MetaData/latitude", {Location}, float_params); + latitude.atts.add("units", std::string("degrees_north")); + + Variable longitude + = out.vars.createWithScales("MetaData/longitude", {Location}, float_params); + longitude.atts.add("units", std::string("degrees_east")); + + Variable sensorAzimuthAngle + = out.vars.createWithScales("MetaData/sensorAzimuthAngle", {Location}, float_params); + sensorAzimuthAngle.atts.add("units", std::string("degrees")); + + Variable sensorViewAngle + = out.vars.createWithScales("MetaData/sensorViewAngle", {Location}, float_params); + sensorViewAngle.atts.add("units", std::string("degrees")); + + Variable sensorZenithAngle + = out.vars.createWithScales("MetaData/sensorZenithAngle", {Location}, float_params); + sensorZenithAngle.atts.add("units", std::string("degrees")); + + Variable solarAzimuthAngle + = out.vars.createWithScales("MetaData/solarAzimuthAngle", {Location}, float_params); + solarAzimuthAngle.atts.add("units", std::string("degrees")); + + Variable solarZenithAngle + = out.vars.createWithScales("MetaData/solarZenithAngle", {Location}, float_params); + solarZenithAngle.atts.add("units", std::string("degrees")); + + Variable lunarAzimuthAngle + = out.vars.createWithScales("MetaData/lunarAzimuthAngle", {Location}, float_params); + lunarAzimuthAngle.atts.add("units", std::string("degrees")); + + Variable lunarZenithAngle + = out.vars.createWithScales("MetaData/lunarZenithAngle", {Location}, float_params); + lunarZenithAngle.atts.add("units", std::string("degrees")); + + // Diagnostic flags + + Variable flagAscDesc = out.vars.createWithScales("MetaData/satelliteAscendingFlag", + {Location}, uint8_params); + flagAscDesc.atts.add("description", std::string("True if in descending portion of orbit.")); + flagAscDesc.atts.add("long_name", std::string("Ascending/Descending flag")); + + Variable flagColdCal = out.vars.createWithScales("MetaData/coldCalibrationFlag", + {Location, Channel}, uint8_params); + flagColdCal.atts.add("description", + std::string("Outlier detection flag for deep space calibration spots.")); + flagColdCal.atts.add("long_name", std::string("Cold Calibration Spot Flag")); + + Variable flagDayNight + = out.vars.createWithScales("MetaData/dayOrNightQualifier", {Location}, uint8_params); + flagDayNight.atts.add("description", + std::string("True if earth is between the sun and spacecraft.")); + flagDayNight.atts.add("long_name", std::string("Day/Night flag")); + + Variable flagManeuver = out.vars.createWithScales("MetaData/spacecraftManeuverFlag", + {Location}, uint8_params); + flagManeuver.atts.add("description", + std::string("True if the spacecraft is in an active maneuver.")); + flagManeuver.atts.add("long_name", std::string("Spacecraft maneuver flag")); + + Variable flagNonOcean + = out.vars.createWithScales("MetaData/nonOceanFlag", {Location}, uint8_params); + flagNonOcean.atts.add( + "description", std::string("0 is ocean, 1 is land, coastline, or undefined. Uses ch 1 spot.")); + flagNonOcean.atts.add("long_name", std::string("Non-ocean Flag")); + + Variable flagPLOrientation = out.vars.createWithScales("MetaData/payloadOrientationFlag", + {Location}, uint8_params); + flagPLOrientation.atts.add("description", + std::string("True if spacecraft is flying payload-first.")); + flagPLOrientation.atts.add("long_name", std::string("Payload Orientation flag")); + + // ObsValue and ObsError + + Variable brightnessTemperature = out.vars.createWithScales( + "ObsValue/brightnessTemperature", {Location, Channel}, float_params); + brightnessTemperature.atts.add("units", std::string("K")); + // ObsError/brightnessTemperature is needed for ioda to read the file. We fill with dummy values since + // we will assign error upstream in later processing. + Variable brightnessTemperatureError = out.vars.createWithScales( + "ObsError/brightnessTemperature", {Location, Channel}, float_params); + brightnessTemperatureError.atts.add("units", std::string("K")); + + // Satellite instrument specs + const size_t numLocs = Location.getDimensions().dimsCur[0]; + const size_t numScans = numLocs / 81; + create_sat_inst_specs(typ, sat_bufr_id, sat_sub_bufr_id, inst_bufr_id, out, Location, Channel, + numScans); + + if (typ == TMS_Type::TIO) create_vars_tio(out, Location, Channel, uint8_params, float_params); + if (typ == TMS_Type::TROPICS) + create_vars_tropics(out, Location, Channel, uint8_params, float_params); +} + +} // namespace tio_converter diff --git a/src/tomorrow_io/creation_tio_sat_inst_specs.cc b/src/tomorrow_io/creation_tio_sat_inst_specs.cc new file mode 100644 index 00000000..a0ae3c30 --- /dev/null +++ b/src/tomorrow_io/creation_tio_sat_inst_specs.cc @@ -0,0 +1,74 @@ +/* + * (C) Copyright 2024 The Tomorrow Companiec, Inc. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ +#include +#include + +#include "product.h" + +namespace tio_converter { +void create_sat_inst_specs(TMS_Type typ, int sat_bufr_id, int sat_sub_bufr_id, int inst_bufr_id, + ioda::ObsGroup &out, ioda::Variable Location, ioda::Variable Channel, + size_t numScans) { + using namespace ioda; + using namespace std; + const size_t numLocs = Location.getDimensions().dimsCur[0]; + + vector sensorCentralFrequency_GHz + = (typ == TMS_Type::TROPICS) + // TROPICS + ? vector{91.655, 114.50, 115.95, 116.65, 117.25, 117.80, + 118.24, 118.58, 184.41, 186.51, 190.31, 204.80} // TIO + : vector{91.65, 118.75, 118.75, 118.75, 118.75, 118.75, + 118.75, 118.75, 184.41, 186.51, 190.31, 204.8}; + + Variable sensorChannelNumber + = out.vars.createWithScales("MetaData/sensorChannelNumber", {Channel}); + Variable sensorCentralFrequency + = out.vars.createWithScales("MetaData/sensorCentralFrequency", {Channel}); + sensorCentralFrequency.atts.add("units", string("Hz")); + + // Custom missing value for bufr id fields. + VariableCreationParameters id_params; + id_params.chunk = true; + id_params.compressWithGZIP(); + id_params.setFillValue(999); + + // satelliteIdentifier + vector vSatelliteIdentifier(numLocs, sat_bufr_id); + Variable satelliteIdentifier + = out.vars.createWithScales("MetaData/satelliteIdentifier", {Location}, id_params); + satelliteIdentifier.write(vSatelliteIdentifier); + + // satelliteSubIdentifier + vector vSatelliteSubIdentifier(numLocs, sat_sub_bufr_id); + Variable satelliteSubIdentifier + = out.vars.createWithScales("MetaData/satelliteSubIdentifier", {Location}, id_params); + satelliteSubIdentifier.write(vSatelliteSubIdentifier); + + // instrumentIdentifier + vector vSatelliteInstrument(numLocs, inst_bufr_id); + Variable satelliteInstrument + = out.vars.createWithScales("MetaData/instrumentIdentifier", {Location}, id_params); + satelliteInstrument.write(vSatelliteInstrument); + + // sensorChannelNumber is a sequence if ints from 1 to number of channels + vector chnum(12); + iota(chnum.begin(), chnum.end(), 1); + sensorChannelNumber.write(chnum); + Channel.write(chnum); + // sensorCentralFrequency + using namespace ioda::udunits; + auto GHzToHz = Units("GHz").getConverterTo(Units("Hz")); + vector cfhz(sensorCentralFrequency_GHz.size()); + GHzToHz->convert(sensorCentralFrequency_GHz.data(), sensorCentralFrequency_GHz.size(), + cfhz.data()); + sensorCentralFrequency.write(cfhz); + + // NOTE (RH): sensorScanPosition is deliberately not written. Users + // should instead use sensorViewAngle and sensorZenithAngle when filtering. +} +} // namespace tio_converter diff --git a/src/tomorrow_io/detect_tms_type.cc b/src/tomorrow_io/detect_tms_type.cc new file mode 100644 index 00000000..2ab6cf8b --- /dev/null +++ b/src/tomorrow_io/detect_tms_type.cc @@ -0,0 +1,57 @@ +/* + * (C) Copyright 2024 The Tomorrow Companiec, Inc. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ +#include +#include +#include + +#include "ioda/Engines/HH.h" +#include "ioda/Group.h" +#include "product.h" +namespace tio_converter { + +void detect_tms(const std::string &input_file, TMS_Type &typ, int &sat_bufr_id, + int &sat_sub_bufr_id, int &inst_bufr_id) { + using namespace ioda; + using namespace std; + typ = TMS_Type::TIO; + sat_bufr_id = 999; + sat_sub_bufr_id = 999; + inst_bufr_id = 999; + + // Query the input file to get the satellite identifier. + Group in = Engines::HH::openFile(input_file, Engines::BackendOpenModes::Read_Only); + + if (in.atts.exists("Source")) { + string source_sat = in.atts["Source"].read(); + if (source_sat.substr(0, 7) == "TROPICS") { + typ = TMS_Type::TROPICS; + inst_bufr_id = 433; // See https://codes.ecmwf.int/odb/satelliteinstrument/. + const map tropics_ids = {{"TROPICS01", 709}, + {"TROPICS03", 228}, + {"TROPICS05", 263}, + {"TROPICS06", 264}, + {"TROPICS07", 284}}; + if (tropics_ids.count(source_sat)) sat_bufr_id = tropics_ids.at(source_sat); + } + } else if (in.atts.exists("Filename")) { + string source_filename = in.atts["Filename"].read(); + if (source_filename.substr(0, 3) == "TMS") { + // NOTE (RH): WMO BUFR values are not yet assigned. Placeholder values + // are used instead. These values were conveyed to Ben Ruston by Jeff Ator, + // and are consistent with usage in the JCSDA converter. + sat_bufr_id = 769; // Placeholder! + inst_bufr_id = 1100; // Purely a placeholder value. + // TMS-S02 is the fourth T.io satellite and the only one used with this converter + // so far. Expect code updates once more instruments launch. + sat_sub_bufr_id = 4; // Placeholder! + } + } else { + throw domain_error("Input file is not a TROPICS or TMS L1B data file."); + } +} + +} // namespace tio_converter diff --git a/src/tomorrow_io/main_tio.cc b/src/tomorrow_io/main_tio.cc new file mode 100644 index 00000000..54affd4b --- /dev/null +++ b/src/tomorrow_io/main_tio.cc @@ -0,0 +1,52 @@ +/* + * (C) Copyright 2024 The Tomorrow Companiec, Inc. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ +#include +#include +#include + +#include "ioda/Engines/HH.h" +#include "ioda/ObsGroup.h" +#include "ioda/Units.h" +#include "product.h" + +void doHelp() { + using namespace std; + cout << "Usage: convert_tio_tms.x input_file_1 [input_file_2 ...] output_file\n"; + exit(1); +} + +int main(int argc, char **argv) { + using namespace std; + using namespace ioda; + using namespace tio_converter; + try { + if (argc < 3) doHelp(); + vector input_files; + for (int i = 1; i + 1 < argc; ++i) input_files.push_back(argv[i]); + if (!input_files.size()) doHelp(); + + TMS_Type typ; + int sat_bufr_id, sat_sub_bufr_id, inst_bufr_id; + detect_tms(input_files[0], typ, sat_bufr_id, sat_sub_bufr_id, inst_bufr_id); + + string sOutFile(argv[argc - 1]); + + ioda::ObsGroup out + = prepare_output_file(typ, sat_bufr_id, sat_sub_bufr_id, inst_bufr_id, input_files, sOutFile); + size_t location_start = 0; + for (const auto &input_file : input_files) { + Group in = Engines::HH::openFile(input_file, Engines::BackendOpenModes::Read_Only); + Converter_Params cp(in, out, location_start); + convert(typ, cp); + location_start += cp.inNumLocs; + } + return 0; + } catch (const std::exception &e) { + std::cerr << e.what() << endl; + return 1; + } +} diff --git a/src/tomorrow_io/product.cc b/src/tomorrow_io/product.cc new file mode 100644 index 00000000..6e9370a6 --- /dev/null +++ b/src/tomorrow_io/product.cc @@ -0,0 +1,65 @@ +/* + * (C) Copyright 2024 The Tomorrow Companiec, Inc. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ +#include "product.h" + +#include +#include +#include + +#include "ioda/Engines/HH.h" + +namespace tio_converter { + +Converter_Params::Converter_Params(const ioda::Group &src, ioda::ObsGroup &dest, int locationOffset) + : in{src}, out{dest} { + // Query the input file for the sizes of key dimensions + this->inNumLocs + = in.vars["scans"].getDimensions().dimsCur[0] * in.vars["spots"].getDimensions().dimsCur[0]; + this->inNumChans = in.vars["channels"].getDimensions().dimsCur[0]; + this->inNumScans = in.vars["scans"].getDimensions().dimsCur[0]; + this->inNumSpots = in.vars["spots"].getDimensions().dimsCur[0]; + + this->locationOffset = locationOffset; +} + +ioda::ObsGroup prepare_output_file(TMS_Type typ, int sat_bufr_id, int sat_sub_bufr_id, + int inst_bufr_id, const std::vector &input_files, + const std::string &output_file) { + // Open each of the input files and get the number of locations within the file. + using namespace ioda; + const size_t numChans = 12; + size_t numLocs = 0; + for (const auto &input_filename : input_files) { + Group in = Engines::HH::openFile(input_filename, Engines::BackendOpenModes::Read_Only); + const size_t inNumScans = in.vars["scans"].getDimensions().dimsCur[0]; + const size_t inNumSpots = in.vars["spots"].getDimensions().dimsCur[0]; + numLocs += inNumScans * inNumSpots; + } + + // Create the output file. + // Make a new observation space and add in appropriate dimensions and variables + NewDimensionScales_t newDims; + // 81 scan positions per line. 400 lines per file is typical. Each file should have multiple chunks + // to avoid excessive open / close operations when constructing the ObsGroup. + const size_t maxNumLocsPerChunk = 10000; + newDims.push_back(NewDimensionScale( + "Location", numLocs, numLocs, (numLocs > maxNumLocsPerChunk) ? maxNumLocsPerChunk : numLocs)); + newDims.push_back(NewDimensionScale("Channel", numChans, numChans, numChans)); + + Group out_group + = Engines::HH::createFile(output_file, Engines::BackendCreateModes::Truncate_If_Exists, + {Engines::HH::HDF5_Version::V18, Engines::HH::HDF5_Version::Latest}); + + ioda::ObsGroup out = ObsGroup::generate(out_group, newDims); + ioda::Variable Location = out.vars["Location"]; + ioda::Variable Channel = out.vars["Channel"]; + create_vars(typ, sat_bufr_id, sat_sub_bufr_id, inst_bufr_id, out, Location, Channel); + + return out; +} + +} // end namespace tio_converter diff --git a/src/tomorrow_io/product.h b/src/tomorrow_io/product.h new file mode 100644 index 00000000..2ab1ab9c --- /dev/null +++ b/src/tomorrow_io/product.h @@ -0,0 +1,76 @@ +#pragma once +/* + * (C) Copyright 2024 The Tomorrow Companiec, Inc. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + */ +#include +#include +#include +#include +#include + +#include "ioda/ObsGroup.h" + +/// @brief The namespace for the Tomorrow.io TMS converter. +namespace tio_converter { + +/// @brief Convenience structure passed to the data copying functions. +struct Converter_Params { + const ioda::Group ∈ ///< The input file + ioda::ObsGroup out; ///< The output file + + ioda::Variable Location; ///< The output file's Location dimension scale. + ioda::Variable Channel; ///< The output file's Channel dimension scale. + int inNumLocs; ///< Number of locations in the input file (numScans * numSpots). + int inNumChans; ///< Number of channels in the input file. + int inNumScans; ///< Number of scan lines in the input file. + int inNumSpots; ///< Number of scan spots per line in the input file. + int locationOffset; ///< An offset value used when concatenating multiple inputs. + + Converter_Params(const ioda::Group &src, ioda::ObsGroup &dest, int locationOffset = 0); +}; + +/// @brief Describes whether the input is a Tio TMS instrument or TROPICS. +enum class TMS_Type { TIO, TROPICS }; + +/// @brief Convenience function to read an input file and infer some basic information. +/// @param input_file is the path to the L1B file to be tested. +/// @param[out] typ indicates whether the data are from Tio or TROPICS. +/// @param[out] sat_bufr_id is the WMO BUFR id for this satellite. +/// @param[out] sat_sub_bufr_id is the sub-satellite id for this satellite. +/// @param[out] inst_bufr_id is the WMO BUFR id for this instrument. +void detect_tms(const std::string &input_file, TMS_Type &typ, int &sat_bufr_id, + int &sat_sub_bufr_id, int &inst_bufr_id); + +/// @brief Create the skeleton of the output IODA file. This will be an ObsGroup +/// with the corrrect number of Locations, plus the correct information regarding +/// satellite / instrument IDs. +/// @param typ indicates whether the data are from Tio or TROPICS. +/// @param sat_bufr_id is the WMO BUFR id for this satellite. +/// @param sat_sub_bufr_id is the sub-satellite id for this satellite. +/// @param inst_bufr_id is the WMO BUFR id for this instrument. +/// @param input_files are the paths of the input L1B files. +/// @param output_file is the path to the output file. +/// @see detect_tms to get the ID fields. +ioda::ObsGroup prepare_output_file(TMS_Type typ, int sat_bufr_id, int sat_sub_bufr_id, + int inst_bufr_id, const std::vector &input_files, + const std::string &output_file); + +/// @brief Function to copy data from an input file to the destination. +void convert(TMS_Type typ, Converter_Params &); + +// Private functions + +/// @brief Function called by prepare_output_file to create the variables in an output file. +void create_vars(TMS_Type typ, int sat_bufr_id, int sat_sub_bufr_id, int inst_bufr_id, + ioda::ObsGroup &out, ioda::Variable Location, ioda::Variable Channel); +/// @brief Function called by prepare_output_file to fill in more variables in the output file. +void create_sat_inst_specs(TMS_Type typ, int sat_bufr_id, int sat_sub_bufr_id, int inst_bufr_id, + ioda::ObsGroup &out, ioda::Variable Location, ioda::Variable Channel, + size_t numScans); +/// @brief Function to set datetimes in the IODA file. Called by "convert". +void convert_datetime(TMS_Type typ, Converter_Params &p); + +} // end namespace tio_converter diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 6b242c1d..37375b21 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -104,6 +104,8 @@ list( APPEND test_input testinput/pandora_site_classification.csv testinput/balloonsonde_sample_20241001T00Z_PT15M.json testinput/tec_groundBased_TENET_20201001T00Z.tec + testinput/TMS02.1B-TB.V00-07.NRT.ST20241115-125911.ET20241115-130757.CT20241119-074232.nc + testinput/TROPICS03.BRTT.L1B.Orbit07819.V05-02.ST20241023-003347.ET20241023-020756.CT20241023-145302.nc ) list( APPEND test_output @@ -200,6 +202,8 @@ list( APPEND test_output testoutput/Pandora57s1_BoulderCO_L2_rnvs3p1-8.nc testoutput/obs.20241001T000000Z_PT15M_balloonsonde.nc4 testoutput/obs.20201001T000000Z_PT1M_tec_groundbased.nc4 + testoutput/tio_converter_tio_tms_l1b_single_V00-07.nc + testoutput/tio_converter_tropics03.nc ) if( iodaconv_gnssro_ENABLED ) @@ -2478,6 +2482,27 @@ ecbuild_add_test( TARGET test_${PROJECT_NAME}_madis_snod -o testrun/madis_snod_2021010100.nc" madis_snod_2021010100.nc ${IODA_CONV_COMP_TOL}) +ecbuild_add_test( TARGET test_${PROJECT_NAME}_tio_conv_tms_tio + TYPE SCRIPT + COMMAND bash + ARGS ${CMAKE_BINARY_DIR}/bin/iodaconv_comp.sh + netcdf + "${CMAKE_BINARY_DIR}/bin/convert_tio.x + testinput/TMS02.1B-TB.V00-07.NRT.ST20241115-125911.ET20241115-130757.CT20241119-074232.nc + testrun/tio_converter_tio_tms_l1b_single_V00-07.nc" + tio_converter_tio_tms_l1b_single_V00-07.nc ${IODA_CONV_COMP_TOL}) + +ecbuild_add_test( TARGET test_${PROJECT_NAME}_tio_conv_tms_tropics + TYPE SCRIPT + COMMAND bash + ARGS ${CMAKE_BINARY_DIR}/bin/iodaconv_comp.sh + netcdf + "${CMAKE_BINARY_DIR}/bin/convert_tio.x + testinput/TROPICS03.BRTT.L1B.Orbit07819.V05-02.ST20241023-003347.ET20241023-020756.CT20241023-145302.nc + testrun/tio_converter_tropics03.nc" + tio_converter_tropics03.nc ${IODA_CONV_COMP_TOL}) + + ####################################################### # Auto-generated tests for ioda file validation # For these tests we check the testoutput directory diff --git a/test/testinput/TMS02.1B-TB.V00-07.NRT.ST20241115-125911.ET20241115-130757.CT20241119-074232.nc b/test/testinput/TMS02.1B-TB.V00-07.NRT.ST20241115-125911.ET20241115-130757.CT20241119-074232.nc new file mode 100644 index 00000000..b8e5e115 --- /dev/null +++ b/test/testinput/TMS02.1B-TB.V00-07.NRT.ST20241115-125911.ET20241115-130757.CT20241119-074232.nc @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5ec83a5ff387ec05b645ddb68a1002cf015a3c59ad92563563cc993b64b580e5 +size 3147982 diff --git a/test/testinput/TROPICS03.BRTT.L1B.Orbit07819.V05-02.ST20241023-003347.ET20241023-020756.CT20241023-145302.nc b/test/testinput/TROPICS03.BRTT.L1B.Orbit07819.V05-02.ST20241023-003347.ET20241023-020756.CT20241023-145302.nc new file mode 100644 index 00000000..a0587eee --- /dev/null +++ b/test/testinput/TROPICS03.BRTT.L1B.Orbit07819.V05-02.ST20241023-003347.ET20241023-020756.CT20241023-145302.nc @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8483c81449aefff5da77e19d473de9917afaa8afd171203ee27f8e254a5f823f +size 57562181 diff --git a/test/testoutput/tio_converter_tio_tms_l1b_single_V00-07.nc b/test/testoutput/tio_converter_tio_tms_l1b_single_V00-07.nc new file mode 100644 index 00000000..5a815ef8 --- /dev/null +++ b/test/testoutput/tio_converter_tio_tms_l1b_single_V00-07.nc @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:a85c2f446e8fb675074ecba93ff4f1a731d0f5bf50c1d2d1b7bae077f394db6a +size 1303287 diff --git a/test/testoutput/tio_converter_tropics03.nc b/test/testoutput/tio_converter_tropics03.nc new file mode 100644 index 00000000..e09fac5c --- /dev/null +++ b/test/testoutput/tio_converter_tropics03.nc @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1ffce01e562810d5f2b40dbc814fb07c06d6e1188f8eb79a357d02738fcca1ce +size 17828749