Spherical sampling

AbstractLatLonSampling <: AbstractTensorProductSphericalSampling

Geophysical / model latitude–longitude grids. No exact band-limited SHT claim.

AbstractReducedGaussianSampling <: AbstractRingSampling

Ring samplings on Gaussian latitudes. Default: OctahedralGaussianSampling.

AbstractRingSampling <: AbstractSphericalSampling

Iso-latitude rings whose longitude count VARIES by ring, so the layout is not a tensor product. A reduced Gaussian grid is the canonical case: latitudes are the Gaussian ones, but each ring carries only as many longitudes as its circumference warrants.

AbstractSpectralQuadratureSampling <: AbstractTensorProductSphericalSampling

Samplings with a known exact quadrature for band-limited spherical harmonics (up to a documented lmax ↔ grid-size relation).

AbstractSphericalSampling

How points are placed on the sphere. Orthogonal to FlowGeometries.Geometry.AbstractSphericalGeometry (the metric / radius) and to grid architecture (StructuredGrid vs unstructured).

AbstractTensorProductSphericalSampling <: AbstractSphericalSampling

Iso-latitude × equispaced-longitude (or arbitrary lon/lat) layouts that fit a 2D StructuredGrid with axes (λ, φ).

ClenshawCurtisSampling <: AbstractClenshawCurtisSampling

Open equiangular Clenshaw–Curtis grid: $θ_i = π(i−1/2)/N_θ$ (no poles), $N_λ = 2N_θ − 1$, $l_{\max} = N_θ − 1$.

This is not the same as Driscoll–Healy: different $N_θ(l_{\max})$, open vs polar-cap nodes, and a different quadrature.

CubedSphereSampling <: AbstractCubedSphereSampling

Inscribed-cube projection: six logically rectangular panels (gnomonic by default). Quasi-uniform alternative to global lat–lon (e.g. FV3).

DriscollHealyEqualSampling <: AbstractDriscollHealySampling

Driscoll–Healy DH1 layout: $N_θ × N_θ$ with $N_θ = 2(l_{\max}+1)$, $Δθ = π/N_θ$, $Δλ = 2π/N_θ$.

DriscollHealySampling <: AbstractDriscollHealySampling

Driscoll–Healy equiangular grid, rectangular (DH2) layout: $N_θ × 2N_θ$ with $N_θ = 2(l_{\max}+1)$, $Δθ = π/N_θ$, $Δλ = π/N_θ$. Includes the north pole band, excludes the south (weight at the north pole is zero).

See also DriscollHealyEqualSampling for the square DH1 $N_θ × N_θ$ layout.

FibonacciSampling(n)

n points on the spherical Fibonacci (golden-spiral) lattice: z_k = (2k+1)/n - 1 with λ_k = 2πk/φ mod 2π, φ the golden ratio.

z advances in equal steps, so the points are spread one per equal-area band, and the golden angle in longitude — the most irrational rotation there is — keeps them from lining up into visible spokes at any n. Quasi-uniform with no polar clustering and no panel seams.

GaussLegendreSampling <: AbstractGaussLegendreSampling

Gauss–Legendre (Gauss–Neumann) latitudes: $μ = \cosθ$ at the $N_θ$ roots of $P_{N_θ}$, with $N_λ = 2N_θ − 1$ equispaced longitudes. Exact for band-limit $l_{\max} = N_θ − 1$.

HEALPixSampling <: AbstractHEALPixSampling

Hierarchical Equal Area isoLatitude Pixelisation (Górski et al.). Equal-area pixels on $4 N_{\mathrm{side}} − 1$ iso-latitude rings; not a tensor-product $N_λ × N_φ$ array ($n_λ$ varies by ring).

IcosahedralSampling <: AbstractIcosahedralSampling

Icosahedral geodesic discretization (triangular faces; hexagonal dual is the Voronoi mesh). Quasi-uniform; unstructured connectivity.

LatLonSampling <: AbstractLatLonSampling

Arbitrary model lat–lon (regional or global, poles optional). No spectral exactness claim.

McEwenWiauxSampling <: AbstractMcEwenWiauxSampling

McEwen–Wiaux (2011) equiangular sampling: $θ_t = π(2t+1)/(2L−1)$ for $t = 0,…,L−1$, $λ_p = 2π p/(2L−1)$ for $p = 0,…,2L−2$, with $L = l_{\max}+1$. Requires asymptotically $∼ 2 L^2$ samples (about half of classical DH).

Nested()

Pixels numbered so that each is subdivided into four contiguous children — a quadtree per base face. The ordering that makes a neighbourhood contiguous. Requires nside to be a power of two.

OctahedralGaussianSampling(N)

Octahedral reduced Gaussian grid with N latitude rings between pole and equator (2N rings in all).

Ring i counted from either pole carries 4i + 16 longitudes — 20 on the ring nearest the pole and four more on each successive ring — which totals 4N(N+9) points. Latitudes are the 2N Gaussian latitudes, so the latitude quadrature is the Gauss–Legendre one.

OpenNodes(), ClosedNodes()

The two equiangular colatitude families: open θᵢ = π(i−½)/N (Clenshaw–Curtis) and closed θᵢ = π(i−1)/N (Driscoll–Healy). Named rather than passed as a θ(i) closure so the sum below can dispatch on which one it is.

ReducedGaussianSampling(nlon_per_ring)

Reduced Gaussian grid with an explicit longitude count per ring, north to south.

The classical reduced grids are published as tables rather than formulas — the count holds constant across blocks of latitudes and jumps between them — so the table is the input. Use OctahedralGaussianSampling for the octahedral rule, which is a formula.

Ring()

Pixels numbered along iso-latitude rings, north to south and east within a ring. The ordering that makes a ring contiguous, so a longitude transform per ring is possible.

RingScheme

Which HEALPix pixel ordering is meant: Ring or Nested.

ScatteredSphericalSampling <: AbstractScatteredSphericalSampling

Arbitrary $(λ, φ)$ point set (NUFFT / NUFSHT paths).

YinYangSampling <: AbstractYinYangSampling

Overset of two low-latitude lat–lon patches (Kageyama–Sato), avoiding polar singularities while keeping structured coordinates on each panel.

_equiangular_sums!(s, family, nlat, nterm) -> s

s[i] = Σ_{k=0}^{nterm-1} sin((2k+1)·θᵢ)/(2k+1) for the given node family.

Evaluated by the angle-addition recurrence sin((2k+1)θ) = 2cos(2θ)·sin((2k−1)θ) − sin((2k−3)θ), seeded with s₋₁ = −sin θ, s₀ = sin θ, so each term costs two multiplies rather than a sin. This is O(nlat·nterm); loading an FFT implementation replaces it with one length-nlat transform (see the AbstractFFTs extension).

_equiangular_weights!(w, family, nlat) -> w

Latitude quadrature weights for an equiangular node set, from the sine-series expansion of the sinθ Jacobian:

wᵢ = (4/N)·sinθᵢ·Σ_{k=0}^{⌊N/2⌋} sin((2k+1)θᵢ)/(2k+1)

family is OpenNodes or ClosedNodes, so one construction serves both. Weights sum to ∫₀^π sinθ dθ = 2 — see latitude_weights for why every sampling here uses that normalization.

_gauss_legendre!(T, μ, w, n)

Core solve. Either output may be nothing, so a caller that wants only the nodes or only the weights needs no scratch vector for the other half.

_gauss_legendre_newton!(T, TW, μ, w, n, m)

Nodes and weights by Newton on Pₙ, evaluated by the Bonnet recurrence from a Tricomi start.

O(n) per root, so O(n²) overall — but it is exact arithmetic converging to eps(TW), which the asymptotic expansion's fixed Float64 coefficient set cannot do. That makes it the right method for wide element types at any n, and for small n, where the expansion is inaccurate and the quadratic cost is microseconds.

_gauss_legendre_μ!(μ, w) -> NamedTuple{(:μ,:w)}

Write n = length(μ)-point Gauss–Legendre nodes/weights on $μ ∈ (-1, 1)$ into the provided buffers, ascending in μ (length(μ) == length(w)).

Newton's method on $Pₙ$, evaluated by the Bonnet recurrence, from a Tricomi starting estimate accurate to $O(n⁻³)$ — so 2–4 iterations reach machine precision. Roots come in $±$ pairs, so only the upper half is solved.

Needs $O(1)$ scratch. Time is still $O(n²)$ — each of the $n/2$ roots costs an $O(n)$ recurrence; an $O(n)$ total needs Bogaert-style asymptotic expansions instead of Newton.

_hp_ang2xyf(nside, θ, ϕ) -> (ix, iy, face)

Face-local coordinates of the pixel containing colatitude θ, longitude ϕ, by the HEALPix projection (Górski et al. 2005). Written with division and remainder rather than shift and mask, so it holds for any nside rather than only a power of two.

_ico_node_id(f, face, i, j, ν, edge_index, nint, face_base) -> Int

Global vertex id of barycentric lattice node (i, j) (with i + j ≤ ν, weights ν-i-j, i, j on the face's corners A, B, C) of face f. Nodes on a corner or a macro-edge resolve to that shared entity's id, which is what makes the two faces meeting at an edge agree without any lookup table.

_icosahedron_base(T) -> NTuple{12,NTuple{3,T}}

The 12 unit vertices of the base icosahedron, in the order _ICOSAHEDRON_FACES indexes them. A tuple, so it costs no allocation. Every raw vertex is a permutation of (0, ±1, ±φ) and so shares the norm √(1+φ²), formed once.

_put_lonlat!(λ, φ, v, p, T)

Normalize p and write vertex v's longitude/latitude. Top-level rather than a closure so nothing is captured.

admits_exact_bandlimited_quadrature(sampling) -> Bool

Whether this sampling's latitude_weights integrate the PRODUCTS that spectral analysis actually forms — two degree-lmax functions, hence degree 2·lmax — exactly at the sampling's own bandlimit.

That is a stronger statement than "the quadrature integrates a single P_l up to lmax", and the distinction decides the answer here. Measured exactness of the weights in this package:

Clenshaw–Curtis's band limit describes what its grid can REPRESENT; its quadrature supports quadrature-based analysis only to lmax ≈ (N−1)/2. Use GaussLegendreSampling when analysis must be exact at the stated band limit.

ang2pix(nside, θ, ϕ; scheme = Ring()) -> Int

The 0-based index of the pixel containing colatitude θ ∈ [0, π] and longitude ϕ.

θ is a COLATITUDE, matching the HEALPix convention throughout this section; use colatitude to convert a geographic latitude.

axes_lengths(sampling, nlat; nlon=nothing) -> NamedTuple{(:nlon,:nlat)}

Buffer lengths required by spherical_axes!: (; nlon, nlat).

bandlimit(sampling, nlat) -> lmax

Maximum spherical-harmonic degree supported by nlat latitude samples under sampling.

Colatitude $θ ∈ [0, π]$ from geographic latitude $φ ∈ [-π/2, π/2]$.

cubed_sphere_points!(λ, φ, panel, n; backend=nothing) -> NamedTuple{(:λ,:φ,:panel)}

Gnomonic cubed-sphere CELL CENTRES into caller-owned buffers of length 6n², plus each point's panel index. Pass panel = nothing when the panel id is not wanted, and it is not computed.

See cubed_sphere_points for the allocating form.

cubed_sphere_points([T = Float64], n; backend=nothing) -> NamedTuple{(:λ,:φ,:panel)}

Gnomonic cubed-sphere CELL CENTRES: 6n² distinct points, plus each point's panel index.

Allocating wrapper around cubed_sphere_points!. Use spherical_points(CubedSphereSampling(), n) when the panel id is not needed.

Geographic latitude $φ$ from colatitude $θ$.

icosahedral_mesh([T = Float64], frequency) -> (; λ, φ, edges, triangles, verts)

Geodesic vertices at frequency ν as both lon/lat (λ, φ) and unit vectors (verts), plus the 10ν²+2 mesh's undirected edges (i,j) with i < j and its 20ν² triangles (i,j,k) (1-based). Vertex numbering is topological — corners, then macro-edge interiors, then face interiors — so it is deterministic and every vertex is generated exactly once.

icosahedral_vertices([T = Float64], frequency=1) -> NamedTuple{(:λ,:φ)}

The 10ν²+2 geodesic vertices as longitude/latitude, without building the mesh topology — several times faster than icosahedral_mesh at large ν, which also returns edges and triangles.

icosahedral_vertices!(λ, φ, frequency = 1) -> NamedTuple{(:λ,:φ)}

Write the 10ν²+2 geodesic vertices' longitude/latitude into the caller's buffers, allocating nothing. The numbering is icosahedral_mesh's topological one — corners, then macro-edge interiors, then face interiors — so each vertex is emitted at its own index with no intermediate array and no lookup.

latitude_weights!(w, sampling, nlat) -> w

Fill preallocated latitude quadrature weights (length(w) == nlat).

latitude_weights!(w, ::AbstractClenshawCurtisSampling, nlat)

Weights for the open nodes θᵢ = π(i−½)/N, from the same sine-series rule as the closed equiangular families.

These integrate a single P_l exactly for l ≤ N−1, which is weaker than bandlimit(ClenshawCurtisSampling(), N) = N−1 suggests: spectral analysis integrates PRODUCTS of two degree-lmax functions, so a quadrature exact to l ≤ N−1 supports quadrature-based analysis only up to lmax ≈ (N−1)/2. The reported band limit describes what the grid can represent, not what this quadrature can integrate. Use GaussLegendreSampling (exact to 2N−1) when analysis must be exact at the stated band limit.

latitude_weights([T = Float64], sampling::AbstractReducedGaussianSampling) -> Vector{T}

Gauss–Legendre weights for the grid's rings, north to south, normalized as everywhere else in this module so that Σw = 2. A full-sphere integral is then Σⱼ wⱼ (2π/nlonⱼ) Σᵢ f, the longitude factor varying by ring because the ring populations do.

latitude_weights([T = Float64], s, nlat) -> Vector{T}
latitude_weights!(w, s, nlat) -> w

Latitude quadrature weights wⱼ for sampling s, normalized so that

Σⱼ wⱼ = ∫₀^π sinθ dθ = 2

for EVERY sampling that provides them. The weights therefore carry the sinθ Jacobian and nothing else; the longitude factor is the caller's, so a full-sphere integral is always

∫ f dΩ ≈ (2π/nlon) · Σⱼ wⱼ Σᵢ f(λᵢ, φⱼ)

regardless of which sampling produced the weights. Not every sampling has them — LatLonSampling has no spectral quadrature at all, and McEwenWiauxSampling's is a different construction.

nlat_for_bandlimit(sampling, lmax) -> nlat

nlon_for_nlat(sampling, nlat) -> nlon

nlon_in_ring(sampling, ring) -> Int
nlon_in_ring(sampling, nlat, ring) -> Int

Longitudes on ONE iso-latitude ring, counted from the north pole, in O(1) and allocating nothing.

The per-ring form of nlon_per_ring, and the one a ring-by-ring loop wants: the table costs an allocation and an O(nrings) build per call, which a loop over rings pays again on every iteration if it asks for it there.

nlon_per_ring(sampling) -> Vector{Int}
nlon_per_ring(sampling, nlat; nlon=nothing) -> Vector{Int}

Longitudes on each iso-latitude ring, north to south.

Defined for every sampling laid out in rings, so a caller walking a map ring by ring — a per-ring longitude transform, a zonal reduction, a row-wise sweep — writes one loop rather than a branch per sampling. A sampling that carries its own size (HEALPixSampling, the reduced Gaussians) answers from itself; a tensor-product one takes the nlat that fixes its shape, as npoints does.

npoints(::AbstractScatteredSphericalSampling, λ, φ) -> Int

The number of points in a scattered set, i.e. length(λ).

npoints(sampling, args...) -> Int

Number of geographic samples produced by spherical_points! / spherical_points.

nrings(sampling) -> Int
nrings(sampling, nlat) -> Int

Number of iso-latitude rings, for any sampling laid out in them — the loop bound that goes with nlon_per_ring.

pix2ang([T = Float64], nside, pix; scheme = Ring()) -> (θ, ϕ)

Colatitude and longitude of pixel pix's centre (0-based index).

pix2vec([T = Float64], nside, pix; scheme = Ring()) -> NTuple{3}

Unit vector to pixel pix's centre.

ring2nest(nside, pix) -> Int
nest2ring(nside, pix) -> Int

Convert a 0-based pixel index between the two orderings. Both need nside to be a power of two, which is what makes the NESTED quadtree exist.

ring_info([T = Float64], nside, ring) -> NamedTuple

What HEALPix ring ring ∈ 1:(4·nside-1) contains, counted from the north pole: startpix (the 0-based RING index of its first pixel, matching ang2pix), ringpix (how many pixels it holds), colatitude, latitude, and shifted — whether its pixel centres are offset half a pixel in ϕ.

Ring width grows 4, 8, … through the polar cap, is 4·nside across the equatorial belt, and shrinks again symmetrically, so this is how to walk a HEALPix map ring by ring without decoding every pixel.

ring_latitudes([T = Float64], sampling) -> Vector{T}

The Gaussian latitudes of a reduced Gaussian grid's rings, north to south.

ring_range(sampling, ring) -> UnitRange{Int}
ring_range(sampling, nlat, ring) -> UnitRange{Int}

The 1-based slice of the flattened point vector holding one iso-latitude ring, north to south — the indices spherical_points writes that ring into.

O(1) for every sampling, so a ring-by-ring pass is a loop over slices carrying no running offset. The octahedral rule's offset is a closed form in the ring index; the tabulated reduced grid carries cumulative counts built with the sampling rather than rescanned per lookup.

spherical_axes([T = Float64], sampling, nlat; nlon=…) -> NamedTuple{(:λ,:φ)}

Allocating wrapper around spherical_axes!.

The element type leads, as it does for zeros and rand, so that it takes part in dispatch: a type supplied as a keyword is invisible to dispatch, and the returned eltype would be unknown to the caller for any call the compiler cannot constant-fold whole.

spherical_axes!(λ, φ, sampling, nlat; nlon=…) -> NamedTuple{(:λ,:φ)}

Fill preallocated longitude / geographic-latitude axes. Requires length(λ) == sz.nlon and length(φ) == sz.nlat for sz = axes_lengths(...).

spherical_points!(λ, φ, sampling, args...) -> NamedTuple{(:λ,:φ)}

Fill preallocated point-coordinate buffers. See npoints for required lengths.

spherical_points!(λ_out, φ_out, ::AbstractScatteredSphericalSampling, λ, φ) -> NamedTuple

Copy a scattered point set into caller-owned buffers.

The source arrays are required: a scattered sampling carries no rule from which points could be generated, so a form taking only the destinations would have nothing to write.

spherical_points!(λ, φ, sampling::AbstractReducedGaussianSampling; scratch = nothing) -> NamedTuple

Ring-by-ring points of a reduced Gaussian grid, north to south, longitudes equispaced within each ring starting at zero. Buffer length is npoints.

The Gaussian latitudes need one value per ring, which cannot overlap the output here (see the note in the body), so scratch — any vector of at least nrings(sampling) elements — makes the fill allocation-free for a caller filling many grids.

spherical_points!(λ, φ, sampling::FibonacciSampling) -> NamedTuple

The n golden-spiral points, allocating nothing.

spherical_points(::AbstractScatteredSphericalSampling, λ, φ) -> NamedTuple{(:λ,:φ)}

A scattered sampling's points are the caller's arrays, so this hands them back. It exists so a scattered set can be driven through the same entry point as a generated one.

spherical_points([T = Float64], sampling, args...) -> NamedTuple{(:λ,:φ)}

Every point of the sampling, flattened, as longitude/latitude vectors. Allocating wrapper around spherical_points!; use npoints to size buffers for the in-place form.

For a tensor-product sampling this is the outer product of its axes, so prefer spherical_axes when the separable form will do.

spherical_quadrature!(λ, φ, w, sampling, nlat; nlon=…) -> NamedTuple{(:λ,:φ,:w)}

Fill preallocated axes and latitude weights together. Buffer lengths are as for spherical_axes!, with length(w) == sz.nlat.

This is the entry point to use whenever both are needed. Gauss–Legendre nodes and weights fall out of a single root solve, but spherical_axes! keeps only the nodes and latitude_weights! only the weights — so calling them in sequence pays the $O(n²)$ solve twice. Every other sampling has independent closed forms for the two, and gets the generic method.

spherical_quadrature([T = Float64], sampling, nlat; nlon=…) -> NamedTuple{(:λ,:φ,:w)}

Allocating wrapper around spherical_quadrature!.

vec2pix(nside, v; scheme = Ring()) -> Int

The 0-based index of the pixel containing direction v, which need not be normalized.

yin_yang_panels!(λyin, φyin, λyang, φyang, nlon, nlat) -> (; yin, yang)

The two Kageyama–Sato panels. yin is a pair of AXES (nlon and nlat long): in its own frame the panel is a separable lat–lon patch. yang is that panel rotated onto the sphere, which is no longer separable in global lon/lat, so it is a pair of nlon × nlat FIELDS — one (λ, φ) per cell.

The shapes differ because the geometry does, not by convention; the argument types say so.

yin_yang_panels([T = Float64], nlon, nlat) -> (; yin, yang)

Allocating wrapper around yin_yang_panels!.