From e96808d14d32152928e96eaa07ac95de735d8d90 Mon Sep 17 00:00:00 2001 From: Dan Paulat Date: Fri, 5 Jan 2024 23:38:17 -0600 Subject: [PATCH] Level 2 data level function implementations --- .../scwx/qt/view/level2_product_view.cpp | 202 +++++++++++++++++- wxdata/include/scwx/wsr88d/wsr88d_types.hpp | 5 + wxdata/source/scwx/wsr88d/wsr88d_types.cpp | 10 + 3 files changed, 208 insertions(+), 9 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 43b487ca..6c89ea6d 100644 --- a/scwx-qt/source/scwx/qt/view/level2_product_view.cpp +++ b/scwx-qt/source/scwx/qt/view/level2_product_view.cpp @@ -468,7 +468,6 @@ void Level2ProductView::ComputeSweep() const uint32_t gates = momentData0->number_of_data_moment_gates(); - auto radialData0 = radarData0->radial_data_block(); auto volumeData0 = radarData0->volume_data_block(); p->latitude_ = volumeData0->latitude(); p->longitude_ = volumeData0->longitude(); @@ -782,24 +781,209 @@ void Level2ProductViewImpl::ComputeCoordinates( std::optional Level2ProductView::GetBinLevel(const common::Coordinate& coordinate) const { - // TODO - Q_UNUSED(coordinate); - return std::nullopt; + auto radarData = p->elevationScan_; + auto dataBlockType = p->dataBlockType_; + + if (radarData == nullptr) + { + return std::nullopt; + } + + auto radarProductManager = radar_product_manager(); + auto radarSite = radarProductManager->radar_site(); + const double radarLatitude = radarSite->latitude(); + const double radarLongitude = radarSite->longitude(); + + // Determine distance and azimuth of coordinate relative to radar location + double s12; // Distance (meters) + double azi1; // Azimuth (degrees) + double azi2; // Unused + util::GeographicLib::DefaultGeodesic().Inverse(radarLatitude, + radarLongitude, + coordinate.latitude_, + coordinate.longitude_, + s12, + azi1, + azi2); + + if (std::isnan(azi1)) + { + // If a problem occurred with the geodesic inverse calculation + return std::nullopt; + } + + // Azimuth is returned as [-180, 180) from the geodesic inverse, we need a + // range of [0, 360) + while (azi1 < 0.0) + { + azi1 += 360.0; + } + + // Find Radial + const std::uint16_t numRadials = + static_cast(radarData->size()); + auto radials = boost::irange(0u, numRadials); + + auto radial = std::find_if( // + std::execution::par_unseq, + radials.begin(), + radials.end(), + [&](std::uint32_t i) + { + bool found = false; + const float startAngle = (*radarData)[i]->azimuth_angle(); + const float nextAngle = + (*radarData)[(i + 1) % numRadials]->azimuth_angle(); + + if (startAngle < nextAngle) + { + if (startAngle <= azi1 && azi1 < nextAngle) + { + found = true; + } + } + else + { + // If the bin crosses 0/360 degrees, special handling is needed + if (startAngle <= azi1 || azi1 < nextAngle) + { + found = true; + } + } + + return found; + }); + + if (radial == radials.end()) + { + // No radial was found (not likely to happen without a gap in data) + return std::nullopt; + } + + // Compute gate interval + auto momentData = (*radarData)[*radial]->moment_data_block(dataBlockType); + const std::uint16_t dataMomentRange = momentData->data_moment_range_raw(); + const std::uint16_t dataMomentInterval = + momentData->data_moment_range_sample_interval_raw(); + const std::uint16_t dataMomentIntervalH = dataMomentInterval / 2; + + // Compute gate size (number of base 250m gates per bin) + const std::uint16_t gateSizeMeters = + static_cast(radarProductManager->gate_size()); + + // Compute gate range [startGate, endGate) + const std::uint16_t startGate = + (dataMomentRange - dataMomentIntervalH) / gateSizeMeters; + const std::uint16_t numberOfDataMomentGates = + momentData->number_of_data_moment_gates(); + + const std::uint16_t gate = s12 / dataMomentInterval - startGate; + + if (gate > numberOfDataMomentGates || gate > common::MAX_DATA_MOMENT_GATES) + { + // Coordinate is beyond radar range + return std::nullopt; + } + + // Compute threshold at which to display an individual bin (minimum of 2) + const std::uint16_t snrThreshold = + std::max(2, momentData->snr_threshold_raw()); + std::uint16_t level; + + if (momentData->data_word_size() == 8) + { + level = + reinterpret_cast(momentData->data_moments())[gate]; + } + else + { + level = + reinterpret_cast(momentData->data_moments())[gate]; + } + + if (level < snrThreshold && level != RANGE_FOLDED) + { + return std::nullopt; + } + + return level; } std::optional Level2ProductView::GetDataLevelCode(std::uint16_t level) const { - // TODO - Q_UNUSED(level); + switch (p->product_) + { + case common::Level2Product::Reflectivity: + case common::Level2Product::Velocity: + case common::Level2Product::SpectrumWidth: + case common::Level2Product::DifferentialReflectivity: + case common::Level2Product::DifferentialPhase: + case common::Level2Product::CorrelationCoefficient: + if (level == RANGE_FOLDED) + { + return wsr88d::DataLevelCode::RangeFolded; + } + break; + + case common::Level2Product::ClutterFilterPowerRemoved: + switch (level) + { + case 0: + return wsr88d::DataLevelCode::ClutterFilterNotApplied; + case 1: + return wsr88d::DataLevelCode::ClutterFilterApplied; + case 2: + return wsr88d::DataLevelCode::DualPolVariablesFiltered; + case 3: + case 4: + case 5: + case 6: + case 7: + return wsr88d::DataLevelCode::Reserved; + default: + break; + } + break; + + default: + break; + } + return std::nullopt; } std::optional Level2ProductView::GetDataValue(std::uint16_t level) const { - // TODO - Q_UNUSED(level); - return std::nullopt; + const float offset = p->momentDataBlock0_->offset(); + const float scale = p->momentDataBlock0_->scale(); + std::uint16_t threshold = std::numeric_limits::max(); + + switch (p->product_) + { + case common::Level2Product::Reflectivity: + case common::Level2Product::Velocity: + case common::Level2Product::SpectrumWidth: + case common::Level2Product::DifferentialReflectivity: + case common::Level2Product::DifferentialPhase: + case common::Level2Product::CorrelationCoefficient: + threshold = 2; + break; + + case common::Level2Product::ClutterFilterPowerRemoved: + threshold = 8; + break; + + default: + break; + } + + if (level < threshold) + { + return std::nullopt; + } + + return (level - offset) / scale; } std::shared_ptr Level2ProductView::Create( diff --git a/wxdata/include/scwx/wsr88d/wsr88d_types.hpp b/wxdata/include/scwx/wsr88d/wsr88d_types.hpp index 428f3765..ae63a155 100644 --- a/wxdata/include/scwx/wsr88d/wsr88d_types.hpp +++ b/wxdata/include/scwx/wsr88d/wsr88d_types.hpp @@ -50,6 +50,11 @@ enum class DataLevelCode Z8, SI, + // Clutter Filter Power Removed + ClutterFilterNotApplied, + ClutterFilterApplied, + DualPolVariablesFiltered, + Unknown }; diff --git a/wxdata/source/scwx/wsr88d/wsr88d_types.cpp b/wxdata/source/scwx/wsr88d/wsr88d_types.cpp index 4ea72757..1cc8747e 100644 --- a/wxdata/source/scwx/wsr88d/wsr88d_types.cpp +++ b/wxdata/source/scwx/wsr88d/wsr88d_types.cpp @@ -50,6 +50,11 @@ static const std::unordered_map dataLevelCodeName_ { {DataLevelCode::Z8, "R(Z) * 0.8"}, {DataLevelCode::SI, "R(Z) * multiplier"}, + // Clutter Filter Power Removed + {DataLevelCode::ClutterFilterNotApplied, "Clutter Filter Not Applied"}, + {DataLevelCode::ClutterFilterApplied, "Clutter Filter Applied"}, + {DataLevelCode::DualPolVariablesFiltered, "Dual Pol Variables Filtered"}, + {DataLevelCode::Unknown, "?"}}; static const std::unordered_map @@ -95,6 +100,11 @@ static const std::unordered_map {DataLevelCode::Z8, "Z8"}, {DataLevelCode::SI, "SI"}, + // Clutter Filter Power Removed + {DataLevelCode::ClutterFilterNotApplied, ""}, + {DataLevelCode::ClutterFilterApplied, ""}, + {DataLevelCode::DualPolVariablesFiltered, ""}, + {DataLevelCode::Unknown, "?"}}; const std::string& GetDataLevelCodeName(DataLevelCode dataLevelCode)