ERS PATC implementation

* PATC parsing and usage implementation. * Refactored and removed remnants of development. Updated SOFTWARE_VERSION field. * Unit tests for ISP adjustment and PATC register. Overflow condition.
cgi-estonia-space · Mar 15, 2024 · 8f8026c · 8f8026c
1 parent 99d3443
commit 8f8026c
Show file tree

Hide file tree

Showing 20 changed files with 429 additions and 36 deletions.
diff --git a/VERSION b/VERSION
@@ -1,5 +1,5 @@
 #define VERSION_MAJOR 0
-#define VERSION_MINOR 2
+#define VERSION_MINOR 3
 #define VERSION_PATCH 0
 
 #define STRINGIZE(s) #s

diff --git a/envisat_format/CMakeLists.txt b/envisat_format/CMakeLists.txt
@@ -11,8 +11,10 @@ list(APPEND ENVISAT_FORMAT_SOURCES
         ${CMAKE_CURRENT_LIST_DIR}/src/envisat_lvl0_parser.cc
         ${CMAKE_CURRENT_LIST_DIR}/src/envisat_lvl1_writer.cc
         ${CMAKE_CURRENT_LIST_DIR}/src/envisat_parse_util.cu
-        ${CMAKE_CURRENT_LIST_DIR}/src/ers_im_parse.cc
         ${CMAKE_CURRENT_LIST_DIR}/src/ers_aux_file.cc
+        ${CMAKE_CURRENT_LIST_DIR}/src/ers_im_parse.cc
+        ${CMAKE_CURRENT_LIST_DIR}/src/ers_sbt.cc
+        ${CMAKE_CURRENT_LIST_DIR}/src/patc.cc
 )
 
 add_library(envisat-format-static STATIC ${ENVISAT_FORMAT_SOURCES})

diff --git a/envisat_format/include/envisat_aux_file.h b/envisat_format/include/envisat_aux_file.h
@@ -10,12 +10,14 @@
 #pragma once
 
 #include <filesystem>
+#include <optional>
 #include <string_view>
 
 #include <boost/date_time/posix_time/posix_time.hpp>
 
 #include "bswap_util.h"
 #include "envisat_aux_xca.h"
+#include "patc.h"
 
 /**
  * Parser for ASA_CON and ASA_INS files, auxiliary files used for metadata extractions
@@ -715,4 +717,5 @@ void FindCONFile(std::string aux_root, boost::posix_time::ptime start, Configura
 
 namespace alus::asar::envformat::aux {
 void GetXca(std::string_view aux_root, boost::posix_time::ptime start, ExternalCalibration& xca, std::string& filename);
+std::optional<Patc> GetTimeCorrelation(std::string_view aux_root);
 }
diff --git a/envisat_format/include/envisat_lvl0_parser.h b/envisat_format/include/envisat_lvl0_parser.h
@@ -38,6 +38,8 @@ struct ASARMetadata {
     std::string orbit_dataset_name;
     std::string acquistion_station;
     std::string processing_station;
+    uint32_t rel_orbit;
+    uint32_t abs_orbit;
     boost::posix_time::ptime sensing_start;
     boost::posix_time::ptime sensing_stop;
     boost::posix_time::ptime first_line_time;

diff --git a/envisat_format/include/ers_env_format.h b/envisat_format/include/ers_env_format.h
@@ -22,6 +22,7 @@ constexpr size_t FEP_ANNOTATION_SIZE_BYTES{20};
 constexpr size_t MDSR_ISP_SIZE_BYTES{11466};  // Adding ISP sensing (12) and FEP annotations (20) one gets sum above.
 constexpr size_t DATA_RECORD_NUMBER_SIZE_BYTES{4};
 constexpr size_t IDHT_HEADER_SIZE_BYTES{10};
+constexpr size_t SBT_OFFSET_BYTES{PROCESSOR_ANNOTATION_ISP_SIZE_BYTES + FEP_ANNOTATION_SIZE_BYTES + 16};
 constexpr size_t AUXILIARY_REPLICA_CALIBRATION_SIZE_BYTES{220};
 constexpr size_t MEASUREMENT_DATA_SIZE_BYTES{11232};
 constexpr size_t MEASUREMENT_DATA_SAMPLE_COUNT{5616};

diff --git a/envisat_format/include/ers_sbt.h b/envisat_format/include/ers_sbt.h
@@ -0,0 +1,14 @@
+
+#pragma once
+
+#include <cstdint>
+
+#include "envisat_types.h"
+#include "patc.h"
+
+namespace alus::asar::envformat::ers {
+
+void RegisterPatc(const aux::Patc& patc);
+void AdjustIspSensingTime(mjd& isp_sensing_time, uint32_t sbt, uint32_t sbt_repeat, bool overflow = false);
+
+}  // namespace alus::asar::envformat::ers
diff --git a/envisat_format/include/patc.h b/envisat_format/include/patc.h
@@ -0,0 +1,23 @@
+
+#pragma once
+
+#include <cstdint>
+#include <filesystem>
+
+namespace alus::asar::envformat::aux {
+
+// Should not be under this module. An exception, as a surprise, there is a file format
+// which is NOT packaged into envisat format, nice.
+struct Patc {
+    uint16_t cyclic_counter;
+    char mission[2];
+    uint32_t orbit_number;
+    int32_t epoch_1950; // Days
+    int32_t milliseconds;
+    uint32_t sbt_counter;
+    uint32_t sbt_period;  // Tick time value in nanoseconds
+};
+
+Patc ParsePatc(const std::filesystem::path& file_loc);
+
+}  // namespace alus::asar::envformat::aux
diff --git a/envisat_format/src/envisat_aux_file.cc b/envisat_format/src/envisat_aux_file.cc
@@ -12,6 +12,8 @@
 #include "alus_log.h"
 #include "asar_aux.h"
 #include "checks.h"
+#include "filesystem_util.h"
+#include "patc.h"
 
 namespace {
 std::string DetermineFilePath(std::string aux_root, boost::posix_time::ptime start, alus::asar::aux::Type t) {
@@ -123,4 +125,26 @@ void GetXca(std::string_view aux_root, boost::posix_time::ptime start, ExternalC
 
     filename = std::filesystem::path(file_path).filename();
 }
+
+std::optional<Patc> GetTimeCorrelation(std::string_view aux_root) {
+    const std::regex patc_reg{R"(PATC.*)"};
+
+    if (!std::filesystem::exists(aux_root)) {
+        throw std::runtime_error("The auxiliary folder - " + std::string(aux_root) + " - does not exist.");
+    }
+
+    const auto patc_listing = util::filesystem::GetFileListingAt(aux_root, patc_reg);
+    if (patc_listing.size() > 1) {
+        throw std::runtime_error("Expected single PATC file in the '" + std::string(aux_root) + "'. There were " +
+                                 std::to_string(patc_listing.size()) + " found.");
+    }
+
+    if (patc_listing.size() == 0) {
+        LOGW << "No PATC files found at " << aux_root;
+        return std::nullopt;
+    }
+
+    return ParsePatc(patc_listing.front());
+}
+
 }  // namespace alus::asar::envformat::aux
diff --git a/envisat_format/src/envisat_lvl0_parser.cc b/envisat_format/src/envisat_lvl0_parser.cc
@@ -125,6 +125,8 @@ void ParseLevel0Header(const std::vector<char>& file_data, ASARMetadata& asar_me
     asar_meta.sensing_stop = StrToPtime(mph.Get("SENSING_STOP"));
     asar_meta.acquistion_station = mph.Get("ACQUISITION_STATION");
     asar_meta.processing_station = mph.Get("PROC_CENTER");
+    asar_meta.rel_orbit = std::stoul(mph.Get("REL_ORBIT"));
+    asar_meta.abs_orbit = std::stoul(mph.Get("ABS_ORBIT"));
 }
 
 RawSampleMeasurements ParseLevel0Packets(const std::vector<char>& file_data, size_t mdsr_offset_bytes,

diff --git a/envisat_format/src/ers_im_parse.cc b/envisat_format/src/ers_im_parse.cc
@@ -21,6 +21,7 @@
 #include "envisat_format_kernels.h"
 #include "envisat_utils.h"
 #include "ers_env_format.h"
+#include "ers_sbt.h"
 #include "parse_util.h"
 
 #define DEBUG_PACKETS 0
@@ -90,6 +91,7 @@ constexpr auto DR_NO_MAX{alus::asar::envformat::parseutil::MaxValueForBits<uint3
 constexpr auto PACKET_COUNTER_MAX{alus::asar::envformat::parseutil::MaxValueForBits<uint8_t, 8>()};
 constexpr uint8_t SUBCOMMUTATION_COUNTER_MAX{48};
 constexpr auto IMAGE_FORMAT_COUNTER_MAX{alus::asar::envformat::parseutil::MaxValueForBits<uint32_t, 24>()};
+constexpr auto SBT_MAX{alus::asar::envformat::parseutil::MaxValueForBits<uint32_t, 32>()};
 
 void InitializeDataRecordNo(const uint8_t* packet_start, uint32_t& dr_no) {
     FetchUint32(packet_start + alus::asar::envformat::ers::highrate::PROCESSOR_ANNOTATION_ISP_SIZE_BYTES +
@@ -102,6 +104,19 @@ void InitializeDataRecordNo(const uint8_t* packet_start, uint32_t& dr_no) {
     }
 }
 
+void FetchSatelliteBinaryTime(const uint8_t* packet_start, uint32_t& sbt) {
+    /*
+         * The update of that binary counter shall occur every 4th PRI.
+         * The time relation to the echo data in the format is as following:
+         * the transfer of the ICU on-board time to the auxiliary memory
+         * occurs t2 before the RF transmit pulse as depicted in
+         * Fig. 4.4.2.4.6-3. The least significant bit is equal to 1/256 sec.
+         *
+         * ER-IS-ESA-GS-0002  - pg. 4.4 - 24
+     */
+    FetchUint32(packet_start + alus::asar::envformat::ers::highrate::SBT_OFFSET_BYTES, sbt);
+}
+
 void InitializeCounters(const uint8_t* packets_start, uint32_t& dr_no, uint8_t& packet_counter,
                         uint8_t& subcommutation_counter, uint32_t& image_format_counter) {
     InitializeDataRecordNo(packets_start, dr_no);
@@ -156,7 +171,7 @@ inline void FetchPriCode(const uint8_t* start_array, uint16_t& var) {
     var |= start_array[1];
 }
 
-// TODO - Use it inside the parsing or metadata calculation.
+// This is deprecated, will be removed along with the initial parsing flow.
 inline void CalculatePrf(uint16_t pri_code, double& pri, double& prf) {
     // ISCE2 uses 210.943006e-9, but
     // ERS-1-Satellite-to-Ground-Segment-Interface-Specification_01.09.1993_v2D_ER-IS-ESA-GS-0001_EECF05
@@ -170,15 +185,16 @@ inline void CalculatePrf(uint16_t pri_code, double& pri, double& prf) {
     }
 }
 
-inline const uint8_t* FetchPacketFrom(const uint8_t* start, alus::asar::envformat::ForecastMeta& meta,
-                                      uint32_t& dr_no) {
+inline const uint8_t* FetchPacketFrom(const uint8_t* start, alus::asar::envformat::ForecastMeta& meta, uint32_t& dr_no,
+                                      uint32_t& sbt) {
     // https://earth.esa.int/eogateway/documents/20142/37627/ERS-products-specification-with-Envisat-format.pdf
     // https://asf.alaska.edu/wp-content/uploads/2019/03/ers_ceos.pdf and ER-IS-ESA-GS-0002 mixed between three.
     auto it = CopyBSwapPOD(meta.isp_sensing_time, start);
     it += 12;  // No MJD for ERS in FEP.
     uint16_t isp_length;
     it = CopyBSwapPOD(isp_length, it);
     InitializeDataRecordNo(start, dr_no);
+    FetchSatelliteBinaryTime(start, sbt);
 
     // ISP size - 1.
     if (isp_length != 11465) {
@@ -229,15 +245,25 @@ std::vector<ForecastMeta> FetchErsL0ImForecastMeta(const std::vector<char>& file
     size_t dr_no_total_missing{};
     uint32_t last_dr_no;
     InitializeDataRecordNo(next_packet, last_dr_no);
+    uint32_t last_sbt;
+    uint32_t sbt_repeat{0};
+    FetchSatelliteBinaryTime(next_packet, last_sbt);
 
     size_t i{};
     for (; i < mdsr.num_dsr; i++) {
         ForecastMeta current_meta;
         uint32_t dr_no;
         current_meta.packet_start_offset_bytes = next_packet - packets_start;
-        next_packet = FetchPacketFrom(next_packet, current_meta, dr_no);
-
-        // LOGD << i << " " << MjdToPtime(current_meta.isp_sensing_time);
+        uint32_t sbt;
+        next_packet = FetchPacketFrom(next_packet, current_meta, dr_no, sbt);
+        const auto sbt_delta = parseutil::CounterGap<uint32_t, SBT_MAX>(last_sbt, sbt);
+        if (sbt_delta == 0) {
+            sbt_repeat++;
+        } else {
+            sbt_repeat = 0;
+        }
+        bool overflow = (sbt < last_sbt ? true : false);
+        ers::AdjustIspSensingTime(current_meta.isp_sensing_time, sbt, sbt_repeat, overflow);
 
         const auto dr_no_gap = parseutil::CounterGap<uint32_t, DR_NO_MAX>(last_dr_no, dr_no);
         // When gap is 0, this is not handled - could be massive rollover or duplicate or the first packet.
@@ -314,8 +340,8 @@ std::vector<ForecastMeta> FetchErsL0ImForecastMeta(const std::vector<char>& file
     packets_before_start = packets_before_requested_start;
     packets_after_stop = packets_after_requested_stop;
 
-    LOGD << "Data record number discontinuities/total missing - " << dr_no_oos << "/"
-         << dr_no_total_missing << " these were inserted as to be parsed/duplicated during prefetch";
+    LOGD << "Data record number discontinuities/total missing - " << dr_no_oos << "/" << dr_no_total_missing
+         << " these were inserted as to be parsed/duplicated during prefetch";
 
     return forecast_meta;
 }
@@ -350,6 +376,8 @@ RawSampleMeasurements ParseErsLevel0ImPackets(const std::vector<char>& file_data
     uint8_t last_subcommutation_counter{};
     uint32_t last_image_format_counter{};
     InitializeCounters(it, last_dr_no, last_packet_counter, last_subcommutation_counter, last_image_format_counter);
+    uint32_t last_sbt{};
+    FetchSatelliteBinaryTime(it, last_sbt);
 #if HANDLE_DIAGNOSTICS
     PacketParseDiagnostics diagnostics{};
 #endif
@@ -361,7 +389,16 @@ RawSampleMeasurements ParseErsLevel0ImPackets(const std::vector<char>& file_data
         it = packets_start + entry.packet_start_offset_bytes;
         EchoMeta echo_meta = {};
 
+        uint32_t sbt;
+        FetchSatelliteBinaryTime(it, sbt);
         it = CopyBSwapPOD(echo_meta.isp_sensing_time, it);
+        bool overflow = (sbt < last_sbt ? true : false);
+        // sbt_repeat is not handled, because no mechanism could be developed.
+        // SBT stays constant for some period of packets due to being updated after every 4th PRI.
+        // When dividing SBT update time period with PRI no sensible period did develop between packets,
+        // hence this is not used right now. There would be minimal range lines skipped/included due to this anyway.
+        ers::AdjustIspSensingTime(echo_meta.isp_sensing_time, sbt, 0, overflow);
+        last_sbt = sbt;
 
         if (echo_meta.isp_sensing_time != entry.isp_sensing_time) {
             throw std::runtime_error(fmt::format(
@@ -415,15 +452,11 @@ RawSampleMeasurements ParseErsLevel0ImPackets(const std::vector<char>& file_data
          * The time relation to the echo data in the format is as following:
          * the transfer of the ICU on-board time to the auxiliary memory
          * occurs t2 before the RF transmit pulse as depicted in
-         * Fig. 4.4.2.4.6-3. The last significant bit is equal to 1/256 sec.
+         * Fig. 4.4.2.4.6-3. The least significant bit is equal to 1/256 sec.
          *
          * ER-IS-ESA-GS-0002  - pg. 4.4 - 24
          */
-        echo_meta.onboard_time = 0;
-        echo_meta.onboard_time |= static_cast<uint64_t>(it[0]) << 24;
-        echo_meta.onboard_time |= static_cast<uint64_t>(it[1]) << 16;
-        echo_meta.onboard_time |= static_cast<uint64_t>(it[2]) << 8;
-        echo_meta.onboard_time |= static_cast<uint64_t>(it[3]) << 0;
+        echo_meta.onboard_time = sbt;
         it += 4;
         /*
          * This word shall define the activity task within the mode of
@@ -984,4 +1017,4 @@ void ParseErsLevel0ImPackets(const std::vector<char>& file_data, const DSD_lvl0&
 #endif
 }
 
-}  // namespace alus::asar::envformat
+}  // namespace alus::asar::envformat
diff --git a/envisat_format/src/ers_sbt.cc b/envisat_format/src/ers_sbt.cc
@@ -0,0 +1,71 @@
+#include "ers_sbt.h"
+
+#include <cstdint>
+#include <stdexcept>
+
+#include <boost/date_time/posix_time/posix_time.hpp>
+
+#include "date_time_util.h"
+#include "envisat_types.h"
+#include "envisat_utils.h"
+#include "parse_util.h"
+#include "patc.h"
+
+namespace {
+
+bool patc_cor_set{false};
+alus::asar::envformat::aux::Patc patc_cor{};
+boost::posix_time::ptime patc_time_point;
+
+constexpr auto SBT_MAX{alus::asar::envformat::parseutil::MaxValueForBits<uint32_t, 32>()};
+
+}  // namespace
+
+namespace alus::asar::envformat::ers {
+
+void RegisterPatc(const aux::Patc& patc) {
+    if (patc.epoch_1950 < 0) {
+        throw std::runtime_error("Supplied PATC has negative epoch for days since 1950.");
+    }
+
+    if (patc.milliseconds < 0) {
+        throw std::runtime_error("Supplied PATC has negative epoch for milliseconds.");
+    }
+
+    const auto epoch = alus::util::date_time::YYYYMMDD("19500101");
+
+    // Days between 19500101 and 19910101
+    if (patc.epoch_1950 < 14975) {
+        const auto err_date = epoch + boost::posix_time::time_duration(patc.epoch_1950 * 24, 0, 0, 0);
+        throw std::runtime_error(
+            "Supplied PATC epoch days has value that is earlier than year 1990. Resulting date - " +
+            to_simple_string(err_date));
+    }
+
+    patc_cor = patc;
+    patc_time_point = epoch + boost::posix_time::time_duration(patc.epoch_1950 * 24, 0, 0, 0);
+    patc_time_point += boost::posix_time::milliseconds(patc.milliseconds);
+    patc_cor_set = true;
+}
+
+void AdjustIspSensingTime(mjd& isp_sensing_time, uint32_t sbt, uint32_t, bool overflow) {
+    if (!patc_cor_set) {
+        return;
+    }
+
+    boost::posix_time::ptime adjusted_ptime;
+    int64_t delta_count{};
+    if (sbt < patc_cor.sbt_counter && !overflow) {
+        delta_count = static_cast<int64_t>(sbt) - patc_cor.sbt_counter;
+    }
+    else {
+        delta_count = parseutil::CounterGap<uint32_t, SBT_MAX>(patc_cor.sbt_counter, sbt);
+    }
+    const int64_t nanoseconds_diff = static_cast<int64_t>(patc_cor.sbt_period) * delta_count;
+    const auto microseconds_val = static_cast<int64_t>(nanoseconds_diff / 1000);
+    adjusted_ptime = patc_time_point + boost::posix_time::microseconds(microseconds_val);
+
+    isp_sensing_time = PtimeToMjd(adjusted_ptime);
+}
+
+}  // namespace alus::asar::envformat::ers