Neighbours and stencils

Use this API to find cells around a cell, build an adjacency table, or compute new values from a neighbourhood. For the edge of an entire region, use Region boundaries. The stencil tutorial works through smoothing, edge detection and graph traversal.

Find neighbours around a cell

neighbors(grid, cell, k) returns cells within k adjacency steps, excluding the centre. ring(grid, cell, k) returns only the cells at step k. Each ring runs counter-clockwise as seen from outside the sphere.

import DiscreteGlobalGrids as DGG

grid = DGG.levelgrid(DGG.HEALPixSystem(), 3)
cell = DGG.cellat(grid, 8.5, 47.4)

(; first_ring = length(DGG.neighbors(grid, cell)),
   second_ring = length(DGG.ring(grid, cell, 2)),
   within_two = length(DGG.neighbors(grid, cell, 2)))

(first_ring = 8, second_ring = 15, within_two = 23)

Pass a cell id to receive cell ids, or a local integer index to receive local indices. On a subset, results include only members of that subset.

DiscreteGlobalGrids.neighbors Function

neighbors(grid::AbstractGrid, c::AbstractCellIndex, k::Integer = 1; connectivity::Connectivity = Vertex())
neighbors(grid::AbstractGrid, p::Int, k::Integer = 1; connectivity::Connectivity = Vertex()) -> AbstractVector{Int}

All cells of grid within k adjacency steps of c, excluding c.

Connectivity

Vertex(), which includes vertex contact, is the default. Edge() requires a shared edge. They coincide where exactly three cells meet at every vertex — the icosahedral hexagon-with-pentagon family (IGeo7, H3) — and differ wherever a vertex carries more, including one pentagonal system: A5's Cairo-style tiling gives 11 vertex against 3 edge neighbours at resolution 1. Quadrilateral grids add corner neighbours under Vertex(): four cells to a vertex in the lattice interior, five where ISEA4R's diamonds meet an icosahedral vertex.

Order

Order is part of the contract, and it is ONE order: every verb in this package that hands back a neighbourhood hands it back counter-clockwise. Results concatenate rings outward:

neighbors(grid, c, k) == vcat(ring(grid, c, 1), ring(grid, c, 2), ..., ring(grid, c, k))

Counter-clockwise, exactly. Take n = cell_centroid(grid, c) as the outward normal and project each ring member's centroid into the tangent plane at n. Read in order, the azimuths of those projections in a right-handed frame (e₁, e₂ = n × e₁) increase and wrap through 2π exactly once. That is the same rotational sense cell_boundary(grid, c) winds in: the package has one handedness, fixed by the boundary contract, and every ring agrees with it. A ring read backwards wraps length - 1 times; an id-sorted one wraps some arbitrary number of times. Neither is this order.

Start. The direction is guaranteed everywhere; the phase is guaranteed only within a system. Ring 1 begins at the system's own documented direction — +s for S2, SW for HEALPix, NW for CopernicusDEM, the development frame's +1 for IGeo7, the smallest-id neighbour for A5 and for the geometric fallback — and rings 2:k of the same call begin on the same spoke as ring 1, the azimuth of ring(grid, c, 1)[1]. So a disc reads as concentric rings all starting in one direction, but which cell that is differs by system and is not a portable fact. Exact azimuth ties break by canonical id. The start is deterministic — a property of the system and the cell alone, the same first member every time the cell is asked, in every idiom, independent of any region or table the ring is read through.

Cells with fewer neighbours yield shorter rings without padding, and an omitted neighbour leaves no gap in the sequence.

The idioms that carry this order

ring, the index forms below, adjacency, member_neighbors, the one-argument neighbors iterator, and the rings mapneighbors and foreachneighbors pass to their callbacks — all of them. Nothing in this package answers a neighbourhood question in ascending id or index.

The verb that is ascending is halo, and it is not a ring: it is a fetch list, ordered so that a read is sequential. It says so where it is documented.

Container

Any ordered, indexable collection with the grid's cell-index eltype. Index forms preserve the id form's container family while replacing its element type with Int; in particular, a fixed-capacity one-ring remains a fixed-capacity one-ring after conversion to indices.

Coverage, and what a subset means

Neighbours outside the grid's coverage are omitted, not padded.

On a subset of a complete level — PartialGrid, CellVector, CellLookup — adjacency and distance are the complete level's, clipped to membership:

ring(sub, c, k) == filter(in(sub), ring(levelgrid(system(sub), level(sub)), c, k))

Distance is therefore measured in the system, never inside the subset. A hole in the subset removes cells; it does not lengthen the path around itself, and a cell reachable only by leaving the subset and coming back keeps the distance the complete level gives it. The two readings coincide at k == 1 and part company from k == 2, which is why this is stated rather than left to the reader.

The filter also pins the rotation: a clipped ring is the complete ring, read from its canonical start, with non-members dropped in place. Its length is the in-set degree, and which absolute slot a surviving member occupied is deliberately not recoverable from the clipped ring — a consumer that needs slot identity reads the complete level's ring.

c outside the subset is an ArgumentError, not a complete-grid answer.

k must be ≥ 0; k == 0 returns an empty collection. See ring for the cells at exactly distance k.

Indices

Given a local index, both verbs answer with in-set local indices in the rotational order above: neighbors(grid, p, k) is neighbors(grid, cellindex(grid, p), k) mapped through localindex, element for element, with non-members dropped. adjacency is this form for a whole region at once.

The result is therefore not sorted. An index list read only by membership does not care, and one read by direction — a gradient, an upwind stencil — cannot be written at all against a sorted list, which is why the index forms carry the same order the id forms do.

source

DiscreteGlobalGrids.ring Function

ring(grid::AbstractGrid, c::AbstractCellIndex, k::Integer; connectivity::Connectivity = Vertex())
ring(grid::AbstractGrid, p::Int, k::Integer; connectivity::Connectivity = Vertex()) -> AbstractVector{Int}

The cells at adjacency distance exactly k from c. ring(grid, c, 0) is c alone. The ordered result satisfies

neighbors(grid, c, k) ==
    reduce(vcat, [ring(grid, c, j) for j in 1:k]; init = eltype(grid)[])

so ring k is the final ordered block of neighbors(grid, c, k). Overrides must preserve this equality.

init is load-bearing, and the splatted vcat(ring(grid, c, 1), ...) is not an equivalent spelling. Rings at different k may arrive in containers of different capacity — a SmallVector{6} at k == 1 beside a heap Vector at k == 2 — and vcat takes similar from its first argument, so concatenating onto the one-ring overflows that one-ring's capacity and throws. Seeding an empty Vector fixes the result type; reduce also keeps the call out of a splat, which is what stops the arity from specialising per k.

ring carries neighbors' order, container, coverage and subset-clipping contracts unchanged — including the local-index form's, which is the same counter-clockwise order read through localindex.

k as a type

Both verbs also accept Val(k), which answers the same cells in the same order and differs only in what the compiler is told. With k in the type, a system's declared maxring folds into a fixed buffer capacity, so the shell is built and returned without reaching the heap:

ring(grid, c, 2)       # Vector, capacity found at run time
ring(grid, c, Val(2))  # SmallVector, capacity found at compile time

Worth reaching for in a focal loop over many cells and not otherwise. The Integer form is never wrong and never slower than it was; Val(k) only removes allocations, and only for a system that has opted in — the rest forward to the Integer form and lose nothing.

source

Choose connectivity

The default Vertex() includes cells touching at a corner. Edge() requires a shared edge. These choices often coincide on hexagonal grids and differ on quadrilateral grids.

DiscreteGlobalGrids.Connectivity Type

abstract type Connectivity

How two cells must meet to count as adjacent. See neighbors.

Concrete singletons: Vertex (Moore, the default) and Edge (von Neumann).

source

DiscreteGlobalGrids.Vertex Type

Vertex() <: Connectivity

Moore connectivity: cells are adjacent when they share at least a vertex. This is the default connectivity.

Vertex() and Edge() coincide where exactly three cells meet at each vertex. Higher-valence vertices add corner-only neighbours. This includes A5's 4-valent corners and ISEA4R's 4-valent lattice corners and 5-valent icosahedral vertices.

source

DiscreteGlobalGrids.Edge Type

Edge() <: Connectivity

von Neumann connectivity: two cells are adjacent only if they share a whole edge. The opt-in restriction of the Vertex() default.

source

DiscreteGlobalGrids.neighborcount Function

neighborcount(grid::AbstractGrid, c::AbstractCellIndex; connectivity::Connectivity = Vertex()) -> Int

Return length(neighbors(grid, c)). Systems may compute structural degrees without constructing the ring. Subsets count only neighbours within the subset, and an out-of-set cell throws as neighbors does.

A subset cell is interior exactly when length(neighbors(sub, c)) == neighborcount(complete, c), so a border scan needs one ring and one count rather than two rings.

source

Build an adjacency table

adjacency(region) caches a neighbour list for each cell. Row i contains local indices for the neighbours of cell i, in ring order. Reuse the table when an algorithm repeatedly traverses the same cells.

The halo keyword controls neighbours outside the region:

Value	Row entries
0 (default)	Only neighbours belonging to the region
1	Indices into a combined region-and-halo buffer
:mark	Complete neighbour slots, with 0 for an outside neighbour

For halo = 1, haloindices(table) and halocells(table) identify the cells in the buffer's halo portion. See Region boundaries for the relationship between that halo and the region's border.

DiscreteGlobalGrids.adjacency Function

adjacency(region; halo = 0, connectivity = Vertex(), threaded = true) -> AdjacencyTable
adjacency(region, hpos::AbstractVector{<:Integer}; connectivity = Vertex(), threaded = true)

The one-ring of every cell of region at once, cached: adj[p] is the ring of in-region local index p as a non-allocating view, in the counter-clockwise order neighbors states. length(adj) is the region size, which is what an entry is compared against to tell a region slot from a halo slot.

The region types are halo's, the complete grid included — adjacency(levelgrid(sys, l)) is a whole level's adjacency.

The three row shapes

• halo = 0 — rings CLIPPED to the region. Members outside it are dropped, the survivors keep their order, and the row length is the in-region degree.

• halo = 1 — rings COMPLETE, addressing a [region; halo] buffer: 1:n names a cell of the region and n + j the j-th cell of halo(region) under the same connectivity. The halo is walked once, here, and the table keeps it — see haloindices.

• halo = :mark — rings COMPLETE, with 0 where a member is outside the region. Slot geometry with no halo to walk, materialise or fetch, which is what a direction codec reads.

halo above 1 throws. A row exists only for an in-region local index, so a wider receptive field is a wider region: adjacency(grow(region, n); halo = 1).

The anchor

Every row is a window onto the canonical ring the COMPLETE level answers — neighbors(levelgrid(system, level), cell) — whose start is a deterministic property of the system and the cell alone. The complete-width shapes preserve SLOT INDICES — slot k of a row is ring member k, in every table of every region containing that cell — so a direction code is a property of the cell, persistable and stable across tables. The clipped shape preserves ORDER but not slots: dropping a member shifts the ones after it, and recovering slot identity from a clipped row is deliberately impossible. Reach for :mark when the slot is the answer.

The second form takes a halo the caller already walked, as strictly ascending complete-level global indices, and builds the halo = 1 table against it. A neighbour in neither half is an ArgumentError naming the cell, not a short row. The list need not be minimal, only ascending and covering.

threaded builds contiguous chunks in separate tasks and accepts Bool or GeometryOps' True()/False(); the result is identical either way.

source

DiscreteGlobalGrids.Engine.AdjacencyTable Type

AdjacencyTable <: AbstractVector

The product of adjacency: t[p] is the one-ring of in-region local index p as a non-allocating view, counter-clockwise, and length(t) is the region size — the quantity an entry is compared against to tell a region slot (1:length(t)) from a halo slot (length(t) + j).

CSR, and the two arrays are public: t.offsets is length(t) + 1 long with row p occupying t.offsets[p] : t.offsets[p+1] - 1, and t.indices is the flat array those slices cut, so a kernel can loop over it directly instead of over rows.

haloindices and halocells give the buffer's second half.

source

DiscreteGlobalGrids.Engine.halocells Function

haloindices(t::AdjacencyTable) -> Vector{Int}
halocells(t::AdjacencyTable) -> Vector{<:AbstractCellIndex}

The second half of the buffer t's entries address, as complete-level global indices or as cell ids: slot length(t) + j of a row is element j of these, which is element j of halo(region) under the table's connectivity.

Both are empty unless the table was built with halo = 1 — the clipped and marked forms name no cell outside the region, so they walk no halo. halocells resolves ids from the indices on each call.

source

DiscreteGlobalGrids.Engine.haloindices Function

haloindices(t::AdjacencyTable) -> Vector{Int}
halocells(t::AdjacencyTable) -> Vector{<:AbstractCellIndex}

The second half of the buffer t's entries address, as complete-level global indices or as cell ids: slot length(t) + j of a row is element j of these, which is element j of halo(region) under the table's connectivity.

Both are empty unless the table was built with halo = 1 — the clipped and marked forms name no cell outside the region, so they walk no halo. halocells resolves ids from the indices on each call.

source

Compute with neighbourhoods

mapneighbors applies a neighborhood kernel to each cell and collects the results. mapneighbors! writes into an existing destination, and foreachneighbors runs a callback without collecting its return values.

All three sweeps accept neighborhood = Disc(k) or Ring(k). The default is Disc(1). A wider chunk halo supplies input context; use neighborhood to choose the callback reach.

For a dimensional array, choose the callback shape with pass:

Mode	Callback	Use it when
Neighbors()	f(cell, neighbors) with indexed cell handles	The callback needs cell identity or controls data access
Values()	f(cell, value, neighbors) with values at one non-spatial position	Each time or band slice needs the same scalar kernel
NeighborSlices()	f(cell, center_slice, neighbor_slices)	One call needs each cell's whole time series or other non-spatial slice

Values() preserves the input dimensions and visits each non-spatial position. NeighborSlices() makes one call per cell and collects one return value per cell. See its reference below for the exact slice and output contracts.

import DiscreteGlobalGrids as DGG
import DimensionalData as DD

grid = DGG.levelgrid(DGG.HEALPixSystem(), 1)
lookup = DGG.CellLookup(grid)
values = [Float64(t + c) for t in 1:3, c in 1:length(lookup)]
cube = DD.DimArray(values, (DD.Dim{:time}(1:3), DGG.Cells(lookup)))

per_time = DGG.mapneighbors((c, x, ns) -> x + sum(ns), cube;
    pass=DGG.Values(), threaded=false)
per_cell = DGG.mapneighbors((c, xs, ns) -> sum(xs) + sum(sum, ns), cube;
    pass=DGG.NeighborSlices(), threaded=false)
@assert size(per_time) == size(cube)
@assert size(per_cell) == (length(lookup),)
@assert parent(per_cell) ≈ vec(sum(parent(per_time); dims=1))
(; per_time=size(per_time), per_cell=size(per_cell))

(per_time = (3, 48), per_cell = (48,))

Use pass = Values() to receive the centre value and neighbouring values. For kernels that also need geometry, indices or multiple variables, see Requesting neighbour fields. For stored data, see chunked sweeps.

The neighborhood keyword chooses which cells each visit hands the kernel. Disc(k) hands it neighbors(grid, cell, k), every cell within k steps; Ring(k) hands it ring(grid, cell, k), the cells at exactly k steps. The default is Disc(1), the one-ring. Every callback form keeps its arity: only the ring argument widens. The same keyword drives the one-argument neighbors iterator and the chunked mapneighbors!, whose halo follows the radius.

cells = DGG.CellVector(grid)
i = DGG.localindex(cells, cell)

within3 = DGG.mapneighbors((c, nbrs) -> [DGG.localindex(h) for h in nbrs], cells;
    neighborhood = DGG.Disc(3))
at3 = DGG.mapneighbors((c, nbrs) -> [DGG.localindex(h) for h in nbrs], cells;
    neighborhood = DGG.Ring(3))

(; disc_is_neighbors = within3[i] == DGG.neighbors(cells, i, 3),
   ring_is_ring = at3[i] == DGG.ring(cells, i, 3),
   one_ring_leads = within3[i][1:length(DGG.neighbors(cells, i))] == DGG.neighbors(cells, i))

(disc_is_neighbors = true, ring_is_ring = true, one_ring_leads = true)

The order is the one neighbors and ring fix: a Disc(k) ring holds ring 1, then ring 2, out to ring k, each counter-clockwise, and a subset drops its non-members in place. The ring boundaries carry no marker, so a kernel that weights by distance sweeps Ring(j) once per j.

One pass at radius k is a different computation from k one-ring passes. The disc pass weights every cell within k steps once; repeated one-ring passes compound, reaching a cell at distance two through every path of length two, so their weights fall off with distance.

mean3(c, x, nbrs) = (x + sum(nbrs)) / (1 + length(nbrs))
values = sin.(eachindex(cells) ./ 7)

radius3 = DGG.mapneighbors(mean3, cells, values; neighborhood = DGG.Disc(3))
threepasses = foldl((v, _) -> DGG.mapneighbors(mean3, cells, v), 1:3; init = values)
maximum(abs, radius3 .- threepasses)

0.411084937291791

Disc(0) is legal and hands every cell an empty ring; Ring(0) throws, because the sweep already hands the centre to the kernel as its first argument. On systems whose winding is CounterClockwise or Clockwise a k ≥ 2 visit is one shell walk, and H3 answers from libh3's own automaton, both on the stack, so a sequential sweep through Disc(3) and Ring(3) allocates nothing there; on CustomOrder or undeclared-winding systems it sorts each shell by centroid azimuth, so a wide sweep there costs more than k one-ring sweeps. adjacency stays a one-ring table; widen the region with grow when a wider table is wanted.

DiscreteGlobalGrids.Engine.mapneighbors Function

mapneighbors(f, A::AbstractDimArray; spatialdim = nothing, pass = Neighbors(),
             order = StorageOrder(), threaded = true, connectivity = Vertex(),
             neighborhood = Disc(1))
mapneighbors(f, A::AbstractDimArray; needs = (Value(a), Centroid()), ...)

Apply a neighborhood callback along a cell dimension. The default dimension is the first cell lookup; spatialdim accepts a DimensionalData dimension selector. A missing or non-cell dimension raises ArgumentError. neighborhood = Disc(k) selects cells within k steps, excluding the center. Ring(k) selects cells at exactly k steps. The default is Disc(1).

Neighbors passes handles and returns one result per cell. Values passes scalar values and preserves all input dimensions. NeighborSlices passes views across other dimensions and returns one result per cell; it requires at least two dimensions. Concrete tuple returns produce one array per component, using the input wrapper and relevant lookups.

needs replaces the pass contract with f(center, rings) and produces one result per cell. Values come from its Value requests, not implicitly from A. Combining needs with a nondefault pass raises ArgumentError.

Stored arrays use chunked execution for Values() or needs with storage order. Index(Local()) still refers to the original cell axis. A permutation order uses the ordinary traversal instead. See Neighbours and stencils for mode examples and Workflow execution details for routing details.

source

mapneighbors(f, A::AbstractDimArray, dims; kw...)
foreachneighbors(f, A::AbstractDimArray, dims; kw...)

The positional spelling of spatialdim = dims: dims names the cell dimension as neighbors reads it. Every other keyword is as in the two-argument forms.

source

mapneighbors(f, cv; order = StorageOrder(), threaded = true,
             connectivity = Vertex(), neighborhood = Disc(1))
mapneighbors(f, cv, data::AbstractVector; ...)
mapneighbors(f, cv; needs = (Value(data), Centroid()), ...)

Apply f to each cell and its clipped neighborhood. cv accepts a CellVector, PartialGrid, or CellLookup.

neighborhood = Disc(k) selects rings 1:k; Ring(k) selects ring k. The default is Disc(1). Disc(0) gives an empty neighborhood; Ring(0) raises ArgumentError. Clipping removes cells outside the subset and preserves order. Disc results contain no ring-boundary markers. Use separate Ring(j) sweeps when a kernel needs each neighbor's distance.

• Without data, f(cell, neighbors) receives indexed cell handles.

• With a same-order vector, f(cell, value, neighbor_values) receives values.

• With needs, f(center, rings) receives one center entry and one neighbor sequence per requested field. rings[j][i] is field j for neighbor i. A positional data vector cannot be combined with needs.

Results follow collection index order. A concrete tuple return produces one output vector per component. Neighbor order follows neighbors; clipping preserves order but does not preserve missing directional slots.

order is StorageOrder or a permutation of local positions. Invalid permutations raise ArgumentError. Threaded callbacks must be order-independent. A callback failure raises NeighborCallbackError with the cell and index. foreachneighbors discards return values and defaults to serial execution.

See Requesting neighbour fields for field requests and Workflow execution details for caching and scheduling contracts.

source

DiscreteGlobalGrids.mapneighbors! Function

mapneighbors!(dest, f, A::AbstractDimArray; neighborhood = Disc(1),
              halo = max(1, radius), chunks = :auto, spatialdim = nothing,
              connectivity = Vertex(), threaded = true)
mapneighbors!(dest, f, A, plan::MapChunkPlan; neighborhood = Disc(1),
              needs = (Value(a), Centroid()), threaded = true)

Apply a neighborhood kernel chunk by chunk and write results into dest. The default callback is f(cell, value, neighbor_values). With needs, it is f(center, rings) and produces one result per cell.

neighborhood accepts Disc(k) or Ring(k). The default halo is max(1, k). A supplied plan must carry at least k rings; a narrower plan raises ArgumentError. A wider halo is allowed.

dest must support range writes along the cell dimension and hold the result shape. It can be an array, a dimensional cube, or writable stored data.

Index(Local()) names positions on the original cube's cell axis, not positions inside a temporary chunk. Results match the whole-axis sweep. threaded controls work within each chunk.

To assign chunks to workers, build one chunkplan before splitting it. See Sweeping a cube along its chunk lines for ownership and Workflow execution details for field translation and routing details.

source

DiscreteGlobalGrids.Engine.foreachneighbors Function

foreachneighbors(f, A::AbstractDimArray; spatialdim = nothing, pass = Neighbors(),
                 order = StorageOrder(), threaded = false,
                 connectivity = Vertex(), neighborhood = Disc(1))
foreachneighbors(f, A::AbstractDimArray; needs = (Value(a), Centroid()), ...)

Call f for each cell and its neighbourhood without collecting results. spatialdim, pass, needs and neighborhood behave as in mapneighbors.

source

foreachneighbors(f, cv; order = StorageOrder(), threaded = false,
                 connectivity = Vertex(), neighborhood = Disc(1))
foreachneighbors(f, cv, data::AbstractVector; ...)
foreachneighbors(f, cv; needs = (Value(data), Centroid()), ...)

Call f for each cell and its clipped neighbourhood, Disc(1) unless neighborhood says otherwise, discarding its return value. The calling forms, neighborhood, needs and order contracts match mapneighbors. Threading is disabled by default; enabling it requires f to be order-independent.

source

DiscreteGlobalGrids.CellLookups.Values Type

Values()

Pass scalar values to f(cell, value, neighbor_values). On an N-D array, mapneighbors runs the stencil independently along the cell dimension for each index of the other dimensions, and keeps A's dimensions.

source

DiscreteGlobalGrids.CellLookups.Neighbors Type

Neighbors()

Pass indexed cell handles to f(cell, neighbors); the handles read data by indexing the array. This is the default for mapneighbors and foreachneighbors on dimensional arrays. mapneighbors returns one result per cell, on the cell dimension.

source

DiscreteGlobalGrids.Engine.Neighborhood Type

Neighborhood

Supertype of the sweep selectors Disc and Ring. A mapneighbors, foreachneighbors, mapneighbors! or one-argument neighbors sweep takes one as its neighborhood keyword and hands the callback that neighbourhood of each visited cell.

source

DiscreteGlobalGrids.Engine.Disc Type

Disc(k)

Sweep selector: the callback's ring argument is neighbors(grid, cell, k), every cell within k adjacency steps of the centre, rings concatenated outward, clipped to the subset. Disc(1) is the default one-ring sweep and Disc(0) visits every cell with an empty ring. k must be non-negative.

source

DiscreteGlobalGrids.Engine.Ring Type

Ring(k)

Sweep selector: the callback's ring argument is ring(grid, cell, k), the cells at exactly k adjacency steps from the centre, clipped to the subset. k must be positive: the sweep already hands the centre to the callback separately, so Ring(0) is refused with a pointer at Disc(0) and Ring(1).

source

DiscreteGlobalGrids.CellLookups.NeighborSlices Type

NeighborSlices()

Pass views to f(cell, slice, neighbor_slices), with the cell dimension removed from each view. mapneighbors returns one result per cell, on the cell dimension. A one-dimensional array is refused — its per-cell slice is a scalar, which is Values.

source

DiscreteGlobalGrids.Engine.StorageOrder Type

StorageOrder()

Traverse cells in storage order. This is the default for mapneighbors and foreachneighbors.

source

DiscreteGlobalGrids.Engine.NeighborCallbackError Type

NeighborCallbackError

A mapneighbors or foreachneighbors callback that threw during a threaded sweep, naming the cell and subset index it was called with. The callback's own exception is err and is shown as the cause.

One callback failure raises one of these: the sweep waits for every chunk, then reports the failure at the lowest index and drops the rest. The sequential path lets the callback's exception through untouched.

source

Neighbours across levels

member_neighbors finds adjacent members of a mixed-level MultiOrderCellSet. Use it when the cells themselves have different resolutions; ordinary neighbors queries a collection at one level.

DiscreteGlobalGrids.Engine.member_neighbors Function

member_neighbors(set::MultiOrderCellSet, c; connectivity = Vertex()) -> Vector

The members of set that share a boundary with the member c — its neighbours in a set whose cells sit at different levels, so a neighbour may be COARSER than c (an ancestor of one of its system-neighbours) or FINER (a cell deep inside one), and this is the one call that finds both.

c must be a member; anything else is an ArgumentError. c itself is never in the result. The order is the package's one order — counter-clockwise about cell_centroid(system(set), c) seen from outside, as neighbors states — measured on each neighbour's own centroid, whatever level it sits at. The start is the smallest-(level, index) neighbour, and exact azimuth ties break the same way, so the answer is deterministic and independent of how the walk found them.

The algorithm

A member's subtree is a region, and two regions of the hierarchy touch exactly when two of their cells touch at a common depth. The set already names one: its REFERENCE LEVEL L, the depth its covering guarantee is stated at and no shallower than any member. So the question becomes a question about level L, and it is answered without expanding anything:

1. walk c's subtree BORDER at L — EdgeCellIterator, O(border) time and O(depth) memory, and the border is where every contact with the outside lives by its own definition;

2. take each border cell's one-ring at L, under the requested connectivity;

3. map each of those cells to the member that contains it, if any. Members are disjoint subtrees, so there is at most one, and on a sorted-subtree system it is found by binary search over the set's own curve keys: the keys are the members' reference-level range starts, ascending and disjoint, so searchsortedlast plus one range test decides it. A neighbour inside c's own subtree maps back to c and is dropped.

O(|border(c, L)| · degree · log|set|) time and O(|answer|) memory beyond the border walk's own O(depth). The border is the square root of the subtree, so a member many levels above L costs the perimeter of its block and never its area.

Exactness, and where it is the hierarchy's answer rather than geometry's

Two members are neighbours here when their level-L cells are, which is geometric boundary sharing exactly where the system's refinement is congruent — HEALPix, S2 and ISEA4R, whose four children tile their parent, so a member's footprint is the union of its level-L descendants'. Under Edge() the same equivalence holds for shared edges, so vertex-only contact is excluded rather than approximated away.

Where children do not tile their parent — IGEO7 and H3 (aperture 7) and A5 — a member's footprint is NOT its descendants' union, and no level-L statement can be a statement about the drawn polygons; see MultiOrderCoverage for the size of that gap. The answer there is the hierarchy's, which is the same relation border is defined by, and it is consistent with every other subtree verb in this package. That is a carve-out about the SYSTEMS, not about this walk.

A5 pays for its missing primitives here too

Without has_sorted_subtrees there are no curve keys to binary search, so the member lookup is a set built per call, O(|set|); and the border iterator materialises the subtree rather than walking it. The answer is unchanged.

source

Neighbour bounds and ordering declarations

Grid implementations use these declarations to describe neighbourhood size and ring orientation. The query functions above handle the resulting traversal.

DiscreteGlobalGrids.maxneighbors Function

maxneighbors(::A5System, connectivity) -> Int

11 under Vertex() and 5 under Edge().

A5 has corner-only neighbours, so Vertex() and Edge() differ.

resolution	Edge()	Vertex()
0	5	5
1	3	11
≥ 2	5	6, 7 or 8

The global bounds are therefore 5 and 11; the latter occurs at level 1.

source

maxneighbors(::H3System, connectivity) -> Int

6, for either connectivity.

H3's cells are hexagons and pentagons, where sharing a vertex and sharing an edge are the same relation — three cells meet at every vertex and any two of them already share an edge — so Vertex() and Edge() coincide, and the bound is the hexagon's six. The twelve pentagons have five.

source

maxneighbors(sys::AbstractHierarchicalGridSystem, connectivity::Connectivity = Vertex()) -> Union{Int,Nothing}
maxneighbors(grid::AbstractGrid, connectivity::Connectivity = Vertex()) -> Union{Int,Nothing}
maxneighbors(cv::CellVector, connectivity::Connectivity = Vertex()) -> Union{Int,Nothing}
maxneighbors(lk::CellLookup, connectivity::Connectivity = Vertex()) -> Union{Int,Nothing}

A static upper bound on the number of connectivity-neighbours of any cell of sys, at any level, or nothing when the system declares no bound.

The grid and collection forms forward through system. They therefore return the containing system's bound for complete grids and subsets alike: clipping a neighbourhood can shorten it, never exceed it. A standalone grid whose system(grid) === nothing returns nothing.

Sizes the neighbourhood family. An Int bound permits the fixed-capacity stack containers behind neighbors and ring on a subset and behind adjacency; the complete-level verbs and the subtree family never ask for it. Defaults to nothing — never a guessed capacity — and the same machinery then buffers one-rings in a heap Vector, allocating once per cell: identical answers, the slow path. Declaring the bound is a speed decision, not a correctness requirement. Individual cells may have fewer neighbours than the bound.

source

maxneighbors(sys, k::Integer, connectivity = Vertex()) -> Union{Int,Nothing}

A static upper bound on length(neighbors(grid, c, k)), or nothing when the system declares no ring law.

Derived, not declared: neighbors(grid, c, k) is the concatenation of rings 1 through k, so this is sum(maxring(sys, j, connectivity) for j in 1:k) and a system that writes maxring gets it for free. A linear ring law gives the quadratic disc bound M * k * (k + 1) / 2 — 3k(k+1) on a hexagonal system.

k == 0 is 0: neighbors(grid, c, 0) is empty, where ring(grid, c, 0) is c alone.

source

maxneighbors(CopernicusDEMSystem{N}(), connectivity) -> Int

36N + 2 under Vertex() and 6 under Edge(). Both bounds are attained.

source

maxneighbors(ISEA4RSystem(), connectivity) -> Int

9 under Vertex(), 4 under Edge().

Interior cells have eight vertex neighbours. At vertices 0 and 11, five diamond corners meet and a cell can have nine. Other valence-3 icosahedron vertices give seven. At level zero every diamond has six vertex neighbours.

source

DiscreteGlobalGrids.maxring Function

maxring(::H3System, k, connectivity) -> Int

6k, as on any hexagonal system: tight at a hexagon, an over-bound at the twelve pentagons, whose rings hold 5k.

source

maxring(sys, k, connectivity = Vertex()) -> Union{Int,Nothing}

A static upper bound on length(ring(grid, c, k)) for any cell of sys at any level, or nothing when the system declares none.

This is where a system writes its ring scaling law. A tiling whose k-ring is a scaled copy of its one-ring has maxring(sys, k) == M * k for the degree M of that turn — 6k on a hexagonal system, 8k on a vertex-connected quad grid, 4k under Edge(). A system whose rings do not grow linearly declares nothing and keeps the heap path.

An override owns every k, including k == 0, which is 1: ring(grid, c, 0) is c alone. The generic method answers k == 0 with 1, k == 1 with maxneighbors, and nothing beyond.

maxneighbors(sys, k, connectivity) sums this, so one method declares both bounds.

source

maxring(::S2System, k, connectivity) -> Int

8k under Vertex() and 4k under Edge(): the quad lattice's k-ring is its one-ring scaled by k, and the cube's 3-valent corners do not exceed it.

source

maxring(::IGeo7System, k, connectivity) -> Int

6k: a hexagon's k-ring is its one-ring scaled by k. Tight — a hexagonal cell attains it — and an over-bound at the twelve pentagons, whose rings hold 5k.

source

DiscreteGlobalGrids.winding Function

winding(sys::AbstractHierarchicalGridSystem, connectivity = Vertex()) -> Winding
winding(grid::AbstractGrid, connectivity = Vertex()) -> Winding

The order one_ring returns a cell's neighbours in, as a trait the engine can read. Defaults to Unordered().

Why it is a trait. neighbors(grid, c, k) and ring(grid, c, k) for k >= 2 are a breadth-first shell walk, and each shell has to come out in the rotational order the two verbs promise. A declared turn lets the walk carry that order outward from the one-rings it is already reading. Without one it has to measure the order instead — a cell_centroid for every cell of every shell, and a sort — which is correct, slower, and the reason an undeclared system pays for k >= 2 what it does.

Declaring a turn is therefore a speed decision, like maxneighbors, and it is checked rather than assumed: test_grid_interface verifies a declared winding against measured azimuth.

The grid form forwards through system; a standalone grid whose system(grid) === nothing is Unordered().

source

DiscreteGlobalGrids.Winding Type

abstract type Winding

The order a system's one_ring arrives in, declared rather than promised in prose. See winding.

What reads it. The shell walk behind neighbors(grid, c, k) and ring(grid, c, k) for k >= 2 has to know the rotational order of each ring. A declared turn (CounterClockwise, Clockwise) lets it propagate that order outward from the one-rings it already has. Otherwise it measures the order geometrically instead — a cell_centroid per cell of every ring, plus a sort.

Concrete singletons: CounterClockwise, Clockwise, CustomOrder and Unordered.

source

DiscreteGlobalGrids.CounterClockwise Type

CounterClockwise() <: Winding

one_ring is one counter-clockwise turn seen from outside the sphere, starting at the system's own start direction. The order neighbors states, and the one every system in this package declares.

source

DiscreteGlobalGrids.Clockwise Type

Clockwise() <: Winding

one_ring is one clockwise turn seen from outside the sphere. A rotational winding like CounterClockwise, read in the other direction: the shell walk reverses it and is otherwise unchanged, so declaring this costs nothing against declaring the counter-clockwise turn.

source

DiscreteGlobalGrids.CustomOrder Type

CustomOrder() <: Winding

one_ring has a deterministic order the shell walk may not carry outward. Callers may rely on it being stable between calls; k >= 2 is measured by azimuth as under Unordered.

Weaker than CounterClockwise, and not the same as having no turn. A5 is the case: its one-rings are counter-clockwise when measured, but its shells are not rotational copies of them — its rings grow 8, 18, 29, 39 rather than linearly — so there is no outward order to propagate and the geometric sort is the answer rather than a fallback. A system whose one-ring is genuinely unsorted wants Unordered instead.

source

DiscreteGlobalGrids.Unordered Type

Unordered() <: Winding

one_ring promises no order at all, not even stability between calls. The weakest declaration, and the safe default for a system that has not stated otherwise. The shell walk measures azimuth, as it does for CustomOrder.

source

Index