nexrad_render/

lib.rs

//! Rendering functions for NEXRAD weather radar data.
//!
//! This crate provides functions to render radar data into visual images. It converts
//! radar moment data (reflectivity, velocity, etc.) into color-mapped images that can
//! be saved to common formats like PNG.
//!
//! # Example
//!
//! ```ignore
//! use nexrad_model::data::Product;
//! use nexrad_render::{render_radials, RenderOptions, nws_reflectivity_scale};
//!
//! let options = RenderOptions::new(800, 800);
//! let image = render_radials(
//!     sweep.radials(),
//!     Product::Reflectivity,
//!     &nws_reflectivity_scale(),
//!     &options,
//! ).unwrap();
//!
//! // Save directly to PNG
//! image.save("radar.png").unwrap();
//! ```
//!
//! # Crate Boundaries
//!
//! This crate provides **visualization and rendering** with the following responsibilities
//! and constraints:
//!
//! ## Responsibilities
//!
//! - Render radar data to images ([`image::RgbaImage`])
//! - Apply color scales to moment data
//! - Handle geometric transformations (polar to Cartesian coordinates)
//! - Consume `nexrad-model` types (Radial, MomentData)
//!
//! ## Constraints
//!
//! - **No data access or network operations**
//! - **No binary parsing or decoding**
//!
//! This crate can be used standalone or through the `nexrad` facade crate (re-exported
//! via the `render` feature, which is enabled by default).

#![forbid(unsafe_code)]
#![deny(clippy::unwrap_used)]
#![deny(clippy::expect_used)]
#![warn(clippy::correctness)]
#![deny(missing_docs)]

pub use image::RgbaImage;
pub use nexrad_model::data::Product;
use nexrad_model::data::{
    CFPMomentValue, CartesianField, DataMoment, GateStatus, MomentValue, Radial, SweepField,
    VerticalField,
};
use nexrad_model::geo::{GeoExtent, RadarCoordinateSystem};
use result::{Error, Result};

mod color;
pub use crate::color::*;

mod render_result;
pub use render_result::{PointQuery, RenderMetadata, RenderResult};

pub mod result;

/// Interpolation method for rendering radar data.
///
/// Controls how pixel values are sampled from the underlying data grid.
/// The default is [`Nearest`](Interpolation::Nearest) for backward compatibility.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum Interpolation {
    /// Nearest-neighbor sampling (fastest, produces blocky/aliased output).
    ///
    /// Each output pixel takes the value of the single closest data gate.
    /// This is the legacy behavior and the default.
    #[default]
    Nearest,
    /// Bilinear interpolation (smoother, anti-aliased output).
    ///
    /// Each output pixel blends the four surrounding data gates using
    /// bilinear weights. Interpolation happens in **data value space**
    /// (the f32 moment value), not in color space, ensuring scientifically
    /// correct results.
    ///
    /// Falls back to nearest-neighbor for any pixel where one or more of
    /// the four neighbors has a non-[`Valid`](GateStatus::Valid) status.
    Bilinear,
}

/// Options for rendering radar radials.
///
/// Use the builder methods to configure rendering options, then pass to
/// [`render_radials`].
///
/// # Example
///
/// ```
/// use nexrad_render::RenderOptions;
///
/// // Render 800x800 with black background (default)
/// let options = RenderOptions::new(800, 800);
///
/// // Render with transparent background for compositing
/// let options = RenderOptions::new(800, 800).transparent();
///
/// // Render with custom background color (RGBA)
/// let options = RenderOptions::new(800, 800).with_background([255, 255, 255, 255]);
/// ```
#[derive(Debug, Clone)]
pub struct RenderOptions {
    pub(crate) size: (usize, usize),
    pub(crate) background: Option<[u8; 4]>,
    pub(crate) extent: Option<GeoExtent>,
    pub(crate) coord_system: Option<RadarCoordinateSystem>,
    pub(crate) interpolation: Interpolation,
}

impl RenderOptions {
    /// Creates new render options with the specified dimensions and black background.
    pub fn new(width: usize, height: usize) -> Self {
        Self {
            size: (width, height),
            background: Some([0, 0, 0, 255]),
            extent: None,
            coord_system: None,
            interpolation: Interpolation::Nearest,
        }
    }

    /// Sets the background to transparent for compositing.
    ///
    /// When rendering with a transparent background, areas without radar data
    /// will be fully transparent, allowing multiple renders to be layered.
    pub fn transparent(mut self) -> Self {
        self.background = None;
        self
    }

    /// Sets a custom background color as RGBA bytes.
    pub fn with_background(mut self, rgba: [u8; 4]) -> Self {
        self.background = Some(rgba);
        self
    }

    /// Sets a geographic extent for the rendered area.
    ///
    /// When set, the image covers exactly this extent. This enables
    /// consistent spatial mapping for side-by-side comparison of multiple renders.
    pub fn with_extent(mut self, extent: GeoExtent) -> Self {
        self.extent = Some(extent);
        self
    }

    /// Sets a radar coordinate system for geographic projection.
    ///
    /// This enables geographic metadata in the [`RenderResult`], including
    /// pixel-to-geo and geo-to-pixel coordinate conversions.
    pub fn with_coord_system(mut self, coord_system: RadarCoordinateSystem) -> Self {
        self.coord_system = Some(coord_system);
        self
    }

    /// Creates render options sized for native resolution of a sweep field.
    ///
    /// Sets both width and height to `gate_count * 2`, which ensures approximately
    /// one pixel per gate at the outer edge of the sweep. This produces the highest
    /// fidelity rendering without wasting pixels.
    pub fn native_for(field: &SweepField) -> Self {
        let size = field.gate_count() * 2;
        Self::new(size, size)
    }

    /// Sets the interpolation method for rendering.
    ///
    /// Default is [`Interpolation::Nearest`]. Use [`Interpolation::Bilinear`]
    /// for smoother output that blends neighboring data gates.
    pub fn with_interpolation(mut self, interpolation: Interpolation) -> Self {
        self.interpolation = interpolation;
        self
    }

    /// Shorthand to enable bilinear interpolation.
    ///
    /// Equivalent to `.with_interpolation(Interpolation::Bilinear)`.
    pub fn bilinear(mut self) -> Self {
        self.interpolation = Interpolation::Bilinear;
        self
    }

    /// Output image dimensions (width, height) in pixels.
    pub fn size(&self) -> (usize, usize) {
        self.size
    }

    /// Background color as RGBA bytes. `None` means transparent (all zeros).
    pub fn background(&self) -> Option<[u8; 4]> {
        self.background
    }

    /// Geographic extent to render. If `None`, auto-computed from data range.
    ///
    /// When set, the image covers exactly this extent, enabling consistent
    /// spatial mapping across multiple renders for side-by-side comparison.
    pub fn extent(&self) -> Option<&GeoExtent> {
        self.extent.as_ref()
    }

    /// Radar coordinate system for geographic projection.
    ///
    /// When provided, the [`RenderResult`] will include geographic metadata
    /// enabling pixel-to-geo and geo-to-pixel coordinate conversions.
    pub fn coord_system(&self) -> Option<&RadarCoordinateSystem> {
        self.coord_system.as_ref()
    }

    /// Interpolation method for pixel sampling.
    ///
    /// Default is [`Interpolation::Nearest`]. Use [`Interpolation::Bilinear`]
    /// for smoother output that blends neighboring data gates.
    pub fn interpolation(&self) -> Interpolation {
        self.interpolation
    }
}

/// Renders radar radials to an RGBA image.
///
/// Converts polar radar data into a Cartesian image representation. Each radial's
/// moment values are mapped to colors using the provided color scale, producing
/// a centered radar image with North at the top.
///
/// # Arguments
///
/// * `radials` - Slice of radials to render (typically from a single sweep)
/// * `product` - The radar product (moment type) to visualize
/// * `scale` - Color scale mapping moment values to colors
/// * `options` - Rendering options (size, background, etc.)
///
/// # Errors
///
/// Returns an error if:
/// - No radials are provided
/// - The requested product is not present in the radials
///
/// # Example
///
/// ```ignore
/// use nexrad_render::{render_radials, Product, RenderOptions, nws_reflectivity_scale};
///
/// let scale = nws_reflectivity_scale();
/// let options = RenderOptions::new(800, 800);
///
/// let image = render_radials(
///     sweep.radials(),
///     Product::Reflectivity,
///     &scale,
///     &options,
/// ).unwrap();
///
/// image.save("radar.png").unwrap();
/// ```
pub fn render_radials(
    radials: &[Radial],
    product: Product,
    scale: &DiscreteColorScale,
    options: &RenderOptions,
) -> Result<RgbaImage> {
    let (width, height) = options.size;
    let mut buffer = vec![0u8; width * height * 4];

    // Fill background
    if let Some(bg) = options.background {
        for pixel in buffer.chunks_exact_mut(4) {
            pixel.copy_from_slice(&bg);
        }
    }

    if radials.is_empty() {
        return Err(Error::NoRadials);
    }

    // Build lookup table for fast color mapping
    let (min_val, max_val) = product.value_range();
    let lut = ColorLookupTable::from_scale(scale, min_val, max_val, 256);

    // Get radar parameters from the first radial
    let first_radial = &radials[0];
    let gate_params = get_gate_params(product, first_radial).ok_or(Error::ProductNotFound)?;
    let first_gate_km = gate_params.first_gate_km;
    let gate_interval_km = gate_params.gate_interval_km;
    let gate_count = gate_params.gate_count;
    let radar_range_km = first_gate_km + gate_count as f64 * gate_interval_km;

    // Pre-extract all moment float values indexed by azimuth for efficient lookup
    let mut radial_data: Vec<(f32, Vec<Option<f32>>)> = Vec::with_capacity(radials.len());
    for radial in radials {
        let azimuth = radial.azimuth_angle_degrees();
        if let Some(values) = get_radial_float_values(product, radial) {
            radial_data.push((azimuth, values));
        }
    }

    // Sort by azimuth for binary search
    radial_data.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));

    // Extract azimuths for binary search
    let sorted_azimuths: Vec<f32> = radial_data.iter().map(|(az, _)| *az).collect();

    // Calculate max azimuth gap threshold based on radial spacing
    // Use 1.5x the expected spacing to allow for minor gaps while rejecting large ones
    let azimuth_spacing = first_radial.azimuth_spacing_degrees();
    let max_azimuth_gap = azimuth_spacing * 1.5;

    let center_x = width as f64 / 2.0;
    let center_y = height as f64 / 2.0;
    let scale_factor = width.max(height) as f64 / 2.0 / radar_range_km;

    // Render each pixel by mapping to radar coordinates
    for y in 0..height {
        let dy = y as f64 - center_y;

        for x in 0..width {
            let dx = x as f64 - center_x;

            // Convert pixel position to distance in km
            let distance_pixels = (dx * dx + dy * dy).sqrt();
            let distance_km = distance_pixels / scale_factor;

            // Skip pixels outside radar coverage
            if distance_km < first_gate_km || distance_km >= radar_range_km {
                continue;
            }

            // Calculate azimuth angle (0 = North, clockwise)
            let azimuth_rad = dx.atan2(-dy);
            let azimuth_deg = (azimuth_rad.to_degrees() + 360.0) % 360.0;

            // Find the closest radial and check if it's within acceptable range
            let (radial_idx, angular_distance) =
                find_closest_radial(&sorted_azimuths, azimuth_deg as f32);

            // Skip pixels where no radial is close enough (partial sweep gaps)
            if angular_distance > max_azimuth_gap {
                continue;
            }

            // Calculate gate index
            let gate_index = ((distance_km - first_gate_km) / gate_interval_km) as usize;
            if gate_index >= gate_count {
                continue;
            }

            // Look up the value and apply color
            let (_, ref values) = radial_data[radial_idx];
            if let Some(Some(value)) = values.get(gate_index) {
                let color = lut.color(*value);
                let pixel_index = (y * width + x) * 4;
                buffer[pixel_index..pixel_index + 4].copy_from_slice(&color);
            }
        }
    }

    // Convert buffer to RgbaImage
    RgbaImage::from_raw(width as u32, height as u32, buffer).ok_or(Error::InvalidDimensions)
}

/// Renders radar radials using the default color scale for the product.
///
/// This is a convenience function that automatically selects an appropriate
/// color scale based on the product type, using standard meteorological conventions.
///
/// # Arguments
///
/// * `radials` - Slice of radials to render (typically from a single sweep)
/// * `product` - The radar product (moment type) to visualize
/// * `options` - Rendering options (size, background, etc.)
///
/// # Errors
///
/// Returns an error if:
/// - No radials are provided
/// - The requested product is not present in the radials
///
/// # Example
///
/// ```ignore
/// use nexrad_render::{render_radials_default, Product, RenderOptions};
///
/// let options = RenderOptions::new(800, 800);
/// let image = render_radials_default(
///     sweep.radials(),
///     Product::Velocity,
///     &options,
/// ).unwrap();
///
/// image.save("velocity.png").unwrap();
/// ```
pub fn render_radials_default(
    radials: &[Radial],
    product: Product,
    options: &RenderOptions,
) -> Result<RgbaImage> {
    let scale = default_scale(product);
    render_radials(radials, product, &scale, options)
}

/// Renders a [`SweepField`] (polar grid) to an image with metadata.
///
/// This is the primary rendering entry point for processed data. It accepts
/// a `SweepField` — the output of data extraction or processing — and produces
/// a [`RenderResult`] that includes both the rendered image and metadata for
/// geographic placement and data inspection.
///
/// # Arguments
///
/// * `field` - The sweep field to render
/// * `color_scale` - Color scale to apply (discrete or continuous)
/// * `options` - Rendering options (size, background, extent, coordinate system)
///
/// # Example
///
/// ```ignore
/// use nexrad_render::{render_sweep, RenderOptions, ColorScale};
/// use nexrad_model::data::{SweepField, Product};
///
/// let field = SweepField::from_radials(sweep.radials(), Product::Reflectivity).unwrap();
/// let scale = ColorScale::from(nexrad_render::nws_reflectivity_scale());
/// let result = render_sweep(&field, &scale, &RenderOptions::new(800, 800))?;
///
/// result.image().save("radar.png").unwrap();
/// let meta = result.metadata();
/// let ring_px = meta.km_to_pixel_distance(100.0);
/// ```
pub fn render_sweep(
    field: &SweepField,
    color_scale: &ColorScale,
    options: &RenderOptions,
) -> Result<RenderResult> {
    let (width, height) = options.size;
    let mut buffer = vec![0u8; width * height * 4];

    // Fill background
    if let Some(bg) = options.background {
        for pixel in buffer.chunks_exact_mut(4) {
            pixel.copy_from_slice(&bg);
        }
    }

    if field.azimuth_count() == 0 {
        return Err(Error::NoRadials);
    }

    let first_gate_km = field.first_gate_range_km();
    let gate_interval_km = field.gate_interval_km();
    let gate_count = field.gate_count();
    let radar_range_km = field.max_range_km();

    // Build LUT from the field's value range or product-typical range
    let (min_val, max_val) = field.value_range().unwrap_or((-32.0, 95.0));
    let lut = ColorLookupTable::from_color_scale(color_scale, min_val, max_val, 256);

    // Calculate max azimuth gap for partial sweep support
    let azimuth_spacing = field.azimuth_spacing_degrees();
    let max_azimuth_gap = azimuth_spacing * 1.5;

    let center_x = width as f64 / 2.0;
    let center_y = height as f64 / 2.0;
    let scale_factor = width.max(height) as f64 / 2.0 / radar_range_km;

    // Render each pixel by mapping to radar coordinates
    let use_bilinear = options.interpolation == Interpolation::Bilinear;

    for y in 0..height {
        let dy = y as f64 - center_y;

        for x in 0..width {
            let dx = x as f64 - center_x;

            let distance_pixels = (dx * dx + dy * dy).sqrt();
            let distance_km = distance_pixels / scale_factor;

            if distance_km < first_gate_km || distance_km >= radar_range_km {
                continue;
            }

            let azimuth_rad = dx.atan2(-dy);
            let azimuth_deg = ((azimuth_rad.to_degrees() + 360.0) % 360.0) as f32;

            // Try bilinear interpolation first, then fall back to nearest-neighbor
            if use_bilinear {
                if let Some(val) =
                    sample_sweep_bilinear(field, azimuth_deg, distance_km, max_azimuth_gap)
                {
                    let color = lut.color(val);
                    let pixel_index = (y * width + x) * 4;
                    buffer[pixel_index..pixel_index + 4].copy_from_slice(&color);
                    continue;
                }
            }

            // Nearest-neighbor path (default, or bilinear fallback)
            let (radial_idx, angular_distance) = find_closest_radial(field.azimuths(), azimuth_deg);

            if angular_distance > max_azimuth_gap {
                continue;
            }

            let gate_index = ((distance_km - first_gate_km) / gate_interval_km) as usize;
            if gate_index >= gate_count {
                continue;
            }

            let (val, status) = field.get(radial_idx, gate_index);
            if status == GateStatus::Valid {
                let color = lut.color(val);
                let pixel_index = (y * width + x) * 4;
                buffer[pixel_index..pixel_index + 4].copy_from_slice(&color);
            }
        }
    }

    let image =
        RgbaImage::from_raw(width as u32, height as u32, buffer).ok_or(Error::InvalidDimensions)?;

    let geo_extent = options.extent.or_else(|| {
        options
            .coord_system
            .as_ref()
            .map(|cs| cs.sweep_extent(radar_range_km))
    });

    let metadata = RenderMetadata {
        width: width as u32,
        height: height as u32,
        center_pixel: (center_x, center_y),
        pixels_per_km: scale_factor,
        max_range_km: radar_range_km,
        elevation_degrees: Some(field.elevation_degrees()),
        geo_extent,
        coord_system: options.coord_system.clone(),
    };

    Ok(RenderResult::new(image, metadata))
}

/// Renders a [`CartesianField`] (geographic grid) to an image with metadata.
///
/// This renders volume-derived products like composite reflectivity, echo tops,
/// and VIL — data that is already projected onto a geographic grid.
///
/// # Arguments
///
/// * `field` - The Cartesian field to render
/// * `color_scale` - Color scale to apply
/// * `options` - Rendering options (size, background)
pub fn render_cartesian(
    field: &CartesianField,
    color_scale: &ColorScale,
    options: &RenderOptions,
) -> Result<RenderResult> {
    let (width, height) = options.size;
    let mut buffer = vec![0u8; width * height * 4];

    if let Some(bg) = options.background {
        for pixel in buffer.chunks_exact_mut(4) {
            pixel.copy_from_slice(&bg);
        }
    }

    let field_width = field.width();
    let field_height = field.height();

    if field_width == 0 || field_height == 0 {
        return Err(Error::NoRadials);
    }

    // Build LUT
    let (min_val, max_val) = field.value_range().unwrap_or((-32.0, 95.0));
    let lut = ColorLookupTable::from_color_scale(color_scale, min_val, max_val, 256);

    // Map output pixels to field cells
    let use_bilinear = options.interpolation == Interpolation::Bilinear;
    let max_row = field_height - 1;
    let max_col = field_width - 1;

    for y in 0..height {
        let row_f = y as f64 / height as f64 * field_height as f64 - 0.5;
        let row_f = row_f.max(0.0);

        for x in 0..width {
            let col_f = x as f64 / width as f64 * field_width as f64 - 0.5;
            let col_f = col_f.max(0.0);

            // Try bilinear interpolation first, then fall back to nearest-neighbor
            if use_bilinear {
                if let Some(val) =
                    sample_grid_bilinear(row_f, col_f, max_row, max_col, |r, c| field.get(r, c))
                {
                    let color = lut.color(val);
                    let pixel_index = (y * width + x) * 4;
                    buffer[pixel_index..pixel_index + 4].copy_from_slice(&color);
                    continue;
                }
            }

            // Nearest-neighbor path (default, or bilinear fallback)
            let field_row = (row_f + 0.5) as usize;
            let field_row = field_row.min(max_row);
            let field_col = (col_f + 0.5) as usize;
            let field_col = field_col.min(max_col);

            let (val, status) = field.get(field_row, field_col);
            if status == GateStatus::Valid {
                let color = lut.color(val);
                let pixel_index = (y * width + x) * 4;
                buffer[pixel_index..pixel_index + 4].copy_from_slice(&color);
            }
        }
    }

    let image =
        RgbaImage::from_raw(width as u32, height as u32, buffer).ok_or(Error::InvalidDimensions)?;

    let metadata = RenderMetadata {
        width: width as u32,
        height: height as u32,
        center_pixel: (width as f64 / 2.0, height as f64 / 2.0),
        pixels_per_km: 0.0, // Not applicable for Cartesian fields
        max_range_km: 0.0,
        elevation_degrees: None,
        geo_extent: Some(*field.extent()),
        coord_system: options.coord_system.clone(),
    };

    Ok(RenderResult::new(image, metadata))
}

/// Renders a [`VerticalField`] (RHI / cross-section display) to an image with metadata.
///
/// This renders vertical cross-sections where the horizontal axis represents
/// distance from the radar and the vertical axis represents altitude.
///
/// # Arguments
///
/// * `field` - The vertical field to render
/// * `color_scale` - Color scale to apply
/// * `options` - Rendering options (size, background)
pub fn render_vertical(
    field: &VerticalField,
    color_scale: &ColorScale,
    options: &RenderOptions,
) -> Result<RenderResult> {
    let (width, height) = options.size;
    let mut buffer = vec![0u8; width * height * 4];

    if let Some(bg) = options.background {
        for pixel in buffer.chunks_exact_mut(4) {
            pixel.copy_from_slice(&bg);
        }
    }

    let field_width = field.width();
    let field_height = field.height();

    if field_width == 0 || field_height == 0 {
        return Err(Error::NoRadials);
    }

    // Build LUT
    let (min_val, max_val) = field.value_range().unwrap_or((-32.0, 95.0));
    let lut = ColorLookupTable::from_color_scale(color_scale, min_val, max_val, 256);

    // Map output pixels to field cells
    let use_bilinear = options.interpolation == Interpolation::Bilinear;
    let max_row = field_height - 1;
    let max_col = field_width - 1;

    for y in 0..height {
        let row_f = y as f64 / height as f64 * field_height as f64 - 0.5;
        let row_f = row_f.max(0.0);

        for x in 0..width {
            let col_f = x as f64 / width as f64 * field_width as f64 - 0.5;
            let col_f = col_f.max(0.0);

            // Try bilinear interpolation first, then fall back to nearest-neighbor
            if use_bilinear {
                if let Some(val) =
                    sample_grid_bilinear(row_f, col_f, max_row, max_col, |r, c| field.get(r, c))
                {
                    let color = lut.color(val);
                    let pixel_index = (y * width + x) * 4;
                    buffer[pixel_index..pixel_index + 4].copy_from_slice(&color);
                    continue;
                }
            }

            // Nearest-neighbor path (default, or bilinear fallback)
            let field_row = (row_f + 0.5) as usize;
            let field_row = field_row.min(max_row);
            let field_col = (col_f + 0.5) as usize;
            let field_col = field_col.min(max_col);

            let (val, status) = field.get(field_row, field_col);
            if status == GateStatus::Valid {
                let color = lut.color(val);
                let pixel_index = (y * width + x) * 4;
                buffer[pixel_index..pixel_index + 4].copy_from_slice(&color);
            }
        }
    }

    let image =
        RgbaImage::from_raw(width as u32, height as u32, buffer).ok_or(Error::InvalidDimensions)?;

    let (d_min, d_max) = field.distance_range_km();
    let max_range = d_max - d_min;

    let metadata = RenderMetadata {
        width: width as u32,
        height: height as u32,
        center_pixel: (0.0, height as f64 / 2.0), // Left edge is range 0
        pixels_per_km: width as f64 / max_range,
        max_range_km: max_range,
        elevation_degrees: None,
        geo_extent: None,
        coord_system: options.coord_system.clone(),
    };

    Ok(RenderResult::new(image, metadata))
}

/// Returns the default color scale for a given product.
///
/// This function selects an appropriate color scale based on the product type,
/// using standard meteorological conventions.
///
/// | Product | Scale |
/// |---------|-------|
/// | Reflectivity | NWS Reflectivity (dBZ) |
/// | Velocity | Divergent Green-Red (-64 to +64 m/s) |
/// | SpectrumWidth | Sequential (0 to 30 m/s) |
/// | DifferentialReflectivity | Divergent (-2 to +6 dB) |
/// | DifferentialPhase | Sequential (0 to 360 deg) |
/// | CorrelationCoefficient | Sequential (0 to 1) |
/// | ClutterFilterPower | Divergent (-20 to +20 dB) |
pub fn default_scale(product: Product) -> DiscreteColorScale {
    match product {
        Product::Reflectivity => nws_reflectivity_scale(),
        Product::Velocity => velocity_scale(),
        Product::SpectrumWidth => spectrum_width_scale(),
        Product::DifferentialReflectivity => differential_reflectivity_scale(),
        Product::DifferentialPhase => differential_phase_scale(),
        Product::CorrelationCoefficient => correlation_coefficient_scale(),
        Product::ClutterFilterPower => clutter_filter_power_scale(),
    }
}

/// Returns the default color scale for a product, wrapped in a [`ColorScale`] enum.
///
/// This is a convenience wrapper around [`default_scale`] that returns the
/// [`ColorScale`] enum directly, avoiding the need to call `.into()` at each call site.
pub fn default_color_scale(product: Product) -> ColorScale {
    default_scale(product).into()
}

/// Find the index in sorted_azimuths closest to the given azimuth and return
/// the angular distance to that radial.
///
/// Returns `(index, angular_distance)` where `angular_distance` is in degrees.
#[inline]
fn find_closest_radial(sorted_azimuths: &[f32], azimuth: f32) -> (usize, f32) {
    let len = sorted_azimuths.len();
    if len == 0 {
        return (0, f32::MAX);
    }

    // Binary search for insertion point
    let pos = sorted_azimuths
        .binary_search_by(|a| a.partial_cmp(&azimuth).unwrap_or(std::cmp::Ordering::Equal))
        .unwrap_or_else(|i| i);

    if pos == 0 {
        // Check if wrapping around (360° boundary) is closer
        let dist_to_first = (sorted_azimuths[0] - azimuth).abs();
        let dist_to_last = 360.0 - sorted_azimuths[len - 1] + azimuth;
        if dist_to_last < dist_to_first {
            return (len - 1, dist_to_last);
        }
        return (0, dist_to_first);
    }

    if pos >= len {
        // Check if wrapping around is closer
        let dist_to_last = (azimuth - sorted_azimuths[len - 1]).abs();
        let dist_to_first = 360.0 - azimuth + sorted_azimuths[0];
        if dist_to_first < dist_to_last {
            return (0, dist_to_first);
        }
        return (len - 1, dist_to_last);
    }

    // Compare distances to neighbors
    let dist_to_prev = (azimuth - sorted_azimuths[pos - 1]).abs();
    let dist_to_curr = (sorted_azimuths[pos] - azimuth).abs();

    if dist_to_prev <= dist_to_curr {
        (pos - 1, dist_to_prev)
    } else {
        (pos, dist_to_curr)
    }
}

/// Find the two radials bracketing the given azimuth and the fractional position between them.
///
/// Returns `Some((lo_idx, hi_idx, frac))` where:
/// - `lo_idx` is the index of the radial at or just below `azimuth`
/// - `hi_idx` is the index of the radial just above `azimuth`
/// - `frac` is in \[0, 1\] — the fractional distance from lo toward hi
///
/// Returns `None` if there are fewer than 2 radials, or if the angular gap between
/// the two bracketing radials exceeds `max_gap` degrees (partial sweep boundary).
#[inline]
fn find_bracketing_radials(
    sorted_azimuths: &[f32],
    azimuth: f32,
    max_gap: f32,
) -> Option<(usize, usize, f32)> {
    let len = sorted_azimuths.len();
    if len < 2 {
        return None;
    }

    let pos = sorted_azimuths
        .binary_search_by(|a| a.partial_cmp(&azimuth).unwrap_or(std::cmp::Ordering::Equal))
        .unwrap_or_else(|i| i);

    // Determine the bracketing pair (lo, hi) such that lo_az <= azimuth <= hi_az
    // with 360-degree wraparound at the edges.
    let (lo_idx, hi_idx) = if pos == 0 || pos >= len {
        // Wraps around: last azimuth → first azimuth
        (len - 1, 0)
    } else {
        (pos - 1, pos)
    };

    let lo_az = sorted_azimuths[lo_idx];
    let hi_az = sorted_azimuths[hi_idx];

    // Angular span with wraparound
    let span = if hi_az >= lo_az {
        hi_az - lo_az
    } else {
        360.0 - lo_az + hi_az
    };

    if span > max_gap || span == 0.0 {
        return None;
    }

    // Fractional position of azimuth between lo and hi
    let offset = if azimuth >= lo_az {
        azimuth - lo_az
    } else {
        360.0 - lo_az + azimuth
    };

    let frac = (offset / span).clamp(0.0, 1.0);
    Some((lo_idx, hi_idx, frac))
}

/// Sample a sweep field using bilinear interpolation in data-value space.
///
/// Returns the interpolated f32 value, or `None` if any of the four neighbors
/// is non-Valid, the pixel falls in a sweep gap, or the coordinates are at the
/// outer edge where interpolation is not possible.
#[inline]
fn sample_sweep_bilinear(
    field: &SweepField,
    azimuth_deg: f32,
    distance_km: f64,
    max_azimuth_gap: f32,
) -> Option<f32> {
    let gate_count = field.gate_count();
    let first_gate_km = field.first_gate_range_km();
    let gate_interval_km = field.gate_interval_km();

    // Azimuth axis: find bracketing radials
    let (az_lo, az_hi, az_frac) =
        find_bracketing_radials(field.azimuths(), azimuth_deg, max_azimuth_gap)?;

    // Range axis: find bracketing gates
    let gate_f = (distance_km - first_gate_km) / gate_interval_km;
    let g_lo = gate_f as usize;
    let g_hi = g_lo + 1;

    if g_hi >= gate_count {
        return None;
    }

    let g_frac = (gate_f - g_lo as f64) as f32;

    // Fetch the four corners and require all Valid
    let (v00, s00) = field.get(az_lo, g_lo);
    let (v01, s01) = field.get(az_lo, g_hi);
    let (v10, s10) = field.get(az_hi, g_lo);
    let (v11, s11) = field.get(az_hi, g_hi);

    if s00 != GateStatus::Valid
        || s01 != GateStatus::Valid
        || s10 != GateStatus::Valid
        || s11 != GateStatus::Valid
    {
        return None;
    }

    // Bilinear blend: lerp along range for each azimuth, then across azimuths
    let v_lo = v00 + (v01 - v00) * g_frac;
    let v_hi = v10 + (v11 - v10) * g_frac;
    Some(v_lo + (v_hi - v_lo) * az_frac)
}

/// Sample a row-major grid using bilinear interpolation in data-value space.
///
/// `row_f` and `col_f` are fractional coordinates into the grid.
/// Returns the interpolated value, or `None` if any of the four neighbors is non-Valid.
#[inline]
fn sample_grid_bilinear(
    row_f: f64,
    col_f: f64,
    max_row: usize,
    max_col: usize,
    get_fn: impl Fn(usize, usize) -> (f32, GateStatus),
) -> Option<f32> {
    let r0 = row_f as usize;
    let c0 = col_f as usize;
    let r1 = (r0 + 1).min(max_row);
    let c1 = (c0 + 1).min(max_col);

    let rf = (row_f - r0 as f64) as f32;
    let cf = (col_f - c0 as f64) as f32;

    let (v00, s00) = get_fn(r0, c0);
    let (v01, s01) = get_fn(r0, c1);
    let (v10, s10) = get_fn(r1, c0);
    let (v11, s11) = get_fn(r1, c1);

    if s00 != GateStatus::Valid
        || s01 != GateStatus::Valid
        || s10 != GateStatus::Valid
        || s11 != GateStatus::Valid
    {
        return None;
    }

    let v_top = v00 + (v01 - v00) * cf;
    let v_bot = v10 + (v11 - v10) * cf;
    Some(v_top + (v_bot - v_top) * rf)
}

/// Gate metadata extracted from a moment data block.
struct GateParams {
    first_gate_km: f64,
    gate_interval_km: f64,
    gate_count: usize,
}

/// Retrieve gate metadata from a radial for the given product.
fn get_gate_params(product: Product, radial: &Radial) -> Option<GateParams> {
    fn extract(m: &impl DataMoment) -> GateParams {
        GateParams {
            first_gate_km: m.first_gate_range_km(),
            gate_interval_km: m.gate_interval_km(),
            gate_count: m.gate_count() as usize,
        }
    }
    if let Some(moment) = product.moment_data(radial) {
        return Some(extract(moment));
    }
    if let Some(cfp) = product.cfp_moment_data(radial) {
        return Some(extract(cfp));
    }
    None
}

/// Extract decoded float values from a radial for the given product.
///
/// Returns `None` for below-threshold, range-folded, and CFP status gates.
fn get_radial_float_values(product: Product, radial: &Radial) -> Option<Vec<Option<f32>>> {
    if let Some(moment) = product.moment_data(radial) {
        return Some(
            moment
                .iter()
                .map(|v| match v {
                    MomentValue::Value(f) => Some(f),
                    _ => None,
                })
                .collect(),
        );
    }

    if let Some(cfp) = product.cfp_moment_data(radial) {
        return Some(
            cfp.iter()
                .map(|v| match v {
                    CFPMomentValue::Value(f) => Some(f),
                    CFPMomentValue::Status(_) => None,
                })
                .collect(),
        );
    }

    None
}