tensor/io/netcdf.go

// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT

package io

import (
	"encoding/binary"
	"maps"
	"math"
	"os"
	"slices"
	"strconv"
	"strings"

	"sourcedock.dev/petrbalvin/tensor/internal/base"
	"sourcedock.dev/petrbalvin/tensor/internal/core"
)

// NetCDF classic I/O. The Network Common Data Form is the archival
// format of climate and ocean science: a self-describing header of
// named dimensions, attributes and variables over big-endian binary
// payloads. This module speaks the classic model, CDF-1 on write and
// CDF-1 or CDF-2 on read. A record dimension (the unlimited first
// axis) interleaves its slabs across the record variables; a file
// carrying one is read, and written back, with the interleaving
// preserved. Variables land the core dtype their classic type code
// carries: NC_BYTE as int8, NC_CHAR as uint8 raw bytes, NC_SHORT as
// int16, NC_INT as int32, and NC_FLOAT and NC_DOUBLE as float64.
// CHAR carries bytes at the array level, never text; whether those
// bytes spell text is the caller's question, not the array's.

// NetCDF type codes as the file stores them.
const (
	ncTypeByte   = 1
	ncTypeChar   = 2
	ncTypeShort  = 3
	ncTypeInt    = 4
	ncTypeFloat  = 5
	ncTypeDouble = 6
)

// Section list tags as the file stores them.
const (
	ncTagDimension = 10
	ncTagVariable  = 11
	ncTagAttribute = 12
)

// NetCDFDim is one named dimension of a NetCDF classic file. A length
// of zero is the record dimension (the unlimited one): only it may be
// zero, it must lead the list, and its extent is the number of records,
// which a record variable carries on its first axis. Reading a file and
// writing it back preserves the record dimension.
type NetCDFDim struct {
	Name   string
	Length int
}

// NetCDFVar is one named variable of a NetCDF classic file. Values
// are row-major with the slowest dimension first, exactly as the file
// stores them, and Dims names the dimensions in that same order.
type NetCDFVar struct {
	Name   string
	Dims   []string
	Values *core.Array
	Attrs  map[string]string
}

// SaveNetCDF writes a NetCDF classic (CDF-1) file: dimensions with
// positive lengths plus at most one record dimension (length zero,
// declared first), float64 as NC_DOUBLE, float32 as NC_FLOAT and int64
// as NC_INT (the classic model's widest integer, so values outside the
// int32 range are an error rather than a silent truncation), and text
// attributes as NC_CHAR. A variable whose first dimension is the record
// dimension is a record variable: it must hold a whole number of
// records, every record variable must agree on that number, and the
// file interleaves one slab of each per record, padded to four bytes,
// the layout the classic model defines. Attribute keys are written in
// sorted order so the same inputs give the same bytes.
func SaveNetCDF(path string, dims []NetCDFDim, vars []NetCDFVar, attrs map[string]string) error {
	const name = "SaveNetCDF"
	recordDim := -1
	seen := make(map[string]bool, len(dims))
	for i, d := range dims {
		if err := ncCheckName(d.Name); err != nil {
			return base.Errf("%s: dimension %d: %w", name, i+1, err)
		}
		if seen[d.Name] {
			return base.Errf("%s: dimension %q is declared twice", name, d.Name)
		}
		seen[d.Name] = true
		if d.Length < 0 {
			return base.Errf("%s: dimension %q has negative length %d", name, d.Name, d.Length)
		}
		if d.Length == 0 {
			// The record dimension: at most one, and it leads the list,
			// as the classic model requires.
			if i != 0 {
				return base.Errf("%s: dimension %q has length 0 but is not the first: the record dimension leads the list",
					name, d.Name)
			}
			recordDim = i
		}
	}
	dimIndex := make(map[string]int, len(dims))
	for i, d := range dims {
		dimIndex[d.Name] = i
	}
	seenVar := make(map[string]bool, len(vars))
	recordVars := make([]bool, len(vars))
	perRecord := make([]int, len(vars))
	numrecs := 0
	numrecsSet := false
	for i := range vars {
		v := &vars[i]
		if err := ncCheckName(v.Name); err != nil {
			return base.Errf("%s: variable %d: %w", name, i+1, err)
		}
		if seenVar[v.Name] {
			return base.Errf("%s: variable %q is declared twice", name, v.Name)
		}
		seenVar[v.Name] = true
		for j, dn := range v.Dims {
			di, ok := dimIndex[dn]
			if !ok {
				return base.Errf("%s: variable %q refers to undeclared dimension %q", name, v.Name, dn)
			}
			if dims[di].Length == 0 && j != 0 {
				return base.Errf("%s: variable %q carries the record dimension %q in position %d: a record variable's record axis leads its dimensions",
					name, v.Name, dn, j)
			}
		}
		if _, _, err := ncExternal(v.Values.Dtype()); err != nil {
			return base.Errf("%s: variable %q: %w", name, v.Name, err)
		}
		if v.Values.Dtype() == core.Int {
			// NC_INT is a signed 32-bit value and the classic model has
			// no wider integer, so a value outside that range is an
			// error rather than a silent truncation.
			for _, x := range v.Values.RawInts() {
				if x < math.MinInt32 || x > math.MaxInt32 {
					return base.Errf("%s: variable %q holds %d, outside the int32 range the classic model stores",
						name, v.Name, x)
				}
			}
		}
		isRecord := recordDim >= 0 && len(v.Dims) > 0 && v.Dims[0] == dims[recordDim].Name
		if isRecord {
			per := 1
			for _, dn := range v.Dims[1:] {
				per *= dims[dimIndex[dn]].Length
			}
			if per <= 0 || v.Values.Len()%per != 0 {
				return base.Errf("%s: record variable %q holds %d values, not a whole number of records of %d",
					name, v.Name, v.Values.Len(), per)
			}
			recs := v.Values.Len() / per
			if numrecsSet && recs != numrecs {
				return base.Errf("%s: record variable %q holds %d records, the others %d",
					name, v.Name, recs, numrecs)
			}
			numrecs, numrecsSet = recs, true
			perRecord[i], recordVars[i] = per, true
			continue
		}
		want := 1
		for _, dn := range v.Dims {
			want *= dims[dimIndex[dn]].Length
		}
		if v.Values.Len() != want {
			return base.Errf("%s: variable %q holds %d values, its dimensions hold %d",
				name, v.Name, v.Values.Len(), want)
		}
	}

	// The image's size follows from the plan before a byte is written:
	// the data section exactly, the header from the names and the
	// attribute texts with a fixed frame each. The buffer is allocated
	// once from that, so the appends below fill it instead of growing
	// and copying the whole file as they go.
	buf := make([]byte, 0, ncImageEstimate(dims, vars, attrs, perRecord, recordVars, numrecs))
	buf = append(buf, 'C', 'D', 'F', 1)
	buf = binary.BigEndian.AppendUint32(buf, uint32(numrecs))
	// dim_list
	if len(dims) == 0 {
		buf = append(buf, 0, 0, 0, 0, 0, 0, 0, 0) // ABSENT
	} else {
		buf = binary.BigEndian.AppendUint32(buf, ncTagDimension)
		buf = binary.BigEndian.AppendUint32(buf, uint32(len(dims)))
		for _, d := range dims {
			buf = ncAppendName(buf, d.Name)
			buf = binary.BigEndian.AppendUint32(buf, uint32(d.Length))
		}
	}

	// gatt_list
	if len(attrs) == 0 {
		buf = append(buf, 0, 0, 0, 0, 0, 0, 0, 0) // ABSENT
	} else {
		buf = binary.BigEndian.AppendUint32(buf, ncTagAttribute)
		buf = binary.BigEndian.AppendUint32(buf, uint32(len(attrs)))
		for _, k := range slices.Sorted(maps.Keys(attrs)) {
			buf = ncAppendName(buf, k)
			buf = binary.BigEndian.AppendUint32(buf, ncTypeChar)
			buf = binary.BigEndian.AppendUint32(buf, uint32(len(attrs[k])))
			buf = append(buf, attrs[k]...)
			buf = ncPad4(buf)
		}
	}

	// var_list. The begin field is the last word of each record, so
	// its position is known as the record is appended; the absolute
	// data offsets follow once the whole header is in place.
	buf = binary.BigEndian.AppendUint32(buf, ncTagVariable)
	buf = binary.BigEndian.AppendUint32(buf, uint32(len(vars)))
	beginAt := make([]int, len(vars))
	for i := range vars {
		v := &vars[i]
		buf = ncAppendName(buf, v.Name)
		buf = binary.BigEndian.AppendUint32(buf, uint32(len(v.Dims)))
		for _, dn := range v.Dims {
			buf = binary.BigEndian.AppendUint32(buf, uint32(dimIndex[dn]))
		}
		if len(v.Attrs) == 0 {
			buf = append(buf, 0, 0, 0, 0, 0, 0, 0, 0) // ABSENT
		} else {
			buf = binary.BigEndian.AppendUint32(buf, ncTagAttribute)
			buf = binary.BigEndian.AppendUint32(buf, uint32(len(v.Attrs)))
			for _, k := range slices.Sorted(maps.Keys(v.Attrs)) {
				buf = ncAppendName(buf, k)
				buf = binary.BigEndian.AppendUint32(buf, ncTypeChar)
				buf = binary.BigEndian.AppendUint32(buf, uint32(len(v.Attrs[k])))
				buf = append(buf, v.Attrs[k]...)
				buf = ncPad4(buf)
			}
		}
		_, size, _ := ncExternal(v.Values.Dtype())
		// A record variable's vsize is one record's worth, padded to a
		// four-byte boundary; a fixed variable's is its whole payload.
		bytes := v.Values.Len() * size
		if recordVars[i] {
			bytes = int(ncPad4Size(int64(perRecord[i]) * int64(size)))
		}
		if uint64(bytes) > math.MaxUint32 {
			return base.Errf("%s: variable %q needs %d bytes, past what a CDF-1 size word can hold",
				name, v.Name, bytes)
		}
		buf = binary.BigEndian.AppendUint32(buf, uint32(ncExternalType(v.Values.Dtype())))
		buf = binary.BigEndian.AppendUint32(buf, uint32(bytes))
		beginAt[i] = len(buf)
		buf = binary.BigEndian.AppendUint32(buf, 0) // patched below
	}

	// data: the fixed variables first, in header order, each padded to
	// a four-byte boundary, then the records. Each record carries the
	// slab of every record variable, again in header order, each slab
	// padded to four bytes, so a record variable's bytes are strided by
	// the whole record size.
	offset := uint32(len(buf))
	for i := range vars {
		if recordVars[i] {
			continue
		}
		v := &vars[i]
		binary.BigEndian.PutUint32(buf[beginAt[i]:], offset)
		start := len(buf)
		payload, err := ncAppendPayload(buf, v.Values, 0, v.Values.Len())
		if err != nil {
			return base.Errf("%s: variable %q: %w", name, v.Name, err)
		}
		buf = payload
		buf = ncPad4(buf)
		offset += uint32(len(buf) - start)
	}
	// A record variable's begin points at its own slab inside the first
	// record, which is what readers (and the model's own arithmetic)
	// expect: record r's slab of that variable then sits a whole record
	// size further on, at begin + r*recordSize.
	slabOffset := offset
	for i := range vars {
		if !recordVars[i] {
			continue
		}
		binary.BigEndian.PutUint32(buf[beginAt[i]:], slabOffset)
		slabOffset += uint32(ncPad4Size(int64(perRecord[i]) * int64(elemSize(vars[i].Values))))
	}
	for rec := range numrecs {
		for i := range vars {
			if !recordVars[i] {
				continue
			}
			start := len(buf)
			lo := rec * perRecord[i]
			payload, err := ncAppendPayload(buf, vars[i].Values, lo, lo+perRecord[i])
			if err != nil {
				return base.Errf("%s: variable %q: %w", name, vars[i].Name, err)
			}
			buf = payload
			for (len(buf)-start)%4 != 0 {
				buf = append(buf, 0)
			}
		}
	}
	// CDF-1 addresses its offsets with a signed 32-bit word, so a file
	// at or past 2 GiB cannot be addressed; refuse rather than wrap the
	// begin fields.
	if uint64(len(buf)) > 1<<31-1 {
		return base.Errf("%s: the file would be %d bytes, past the 2 GiB a CDF-1 offset can address", name, len(buf))
	}
	return os.WriteFile(path, buf, 0o644)
}

// ncImageEstimate bounds the byte size of the file image SaveNetCDF
// builds. The data section is exact: each payload's four-byte-padded
// extent, with a record variable's slab repeated once per record. The
// header is charged a fixed frame per dimension, variable, dimension
// reference and attribute plus twice the bytes of every name and
// attribute text, which is past what the tagged lists cost. The sum is
// kept in int64 and clamped to the 2 GiB a CDF-1 offset can address,
// so no declaration can wrap it into a small or negative capacity.
func ncImageEstimate(dims []NetCDFDim, vars []NetCDFVar, attrs map[string]string, perRecord []int, recordVars []bool, numrecs int) int {
	const budget = int64(1) << 31
	total := int64(2048)
	add := func(n int64) {
		if total += n; total > budget {
			total = budget
		}
	}
	nameBytes := func(k, v string) int64 { return int64(2*(len(k)+len(v))) + 64 }
	for i := range vars {
		size := int64(elemSize(vars[i].Values))
		if recordVars[i] {
			add(ncPad4Size(int64(perRecord[i])*size) * int64(numrecs))
		} else {
			add(ncPad4Size(int64(vars[i].Values.Len()) * size))
		}
	}
	for _, d := range dims {
		add(nameBytes(d.Name, ""))
	}
	for k, v := range attrs {
		add(nameBytes(k, v))
	}
	for i := range vars {
		add(nameBytes(vars[i].Name, "") + 4*int64(len(vars[i].Dims)) + 128)
		for k, v := range vars[i].Attrs {
			add(nameBytes(k, v))
		}
	}
	return int(total)
}

// elemSize returns the external width of an array's dtype; ncExternal
// has already accepted it by the time the writer calls this.
func elemSize(a *core.Array) int {
	_, size, err := ncExternal(a.Dtype())
	if err != nil {
		return 0
	}
	return size
}

// ncAppendPayload appends arr's elements [start, end) in the external
// format of its dtype: the slot is reserved once and each element is
// stored into it, so no element pays for the append's own capacity
// check. Only the dtypes the classic model stores are accepted: a
// narrower array carries no integer payload for the int branch to
// read, so any other dtype is a named error rather than a bounds
// failure or a silent cast. SaveNetCDF's validation refuses the same
// dtypes first, so the error is a defence in depth for direct callers.
func ncAppendPayload(buf []byte, arr *core.Array, start, end int) ([]byte, error) {
	switch arr.Dtype() {
	case core.Float:
		vals := arr.RawFloats()[start:end]
		at := len(buf)
		buf = append(buf, make([]byte, 8*len(vals))...)
		for i, x := range vals {
			binary.BigEndian.PutUint64(buf[at+i*8:], math.Float64bits(x))
		}
	case core.Float32:
		vals := arr.RawFloat32s()[start:end]
		at := len(buf)
		buf = append(buf, make([]byte, 4*len(vals))...)
		for i, x := range vals {
			binary.BigEndian.PutUint32(buf[at+i*4:], math.Float32bits(x))
		}
	case core.Int:
		vals := arr.RawInts()[start:end]
		at := len(buf)
		buf = append(buf, make([]byte, 4*len(vals))...)
		for i, x := range vals {
			binary.BigEndian.PutUint32(buf[at+i*4:], uint32(int32(x)))
		}
	default:
		return nil, base.Errf("SaveNetCDF: cannot store dtype %s in a NetCDF classic file; convert with Astype", arr.Dtype())
	}
	return buf, nil
}

// LoadNetCDF reads a NetCDF classic file (CDF-1 and CDF-2), returning
// its dimensions, its variables and its global attributes. Every
// variable lands the core dtype its classic type code carries:
// NC_BYTE as int8 (signed in the classic model), NC_CHAR as uint8 raw
// bytes (CHAR carries bytes at the array level, never text), NC_SHORT
// as int16, NC_INT as int32, and NC_FLOAT and NC_DOUBLE as float64,
// where every classic value is exact. A type code beyond the classic
// six is refused by name: the unsigned and 64-bit codes exist only in
// formats this module does not speak. A variable that lands a narrow
// dtype here is refused by SaveNetCDF, whose writer stores float64,
// float32 and int64 arrays only; convert with Astype first. Variable
// and dimension attributes are returned on the variables themselves;
// only the global attributes come back in the map. A record dimension
// comes back with length zero, and each record variable's first axis
// carries the record count the file declares.
func LoadNetCDF(path string) ([]NetCDFDim, []NetCDFVar, map[string]string, error) {
	data, err := os.ReadFile(path)
	if err != nil {
		return nil, nil, nil, base.Errf("LoadNetCDF: %w", err)
	}
	return parseNetCDF(data)
}

// maxIndex is the largest int this platform holds: a declared
// dimension length beyond it cannot appear in the shapes the core
// builds.
const maxIndex = int64(^uint(0) >> 1)

// ncReader walks the file bytes with bounds checks, recording the
// first violation so the parse can read straight through a structure
// and report once at the end.
type ncReader struct {
	data []byte
	pos  int
	err  error
}

func (r *ncReader) fail(format string, args ...any) {
	if r.err == nil {
		r.err = base.Errf("LoadNetCDF: "+format, args...)
	}
}

func (r *ncReader) take(n int) []byte { return r.take64(int64(n)) }

// take64 is take with the arithmetic in int64, so a length read from
// the header cannot wrap on a 32-bit platform before the bounds check
// has seen it.
func (r *ncReader) take64(n int64) []byte {
	if r.err != nil {
		return nil
	}
	if n < 0 || n > int64(len(r.data)-r.pos) {
		r.fail("file ends inside the header at byte %d", r.pos)
		return nil
	}
	b := r.data[r.pos : r.pos+int(n)]
	r.pos += int(n)
	return b
}

// boundedCount caps a record count declared in the header by the bytes
// that could possibly hold it: every record of a list is at least
// minBytes long, so n of them cannot fit into what is left of the
// file. A header that lies about its sizes fails here instead of
// making the reader allocate the claimed amount, which for a uint32
// count is tens of gigabytes from a file of a few bytes.
func (r *ncReader) boundedCount(n uint32, minBytes int64, what string) int {
	if r.err != nil {
		return 0
	}
	left := int64(len(r.data) - r.pos)
	if int64(n) > left/minBytes {
		r.fail("the header declares %d %s, more than the %d remaining bytes can hold", n, what, left)
		return 0
	}
	return int(n)
}

func (r *ncReader) u32() uint32 {
	b := r.take(4)
	if b == nil {
		return 0
	}
	return binary.BigEndian.Uint32(b)
}

func (r *ncReader) u64() uint64 {
	b := r.take(8)
	if b == nil {
		return 0
	}
	return binary.BigEndian.Uint64(b)
}

// name reads a length-prefixed, 4-byte-padded name.
func (r *ncReader) name() string {
	n := r.u32()
	b := r.take64(int64(n) + int64(ncPadLen(int(n))))
	if b == nil {
		return ""
	}
	return string(b[:n])
}

func parseNetCDF(data []byte) ([]NetCDFDim, []NetCDFVar, map[string]string, error) {
	r := &ncReader{data: data}
	magic := r.take(4)
	if r.err != nil {
		return nil, nil, nil, r.err
	}
	if string(magic[:3]) != "CDF" {
		return nil, nil, nil, base.Errf("LoadNetCDF: not a NetCDF file, the magic is %q", magic[:3])
	}
	wide := false
	switch magic[3] {
	case 1:
	case 2:
		wide = true // 64-bit offsets; numrecs stays 32-bit
	default:
		return nil, nil, nil, base.Errf("LoadNetCDF: unsupported NetCDF version %d", magic[3])
	}
	numrecs := int64(r.u32())

	dims := r.dimList()
	attrs := r.attList()
	vars := r.varList(wide, dims)
	if r.err != nil {
		return nil, nil, nil, r.err
	}

	// The record variables' slabs interleave, so a record variable's
	// bytes are not contiguous: the span check and the gather both need
	// the size of a whole record, which is the sum of every record
	// variable's slab padded to a four-byte boundary. Each factor is
	// bounded per variable here, exactly as the read below bounds it
	// before it multiplies: without that, one hostile declaration
	// inflates the record size for every other record variable, ones
	// already validated included, and the wrapped span that follows
	// walks the reader off the end of the file.
	slabs := make([]int64, len(vars))
	recordSize := int64(0)
	for i, v := range vars {
		if !isRecordVar(dims, v.dimIDs) {
			continue
		}
		width, ok := ncTypeWidth(v.ncType)
		if !ok {
			continue // the read below reports the unknown type
		}
		per := int64(1)
		for _, did := range v.dimIDs[1:] {
			l := int64(dims[did].Length)
			if l <= 0 || per > int64(len(r.data))/l {
				return nil, nil, nil, base.Errf("LoadNetCDF: variable %q declares more elements than the file holds", v.name)
			}
			per *= l
		}
		slabs[i] = ncPad4Size(per * width)
		if slabs[i] > math.MaxInt64-recordSize {
			return nil, nil, nil, base.Errf("LoadNetCDF: the record variables declare more bytes per record than can be addressed")
		}
		recordSize += slabs[i]
	}

	out := make([]NetCDFVar, len(vars))
	for i, v := range vars {
		isRecord := isRecordVar(dims, v.dimIDs)
		// The element count is a product of header fields, so it can
		// wrap: a variable with more elements than the file has bytes
		// cannot be backed by it, and the division keeps the product
		// itself inside int64 while the check runs. The record axis
		// counts once here; its extent is numrecs, applied below.
		count := int64(1)
		first := 0
		if isRecord {
			first = 1
		}
		for _, did := range v.dimIDs[first:] {
			l := int64(dims[did].Length)
			if l <= 0 || count > int64(len(r.data))/l {
				return nil, nil, nil, base.Errf("LoadNetCDF: variable %q declares more elements than the file holds", v.name)
			}
			count *= l
		}
		perRecord := count
		if isRecord {
			if numrecs > int64(len(r.data)) || (perRecord > 0 && numrecs > int64(len(r.data))/perRecord) {
				return nil, nil, nil, base.Errf("LoadNetCDF: variable %q declares more elements than the file holds", v.name)
			}
			count = perRecord * numrecs
		}
		names := make([]string, len(v.dimIDs))
		shape := make([]int, len(v.dimIDs))
		for j, did := range v.dimIDs {
			names[j] = dims[did].Name
			shape[j] = dims[did].Length
			if isRecord && j == 0 {
				shape[j] = int(numrecs)
			}
		}
		// A rank-0 (scalar) variable has no shape at all, and an Array
		// always carries at least one dimension, so the scalar comes
		// back as a one-element vector with no dimension names.
		if len(shape) == 0 {
			shape = []int{1}
		}
		// The payload is allocated once in the dtype the type code
		// lands and the decode fills it in place, so no element is
		// copied twice and no integer value is rounded by a widening.
		landing, ok := ncLandingDtype(v.ncType)
		if !ok {
			return nil, nil, nil, base.Errf("LoadNetCDF: variable %q has unknown type %d", v.name, v.ncType)
		}
		arr := core.New(landing, shape...)
		if arr == nil {
			return nil, nil, nil, base.Errf("LoadNetCDF: variable %q: shape %v cannot build an array", v.name, shape)
		}
		if isRecord {
			r.valuesRecord(arr, v.begin, v.ncType, perRecord, numrecs, recordSize)
		} else {
			r.values(arr, v.begin, v.ncType, count)
		}
		if r.err != nil {
			return nil, nil, nil, r.err
		}
		out[i] = NetCDFVar{Name: v.name, Dims: names, Values: arr, Attrs: v.attrs}
	}
	return dims, out, attrs, nil
}

// ncVar is the parsed header record of one variable.
type ncVar struct {
	name   string
	dimIDs []int
	attrs  map[string]string
	ncType uint32
	begin  int64
}

func (r *ncReader) dimList() []NetCDFDim {
	if r.err != nil {
		return nil
	}
	if r.u32() == 0 {
		r.u32() // ABSENT is two zero words
		return nil
	}
	n := r.boundedCount(r.u32(), 8, "dimensions")
	dims := make([]NetCDFDim, 0, n)
	seen := map[string]bool{}
	for range n {
		name := r.name()
		length := r.u32()
		if r.err != nil {
			return nil
		}
		// Names must be unique: a variable refers to dimensions by
		// index, and with a duplicated name its own name list would
		// resolve differently by position and by name, so any answer
		// would be a silent guess.
		if seen[name] {
			r.fail("dimension %q is declared twice", name)
			return nil
		}
		seen[name] = true
		if length == 0 {
			// The record dimension: exactly one, and only the first.
			// Its extent is the number of records, which each record
			// variable carries on its leading axis.
			if len(dims) != 0 {
				r.fail("dimension %q has length 0 but is not the first: only one record dimension exists, and it leads the list", name)
				return nil
			}
			dims = append(dims, NetCDFDim{Name: name, Length: 0})
			continue
		}
		// A declared length must fit an int on every platform the
		// library targets. Whether the file actually holds the data is
		// decided per variable, where the elements are read.
		if int64(length) > maxIndex {
			r.fail("dimension %q is declared %d long, beyond this platform's index range", name, length)
			return nil
		}
		dims = append(dims, NetCDFDim{Name: name, Length: int(length)})
	}
	return dims
}

func (r *ncReader) attList() map[string]string {
	if r.err != nil {
		return nil
	}
	if r.u32() == 0 {
		r.u32() // ABSENT
		return nil
	}
	n := r.boundedCount(r.u32(), 12, "global attributes")
	attrs := make(map[string]string, n)
	seen := map[string]bool{}
	for range n {
		name := r.name()
		ncType := r.u32()
		count := r.boundedCount(r.u32(), 1, "attribute elements")
		if r.err != nil {
			return nil
		}
		if seen[name] {
			r.fail("attribute %q is declared twice", name)
			return nil
		}
		seen[name] = true
		attrs[name] = r.attValue(ncType, count)
		if r.err != nil {
			return nil
		}
	}
	return attrs
}

// attValue reads one attribute payload. Text comes back verbatim;
// numeric types, which the classic model allows here, are formatted as
// decimal so callers see the value the file carries either way.
func (r *ncReader) attValue(ncType uint32, count int) string {
	if r.err != nil {
		return ""
	}
	switch ncType {
	case ncTypeChar:
		b := r.take(count + ncPadLen(count))
		if b == nil {
			return ""
		}
		return string(b[:count])
	case ncTypeByte:
		b := r.take(count + ncPadLen(count))
		if b == nil {
			return ""
		}
		parts := make([]string, count)
		for i := range count {
			parts[i] = strconv.FormatInt(int64(int8(b[i])), 10)
		}
		return strings.Join(parts, ", ")
	case ncTypeShort:
		b := r.take(count*2 + ncPadLen(count*2))
		if b == nil {
			return ""
		}
		parts := make([]string, count)
		for i := range count {
			parts[i] = strconv.FormatInt(int64(int16(binary.BigEndian.Uint16(b[i*2:]))), 10)
		}
		return strings.Join(parts, ", ")
	case ncTypeInt:
		b := r.take(count * 4)
		if b == nil {
			return ""
		}
		parts := make([]string, count)
		for i := range count {
			parts[i] = strconv.FormatInt(int64(int32(binary.BigEndian.Uint32(b[i*4:]))), 10)
		}
		return strings.Join(parts, ", ")
	case ncTypeFloat:
		b := r.take(count * 4)
		if b == nil {
			return ""
		}
		parts := make([]string, count)
		for i := range count {
			parts[i] = strconv.FormatFloat(float64(math.Float32frombits(binary.BigEndian.Uint32(b[i*4:]))), 'g', -1, 32)
		}
		return strings.Join(parts, ", ")
	case ncTypeDouble:
		b := r.take(count * 8)
		if b == nil {
			return ""
		}
		parts := make([]string, count)
		for i := range count {
			parts[i] = strconv.FormatFloat(math.Float64frombits(binary.BigEndian.Uint64(b[i*8:])), 'g', -1, 64)
		}
		return strings.Join(parts, ", ")
	default:
		r.fail("unknown attribute type %d", ncType)
		return ""
	}
}

func (r *ncReader) varList(wide bool, dims []NetCDFDim) []ncVar {
	if r.err != nil {
		return nil
	}
	if r.u32() == 0 {
		r.u32() // ABSENT
		return nil
	}
	n := r.boundedCount(r.u32(), 28, "variables")
	vars := make([]ncVar, 0, n)
	seenVar := map[string]bool{}
	for range n {
		v := ncVar{name: r.name()}
		if seenVar[v.name] {
			r.fail("variable %q is declared twice", v.name)
			return nil
		}
		seenVar[v.name] = true
		rank := int(r.u32())
		if rank < 0 || int64(rank) > int64(len(r.data)-r.pos)/4 {
			r.fail("variable %q declares %d dimensions, more than the file can hold", v.name, rank)
			return nil
		}
		v.dimIDs = make([]int, rank)
		for j := range rank {
			id := int(r.u32())
			if id < 0 || id >= len(dims) {
				r.fail("variable %q refers to dimension id %d of %d", v.name, id, len(dims))
				return nil
			}
			v.dimIDs[j] = id
		}
		v.attrs = r.attList()
		v.ncType = r.u32()
		r.u32() // vsize: redundant, the dimensions carry the count
		if wide {
			v.begin = int64(r.u64())
		} else {
			v.begin = int64(r.u32())
		}
		if r.err != nil {
			return nil
		}
		vars = append(vars, v)
	}
	return vars
}

// isRecordVar reports whether the variable's leading dimension is the
// record dimension, which is what makes its slabs interleave.
func isRecordVar(dims []NetCDFDim, dimIDs []int) bool {
	return len(dimIDs) > 0 && dims[dimIDs[0]].Length == 0
}

// ncPad4Size rounds a byte count up to a four-byte boundary, the
// padding the classic format gives every variable slab.
func ncPad4Size(n int64) int64 { return (n + 3) &^ 3 }

// ncTypeWidth returns the stored size of a classic type code. CHAR
// carries one byte per element, which the array level lands as uint8
// raw bytes rather than as text.
func ncTypeWidth(ncType uint32) (int64, bool) {
	switch ncType {
	case ncTypeByte, ncTypeChar:
		return 1, true
	case ncTypeShort:
		return 2, true
	case ncTypeInt, ncTypeFloat:
		return 4, true
	case ncTypeDouble:
		return 8, true
	}
	return 0, false
}

// ncLandingDtype maps a classic type code onto the core dtype its
// values land in. NC_BYTE is signed in the classic model; NC_CHAR
// lands uint8 raw bytes, which say nothing about text at the array
// level. NC_FLOAT and NC_DOUBLE keep their float64 landing, where
// every classic value is exact. Codes 7 to 11 exist only in the
// unsigned and 64-bit extensions of the classic formats this module
// speaks, and stay refused by name: no supported writer emits them.
func ncLandingDtype(ncType uint32) (core.Dtype, bool) {
	switch ncType {
	case ncTypeByte:
		return core.Int8, true
	case ncTypeChar:
		return core.Uint8, true
	case ncTypeShort:
		return core.Int16, true
	case ncTypeInt:
		return core.Int32, true
	case ncTypeFloat, ncTypeDouble:
		return core.Float, true
	}
	return 0, false
}

// valuesRecord reads a record variable into arr's payload: each of the
// numrecs records carries its own slab of perRecord elements at the
// variable's begin offset plus a whole multiple of recordSize, because
// the other record variables' slabs sit between them. The padding
// inside a slab is skipped, never decoded.
func (r *ncReader) valuesRecord(arr *core.Array, begin int64, ncType uint32, perRecord, numrecs, recordSize int64) {
	if r.err != nil {
		return
	}
	width, ok := ncTypeWidth(ncType)
	if !ok {
		r.fail("variable has unknown type %d", ncType)
		return
	}
	if perRecord < 0 || numrecs < 0 || recordSize < perRecord*width {
		r.fail("record variable declares %d elements per record over %d records", perRecord, numrecs)
		return
	}
	// Every slab must lie inside the file: with records present the last
	// record's slab ends at begin + (numrecs-1)*recordSize +
	// perRecord*width, and with none there is nothing to read at all.
	// The record count is divided out before it is multiplied: a hostile
	// record size makes the product wrap, and a wrapped span passes a
	// check that compares the product itself.
	failSpan := func() {
		r.fail("variable data at offset %d (%d records of %d bytes) runs past the end of the %d-byte file",
			begin, numrecs, recordSize, len(r.data))
	}
	span := int64(0)
	if numrecs > 0 {
		span = perRecord * width
		fits := span <= int64(len(r.data))
		if numrecs > 1 {
			fits = fits && recordSize <= (int64(len(r.data))-span)/(numrecs-1)
		}
		if !fits {
			failSpan()
			return
		}
		span += (numrecs - 1) * recordSize
	}
	if begin < 0 || begin > int64(len(r.data))-span {
		failSpan()
		return
	}
	for rec := range numrecs {
		slab := begin + rec*recordSize
		// The window, not the tail: decodeCells fills [lo, hi) cells.
		lo := int(rec * perRecord)
		if err := decodeCells(arr, r.data[slab:slab+perRecord*width], ncType, lo, lo+int(perRecord)); err != nil {
			r.fail("%v", err)
			return
		}
	}
}

// values reads count elements of the given type from the absolute file
// offset into arr's payload, in the dtype the type code lands: every
// classic integer type fits its own narrow dtype exactly, and the
// float codes keep their float64 landing.
func (r *ncReader) values(arr *core.Array, begin int64, ncType uint32, count int64) {
	if r.err != nil {
		return
	}
	width, ok := ncTypeWidth(ncType)
	if !ok {
		r.fail("variable has unknown type %d", ncType)
		return
	}
	// The declared count is a product of header fields, so it can be
	// negative (a wrapped multiply) or larger than the address space.
	// The comparison is done in int64 against the bytes actually left
	// in the file, which rejects both without ever forming a size the
	// allocator would have to refuse.
	span := count * width
	if count < 0 || span < 0 || begin < 0 || begin > int64(len(r.data))-span {
		r.fail("variable data at offset %d (%d %d-byte elements) runs past the end of the %d-byte file",
			begin, count, width, len(r.data))
		return
	}
	if err := decodeCells(arr, r.data[begin:begin+span], ncType, 0, int(count)); err != nil {
		r.fail("%v", err)
	}
}

// decodeCells decodes one contiguous run of stored cells into the
// payload window [lo, hi) of arr, which the caller allocated in the
// dtype ncLandingDtype maps ncType onto; the raw run and the window
// must agree in element count for the type. A type code without a case
// is refused: the classic formats carry only the codes
// ncLandingDtype maps, and an unhandled one left the payload silently
// zero once, which is the worst answer a reader can give.
func decodeCells(arr *core.Array, raw []byte, ncType uint32, lo, hi int) error {
	switch ncType {
	case ncTypeByte:
		dst := arr.RawInt8s()[lo:hi]
		for i := range dst {
			dst[i] = int8(raw[i])
		}
	case ncTypeChar:
		copy(arr.RawUint8s()[lo:hi], raw[:hi-lo])
	case ncTypeShort:
		dst := arr.RawInt16s()[lo:hi]
		for i := range dst {
			dst[i] = int16(binary.BigEndian.Uint16(raw[i*2:]))
		}
	case ncTypeInt:
		dst := arr.RawInt32s()[lo:hi]
		for i := range dst {
			dst[i] = int32(binary.BigEndian.Uint32(raw[i*4:]))
		}
	case ncTypeFloat:
		dst := arr.RawFloats()[lo:hi]
		for i := range dst {
			dst[i] = float64(math.Float32frombits(binary.BigEndian.Uint32(raw[i*4:])))
		}
	case ncTypeDouble:
		dst := arr.RawFloats()[lo:hi]
		for i := range dst {
			dst[i] = math.Float64frombits(binary.BigEndian.Uint64(raw[i*8:]))
		}
	default:
		return base.Errf("variable has unknown type %d", ncType)
	}
	return nil
}

// ncExternal maps a dtype onto the NetCDF type code and the external
// size the file stores it in.
func ncExternal(dt core.Dtype) (code uint32, size int, err error) {
	switch dt {
	case core.Float:
		return ncTypeDouble, 8, nil
	case core.Float32:
		return ncTypeFloat, 4, nil
	case core.Int:
		return ncTypeInt, 4, nil
	default:
		return 0, 0, base.Errf("the classic model stores float64, float32 and int arrays, got dtype %s", dt)
	}
}

func ncExternalType(dt core.Dtype) uint32 {
	code, _, _ := ncExternal(dt)
	return code
}

// ncCheckName enforces the traditional NetCDF name grammar: the first
// character alphanumeric or '_', the rest alphanumeric or one of
// '_.@+-', which every classic reader accepts without escaping.
func ncCheckName(s string) error {
	if s == "" {
		return base.Errf("names must not be empty")
	}
	for i := range s {
		c := s[i]
		ok := c == '_' || c >= '0' && c <= '9' || c >= 'A' && c <= 'Z' || c >= 'a' && c <= 'z'
		if i > 0 {
			ok = ok || c == '.' || c == '@' || c == '+' || c == '-'
		}
		if !ok {
			return base.Errf("name %q contains the character %q, which the traditional name grammar forbids", s, c)
		}
	}
	return nil
}

// ncAppendName writes a length-prefixed name padded to 4 bytes.
func ncAppendName(buf []byte, s string) []byte {
	buf = binary.BigEndian.AppendUint32(buf, uint32(len(s)))
	buf = append(buf, s...)
	for range ncPadLen(len(s)) {
		buf = append(buf, 0)
	}
	return buf
}

// ncPadLen is the number of bytes needed to reach the next 4-byte
// boundary.
func ncPadLen(n int) int { return (4 - n%4) % 4 }

// ncPad4 appends null bytes up to the next 4-byte boundary. The
// header pads with zeros; the payload types this module writes are all
// multiples of four bytes wide, so no fill value is ever needed.
func ncPad4(buf []byte) []byte {
	for range ncPadLen(len(buf)) {
		buf = append(buf, 0)
	}
	return buf
}