Merge pull request #279 from dpaulat/hotfix/missing-level2-radials

Handle Missing Level 2 Radials

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

@@ -106,14 +106,17 @@ public:
       threadPool_.join();
    };
 
-   void
+   void ComputeCoordinates(
-   ComputeCoordinates(std::shared_ptr<wsr88d::rda::ElevationScan> radarData);
+      const std::shared_ptr<wsr88d::rda::ElevationScan>& radarData);
 
    void SetProduct(const std::string& productName);
    void SetProduct(common::Level2Product product);
    void UpdateOtherUnits(const std::string& name);
    void UpdateSpeedUnits(const std::string& name);
 
+   static bool IsRadarDataIncomplete(
+      const std::shared_ptr<const wsr88d::rda::ElevationScan>& radarData);
+
    Level2ProductView* self_;
 
    boost::asio::thread_pool threadPool_ {1u};
@@ -536,7 +539,22 @@ void Level2ProductView::ComputeSweep()
       return;
    }
 
-   const size_t radials = radarData->size();
+   std::size_t radials       = radarData->crbegin()->first + 1;
+   std::size_t vertexRadials = radials;
+
+   // When there is missing data, insert another empty vertex radial at the end
+   // to avoid stretching
+   const bool isRadarDataIncomplete =
+      Level2ProductViewImpl::IsRadarDataIncomplete(radarData);
+   if (isRadarDataIncomplete)
+   {
+      ++vertexRadials;
+   }
+
+   // Limit radials
+   radials = std::min<std::size_t>(radials, common::MAX_0_5_DEGREE_RADIALS);
+   vertexRadials =
+      std::min<std::size_t>(vertexRadials, common::MAX_0_5_DEGREE_RADIALS);
 
    p->ComputeCoordinates(radarData);
 
@@ -574,7 +592,8 @@ void Level2ProductView::ComputeSweep()
    std::vector<float>& vertices = p->vertices_;
    size_t              vIndex   = 0;
    vertices.clear();
-   vertices.resize(radials * gates * VERTICES_PER_BIN * VALUES_PER_VERTEX);
+   vertices.resize(vertexRadials * gates * VERTICES_PER_BIN *
+                   VALUES_PER_VERTEX);
 
    // Setup data moment vector
    std::vector<uint8_t>&  dataMoments8  = p->dataMoments8_;
@@ -721,15 +740,16 @@ void Level2ProductView::ComputeSweep()
          {
             const std::uint16_t baseCoord = gate - 1;
 
-            std::size_t offset1 = ((startRadial + radial) % radials *
+            std::size_t offset1 = ((startRadial + radial) % vertexRadials *
                                       common::MAX_DATA_MOMENT_GATES +
                                    baseCoord) *
                                   2;
             std::size_t offset2 = offset1 + gateSize * 2;
-            std::size_t offset3 = (((startRadial + radial + 1) % radials) *
+            std::size_t offset3 =
-                                      common::MAX_DATA_MOMENT_GATES +
+               (((startRadial + radial + 1) % vertexRadials) *
-                                   baseCoord) *
+                   common::MAX_DATA_MOMENT_GATES +
-                                  2;
+                baseCoord) *
+               2;
             std::size_t offset4 = offset3 + gateSize * 2;
 
             vertices[vIndex++] = coordinates[offset1];
@@ -756,14 +776,15 @@ void Level2ProductView::ComputeSweep()
          {
             const std::uint16_t baseCoord = gate;
 
-            std::size_t offset1 = ((startRadial + radial) % radials *
+            std::size_t offset1 = ((startRadial + radial) % vertexRadials *
-                                      common::MAX_DATA_MOMENT_GATES +
-                                   baseCoord) *
-                                  2;
-            std::size_t offset2 = (((startRadial + radial + 1) % radials) *
                                       common::MAX_DATA_MOMENT_GATES +
                                    baseCoord) *
                                   2;
+            std::size_t offset2 =
+               (((startRadial + radial + 1) % vertexRadials) *
+                   common::MAX_DATA_MOMENT_GATES +
+                baseCoord) *
+               2;
 
             vertices[vIndex++] = p->latitude_;
             vertices[vIndex++] = p->longitude_;
@@ -807,7 +828,7 @@ void Level2ProductView::ComputeSweep()
 }
 
 void Level2ProductViewImpl::ComputeCoordinates(
-   std::shared_ptr<wsr88d::rda::ElevationScan> radarData)
+   const std::shared_ptr<wsr88d::rda::ElevationScan>& radarData)
 {
    logger_->debug("ComputeCoordinates()");
 
@@ -828,52 +849,122 @@ void Level2ProductViewImpl::ComputeCoordinates(
    auto& radarData0  = (*radarData)[0];
    auto  momentData0 = radarData0->moment_data_block(dataBlockType_);
 
-   const std::uint16_t numRadials =
+   std::uint16_t numRadials =
-      static_cast<std::uint16_t>(radarData->size());
+      static_cast<std::uint16_t>(radarData->crbegin()->first + 1);
    const std::uint16_t numRangeBins =
       std::max(momentData0->number_of_data_moment_gates() + 1u,
                common::MAX_DATA_MOMENT_GATES);
 
+   // Add an extra radial when incomplete data exists
+   if (IsRadarDataIncomplete(radarData))
+   {
+      ++numRadials;
+   }
+
+   // Limit radials
+   numRadials =
+      std::min<std::uint16_t>(numRadials, common::MAX_0_5_DEGREE_RADIALS);
+
    auto radials = boost::irange<std::uint32_t>(0u, numRadials);
    auto gates   = boost::irange<std::uint32_t>(0u, numRangeBins);
 
-   std::for_each(std::execution::par_unseq,
+   std::for_each(
-                 radials.begin(),
+      std::execution::par_unseq,
-                 radials.end(),
+      radials.begin(),
-                 [&](std::uint32_t radial)
+      radials.end(),
-                 {
+      [&](std::uint32_t radial)
-                    const units::degrees<float> angle =
+      {
-                       (*radarData)[radial]->azimuth_angle();
+         units::degrees<float> angle {};
 
-                    std::for_each(std::execution::par_unseq,
+         auto radialData = radarData->find(radial);
-                                  gates.begin(),
+         if (radialData != radarData->cend())
-                                  gates.end(),
+         {
-                                  [&](std::uint32_t gate)
+            angle = radialData->second->azimuth_angle();
-                                  {
+         }
-                                     const std::uint32_t radialGate =
+         else
-                                        radial * common::MAX_DATA_MOMENT_GATES +
+         {
-                                        gate;
+            auto prevRadial1 = radarData->find(
-                                     const float range = (gate + 1) * gateSize;
+               (radial >= 1) ? radial - 1 : numRadials - (1 - radial));
-                                     const std::size_t offset = radialGate * 2;
+            auto prevRadial2 = radarData->find(
+               (radial >= 2) ? radial - 2 : numRadials - (2 - radial));
 
-                                     double latitude;
+            if (prevRadial1 != radarData->cend() &&
-                                     double longitude;
+                prevRadial2 != radarData->cend())
+            {
+               const units::degrees<float> prevAngle1 =
+                  prevRadial1->second->azimuth_angle();
+               const units::degrees<float> prevAngle2 =
+                  prevRadial2->second->azimuth_angle();
 
-                                     geodesic.Direct(radarLatitude,
+               // No wrapping required since angle is only used for geodesic
-                                                     radarLongitude,
+               // calculation
-                                                     angle.value(),
+               const units::degrees<float> deltaAngle = prevAngle1 - prevAngle2;
-                                                     range,
-                                                     latitude,
-                                                     longitude);
 
-                                     coordinates_[offset]     = latitude;
+               angle = prevAngle1 + deltaAngle;
-                                     coordinates_[offset + 1] = longitude;
+            }
-                                  });
+            else if (prevRadial1 != radarData->cend())
-                 });
+            {
+               const units::degrees<float> prevAngle1 =
+                  prevRadial1->second->azimuth_angle();
+
+               // Assume a half degree delta if there aren't enough angles
+               // to determine a delta angle
+               constexpr units::degrees<float> deltaAngle {0.5f};
+
+               angle = prevAngle1 + deltaAngle;
+            }
+            else
+            {
+               // Not enough angles present to determine an angle
+               return;
+            }
+         }
+
+         std::for_each(std::execution::par_unseq,
+                       gates.begin(),
+                       gates.end(),
+                       [&](std::uint32_t gate)
+                       {
+                          const std::uint32_t radialGate =
+                             radial * common::MAX_DATA_MOMENT_GATES + gate;
+                          const float       range  = (gate + 1) * gateSize;
+                          const std::size_t offset = radialGate * 2;
+
+                          double latitude;
+                          double longitude;
+
+                          geodesic.Direct(radarLatitude,
+                                          radarLongitude,
+                                          angle.value(),
+                                          range,
+                                          latitude,
+                                          longitude);
+
+                          coordinates_[offset]     = latitude;
+                          coordinates_[offset + 1] = longitude;
+                       });
+      });
    timer.stop();
    logger_->debug("Coordinates calculated in {}", timer.format(6, "%ws"));
 }
 
+bool Level2ProductViewImpl::IsRadarDataIncomplete(
+   const std::shared_ptr<const wsr88d::rda::ElevationScan>& radarData)
+{
+   // Assume the data is incomplete when the delta between the first and last
+   // angles is greater than 2.5 degrees.
+   constexpr units::degrees<float> kIncompleteDataAngleThreshold_ {2.5};
+
+   const units::degrees<float> firstAngle =
+      radarData->cbegin()->second->azimuth_angle();
+   const units::degrees<float> lastAngle =
+      radarData->crbegin()->second->azimuth_angle();
+   const units::degrees<float> angleDelta =
+      common::GetAngleDelta(firstAngle, lastAngle);
+
+   return angleDelta > kIncompleteDataAngleThreshold_;
+}
+
 std::optional<std::uint16_t>
 Level2ProductView::GetBinLevel(const common::Coordinate& coordinate) const
 {
@@ -916,8 +1007,19 @@ Level2ProductView::GetBinLevel(const common::Coordinate& coordinate) const
    }
 
    // Find Radial
-   const std::uint16_t numRadials =
+   std::uint16_t numRadials =
-      static_cast<std::uint16_t>(radarData->size());
+      static_cast<std::uint16_t>(radarData->crbegin()->first + 1);
+
+   // Add an extra radial when incomplete data exists
+   if (Level2ProductViewImpl::IsRadarDataIncomplete(radarData))
+   {
+      ++numRadials;
+   }
+
+   // Limit radials
+   numRadials =
+      std::min<std::uint16_t>(numRadials, common::MAX_0_5_DEGREE_RADIALS);
+
    auto radials = boost::irange<std::uint32_t>(0u, numRadials);
 
    auto radial = std::find_if( //
@@ -926,25 +1028,59 @@ Level2ProductView::GetBinLevel(const common::Coordinate& coordinate) const
       radials.end(),
       [&](std::uint32_t i)
       {
-         bool                        found = false;
+         bool hasNextAngle = false;
-         const units::degrees<float> startAngle =
+         bool found        = false;
-            (*radarData)[i]->azimuth_angle();
-         const units::degrees<float> nextAngle =
-            (*radarData)[(i + 1) % numRadials]->azimuth_angle();
 
-         if (startAngle < nextAngle)
+         units::degrees<float> startAngle {};
+         units::degrees<float> nextAngle {};
+
+         auto radialData = radarData->find(i);
+         if (radialData != radarData->cend())
          {
-            if (startAngle.value() <= azi1 && azi1 < nextAngle.value())
+            startAngle = radialData->second->azimuth_angle();
+
+            auto nextRadial = radarData->find((i + 1) % numRadials);
+            if (nextRadial != radarData->cend())
             {
-               found = true;
+               nextAngle    = nextRadial->second->azimuth_angle();
+               hasNextAngle = true;
+            }
+            else
+            {
+               // Next angle is not available, interpolate
+               auto prevRadial =
+                  radarData->find((i >= 1) ? i - 1 : numRadials - (1 - i));
+
+               if (prevRadial != radarData->cend())
+               {
+                  const units::degrees<float> prevAngle =
+                     prevRadial->second->azimuth_angle();
+
+                  const units::degrees<float> deltaAngle =
+                     common::GetAngleDelta(startAngle, prevAngle);
+
+                  nextAngle    = startAngle + deltaAngle;
+                  hasNextAngle = true;
+               }
             }
          }
-         else
+
+         if (hasNextAngle)
          {
-            // If the bin crosses 0/360 degrees, special handling is needed
+            if (startAngle < nextAngle)
-            if (startAngle.value() <= azi1 || azi1 < nextAngle.value())
             {
-               found = true;
+               if (startAngle.value() <= azi1 && azi1 < nextAngle.value())
+               {
+                  found = true;
+               }
+            }
+            else
+            {
+               // If the bin crosses 0/360 degrees, special handling is needed
+               if (startAngle.value() <= azi1 || azi1 < nextAngle.value())
+               {
+                  found = true;
+               }
             }
          }
 

wxdata/include/scwx/common/geographic.hpp

@@ -3,6 +3,8 @@
 #include <string>
 #include <vector>
 
+#include <units/angle.h>
+
 namespace scwx
 {
 namespace common
@@ -46,6 +48,17 @@ enum class DistanceType
    Miles
 };
 
+/**
+ * Calculate the absolute angle delta between two angles.
+ *
+ * @param [in] angle1 First angle
+ * @param [in] angle2 Second angle
+ *
+ * @return Absolute angle delta normalized to [0, 360)
+ */
+units::degrees<float> GetAngleDelta(units::degrees<float> angle1,
+                                    units::degrees<float> angle2);
+
 /**
  * Calculate the geographic midpoint of a set of coordinates. Uses Method A
  * described at http://www.geomidpoint.com/calculation.html.

wxdata/source/scwx/common/geographic.cpp

@@ -14,6 +14,34 @@ static std::string GetDegreeString(double             degrees,
                                    DegreeStringType   type,
                                    const std::string& suffix);
 
+units::degrees<float> GetAngleDelta(units::degrees<float> angle1,
+                                    units::degrees<float> angle2)
+{
+   // Normalize angles to [0, 360)
+   while (angle1.value() < 0.0f)
+   {
+      angle1 += units::degrees<float> {360.0f};
+   }
+   while (angle2.value() < 0.0f)
+   {
+      angle2 += units::degrees<float> {360.0f};
+   }
+   angle1 = units::degrees<float> {std::fmod(angle1.value(), 360.f)};
+   angle2 = units::degrees<float> {std::fmod(angle2.value(), 360.f)};
+
+   // Calculate the absolute difference
+   auto delta = angle1 - angle2;
+   if (delta < units::degrees<float> {0.0f})
+   {
+      delta *= -1.0f;
+   }
+
+   // Account for wrapping
+   delta = std::min(delta, units::degrees<float> {360.0f} - delta);
+
+   return delta;
+}
+
 Coordinate GetCentroid(const std::vector<Coordinate>& coordinates)
 {
    double x = 0.0;