Load level 2 radar data for rendering

scwx-qt/scwx-qt.cmake

@@ -110,6 +110,10 @@ qt_add_executable(scwx-qt
 
 qt6_create_translation_scwx(QM_FILES ${scwx-qt_SOURCE_DIR} ${TS_FILES})
 
+if (WIN32)
+    target_compile_definitions(scwx-qt PRIVATE WIN32_LEAN_AND_MEAN)
+endif()
+
 target_include_directories(scwx-qt PRIVATE ${scwx-qt_SOURCE_DIR}/source
                                            ${MBGL_INCLUDE_DIR})
 

scwx-qt/source/scwx/qt/manager/radar_manager.cpp

@@ -1,6 +1,7 @@
 #include <scwx/qt/manager/radar_manager.hpp>
 #include <scwx/common/constants.hpp>
 
+#include <deque>
 #include <execution>
 
 #include <boost/log/trivial.hpp>
@@ -27,6 +28,9 @@ static constexpr uint32_t NUM_COORIDNATES_0_5_DEGREE =
 static constexpr uint32_t NUM_COORIDNATES_1_DEGREE =
    NUM_RADIAL_GATES_1_DEGREE * 2;
 
+// TODO: Configure this in settings for radar loop
+static constexpr size_t MAX_LEVEL2_FILES = 6;
+
 class RadarManagerImpl
 {
 public:
@@ -35,6 +39,8 @@ public:
 
    std::vector<float> coordinates0_5Degree_;
    std::vector<float> coordinates1Degree_;
+
+   std::deque<std::shared_ptr<wsr88d::Ar2vFile>> level2Data_;
 };
 
 RadarManager::RadarManager() : p(std::make_unique<RadarManagerImpl>()) {}
@@ -55,28 +61,43 @@ RadarManager::coordinates(common::RadialSize radialSize) const
    throw std::exception("Invalid radial size");
 }
 
+std::shared_ptr<const wsr88d::Ar2vFile> RadarManager::level2_data() const
+{
+   std::shared_ptr<const wsr88d::Ar2vFile> level2Data = nullptr;
+
+   if (p->level2Data_.size() > 0)
+   {
+      level2Data = p->level2Data_.back();
+   }
+
+   return level2Data;
+}
+
 void RadarManager::Initialize()
 {
    BOOST_LOG_TRIVIAL(debug) << logPrefix_ << "Initialize()";
 
    boost::timer::cpu_timer timer;
 
-   // Calculate coordinates
+   GeographicLib::Geodesic geodesic(GeographicLib::Constants::WGS84_a(),
+                                    GeographicLib::Constants::WGS84_f());
+
+   // TODO: This should be retrieved from configuration
+   const QMapbox::Coordinate radar(38.6986, -90.6828);
+
+   // Calculate half degree azimuth coordinates
    timer.start();
    std::vector<float>& coordinates0_5Degree = p->coordinates0_5Degree_;
 
    coordinates0_5Degree.resize(NUM_COORIDNATES_0_5_DEGREE);
 
-   GeographicLib::Geodesic geodesic(GeographicLib::Constants::WGS84_a(),
+   auto radialGates0_5Degree =
-                                    GeographicLib::Constants::WGS84_f());
+      boost::irange<uint32_t>(0, NUM_RADIAL_GATES_0_5_DEGREE);
-
-   const QMapbox::Coordinate radar(38.6986, -90.6828);
-   auto radialGates = boost::irange<uint32_t>(0, NUM_RADIAL_GATES_0_5_DEGREE);
 
    std::for_each(
       std::execution::par_unseq,
-      radialGates.begin(),
+      radialGates0_5Degree.begin(),
-      radialGates.end(),
+      radialGates0_5Degree.end(),
       [&](uint32_t radialGate) {
          const uint16_t gate =
             static_cast<uint16_t>(radialGate % common::MAX_DATA_MOMENT_GATES);
@@ -100,6 +121,62 @@ void RadarManager::Initialize()
    BOOST_LOG_TRIVIAL(debug)
       << logPrefix_ << "Coordinates (0.5 degree) calculated in "
       << timer.format(6, "%ws");
+
+   // Calculate 1 degree azimuth coordinates
+   timer.start();
+   std::vector<float>& coordinates1Degree = p->coordinates1Degree_;
+
+   coordinates1Degree.resize(NUM_COORIDNATES_1_DEGREE);
+
+   auto radialGates1Degree =
+      boost::irange<uint32_t>(0, NUM_RADIAL_GATES_1_DEGREE);
+
+   std::for_each(
+      std::execution::par_unseq,
+      radialGates1Degree.begin(),
+      radialGates1Degree.end(),
+      [&](uint32_t radialGate) {
+         const uint16_t gate =
+            static_cast<uint16_t>(radialGate % common::MAX_DATA_MOMENT_GATES);
+         const uint16_t radial =
+            static_cast<uint16_t>(radialGate / common::MAX_DATA_MOMENT_GATES);
+
+         const float  angle  = radial * 1.0f - 0.5f; // 1 degree radial
+         const float  range  = (gate + 1) * 250.0f;  // 0.25km gate size
+         const size_t offset = radialGate * 2;
+
+         double latitude;
+         double longitude;
+
+         geodesic.Direct(
+            radar.first, radar.second, angle, range, latitude, longitude);
+
+         coordinates1Degree[offset]     = latitude;
+         coordinates1Degree[offset + 1] = longitude;
+      });
+   timer.stop();
+   BOOST_LOG_TRIVIAL(debug)
+      << logPrefix_ << "Coordinates (1 degree) calculated in "
+      << timer.format(6, "%ws");
+}
+
+void RadarManager::LoadLevel2Data(const std::string& filename)
+{
+   std::shared_ptr<wsr88d::Ar2vFile> ar2vFile =
+      std::make_shared<wsr88d::Ar2vFile>();
+
+   bool success = ar2vFile->LoadFile(filename);
+   if (!success)
+   {
+      return;
+   }
+
+   // TODO: Sort and index these
+   if (p->level2Data_.size() >= MAX_LEVEL2_FILES - 1)
+   {
+      p->level2Data_.pop_front();
+   }
+   p->level2Data_.push_back(ar2vFile);
 }
 
 } // namespace manager

scwx-qt/source/scwx/qt/manager/radar_manager.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <scwx/common/types.hpp>
+#include <scwx/wsr88d/ar2v_file.hpp>
 
 #include <memory>
 #include <vector>
@@ -28,7 +29,11 @@ public:
 
    const std::vector<float>& coordinates(common::RadialSize radialSize) const;
 
+   // TODO: Improve this interface
+   std::shared_ptr<const wsr88d::Ar2vFile> level2_data() const;
+
    void Initialize();
+   void LoadLevel2Data(const std::string& filename);
 
 private:
    std::unique_ptr<RadarManagerImpl> p;

scwx-qt/source/scwx/qt/map/map_widget.cpp

@@ -36,7 +36,12 @@ class MapWidgetImpl
 {
 public:
    explicit MapWidgetImpl(const QMapboxGLSettings& settings) :
-       gl_(), settings_(settings), map_(), lastPos_(), frameDraws_(0)
+       gl_(),
+       settings_(settings),
+       map_(),
+       radarManager_ {std::make_shared<manager::RadarManager>()},
+       lastPos_(),
+       frameDraws_(0)
    {
    }
    ~MapWidgetImpl() = default;
@@ -46,6 +51,8 @@ public:
    QMapboxGLSettings          settings_;
    std::shared_ptr<QMapboxGL> map_;
 
+   std::shared_ptr<manager::RadarManager> radarManager_;
+
    QPointF lastPos_;
 
    uint64_t frameDraws_;
@@ -55,6 +62,13 @@ MapWidget::MapWidget(const QMapboxGLSettings& settings) :
     p(std::make_unique<MapWidgetImpl>(settings))
 {
    setFocusPolicy(Qt::StrongFocus);
+
+   p->radarManager_->Initialize();
+   QString ar2vFile = qgetenv("AR2V_FILE");
+   if (!ar2vFile.isEmpty())
+   {
+      p->radarManager_->LoadLevel2Data(ar2vFile.toUtf8().constData());
+   }
 }
 
 MapWidget::~MapWidget()
@@ -87,12 +101,9 @@ void MapWidget::changeStyle()
 
 void MapWidget::AddLayers()
 {
-   std::shared_ptr<manager::RadarManager> radarManager =
-      std::make_shared<manager::RadarManager>();
    std::shared_ptr<view::RadarView> radarView =
-      std::make_shared<view::RadarView>(radarManager, p->map_);
+      std::make_shared<view::RadarView>(p->radarManager_, p->map_);
 
-   radarManager->Initialize();
    radarView->Initialize();
 
    // QMapboxGL::addCustomLayer will take ownership of the QScopedPointer

scwx-qt/source/scwx/qt/view/radar_view.cpp

@@ -63,40 +63,112 @@ void RadarView::Initialize()
 
    boost::timer::cpu_timer timer;
 
+   // TODO: Pick this based on radar data
    const std::vector<float>& coordinates =
       p->radarManager_->coordinates(common::RadialSize::_0_5Degree);
 
+   std::shared_ptr<const wsr88d::Ar2vFile> level2Data =
+      p->radarManager_->level2_data();
+   if (level2Data == nullptr)
+   {
+      return;
+   }
+
+   // TODO: Pick these based on view settings
+   auto                       radarData = level2Data->radar_data()[0];
+   wsr88d::rda::DataBlockType blockType = wsr88d::rda::DataBlockType::MomentRef;
+
    // Calculate vertices
    timer.start();
+
+   auto momentData0 = radarData[0]->moment_data_block(blockType);
+
    std::vector<float>& vertices = p->vertices_;
-   const uint32_t      radials  = common::MAX_RADIALS;
+   const size_t        radials  = radarData.size();
-   const uint32_t      gates    = common::MAX_DATA_MOMENT_GATES;
+   const uint32_t      gates    = momentData0->number_of_data_moment_gates();
    vertices.clear();
    vertices.resize(radials * gates * VERTICES_PER_BIN * VALUES_PER_VERTEX);
    size_t index = 0;
 
-   for (uint16_t radial = 0; radial < 720; ++radial)
+   // Compute threshold at which to display an individual bin
-   {
+   const float    scale  = momentData0->scale();
-      const float dataMomentRange     = 2.125f * 1000.0f;
+   const float    offset = momentData0->offset();
-      const float dataMomentInterval  = 0.25f * 1000.0f;
+   const uint16_t snrThreshold =
-      const float dataMomentIntervalH = dataMomentInterval * 0.5f;
+      std::lroundf(momentData0->snr_threshold_raw() * scale / 10 + offset);
-      const float snrThreshold        = 2.0f;
 
-      const uint16_t startGate               = 7;
+   // Azimuth resolution spacing:
-      const uint16_t numberOfDataMomentGates = 1832;
+   //   1 = 0.5 degrees
+   //   2 = 1.0 degrees
+   const float radialMultiplier =
+      2.0f /
+      std::clamp<int8_t>(radarData[0]->azimuth_resolution_spacing(), 1, 2);
+
+   const float    startAngle  = radarData[0]->azimuth_angle();
+   const uint16_t startRadial = std::lroundf(startAngle * radialMultiplier);
+
+   for (uint16_t radial = 0; radial < radials; ++radial)
+   {
+      auto radialData = radarData[radial];
+      auto momentData = radarData[radial]->moment_data_block(blockType);
+
+      // Compute gate interval
+      const uint16_t dataMomentRange = momentData->data_moment_range_raw();
+      const uint16_t dataMomentInterval =
+         momentData->data_moment_range_sample_interval_raw();
+      const uint16_t dataMomentIntervalH = dataMomentInterval / 2;
+
+      // Compute gate size (number of base 250m gates per bin)
+      const uint16_t gateSize = std::max<uint16_t>(1, dataMomentInterval / 250);
+
+      // Compute gate range [startGate, endGate)
+      const uint16_t startGate = (dataMomentRange - dataMomentIntervalH) / 250;
+      const uint16_t numberOfDataMomentGates =
+         std::min<uint16_t>(momentData->number_of_data_moment_gates(),
+                            static_cast<uint16_t>(gates));
       const uint16_t endGate =
-         std::min<uint16_t>(numberOfDataMomentGates + startGate,
+         std::min<uint16_t>(startGate + numberOfDataMomentGates * gateSize,
-                            common::MAX_DATA_MOMENT_GATES - 1);
+                            common::MAX_DATA_MOMENT_GATES);
 
-      for (uint16_t gate = startGate; gate < endGate; ++gate)
+      const uint8_t*  dataMoments8  = nullptr;
+      const uint16_t* dataMoments16 = nullptr;
+
+      if (momentData->data_word_size() == 8)
       {
-         size_t offset1 = (radial * common::MAX_DATA_MOMENT_GATES + gate) * 2;
+         dataMoments8 =
-         size_t offset2 = offset1 + 2;
+            reinterpret_cast<const uint8_t*>(momentData->data_moments());
-         size_t offset3 = (((radial + 1) % common::MAX_RADIALS) *
+      }
+      else
+      {
+         dataMoments16 =
+            reinterpret_cast<const uint16_t*>(momentData->data_moments());
+      }
+
+      for (uint16_t gate = startGate, i = 0; gate + gateSize <= endGate;
+           gate += gateSize, ++i)
+      {
+         uint16_t dataValue =
+            (dataMoments8 != nullptr) ? dataMoments8[i] : dataMoments16[i];
+
+         if (dataValue < snrThreshold)
+         {
+            continue;
+         }
+
+         if (gate > 0)
+         {
+            const uint16_t baseCoord = gate - 1;
+
+            size_t offset1 = ((startRadial + radial) % common::MAX_RADIALS *
                                  common::MAX_DATA_MOMENT_GATES +
-                           gate) *
+                              baseCoord) *
                              2;
-         size_t offset4 = offset3 + 2;
+            size_t offset2 = offset1 + gateSize * 2;
+            size_t offset3 =
+               (((startRadial + radial + 1) % common::MAX_RADIALS) *
+                   common::MAX_DATA_MOMENT_GATES +
+                baseCoord) *
+               2;
+            size_t offset4 = offset3 + gateSize * 2;
 
             vertices[index++] = coordinates[offset1];
             vertices[index++] = coordinates[offset1 + 1];
@@ -116,7 +188,33 @@ void RadarView::Initialize()
             vertices[index++] = coordinates[offset2];
             vertices[index++] = coordinates[offset2 + 1];
          }
+         else
+         {
+            const uint16_t baseCoord = gate;
+
+            size_t offset1 = ((startRadial + radial) % common::MAX_RADIALS *
+                                 common::MAX_DATA_MOMENT_GATES +
+                              baseCoord) *
+                             2;
+            size_t offset2 =
+               (((startRadial + radial + 1) % common::MAX_RADIALS) *
+                   common::MAX_DATA_MOMENT_GATES +
+                baseCoord) *
+               2;
+
+            // TODO: Radar location
+            vertices[index++] = 38.6986f;
+            vertices[index++] = -90.6828f;
+
+            vertices[index++] = coordinates[offset1];
+            vertices[index++] = coordinates[offset1 + 1];
+
+            vertices[index++] = coordinates[offset2];
+            vertices[index++] = coordinates[offset2 + 1];
          }
+      }
+   }
+   vertices.resize(index);
    timer.stop();
    BOOST_LOG_TRIVIAL(debug)
       << logPrefix_ << "Vertices calculated in " << timer.format(6, "%ws");