From a5484977671da9c2a4eeb9756d93c65c2d6a0649 Mon Sep 17 00:00:00 2001 From: Dan Paulat Date: Tue, 1 Oct 2024 06:31:02 -0500 Subject: [PATCH 1/4] Add partial handling for missing level 2 radials - Need to handle additional cases (1999-05-03 @ KTLX is a good sample) - Still crashes on getting bin level for hover text --- .../scwx/qt/view/level2_product_view.cpp | 106 +++++++++++++----- 1 file changed, 75 insertions(+), 31 deletions(-) diff --git a/scwx-qt/source/scwx/qt/view/level2_product_view.cpp b/scwx-qt/source/scwx/qt/view/level2_product_view.cpp index 47558d81..a952b863 100644 --- a/scwx-qt/source/scwx/qt/view/level2_product_view.cpp +++ b/scwx-qt/source/scwx/qt/view/level2_product_view.cpp @@ -536,7 +536,7 @@ void Level2ProductView::ComputeSweep() return; } - const size_t radials = radarData->size(); + const std::size_t radials = radarData->crbegin()->first; p->ComputeCoordinates(radarData); @@ -829,7 +829,7 @@ void Level2ProductViewImpl::ComputeCoordinates( auto momentData0 = radarData0->moment_data_block(dataBlockType_); const std::uint16_t numRadials = - static_cast(radarData->size()); + static_cast(radarData->crbegin()->first + 1); const std::uint16_t numRangeBins = std::max(momentData0->number_of_data_moment_gates() + 1u, common::MAX_DATA_MOMENT_GATES); @@ -837,39 +837,83 @@ void Level2ProductViewImpl::ComputeCoordinates( auto radials = boost::irange(0u, numRadials); auto gates = boost::irange(0u, numRangeBins); - std::for_each(std::execution::par_unseq, - radials.begin(), - radials.end(), - [&](std::uint32_t radial) - { - const units::degrees angle = - (*radarData)[radial]->azimuth_angle(); + std::for_each( + std::execution::par_unseq, + radials.begin(), + radials.end(), + [&](std::uint32_t radial) + { + units::degrees angle {}; - 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; + auto radialData = radarData->find(radial); + if (radialData != radarData->cend()) + { + angle = radialData->second->azimuth_angle(); + } + else + { + auto prevRadial1 = radarData->find( + (radial >= 1) ? radial - 1 : numRadials - (1 - radial)); + auto prevRadial2 = radarData->find( + (radial >= 2) ? radial - 2 : numRadials - (2 - radial)); - double latitude; - double longitude; + if (prevRadial1 != radarData->cend() && + prevRadial2 != radarData->cend()) + { + const units::degrees prevAngle1 = + prevRadial1->second->azimuth_angle(); + const units::degrees prevAngle2 = + prevRadial2->second->azimuth_angle(); - geodesic.Direct(radarLatitude, - radarLongitude, - angle.value(), - range, - latitude, - longitude); + // No wrapping required since angle is only used for geodesic + // calculation + const units::degrees deltaAngle = prevAngle1 - prevAngle2; - coordinates_[offset] = latitude; - coordinates_[offset + 1] = longitude; - }); - }); + angle = prevAngle1 + deltaAngle; + } + else if (prevRadial1 != radarData->cend()) + { + const units::degrees prevAngle1 = + prevRadial1->second->azimuth_angle(); + + // Assume a half degree delta if there aren't enough angles + // to determine a delta angle + constexpr units::degrees deltaAngle = + units::degrees {0.5}; + + 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")); } From fe4a324a04d005aaab77aff485434438b3a1311b Mon Sep 17 00:00:00 2001 From: Dan Paulat Date: Wed, 2 Oct 2024 05:56:27 -0500 Subject: [PATCH 2/4] Use an extra vertex radial with missing data to prevent stretching --- .../scwx/qt/view/level2_product_view.cpp | 56 ++++++++++++++++--- wxdata/include/scwx/common/geographic.hpp | 13 +++++ wxdata/source/scwx/common/geographic.cpp | 28 ++++++++++ 3 files changed, 89 insertions(+), 8 deletions(-) diff --git a/scwx-qt/source/scwx/qt/view/level2_product_view.cpp b/scwx-qt/source/scwx/qt/view/level2_product_view.cpp index a952b863..15112410 100644 --- a/scwx-qt/source/scwx/qt/view/level2_product_view.cpp +++ b/scwx-qt/source/scwx/qt/view/level2_product_view.cpp @@ -106,14 +106,17 @@ public: threadPool_.join(); }; - void - ComputeCoordinates(std::shared_ptr radarData); + void ComputeCoordinates( + const std::shared_ptr& 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& radarData); + Level2ProductView* self_; boost::asio::thread_pool threadPool_ {1u}; @@ -536,7 +539,17 @@ void Level2ProductView::ComputeSweep() return; } - const std::size_t radials = radarData->crbegin()->first; + const 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; + } p->ComputeCoordinates(radarData); @@ -574,7 +587,8 @@ void Level2ProductView::ComputeSweep() std::vector& 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& dataMoments8 = p->dataMoments8_; @@ -807,7 +821,7 @@ void Level2ProductView::ComputeSweep() } void Level2ProductViewImpl::ComputeCoordinates( - std::shared_ptr radarData) + const std::shared_ptr& radarData) { logger_->debug("ComputeCoordinates()"); @@ -828,12 +842,22 @@ 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(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(numRadials, common::MAX_0_5_DEGREE_RADIALS); + auto radials = boost::irange(0u, numRadials); auto gates = boost::irange(0u, numRangeBins); @@ -878,8 +902,7 @@ void Level2ProductViewImpl::ComputeCoordinates( // Assume a half degree delta if there aren't enough angles // to determine a delta angle - constexpr units::degrees deltaAngle = - units::degrees {0.5}; + constexpr units::degrees deltaAngle {0.5f}; angle = prevAngle1 + deltaAngle; } @@ -918,6 +941,23 @@ void Level2ProductViewImpl::ComputeCoordinates( logger_->debug("Coordinates calculated in {}", timer.format(6, "%ws")); } +bool Level2ProductViewImpl::IsRadarDataIncomplete( + const std::shared_ptr& radarData) +{ + // Assume the data is incomplete when the delta between the first and last + // angles is greater than 2.5 degrees. + constexpr units::degrees kIncompleteDataAngleThreshold_ {2.5}; + + const units::degrees firstAngle = + radarData->cbegin()->second->azimuth_angle(); + const units::degrees lastAngle = + radarData->crbegin()->second->azimuth_angle(); + const units::degrees angleDelta = + common::GetAngleDelta(firstAngle, lastAngle); + + return angleDelta > kIncompleteDataAngleThreshold_; +} + std::optional Level2ProductView::GetBinLevel(const common::Coordinate& coordinate) const { diff --git a/wxdata/include/scwx/common/geographic.hpp b/wxdata/include/scwx/common/geographic.hpp index 8945db17..8b234fd2 100644 --- a/wxdata/include/scwx/common/geographic.hpp +++ b/wxdata/include/scwx/common/geographic.hpp @@ -3,6 +3,8 @@ #include #include +#include + 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 GetAngleDelta(units::degrees angle1, + units::degrees angle2); + /** * Calculate the geographic midpoint of a set of coordinates. Uses Method A * described at http://www.geomidpoint.com/calculation.html. diff --git a/wxdata/source/scwx/common/geographic.cpp b/wxdata/source/scwx/common/geographic.cpp index bf7d1a2c..e9a494d3 100644 --- a/wxdata/source/scwx/common/geographic.cpp +++ b/wxdata/source/scwx/common/geographic.cpp @@ -14,6 +14,34 @@ static std::string GetDegreeString(double degrees, DegreeStringType type, const std::string& suffix); +units::degrees GetAngleDelta(units::degrees angle1, + units::degrees angle2) +{ + // Normalize angles to [0, 360) + while (angle1.value() < 0.0f) + { + angle1 += units::degrees {360.0f}; + } + while (angle2.value() < 0.0f) + { + angle2 += units::degrees {360.0f}; + } + angle1 = units::degrees {std::fmod(angle1.value(), 360.f)}; + angle2 = units::degrees {std::fmod(angle2.value(), 360.f)}; + + // Calculate the absolute difference + auto delta = angle1 - angle2; + if (delta < units::degrees {0.0f}) + { + delta *= -1.0f; + } + + // Account for wrapping + delta = std::min(delta, units::degrees {360.0f} - delta); + + return delta; +} + Coordinate GetCentroid(const std::vector& coordinates) { double x = 0.0; From 41c112538842567945436ad77e69d32db581387a Mon Sep 17 00:00:00 2001 From: Dan Paulat Date: Thu, 3 Oct 2024 05:44:40 -0500 Subject: [PATCH 3/4] Fix usage of vertexRadials in ComputeSweep --- .../scwx/qt/view/level2_product_view.cpp | 31 ++++++++++++------- 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/scwx-qt/source/scwx/qt/view/level2_product_view.cpp b/scwx-qt/source/scwx/qt/view/level2_product_view.cpp index 15112410..a2c1b42e 100644 --- a/scwx-qt/source/scwx/qt/view/level2_product_view.cpp +++ b/scwx-qt/source/scwx/qt/view/level2_product_view.cpp @@ -539,8 +539,8 @@ void Level2ProductView::ComputeSweep() return; } - const std::size_t radials = radarData->crbegin()->first + 1; - std::size_t vertexRadials = radials; + 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 @@ -551,6 +551,11 @@ void Level2ProductView::ComputeSweep() ++vertexRadials; } + // Limit radials + radials = std::min(radials, common::MAX_0_5_DEGREE_RADIALS); + vertexRadials = + std::min(vertexRadials, common::MAX_0_5_DEGREE_RADIALS); + p->ComputeCoordinates(radarData); const std::vector& coordinates = p->coordinates_; @@ -735,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) * - common::MAX_DATA_MOMENT_GATES + - baseCoord) * - 2; + std::size_t offset3 = + (((startRadial + radial + 1) % vertexRadials) * + common::MAX_DATA_MOMENT_GATES + + baseCoord) * + 2; std::size_t offset4 = offset3 + gateSize * 2; vertices[vIndex++] = coordinates[offset1]; @@ -770,14 +776,15 @@ void Level2ProductView::ComputeSweep() { const std::uint16_t baseCoord = gate; - std::size_t offset1 = ((startRadial + radial) % radials * - common::MAX_DATA_MOMENT_GATES + - baseCoord) * - 2; - std::size_t offset2 = (((startRadial + radial + 1) % radials) * + std::size_t offset1 = ((startRadial + radial) % vertexRadials * 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_; From 1436b7bba64ca1f82552a02fe8fc79d361234ee6 Mon Sep 17 00:00:00 2001 From: Dan Paulat Date: Thu, 3 Oct 2024 05:54:44 -0500 Subject: [PATCH 4/4] Handle missing level 2 radials when getting bin data --- .../scwx/qt/view/level2_product_view.cpp | 73 +++++++++++++++---- 1 file changed, 59 insertions(+), 14 deletions(-) diff --git a/scwx-qt/source/scwx/qt/view/level2_product_view.cpp b/scwx-qt/source/scwx/qt/view/level2_product_view.cpp index a2c1b42e..0938f614 100644 --- a/scwx-qt/source/scwx/qt/view/level2_product_view.cpp +++ b/scwx-qt/source/scwx/qt/view/level2_product_view.cpp @@ -1007,8 +1007,19 @@ Level2ProductView::GetBinLevel(const common::Coordinate& coordinate) const } // Find Radial - const std::uint16_t numRadials = - static_cast(radarData->size()); + std::uint16_t numRadials = + static_cast(radarData->crbegin()->first + 1); + + // Add an extra radial when incomplete data exists + if (Level2ProductViewImpl::IsRadarDataIncomplete(radarData)) + { + ++numRadials; + } + + // Limit radials + numRadials = + std::min(numRadials, common::MAX_0_5_DEGREE_RADIALS); + auto radials = boost::irange(0u, numRadials); auto radial = std::find_if( // @@ -1017,25 +1028,59 @@ Level2ProductView::GetBinLevel(const common::Coordinate& coordinate) const radials.end(), [&](std::uint32_t i) { - bool found = false; - const units::degrees startAngle = - (*radarData)[i]->azimuth_angle(); - const units::degrees nextAngle = - (*radarData)[(i + 1) % numRadials]->azimuth_angle(); + bool hasNextAngle = false; + bool found = false; - if (startAngle < nextAngle) + units::degrees startAngle {}; + units::degrees 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 prevAngle = + prevRadial->second->azimuth_angle(); + + const units::degrees 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.value() <= azi1 || azi1 < nextAngle.value()) + if (startAngle < nextAngle) { - 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; + } } }