grib_reader/

lib.rs

//! Pure-Rust GRIB file reader.
//!
//! The reader supports GRIB1 and GRIB2 regular geographic grids, GRIB2 rotated
//! latitude/longitude, regular Gaussian, and common projected grids, simple and
//! complex and IEEE packing, multi-field bitmap reuse, bounded allocation, and
//! optional image-backed and CCSDS packing codecs.
//!
//! Missing grid points from bitmaps and embedded complex-packing codes decode
//! as `NaN`. Complex-packing substitute metadata remains available through
//! [`ComplexPackingParams::primary_missing_substitute`] and
//! [`ComplexPackingParams::secondary_missing_substitute`] for callers that
//! need producer-specific replacement values.
//!
//! # Example
//!
//! ```no_run
//! use grib_reader::GribFile;
//!
//! let file = GribFile::open("gfs.grib2")?;
//! println!("messages: {}", file.message_count());
//!
//! for msg in file.messages() {
//!     println!(
//!         "  {} {:?} {:04}-{:02}-{:02}T{:02}:{:02}:{:02}Z",
//!         msg.parameter_name(),
//!         msg.grid_shape(),
//!         msg.reference_time().year,
//!         msg.reference_time().month,
//!         msg.reference_time().day,
//!         msg.reference_time().hour,
//!         msg.reference_time().minute,
//!         msg.reference_time().second,
//!     );
//! }
//!
//! let data = file.message(0)?.read_data_as_f64()?;
//! println!("shape: {:?}", data.shape());
//! # Ok::<(), grib_reader::Error>(())
//! ```

pub mod data;
pub mod error;
pub mod grib1;
pub mod grid;
pub mod indicator;
pub mod metadata;
pub mod parameter;
pub mod product;
pub mod sections;

pub use data::{
    CcsdsBlockSize, CcsdsFlags, CcsdsPackingParams, ComplexPackingParams, DataRepresentation,
    DecodeSample, Grib1SecondOrderLayout, Grib1SecondOrderPackingParams, IeeePackingParams,
    IeeePrecision, Jpeg2000PackingParams, PngPackingParams, RunLengthPackingParams,
    ScaledPackingParams, SpatialDifferencingParams,
};
pub use error::{Error, Result};
pub use grib1::{BinaryDataSection, GridDescription, ProductDefinition as Grib1ProductDefinition};
pub use grid::{
    AlbersEqualAreaGrid, GridDefinition, LambertConformalGrid, LatLonGrid, MercatorGrid,
    PolarStereographicGrid, ProjectedGridCore, RegularGaussianGrid, RotatedLatLonGrid,
};
pub use metadata::{ForecastTimeUnit, Parameter, ParameterTableSource, ReferenceTime};
pub use parameter::{
    LocalParameterEntry, LocalParameterTable, OwnedLocalParameterEntry, BUILTIN_LOCAL_PARAMETERS,
    LOCAL_PARAMETER_TABLE_CSV_HEADER,
};
pub use product::{
    AnalysisOrForecastTemplate, DerivedForecastTemplate, DerivedStatisticalProcessTemplate,
    EnsembleStatisticalProcessTemplate, FixedSurface, Identification,
    IndividualEnsembleForecastTemplate, PercentileForecastTemplate,
    PercentileStatisticalProcessTemplate, ProbabilityForecastTemplate, ProbabilityLimit,
    ProbabilityStatisticalProcessTemplate, ProbabilityType, ProductDefinition,
    ProductDefinitionTemplate, ScaledValue, SpatialProcessTemplate, StatisticalInterval,
    StatisticalProcessTemplate, StatisticalTimeRange,
};

use std::io::Read;
use std::path::Path;
use std::sync::Arc;

use grib_core::bit::read_bit;
use grib_core::{ensure_limit, filled_vec};
use memchr::memmem::Finder;
use memmap2::Mmap;
use ndarray::{ArrayD, IxDyn};

use crate::data::{
    bitmap_payload as grib2_bitmap_payload, count_bitmap_present_points,
    decode_field_into_with_group_lengths, decode_payload_into_with_group_lengths,
};
use crate::indicator::Indicator;
use crate::parameter::lookup_parameter_with_local_entries;
use crate::sections::{index_fields, FieldSections, SectionRef};

#[cfg(feature = "rayon")]
use rayon::prelude::*;

const GRIB_MAGIC: &[u8; 4] = b"GRIB";
pub const DEFAULT_MAX_DECODED_POINTS: usize = 25_000_000;
pub const DEFAULT_MAX_AXIS_POINTS: usize = 1_000_000;

#[derive(Debug, Clone, Copy)]
struct ReaderOptions {
    strict: bool,
    max_decoded_points: Option<usize>,
    max_axis_points: Option<usize>,
}

impl Default for ReaderOptions {
    fn default() -> Self {
        Self {
            strict: true,
            max_decoded_points: Some(DEFAULT_MAX_DECODED_POINTS),
            max_axis_points: Some(DEFAULT_MAX_AXIS_POINTS),
        }
    }
}

/// Builder for opening GRIB data with explicit scanning, allocation, and
/// center-defined metadata policies.
#[derive(Debug, Clone, Copy, Default)]
pub struct GribFileBuilder<'a> {
    options: ReaderOptions,
    local_parameters: &'a [LocalParameterEntry<'a>],
    predefined_bitmaps: &'a [PredefinedBitmap<'a>],
}

impl<'a> GribFileBuilder<'a> {
    /// Controls whether malformed `GRIB` candidates abort scanning.
    ///
    /// Strict mode is enabled by default. Tolerant mode skips candidates with
    /// invalid framing or unsupported edition numbers, while still reporting
    /// semantic errors in well-framed GRIB1/GRIB2 messages.
    pub fn strict(mut self, strict: bool) -> Self {
        self.options.strict = strict;
        self
    }

    /// Sets the maximum number of points decoded into a field buffer.
    pub fn max_decoded_points(mut self, max: usize) -> Self {
        self.options.max_decoded_points = Some(max);
        self
    }

    /// Disables the decoded-field point limit.
    pub fn without_decoded_point_limit(mut self) -> Self {
        self.options.max_decoded_points = None;
        self
    }

    /// Sets the maximum number of coordinates allocated by an axis helper.
    pub fn max_axis_points(mut self, max: usize) -> Self {
        self.options.max_axis_points = Some(max);
        self
    }

    /// Disables the coordinate-axis point limit.
    pub fn without_axis_point_limit(mut self) -> Self {
        self.options.max_axis_points = None;
        self
    }

    /// Disables both decoded-field and coordinate-axis limits.
    pub fn without_limits(mut self) -> Self {
        self.options.max_decoded_points = None;
        self.options.max_axis_points = None;
        self
    }

    /// Supplies GRIB2 local parameter definitions used during indexing.
    pub fn local_parameters(mut self, entries: &'a [LocalParameterEntry<'a>]) -> Self {
        self.local_parameters = entries;
        self
    }

    /// Supplies center-defined GRIB1 predefined bitmaps.
    pub fn grib1_predefined_bitmaps(mut self, bitmaps: &'a [PredefinedBitmap<'a>]) -> Self {
        self.predefined_bitmaps = bitmaps;
        self
    }

    /// Opens a file through a read-only memory mapping.
    pub fn open<P: AsRef<Path>>(self, path: P) -> Result<GribFile> {
        let file = std::fs::File::open(path.as_ref())
            .map_err(|err| Error::Io(err, path.as_ref().display().to_string()))?;
        let mmap = unsafe { Mmap::map(&file) }
            .map_err(|err| Error::Io(err, path.as_ref().display().to_string()))?;
        GribFile::from_data(
            GribData::Mmap(mmap),
            self.options,
            self.local_parameters,
            self.predefined_bitmaps,
        )
    }

    /// Opens an owned in-memory GRIB buffer without copying it.
    pub fn from_bytes(self, data: Vec<u8>) -> Result<GribFile> {
        GribFile::from_data(
            GribData::Bytes(data),
            self.options,
            self.local_parameters,
            self.predefined_bitmaps,
        )
    }

    /// Reads GRIB data from an arbitrary byte stream.
    pub fn from_reader<R: Read>(self, mut reader: R) -> Result<GribFile> {
        let mut data = Vec::new();
        let mut chunk = [0u8; 64 * 1024];
        loop {
            let read = match reader.read(&mut chunk) {
                Ok(read) => read,
                Err(err) if err.kind() == std::io::ErrorKind::Interrupted => continue,
                Err(err) => return Err(Error::Io(err, "GRIB input stream".into())),
            };
            if read == 0 {
                break;
            }
            data.try_reserve(read)
                .map_err(|error| Error::allocation("GRIB input bytes", read, error))?;
            data.extend_from_slice(&chunk[..read]);
        }
        self.from_bytes(data)
    }
}

/// Caller-supplied GRIB1 predefined bitmap.
///
/// GRIB1 predefined bitmap table references are center-defined, not globally
/// standardized. `bitmap` is the raw bitmap payload: one bit per grid point in
/// GRIB scan order, with `1` meaning present and `0` meaning missing.
#[derive(Debug, Clone, Copy)]
pub struct PredefinedBitmap<'a> {
    pub center_id: u16,
    pub subcenter_id: Option<u16>,
    pub table_reference: u16,
    pub bitmap: &'a [u8],
}

/// A parsed GRIB field.
#[derive(Debug, Clone)]
pub struct MessageMetadata {
    pub edition: u8,
    pub center_id: u16,
    pub subcenter_id: u16,
    pub discipline: Option<u8>,
    pub reference_time: ReferenceTime,
    pub parameter: Parameter,
    pub grid: GridDefinition,
    pub data_representation: DataRepresentation,
    pub forecast_time_unit: Option<u8>,
    pub forecast_time: Option<i64>,
    pub message_offset: u64,
    pub message_length: u64,
    pub field_index_in_message: usize,
    grib1_product: Option<grib1::ProductDefinition>,
    grib2_identification: Option<Identification>,
    grib2_product: Option<ProductDefinition>,
}

/// A GRIB file containing one or more logical fields.
pub struct GribFile {
    data: GribData,
    messages: Vec<MessageIndex>,
    options: ReaderOptions,
}

#[derive(Clone)]
struct MessageIndex {
    offset: usize,
    length: usize,
    metadata: MessageMetadata,
    decode_plan: DecodePlan,
    row_group_lengths: Option<Arc<[usize]>>,
}

#[derive(Clone)]
enum DecodePlan {
    Grib1 {
        bitmap: Option<Grib1BitmapSource>,
        data: SectionRef,
    },
    Grib2(FieldSections),
}

#[derive(Clone)]
enum Grib1BitmapSource {
    Explicit(SectionRef),
    Predefined(Arc<[u8]>),
}

enum GribData {
    Mmap(Mmap),
    Bytes(Vec<u8>),
}

impl GribData {
    fn as_bytes(&self) -> &[u8] {
        match self {
            GribData::Mmap(m) => m,
            GribData::Bytes(b) => b,
        }
    }
}

impl GribFile {
    /// Starts a configurable GRIB open operation.
    pub fn builder<'a>() -> GribFileBuilder<'a> {
        GribFileBuilder::default()
    }

    /// Open a GRIB file from disk using memory-mapped I/O.
    pub fn open<P: AsRef<Path>>(path: P) -> Result<Self> {
        Self::builder().open(path)
    }

    /// Open a GRIB file from an owned byte buffer.
    pub fn from_bytes(data: Vec<u8>) -> Result<Self> {
        Self::builder().from_bytes(data)
    }

    fn from_data(
        data: GribData,
        options: ReaderOptions,
        local_parameters: &[LocalParameterEntry<'_>],
        predefined_bitmaps: &[PredefinedBitmap<'_>],
    ) -> Result<Self> {
        let messages = scan_messages(
            data.as_bytes(),
            options,
            local_parameters,
            predefined_bitmaps,
        )?;
        if messages.is_empty() {
            return Err(Error::NoMessages);
        }
        Ok(Self {
            data,
            messages,
            options,
        })
    }

    /// Returns the number of logical fields in the file.
    pub fn message_count(&self) -> usize {
        self.messages.len()
    }

    /// Access a field by index.
    pub fn message(&self, index: usize) -> Result<Message<'_>> {
        if index >= self.messages.len() {
            return Err(Error::MessageNotFound(index));
        }
        Ok(self.message_at(index))
    }

    fn message_at(&self, index: usize) -> Message<'_> {
        let record = &self.messages[index];
        let bytes = &self.data.as_bytes()[record.offset..record.offset + record.length];
        Message {
            index,
            bytes,
            metadata: &record.metadata,
            decode_plan: &record.decode_plan,
            row_group_lengths: record.row_group_lengths.as_deref(),
            options: self.options,
        }
    }

    /// Iterate over all fields.
    pub fn messages(&self) -> MessageIter<'_> {
        MessageIter {
            file: self,
            index: 0,
        }
    }

    /// Decode every field in the file.
    pub fn read_all_data_as_f64(&self) -> Result<Vec<ArrayD<f64>>> {
        #[cfg(feature = "rayon")]
        {
            (0..self.message_count())
                .into_par_iter()
                .map(|index| self.message_at(index).read_data_as_f64())
                .collect()
        }
        #[cfg(not(feature = "rayon"))]
        {
            (0..self.message_count())
                .map(|index| self.message_at(index).read_data_as_f64())
                .collect()
        }
    }

    /// Decode every field in the file as `f32`.
    pub fn read_all_data_as_f32(&self) -> Result<Vec<ArrayD<f32>>> {
        #[cfg(feature = "rayon")]
        {
            (0..self.message_count())
                .into_par_iter()
                .map(|index| self.message_at(index).read_data_as_f32())
                .collect()
        }
        #[cfg(not(feature = "rayon"))]
        {
            (0..self.message_count())
                .map(|index| self.message_at(index).read_data_as_f32())
                .collect()
        }
    }
}

/// A single logical GRIB field.
pub struct Message<'a> {
    bytes: &'a [u8],
    metadata: &'a MessageMetadata,
    decode_plan: &'a DecodePlan,
    row_group_lengths: Option<&'a [usize]>,
    options: ReaderOptions,
    index: usize,
}

impl<'a> Message<'a> {
    pub fn edition(&self) -> u8 {
        self.metadata.edition
    }

    pub fn index(&self) -> usize {
        self.index
    }

    pub fn metadata(&self) -> &MessageMetadata {
        self.metadata
    }

    pub fn reference_time(&self) -> &ReferenceTime {
        &self.metadata.reference_time
    }

    pub fn parameter(&self) -> &Parameter {
        &self.metadata.parameter
    }

    pub fn center_id(&self) -> u16 {
        self.metadata.center_id
    }

    pub fn subcenter_id(&self) -> u16 {
        self.metadata.subcenter_id
    }

    pub fn identification(&self) -> Option<&Identification> {
        self.metadata.grib2_identification.as_ref()
    }

    pub fn product_definition(&self) -> Option<&ProductDefinition> {
        self.metadata.grib2_product.as_ref()
    }

    pub fn grib1_product_definition(&self) -> Option<&grib1::ProductDefinition> {
        self.metadata.grib1_product.as_ref()
    }

    /// Returns the complete encoded GRIB2 Section 2, including its five-byte
    /// length and section-number header, when present for this field.
    pub fn local_use_section(&self) -> Option<&'a [u8]> {
        match self.decode_plan {
            DecodePlan::Grib2(field) => field
                .local_use
                .map(|section| section_bytes(self.bytes, section)),
            DecodePlan::Grib1 { .. } => None,
        }
    }

    pub fn grid_definition(&self) -> &GridDefinition {
        &self.metadata.grid
    }

    pub fn parameter_name(&self) -> &str {
        self.metadata.parameter.short_name.as_ref()
    }

    pub fn parameter_description(&self) -> &str {
        self.metadata.parameter.description.as_ref()
    }

    pub fn forecast_time_unit(&self) -> Option<u8> {
        self.metadata.forecast_time_unit
    }

    pub fn forecast_time_unit_kind(&self) -> Option<ForecastTimeUnit> {
        let unit = self.metadata.forecast_time_unit?;
        ForecastTimeUnit::from_edition_and_code(self.metadata.edition, unit)
    }

    /// Forecast offset in the units returned by [`Message::forecast_time_unit`].
    ///
    /// GRIB2 offsets are signed; GRIB1 offsets are widened to the same type.
    pub fn forecast_time(&self) -> Option<i64> {
        self.metadata.forecast_time
    }

    pub fn valid_time(&self) -> Option<ReferenceTime> {
        if let Some(end) = self
            .metadata
            .grib2_product
            .as_ref()
            .and_then(ProductDefinition::end_of_overall_time_interval)
        {
            return Some(end);
        }

        let unit = self.metadata.forecast_time_unit?;
        let lead = match &self.metadata.grib1_product {
            Some(product) => i64::from(product.valid_time_offset()?),
            None => self.metadata.forecast_time?,
        };
        self.metadata
            .reference_time
            .checked_add_forecast_time_by_edition(self.metadata.edition, unit, lead)
    }

    pub fn grid_shape(&self) -> (usize, usize) {
        self.metadata.grid.shape()
    }

    pub fn latitudes(&self) -> Result<Option<Vec<f64>>> {
        self.metadata
            .grid
            .latitudes_with_limit(self.options.max_axis_points)
    }

    pub fn longitudes(&self) -> Result<Option<Vec<f64>>> {
        self.metadata
            .grid
            .longitudes_with_limit(self.options.max_axis_points)
    }

    /// Latitude axis in the rotated coordinate system for Template 3.1.
    pub fn rotated_latitudes(&self) -> Result<Option<Vec<f64>>> {
        self.metadata
            .grid
            .as_rotated_lat_lon()
            .map(|grid| grid.rotated_latitudes_with_limit(self.options.max_axis_points))
            .transpose()
    }

    /// Longitude axis in the rotated coordinate system for Template 3.1.
    pub fn rotated_longitudes(&self) -> Result<Option<Vec<f64>>> {
        self.metadata
            .grid
            .as_rotated_lat_lon()
            .map(|grid| grid.rotated_longitudes_with_limit(self.options.max_axis_points))
            .transpose()
    }

    pub fn projected_x_coordinates(&self) -> Result<Option<Vec<f64>>> {
        self.metadata
            .grid
            .projected_x_coordinates_with_limit(self.options.max_axis_points)
    }

    pub fn projected_y_coordinates(&self) -> Result<Option<Vec<f64>>> {
        self.metadata
            .grid
            .projected_y_coordinates_with_limit(self.options.max_axis_points)
    }

    /// Decodes into row-major reader order.
    ///
    /// Alternating scan rows are normalized, but the GRIB i/j scan directions
    /// are preserved. Use [`Message::decode_into_north_up`] when deterministic
    /// +i columns and -j rows are required.
    /// Bitmap omissions and embedded complex-packing missing codes are written
    /// as `NaN`; substitute metadata is not applied automatically.
    pub fn decode_into<T: DecodeSample>(&self, out: &mut [T]) -> Result<()> {
        self.decode_scan_order_into(out)?;
        self.metadata.grid.reorder_for_ndarray_in_place(out)
    }

    /// Decodes into deterministic row-major orientation: +i columns and -j
    /// rows (north-to-south for conventional geographic grids).
    ///
    /// Bitmap omissions and embedded complex-packing missing codes are written
    /// as `NaN`; substitute metadata is not applied automatically.
    pub fn decode_into_north_up<T: DecodeSample>(&self, out: &mut [T]) -> Result<()> {
        self.decode_scan_order_into(out)?;
        self.metadata.grid.normalize_north_up_in_place(out)
    }

    fn decode_scan_order_into<T: DecodeSample>(&self, out: &mut [T]) -> Result<()> {
        self.metadata.grid.validate_supported_scan_order()?;
        let num_grid_points = self.checked_num_points()?;
        if out.len() != num_grid_points {
            return Err(Error::DataLengthMismatch {
                expected: num_grid_points,
                actual: out.len(),
            });
        }

        match &self.decode_plan {
            DecodePlan::Grib2(field) => {
                let data_section = section_bytes(self.bytes, field.data);
                let bitmap_section = match field.bitmap {
                    Some(section) => grib2_bitmap_payload(section_bytes(self.bytes, section))?,
                    None => None,
                };
                decode_field_into_with_group_lengths(
                    data_section,
                    &self.metadata.data_representation,
                    bitmap_section,
                    num_grid_points,
                    out,
                    self.row_group_lengths,
                )?
            }
            DecodePlan::Grib1 { bitmap, data } => {
                let bitmap_section = match bitmap {
                    Some(Grib1BitmapSource::Explicit(section)) => {
                        grib1::bitmap_payload(section_bytes(self.bytes, *section))?
                    }
                    Some(Grib1BitmapSource::Predefined(payload)) => Some(payload.as_ref()),
                    None => None,
                };
                let data_section = section_bytes(self.bytes, *data);
                if data_section.len() < 11 {
                    return Err(Error::InvalidSection {
                        section: 4,
                        reason: format!("expected at least 11 bytes, got {}", data_section.len()),
                    });
                }
                decode_payload_into_with_group_lengths(
                    &data_section[11..],
                    &self.metadata.data_representation,
                    bitmap_section,
                    num_grid_points,
                    out,
                    self.row_group_lengths,
                )?
            }
        }

        Ok(())
    }

    /// Returns a flat `f64` field in row-major reader order.
    ///
    /// Alternating rows are normalized, while encoded i/j directions are
    /// preserved. Use [`Message::read_flat_data_north_up_as_f64`] for +i
    /// columns and -j rows.
    pub fn read_flat_data_as_f64(&self) -> Result<Vec<f64>> {
        let num_grid_points = self.checked_num_points()?;
        let mut decoded = filled_vec(num_grid_points, 0.0_f64, "decoded f64 field")?;
        self.decode_into(&mut decoded)?;
        Ok(decoded)
    }

    /// Returns a flat `f32` field in row-major reader order.
    ///
    /// Alternating rows are normalized, while encoded i/j directions are
    /// preserved. Use [`Message::read_flat_data_north_up_as_f32`] for +i
    /// columns and -j rows.
    pub fn read_flat_data_as_f32(&self) -> Result<Vec<f32>> {
        let num_grid_points = self.checked_num_points()?;
        let mut decoded = filled_vec(num_grid_points, 0.0_f32, "decoded f32 field")?;
        self.decode_into(&mut decoded)?;
        Ok(decoded)
    }

    /// Returns a flat `f64` field with +i columns and -j rows.
    pub fn read_flat_data_north_up_as_f64(&self) -> Result<Vec<f64>> {
        let num_grid_points = self.checked_num_points()?;
        let mut decoded = filled_vec(num_grid_points, 0.0_f64, "decoded f64 field")?;
        self.decode_into_north_up(&mut decoded)?;
        Ok(decoded)
    }

    /// Returns a flat `f32` field with +i columns and -j rows.
    pub fn read_flat_data_north_up_as_f32(&self) -> Result<Vec<f32>> {
        let num_grid_points = self.checked_num_points()?;
        let mut decoded = filled_vec(num_grid_points, 0.0_f32, "decoded f32 field")?;
        self.decode_into_north_up(&mut decoded)?;
        Ok(decoded)
    }

    /// Returns an `f64` ndarray in row-major reader order.
    ///
    /// Alternating rows are normalized, while encoded i/j directions are
    /// preserved. Use [`Message::read_data_north_up_as_f64`] for deterministic
    /// +i columns and -j rows.
    pub fn read_data_as_f64(&self) -> Result<ArrayD<f64>> {
        let ordered = self.read_flat_data_as_f64()?;
        decoded_array(self.metadata.grid.ndarray_shape(), ordered)
    }

    /// Returns an `f32` ndarray in row-major reader order.
    ///
    /// Alternating rows are normalized, while encoded i/j directions are
    /// preserved. Use [`Message::read_data_north_up_as_f32`] for deterministic
    /// +i columns and -j rows.
    pub fn read_data_as_f32(&self) -> Result<ArrayD<f32>> {
        let ordered = self.read_flat_data_as_f32()?;
        decoded_array(self.metadata.grid.ndarray_shape(), ordered)
    }

    /// Returns an `f64` ndarray with +i columns and -j rows.
    pub fn read_data_north_up_as_f64(&self) -> Result<ArrayD<f64>> {
        let ordered = self.read_flat_data_north_up_as_f64()?;
        decoded_array(self.metadata.grid.ndarray_shape(), ordered)
    }

    /// Returns an `f32` ndarray with +i columns and -j rows.
    pub fn read_data_north_up_as_f32(&self) -> Result<ArrayD<f32>> {
        let ordered = self.read_flat_data_north_up_as_f32()?;
        decoded_array(self.metadata.grid.ndarray_shape(), ordered)
    }

    pub fn raw_bytes(&self) -> &[u8] {
        self.bytes
    }

    fn checked_num_points(&self) -> Result<usize> {
        let num_grid_points = self.metadata.grid.checked_num_points()?;
        ensure_limit(
            "decoded grid points",
            num_grid_points,
            self.options.max_decoded_points,
        )?;
        Ok(num_grid_points)
    }
}

fn decoded_array<T>(shape: Vec<usize>, values: Vec<T>) -> Result<ArrayD<T>> {
    ArrayD::from_shape_vec(IxDyn(&shape), values).map_err(|error| Error::InvalidSection {
        section: 3,
        reason: format!("grid shape cannot represent decoded field: {error}"),
    })
}

/// Iterator over fields in a GRIB file.
pub struct MessageIter<'a> {
    file: &'a GribFile,
    index: usize,
}

impl<'a> Iterator for MessageIter<'a> {
    type Item = Message<'a>;

    fn next(&mut self) -> Option<Self::Item> {
        if self.index >= self.file.message_count() {
            return None;
        }
        let message = self.file.message_at(self.index);
        self.index += 1;
        Some(message)
    }
}

fn section_bytes(msg_bytes: &[u8], section: SectionRef) -> &[u8] {
    &msg_bytes[section.offset..section.offset + section.length]
}

fn scan_messages(
    data: &[u8],
    options: ReaderOptions,
    local_parameters: &[LocalParameterEntry<'_>],
    predefined_bitmaps: &[PredefinedBitmap<'_>],
) -> Result<Vec<MessageIndex>> {
    let mut messages = Vec::new();
    let mut pos = 0usize;
    let finder = Finder::new(GRIB_MAGIC);

    while data.len().saturating_sub(pos) >= 8 {
        let Some(relative_pos) = finder.find(&data[pos..]) else {
            break;
        };
        pos = pos
            .checked_add(relative_pos)
            .ok_or(Error::ArithmeticOverflow {
                operation: "locating a GRIB signature",
            })?;
        if data.len().saturating_sub(pos) < 8 {
            break;
        }

        let (indicator, next_pos) = match locate_message(data, pos) {
            Ok(located) => located,
            Err(err) if !options.strict && is_recoverable_candidate_error(&err) => {
                pos = pos
                    .checked_add(GRIB_MAGIC.len())
                    .ok_or(Error::ArithmeticOverflow {
                        operation: "advancing past a rejected GRIB signature",
                    })?;
                continue;
            }
            Err(err) => return Err(err),
        };

        let message_bytes = &data[pos..next_pos];
        let indexed = match indicator.edition {
            1 => index_grib1_message(message_bytes, pos, predefined_bitmaps)?,
            2 => index_grib2_message(message_bytes, pos, &indicator, local_parameters)?,
            other => return Err(Error::UnsupportedEdition(other)),
        };

        messages.extend(indexed);
        pos = next_pos;
    }

    Ok(messages)
}

fn is_recoverable_candidate_error(err: &Error) -> bool {
    matches!(
        err,
        Error::InvalidMessage(_) | Error::Truncated { .. } | Error::UnsupportedEdition(_)
    )
}

fn locate_message(data: &[u8], pos: usize) -> Result<(Indicator, usize)> {
    let edition = Indicator::edition(&data[pos..]).ok_or_else(|| {
        Error::InvalidMessage(format!("failed to parse indicator at byte offset {pos}"))
    })?;
    if !matches!(edition, 1 | 2) {
        return Err(Error::UnsupportedEdition(edition));
    }

    let indicator = Indicator::parse(&data[pos..]).ok_or_else(|| {
        Error::InvalidMessage(format!("failed to parse indicator at byte offset {pos}"))
    })?;
    let length = usize::try_from(indicator.total_length).map_err(|_| {
        Error::InvalidMessage(format!(
            "message at byte offset {pos} reports a length that exceeds usize"
        ))
    })?;
    if length < 12 {
        return Err(Error::InvalidMessage(format!(
            "message at byte offset {pos} reports impossible length {length}"
        )));
    }
    let end = pos
        .checked_add(length)
        .ok_or_else(|| Error::InvalidMessage("message length overflow".into()))?;
    if end > data.len() {
        return Err(Error::Truncated { offset: end as u64 });
    }
    if &data[end - 4..end] != b"7777" {
        return Err(Error::InvalidMessage(format!(
            "message at byte offset {pos} does not end with 7777"
        )));
    }

    Ok((indicator, end))
}

fn index_grib1_message(
    message_bytes: &[u8],
    offset: usize,
    predefined_bitmaps: &[PredefinedBitmap<'_>],
) -> Result<Vec<MessageIndex>> {
    let sections = grib1::parse_message_sections(message_bytes)?;
    let grid_ref = sections.grid.ok_or_else(|| {
        Error::InvalidSectionOrder(
            "GRIB1 decoding requires an explicit grid definition section".into(),
        )
    })?;
    let grid_description = GridDescription::parse(section_bytes(message_bytes, grid_ref))?;
    let grid = grid_description.grid;
    let num_grid_points = grid.checked_num_points()?;

    let bitmap = match sections.bitmap {
        Some(section) => {
            let section_bytes = section_bytes(message_bytes, section);
            let table_reference = grib1::bitmap_table_reference(section_bytes)?;
            if table_reference == 0 {
                Some(Grib1BitmapSource::Explicit(section))
            } else {
                let payload = resolve_predefined_bitmap(
                    predefined_bitmaps,
                    sections.product.center_id as u16,
                    sections.product.subcenter_id as u16,
                    table_reference,
                )?;
                Some(Grib1BitmapSource::Predefined(Arc::from(payload)))
            }
        }
        None => None,
    };

    let bitmap_present_count = match &bitmap {
        Some(Grib1BitmapSource::Explicit(section)) => {
            let bitmap_payload = grib1::bitmap_payload(section_bytes(message_bytes, *section))?
                .ok_or(Error::MissingBitmap)?;
            count_bitmap_present_points(bitmap_payload, num_grid_points)?
        }
        Some(Grib1BitmapSource::Predefined(payload)) => {
            count_bitmap_present_points(payload, num_grid_points)?
        }
        None => num_grid_points,
    };

    let (_bds, data_representation) = BinaryDataSection::parse(
        section_bytes(message_bytes, sections.data),
        sections.product.decimal_scale,
        bitmap_present_count,
    )?;
    let bitmap_payload = match &bitmap {
        Some(Grib1BitmapSource::Explicit(section)) => {
            grib1::bitmap_payload(section_bytes(message_bytes, *section))?
        }
        Some(Grib1BitmapSource::Predefined(payload)) => Some(payload.as_ref()),
        None => None,
    };
    let row_group_lengths =
        row_group_lengths_for_representation(&grid, bitmap_payload, &data_representation)?;
    let parameter = sections.product.parameter();

    Ok(vec![MessageIndex {
        offset,
        length: message_bytes.len(),
        decode_plan: DecodePlan::Grib1 {
            bitmap,
            data: sections.data,
        },
        row_group_lengths,
        metadata: MessageMetadata {
            edition: 1,
            center_id: sections.product.center_id as u16,
            subcenter_id: sections.product.subcenter_id as u16,
            discipline: None,
            reference_time: sections.product.reference_time,
            parameter,
            grid,
            data_representation,
            forecast_time_unit: Some(sections.product.forecast_time_unit),
            forecast_time: sections.product.forecast_time().map(i64::from),
            message_offset: offset as u64,
            message_length: message_bytes.len() as u64,
            field_index_in_message: 0,
            grib1_product: Some(sections.product),
            grib2_identification: None,
            grib2_product: None,
        },
    }])
}

fn resolve_predefined_bitmap<'a>(
    predefined_bitmaps: &[PredefinedBitmap<'a>],
    center_id: u16,
    subcenter_id: u16,
    table_reference: u16,
) -> Result<&'a [u8]> {
    predefined_bitmaps
        .iter()
        .find(|entry| {
            entry.center_id == center_id
                && entry.subcenter_id == Some(subcenter_id)
                && entry.table_reference == table_reference
        })
        .or_else(|| {
            predefined_bitmaps.iter().find(|entry| {
                entry.center_id == center_id
                    && entry.subcenter_id.is_none()
                    && entry.table_reference == table_reference
            })
        })
        .map(|entry| entry.bitmap)
        .ok_or(Error::UnsupportedBitmapIndicator(table_reference))
}

fn row_group_lengths_for_representation(
    grid: &GridDefinition,
    bitmap_payload: Option<&[u8]>,
    representation: &DataRepresentation,
) -> Result<Option<Arc<[usize]>>> {
    let (expected_groups, expected_values, section) = match representation {
        DataRepresentation::Grib1SecondOrderPacking(params)
            if params.layout() == Grib1SecondOrderLayout::RowByRow =>
        {
            (params.group_count(), params.packing().encoded_values, 4)
        }
        DataRepresentation::ComplexPacking(params) if params.group_splitting_method == 0 => {
            (params.num_groups, params.encoded_values, 5)
        }
        _ => return Ok(None),
    };

    let num_grid_points = grid.checked_num_points()?;
    let (points_per_line, number_of_lines) = grid.scan_line_shape()?;
    let shape_points =
        points_per_line
            .checked_mul(number_of_lines)
            .ok_or(Error::ArithmeticOverflow {
                operation: "computing scan-line grid point count",
            })?;
    if shape_points != num_grid_points {
        return Err(Error::InvalidSection {
            section: 3,
            reason: format!(
                "scan-line shape {points_per_line}x{number_of_lines} contains {shape_points} points, but the grid declares {num_grid_points}"
            ),
        });
    }
    if number_of_lines != expected_groups {
        return Err(Error::InvalidSection {
            section,
            reason: format!(
                "row-wise packing declares {expected_groups} groups, but the grid has {number_of_lines} scan lines"
            ),
        });
    }

    let mut lengths = Vec::new();
    lengths
        .try_reserve_exact(number_of_lines)
        .map_err(|error| Error::allocation("row-wise group lengths", number_of_lines, error))?;
    match bitmap_payload {
        Some(bitmap) => {
            count_bitmap_present_points(bitmap, num_grid_points)?;
            for line_index in 0..number_of_lines {
                let line_start =
                    line_index
                        .checked_mul(points_per_line)
                        .ok_or(Error::ArithmeticOverflow {
                            operation: "computing bitmap scan-line offset",
                        })?;
                let mut present = 0usize;
                for point_index in 0..points_per_line {
                    if read_bit(bitmap, line_start + point_index)? {
                        present += 1;
                    }
                }
                lengths.push(present);
            }
        }
        None => lengths.resize(number_of_lines, points_per_line),
    }

    let actual_values = lengths.iter().try_fold(0usize, |total, length| {
        total.checked_add(*length).ok_or(Error::ArithmeticOverflow {
            operation: "summing row-wise group lengths",
        })
    })?;
    if actual_values != expected_values {
        return Err(Error::DataLengthMismatch {
            expected: expected_values,
            actual: actual_values,
        });
    }

    Ok(Some(Arc::from(lengths)))
}

fn validate_grib2_grid_points(grid: &GridDefinition) -> Result<()> {
    let Some(declared) = grid.declared_num_points() else {
        return Ok(());
    };
    let projected = grid.shape_num_points()?;
    if declared != projected {
        return Err(Error::InvalidSection {
            section: 3,
            reason: format!(
                "grid template 3.{} declares {declared} points but Nx*Ny is {projected}",
                grid.template_number()
            ),
        });
    }
    Ok(())
}

fn index_grib2_message(
    message_bytes: &[u8],
    offset: usize,
    indicator: &Indicator,
    local_parameters: &[LocalParameterEntry<'_>],
) -> Result<Vec<MessageIndex>> {
    let fields = index_fields(message_bytes)?;
    let mut messages = Vec::with_capacity(fields.len());

    for (field_index_in_message, field_sections) in fields.into_iter().enumerate() {
        let identification =
            Identification::parse(section_bytes(message_bytes, field_sections.identification))?;
        let grid_section = section_bytes(message_bytes, field_sections.grid);
        let grid = GridDefinition::parse(grid_section)?;
        validate_grib2_grid_points(&grid)?;
        let product =
            ProductDefinition::parse(section_bytes(message_bytes, field_sections.product))?;
        let data_representation = DataRepresentation::parse(section_bytes(
            message_bytes,
            field_sections.data_representation,
        ))?;
        let bitmap_payload = match field_sections.bitmap {
            Some(section) => grib2_bitmap_payload(section_bytes(message_bytes, section))?,
            None => None,
        };
        let row_group_lengths =
            row_group_lengths_for_representation(&grid, bitmap_payload, &data_representation)?;
        let parameter = lookup_parameter_with_local_entries(
            indicator.discipline,
            product.parameter_category,
            product.parameter_number,
            identification.center_id,
            identification.subcenter_id,
            identification.local_table_version,
            local_parameters,
        );

        messages.push(MessageIndex {
            offset,
            length: message_bytes.len(),
            decode_plan: DecodePlan::Grib2(field_sections),
            row_group_lengths,
            metadata: MessageMetadata {
                edition: indicator.edition,
                center_id: identification.center_id,
                subcenter_id: identification.subcenter_id,
                discipline: Some(indicator.discipline),
                reference_time: ReferenceTime {
                    year: identification.reference_year,
                    month: identification.reference_month,
                    day: identification.reference_day,
                    hour: identification.reference_hour,
                    minute: identification.reference_minute,
                    second: identification.reference_second,
                },
                parameter,
                grid,
                data_representation,
                forecast_time_unit: product.forecast_time_unit(),
                forecast_time: product.forecast_time().map(i64::from),
                message_offset: offset as u64,
                message_length: message_bytes.len() as u64,
                field_index_in_message,
                grib1_product: None,
                grib2_identification: Some(identification),
                grib2_product: Some(product),
            },
        });
    }

    Ok(messages)
}

#[cfg(test)]
mod tests {
    use super::*;

    fn grib_i32_bytes(value: i32) -> [u8; 4] {
        if value >= 0 {
            (value as u32).to_be_bytes()
        } else {
            ((-value) as u32 | 0x8000_0000).to_be_bytes()
        }
    }

    fn build_indicator(total_len: usize, discipline: u8) -> Vec<u8> {
        let mut indicator = Vec::with_capacity(16);
        indicator.extend_from_slice(b"GRIB");
        indicator.extend_from_slice(&[0, 0]);
        indicator.push(discipline);
        indicator.push(2);
        indicator.extend_from_slice(&(total_len as u64).to_be_bytes());
        indicator
    }

    fn build_identification() -> Vec<u8> {
        let mut section = vec![0u8; 21];
        section[..4].copy_from_slice(&(21u32).to_be_bytes());
        section[4] = 1;
        section[5..7].copy_from_slice(&7u16.to_be_bytes());
        section[7..9].copy_from_slice(&0u16.to_be_bytes());
        section[9] = 35;
        section[10] = 1;
        section[11] = 1;
        section[12..14].copy_from_slice(&2026u16.to_be_bytes());
        section[14] = 3;
        section[15] = 20;
        section[16] = 12;
        section[17] = 0;
        section[18] = 0;
        section[19] = 0;
        section[20] = 1;
        section
    }

    fn build_local_use(payload: &[u8]) -> Vec<u8> {
        let mut section = vec![0u8; payload.len() + 5];
        section[..4].copy_from_slice(&((payload.len() + 5) as u32).to_be_bytes());
        section[4] = 2;
        section[5..].copy_from_slice(payload);
        section
    }

    fn build_grid(ni: u32, nj: u32, scanning_mode: u8) -> Vec<u8> {
        let mut section = vec![0u8; 72];
        section[..4].copy_from_slice(&(72u32).to_be_bytes());
        section[4] = 3;
        section[6..10].copy_from_slice(&(ni * nj).to_be_bytes());
        section[12..14].copy_from_slice(&0u16.to_be_bytes());
        section[30..34].copy_from_slice(&ni.to_be_bytes());
        section[34..38].copy_from_slice(&nj.to_be_bytes());
        section[46..50].copy_from_slice(&grib_i32_bytes(50_000_000));
        section[50..54].copy_from_slice(&grib_i32_bytes(-120_000_000));
        section[55..59].copy_from_slice(&grib_i32_bytes(49_000_000));
        section[59..63].copy_from_slice(&grib_i32_bytes(-119_000_000));
        section[63..67].copy_from_slice(&1_000_000u32.to_be_bytes());
        section[67..71].copy_from_slice(&1_000_000u32.to_be_bytes());
        section[71] = scanning_mode;
        section
    }

    fn build_rotated_grid() -> Vec<u8> {
        let mut section = vec![0u8; 84];
        section[..4].copy_from_slice(&84u32.to_be_bytes());
        section[4] = 3;
        section[6..10].copy_from_slice(&6u32.to_be_bytes());
        section[12..14].copy_from_slice(&1u16.to_be_bytes());
        section[30..34].copy_from_slice(&3u32.to_be_bytes());
        section[34..38].copy_from_slice(&2u32.to_be_bytes());
        section[46..50].copy_from_slice(&grib_i32_bytes(-10_000_000));
        section[50..54].copy_from_slice(&grib_i32_bytes(-20_000_000));
        section[55..59].copy_from_slice(&grib_i32_bytes(0));
        section[59..63].copy_from_slice(&grib_i32_bytes(20_000_000));
        section[63..67].copy_from_slice(&20_000_000u32.to_be_bytes());
        section[67..71].copy_from_slice(&10_000_000u32.to_be_bytes());
        section[71] = 0b0100_0000;
        section[72..76].copy_from_slice(&grib_i32_bytes(-30_000_000));
        section[76..80].copy_from_slice(&10_000_000u32.to_be_bytes());
        section[80..84].copy_from_slice(&15.5f32.to_be_bytes());
        section
    }

    fn build_regular_gaussian_grid() -> Vec<u8> {
        let mut section = vec![0u8; 72];
        section[..4].copy_from_slice(&72u32.to_be_bytes());
        section[4] = 3;
        section[6..10].copy_from_slice(&16u32.to_be_bytes());
        section[12..14].copy_from_slice(&40u16.to_be_bytes());
        section[30..34].copy_from_slice(&4u32.to_be_bytes());
        section[34..38].copy_from_slice(&4u32.to_be_bytes());
        section[46..50].copy_from_slice(&grib_i32_bytes(59_444_408));
        section[50..54].copy_from_slice(&grib_i32_bytes(0));
        section[55..59].copy_from_slice(&grib_i32_bytes(-59_444_408));
        section[59..63].copy_from_slice(&grib_i32_bytes(270_000_000));
        section[63..67].copy_from_slice(&90_000_000u32.to_be_bytes());
        section[67..71].copy_from_slice(&2u32.to_be_bytes());
        section
    }

    fn build_product(parameter_category: u8, parameter_number: u8) -> Vec<u8> {
        let mut section = vec![0u8; 34];
        section[..4].copy_from_slice(&(34u32).to_be_bytes());
        section[4] = 4;
        section[7..9].copy_from_slice(&0u16.to_be_bytes());
        section[9] = parameter_category;
        section[10] = parameter_number;
        section[11] = 2;
        section[17] = 1;
        section[18..22].copy_from_slice(&0u32.to_be_bytes());
        section[22] = 103;
        section[23] = 0;
        section[24..28].copy_from_slice(&850u32.to_be_bytes());
        section[28] = 255;
        section
    }

    fn build_simple_representation(encoded_values: usize, bits_per_value: u8) -> Vec<u8> {
        let mut section = vec![0u8; 21];
        section[..4].copy_from_slice(&(21u32).to_be_bytes());
        section[4] = 5;
        section[5..9].copy_from_slice(&(encoded_values as u32).to_be_bytes());
        section[9..11].copy_from_slice(&0u16.to_be_bytes());
        section[11..15].copy_from_slice(&0f32.to_be_bytes());
        section[19] = bits_per_value;
        section[20] = 0;
        section
    }

    fn pack_u8_values(values: &[u8]) -> Vec<u8> {
        values.to_vec()
    }

    fn build_bitmap(bits: &[bool]) -> Vec<u8> {
        let payload_len = bits.len().div_ceil(8) + 1;
        let mut section = vec![0u8; payload_len + 5];
        section[..4].copy_from_slice(&((payload_len + 5) as u32).to_be_bytes());
        section[4] = 6;
        section[5] = 0;
        for (index, bit) in bits.iter().copied().enumerate() {
            if bit {
                section[6 + index / 8] |= 1 << (7 - (index % 8));
            }
        }
        section
    }

    fn build_reused_bitmap() -> Vec<u8> {
        let mut section = vec![0u8; 6];
        section[..4].copy_from_slice(&6u32.to_be_bytes());
        section[4] = 6;
        section[5] = 254;
        section
    }

    fn build_data(payload: &[u8]) -> Vec<u8> {
        let mut section = vec![0u8; payload.len() + 5];
        section[..4].copy_from_slice(&((payload.len() + 5) as u32).to_be_bytes());
        section[4] = 7;
        section[5..].copy_from_slice(payload);
        section
    }

    fn assemble_grib2_message(sections: &[Vec<u8>]) -> Vec<u8> {
        let total_len = 16 + sections.iter().map(|section| section.len()).sum::<usize>() + 4;
        let mut message = build_indicator(total_len, 0);
        for section in sections {
            message.extend_from_slice(section);
        }
        message.extend_from_slice(b"7777");
        message
    }

    fn build_grib1_message(values: &[u8]) -> Vec<u8> {
        let mut pds = vec![0u8; 28];
        pds[..3].copy_from_slice(&[0, 0, 28]);
        pds[3] = 2;
        pds[4] = 7;
        pds[5] = 255;
        pds[6] = 0;
        pds[7] = 0b1000_0000;
        pds[8] = 11;
        pds[9] = 100;
        pds[10..12].copy_from_slice(&850u16.to_be_bytes());
        pds[12] = 26;
        pds[13] = 3;
        pds[14] = 20;
        pds[15] = 12;
        pds[16] = 0;
        pds[17] = 1;
        pds[18] = 0;
        pds[19] = 0;
        pds[20] = 0;
        pds[24] = 21;
        pds[25] = 0;

        let mut gds = vec![0u8; 32];
        gds[..3].copy_from_slice(&[0, 0, 32]);
        gds[5] = 0;
        gds[6..8].copy_from_slice(&2u16.to_be_bytes());
        gds[8..10].copy_from_slice(&2u16.to_be_bytes());
        gds[10..13].copy_from_slice(&[0x01, 0x4d, 0x50]);
        gds[13..16].copy_from_slice(&[0x81, 0xd4, 0xc0]);
        gds[16] = 0x80;
        gds[17..20].copy_from_slice(&[0x01, 0x49, 0x68]);
        gds[20..23].copy_from_slice(&[0x81, 0xd0, 0xd8]);
        gds[23..25].copy_from_slice(&1000u16.to_be_bytes());
        gds[25..27].copy_from_slice(&1000u16.to_be_bytes());

        let mut bds = vec![0u8; 11 + values.len()];
        let len = bds.len() as u32;
        bds[..3].copy_from_slice(&[(len >> 16) as u8, (len >> 8) as u8, len as u8]);
        bds[3] = 0;
        bds[10] = 8;
        bds[11..].copy_from_slice(values);

        let total_len = 8 + pds.len() + gds.len() + bds.len() + 4;
        let mut message = Vec::new();
        message.extend_from_slice(b"GRIB");
        message.extend_from_slice(&[
            ((total_len >> 16) & 0xff) as u8,
            ((total_len >> 8) & 0xff) as u8,
            (total_len & 0xff) as u8,
            1,
        ]);
        message.extend_from_slice(&pds);
        message.extend_from_slice(&gds);
        message.extend_from_slice(&bds);
        message.extend_from_slice(b"7777");
        message
    }

    fn build_grib1_message_with_predefined_bitmap(values: &[u8], table_reference: u16) -> Vec<u8> {
        let mut message = build_grib1_message(values);
        message[8 + 7] = 0b1100_0000;

        let mut bitmap = [0u8; 6];
        bitmap[..3].copy_from_slice(&[0, 0, 6]);
        bitmap[4..6].copy_from_slice(&table_reference.to_be_bytes());

        let bitmap_offset = 8 + 28 + 32;
        message.splice(bitmap_offset..bitmap_offset, bitmap);
        let total_len = message.len();
        message[4] = ((total_len >> 16) & 0xff) as u8;
        message[5] = ((total_len >> 8) & 0xff) as u8;
        message[6] = (total_len & 0xff) as u8;
        message
    }

    fn build_grib1_second_order_message(layout: Grib1SecondOrderLayout) -> Vec<u8> {
        let mut message = build_grib1_message(&[0; 6]);
        let grid_offset = 8 + 28;
        message[grid_offset + 6..grid_offset + 8].copy_from_slice(&3u16.to_be_bytes());
        message[grid_offset + 8..grid_offset + 10].copy_from_slice(&2u16.to_be_bytes());

        let (length, extended_flags, n1, n2, widths, bitmap_offset) = match layout {
            Grib1SecondOrderLayout::RowByRow => (27usize, 0x10, 24u16, 26u16, &[2, 2][..], None),
            Grib1SecondOrderLayout::ConstantWidth => (27, 0x20, 24, 26, &[2][..], Some(22usize)),
            Grib1SecondOrderLayout::VariableWidth => (28, 0x30, 25, 27, &[2, 2][..], Some(23usize)),
            _ => unreachable!(),
        };
        let mut bds = vec![0u8; length];
        bds[..3].copy_from_slice(&[
            ((length >> 16) & 0xff) as u8,
            ((length >> 8) & 0xff) as u8,
            (length & 0xff) as u8,
        ]);
        bds[3] = 0x44;
        bds[10] = 5;
        bds[11..13].copy_from_slice(&n1.to_be_bytes());
        bds[13] = extended_flags;
        bds[14..16].copy_from_slice(&n2.to_be_bytes());
        bds[16..18].copy_from_slice(&2u16.to_be_bytes());
        bds[18..20].copy_from_slice(&6u16.to_be_bytes());
        bds[21..21 + widths.len()].copy_from_slice(widths);
        if let Some(offset) = bitmap_offset {
            bds[offset] = 0b1001_0000;
        }
        let first_order_index = usize::from(n1) - 1;
        bds[first_order_index..first_order_index + 2].copy_from_slice(&[0x55, 0x00]);
        let second_order_index = usize::from(n2) - 1;
        bds[second_order_index..second_order_index + 2].copy_from_slice(&[0x18, 0xb0]);

        let data_offset = 8 + 28 + 32;
        message.truncate(data_offset);
        message.extend_from_slice(&bds);
        message.extend_from_slice(b"7777");
        let total_len = message.len();
        message[4] = ((total_len >> 16) & 0xff) as u8;
        message[5] = ((total_len >> 8) & 0xff) as u8;
        message[6] = (total_len & 0xff) as u8;
        message
    }

    fn build_grib1_row_by_row_message_with_primary_bitmap() -> Vec<u8> {
        let mut message = build_grib1_second_order_message(Grib1SecondOrderLayout::RowByRow);
        let data_offset = 8 + 28 + 32;
        let mut bds = vec![0u8; 26];
        bds[..3].copy_from_slice(&[0, 0, 26]);
        bds[3] = 0x44;
        bds[10] = 5;
        bds[11..13].copy_from_slice(&24u16.to_be_bytes());
        bds[13] = 0x10;
        bds[14..16].copy_from_slice(&26u16.to_be_bytes());
        bds[16..18].copy_from_slice(&2u16.to_be_bytes());
        bds[18..20].copy_from_slice(&2u16.to_be_bytes());
        bds[21..23].copy_from_slice(&[2, 0]);
        bds[23..25].copy_from_slice(&[0x50, 0x00]);
        bds[25] = 0x20;

        message.truncate(data_offset);
        message[8 + 7] = 0b1100_0000;
        message.extend_from_slice(&[0, 0, 7, 2, 0, 0, 0b1010_0000]);
        message.extend_from_slice(&bds);
        message.extend_from_slice(b"7777");
        let total_len = message.len();
        message[4] = ((total_len >> 16) & 0xff) as u8;
        message[5] = ((total_len >> 8) & 0xff) as u8;
        message[6] = (total_len & 0xff) as u8;
        message
    }

    #[test]
    fn scans_single_grib2_message() {
        let message = assemble_grib2_message(&[
            build_identification(),
            build_grid(2, 2, 0),
            build_product(0, 0),
            build_simple_representation(4, 8),
            build_data(&pack_u8_values(&[1, 2, 3, 4])),
        ]);
        let messages = scan_messages(&message, ReaderOptions::default(), &[], &[]).unwrap();
        assert_eq!(messages.len(), 1);
        assert_eq!(messages[0].metadata.parameter.short_name, "TMP");
    }

    #[test]
    fn scans_across_large_non_grib_prefix() {
        let message = assemble_grib2_message(&[
            build_identification(),
            build_grid(2, 2, 0),
            build_product(0, 0),
            build_simple_representation(4, 8),
            build_data(&pack_u8_values(&[1, 2, 3, 4])),
        ]);
        let mut input = vec![0xa5; 1024 * 1024];
        input.extend_from_slice(&message);

        let messages = scan_messages(&input, ReaderOptions::default(), &[], &[]).unwrap();
        assert_eq!(messages.len(), 1);
        assert_eq!(messages[0].offset, 1024 * 1024);
    }

    #[test]
    fn decodes_simple_grib2_message_to_ndarray() {
        let message = assemble_grib2_message(&[
            build_identification(),
            build_grid(2, 2, 0),
            build_product(0, 0),
            build_simple_representation(4, 8),
            build_data(&pack_u8_values(&[1, 2, 3, 4])),
        ]);
        let file = GribFile::from_bytes(message).unwrap();
        let array = file.message(0).unwrap().read_data_as_f64().unwrap();
        assert_eq!(array.shape(), &[2, 2]);
        assert_eq!(
            array.iter().copied().collect::<Vec<_>>(),
            vec![1.0, 2.0, 3.0, 4.0]
        );
    }

    #[test]
    fn decodes_rotated_latlon_grid_and_exposes_rotated_axes() {
        let message = assemble_grib2_message(&[
            build_identification(),
            build_rotated_grid(),
            build_product(0, 0),
            build_simple_representation(6, 8),
            build_data(&pack_u8_values(&[1, 2, 3, 4, 5, 6])),
        ]);
        let file = GribFile::from_bytes(message).unwrap();
        let field = file.message(0).unwrap();

        assert_eq!(field.grid_shape(), (3, 2));
        assert_eq!(field.latitudes().unwrap(), None);
        assert_eq!(
            field.rotated_latitudes().unwrap().unwrap(),
            vec![-10.0, 0.0]
        );
        assert_eq!(
            field.rotated_longitudes().unwrap().unwrap(),
            vec![-20.0, 0.0, 20.0]
        );
        assert_eq!(
            field
                .read_data_as_f64()
                .unwrap()
                .iter()
                .copied()
                .collect::<Vec<_>>(),
            vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
        );
    }

    #[test]
    fn decodes_regular_gaussian_grid_and_computes_axes() {
        let values: Vec<_> = (0u8..16).collect();
        let message = assemble_grib2_message(&[
            build_identification(),
            build_regular_gaussian_grid(),
            build_product(0, 0),
            build_simple_representation(values.len(), 8),
            build_data(&pack_u8_values(&values)),
        ]);
        let file = GribFile::from_bytes(message).unwrap();
        let field = file.message(0).unwrap();

        assert_eq!(field.grid_shape(), (4, 4));
        assert_eq!(
            field.longitudes().unwrap().unwrap(),
            vec![0.0, 90.0, 180.0, 270.0]
        );
        let latitudes = field.latitudes().unwrap().unwrap();
        assert!((latitudes[0] - 59.444_408_289_166_77).abs() < 1e-12);
        assert!((latitudes[3] + 59.444_408_289_166_77).abs() < 1e-12);
        assert_eq!(field.read_data_as_f64().unwrap().shape(), &[4, 4]);
    }

    #[test]
    fn applies_bitmap_to_missing_values() {
        let message = assemble_grib2_message(&[
            build_identification(),
            build_grid(2, 2, 0),
            build_product(0, 1),
            build_simple_representation(3, 8),
            build_bitmap(&[true, false, true, true]),
            build_data(&pack_u8_values(&[10, 20, 30])),
        ]);
        let file = GribFile::from_bytes(message).unwrap();
        let array = file.message(0).unwrap().read_data_as_f64().unwrap();
        let values = array.iter().copied().collect::<Vec<_>>();
        assert_eq!(values[0], 10.0);
        assert!(values[1].is_nan());
        assert_eq!(values[2], 20.0);
        assert_eq!(values[3], 30.0);
    }

    #[test]
    fn indexes_multiple_fields_in_one_grib2_message() {
        let message = assemble_grib2_message(&[
            build_identification(),
            build_grid(2, 2, 0),
            build_product(0, 0),
            build_simple_representation(4, 8),
            build_data(&pack_u8_values(&[1, 2, 3, 4])),
            build_product(0, 2),
            build_simple_representation(4, 8),
            build_data(&pack_u8_values(&[5, 6, 7, 8])),
        ]);
        let file = GribFile::from_bytes(message).unwrap();
        assert_eq!(file.message_count(), 2);
        assert_eq!(file.message(0).unwrap().parameter_name(), "TMP");
        assert_eq!(file.message(1).unwrap().parameter_name(), "POT");
    }

    #[test]
    fn exposes_local_use_section_for_each_following_field() {
        let local_use = build_local_use(&[0xde, 0xad, 0xbe, 0xef]);
        let message = assemble_grib2_message(&[
            build_identification(),
            local_use.clone(),
            build_grid(2, 2, 0),
            build_product(0, 0),
            build_simple_representation(4, 8),
            build_data(&[1, 2, 3, 4]),
            build_product(0, 2),
            build_simple_representation(4, 8),
            build_data(&[5, 6, 7, 8]),
        ]);

        let file = GribFile::from_bytes(message).unwrap();
        assert_eq!(
            file.message(0).unwrap().local_use_section(),
            Some(local_use.as_slice())
        );
        assert_eq!(
            file.message(1).unwrap().local_use_section(),
            Some(local_use.as_slice())
        );
    }

    #[test]
    fn decodes_fields_that_reuse_a_previous_bitmap() {
        let message = assemble_grib2_message(&[
            build_identification(),
            build_grid(2, 2, 0),
            build_product(0, 0),
            build_simple_representation(2, 8),
            build_bitmap(&[true, false, true, false]),
            build_data(&[1, 2]),
            build_product(0, 2),
            build_simple_representation(2, 8),
            build_reused_bitmap(),
            build_data(&[3, 4]),
        ]);

        let file = GribFile::from_bytes(message).unwrap();
        let first = file.message(0).unwrap().read_flat_data_as_f64().unwrap();
        let second = file.message(1).unwrap().read_flat_data_as_f64().unwrap();
        assert_eq!([first[0], first[2]], [1.0, 2.0]);
        assert!(first[1].is_nan() && first[3].is_nan());
        assert_eq!([second[0], second[2]], [3.0, 4.0]);
        assert!(second[1].is_nan() && second[3].is_nan());
    }

    #[test]
    fn decodes_north_up_without_changing_default_reader_order() {
        let message = assemble_grib2_message(&[
            build_identification(),
            build_grid(2, 2, 0b1100_0000),
            build_product(0, 0),
            build_simple_representation(4, 8),
            build_data(&[1, 2, 3, 4]),
        ]);
        let file = GribFile::from_bytes(message).unwrap();
        let field = file.message(0).unwrap();

        assert_eq!(
            field.read_flat_data_as_f64().unwrap(),
            vec![1.0, 2.0, 3.0, 4.0]
        );
        assert_eq!(
            field.read_flat_data_north_up_as_f64().unwrap(),
            vec![4.0, 3.0, 2.0, 1.0]
        );
    }

    #[test]
    fn decodes_simple_grib1_message_to_ndarray() {
        let message = build_grib1_message(&[1, 2, 3, 4]);
        let file = GribFile::from_bytes(message).unwrap();
        let field = file.message(0).unwrap();
        assert_eq!(field.edition(), 1);
        assert_eq!(field.parameter_name(), "TMP");
        assert!(field.identification().is_none());
        assert!(field.grib1_product_definition().is_some());
        let array = field.read_data_as_f64().unwrap();
        assert_eq!(array.shape(), &[2, 2]);
        assert_eq!(
            array.iter().copied().collect::<Vec<_>>(),
            vec![1.0, 2.0, 3.0, 4.0]
        );
    }

    #[test]
    fn decodes_standard_grib1_second_order_layouts() {
        for layout in [
            Grib1SecondOrderLayout::RowByRow,
            Grib1SecondOrderLayout::ConstantWidth,
            Grib1SecondOrderLayout::VariableWidth,
        ] {
            let opened = GribFile::from_bytes(build_grib1_second_order_message(layout)).unwrap();
            let field = opened.message(0).unwrap();
            assert_eq!(field.grid_shape(), (3, 2));
            assert_eq!(
                field.read_flat_data_as_f64().unwrap(),
                vec![10.0, 11.0, 12.0, 20.0, 22.0, 23.0],
                "layout {layout:?}"
            );
            let DataRepresentation::Grib1SecondOrderPacking(params) =
                &field.metadata().data_representation
            else {
                panic!("expected GRIB1 second-order metadata");
            };
            assert_eq!(params.layout(), layout);
        }
    }

    #[test]
    fn decodes_grib1_row_by_row_groups_reduced_by_primary_bitmap() {
        let opened =
            GribFile::from_bytes(build_grib1_row_by_row_message_with_primary_bitmap()).unwrap();
        let values = opened.message(0).unwrap().read_flat_data_as_f64().unwrap();
        assert_eq!(values[0], 10.0);
        assert!(values[1].is_nan());
        assert_eq!(values[2], 12.0);
        assert!(values[3..].iter().all(|value| value.is_nan()));
    }

    #[test]
    fn decodes_grib1_predefined_bitmap_when_supplied() {
        let message = build_grib1_message_with_predefined_bitmap(&[9, 7], 300);
        let bitmap = [0b1010_0000];
        let predefined = [PredefinedBitmap {
            center_id: 7,
            subcenter_id: Some(0),
            table_reference: 300,
            bitmap: &bitmap,
        }];

        let file = GribFile::builder()
            .grib1_predefined_bitmaps(&predefined)
            .from_bytes(message)
            .unwrap();
        let decoded = file.message(0).unwrap().read_flat_data_as_f64().unwrap();
        assert_eq!(decoded[0], 9.0);
        assert!(decoded[1].is_nan());
        assert_eq!(decoded[2], 7.0);
        assert!(decoded[3].is_nan());
    }

    #[test]
    fn rejects_grib1_predefined_bitmap_without_supplied_table() {
        let message = build_grib1_message_with_predefined_bitmap(&[9, 7], 300);
        let err = match GribFile::from_bytes(message) {
            Ok(_) => panic!("expected unsupported predefined bitmap"),
            Err(err) => err,
        };

        assert!(matches!(err, Error::UnsupportedBitmapIndicator(300)));
    }
}