Requesting neighbour fields

Neighbourhood kernels commonly use fields alongside cell ids: elevations, indices, or centroids. A request declares those inputs once, and the sweep streams them with each clipped ring. For example, needs = (Value(dem), Centroid()) supplies values and geometry without callback lookups.

Each entry of needs is an AbstractNeed. There are four of them, and one — Index — also names the index space it answers in:

need	the visited cell	a ring slot
Cell()	its cell id	the neighbour's cell id
Index(Local())	its index in the collection	the neighbour's
Index(Global())	its globalindex	the neighbour's
Index(T)	its id reindexed to scheme T	the neighbour's
Value(a)	a at its index	a at the neighbour's index
Centroid()	its centroid, on the unit sphere	the neighbour's

Value accepts any vector laid out against the collection, and needs may carry as many fields as the kernel reads. A single field is what Values passes; a request extends that form to multiple fields and geometry. needs is a keyword on mapneighbors and foreachneighbors, alongside order, threaded and connectivity, over a CellVector, a PartialGrid, a CellLookup or a dimensional array.

For the basic kernel forms, see Neighbours and stencils.

The callback

f(center, rings) receives tuples with one entry per need, in the order listed by needs. center holds the visited cell's answers and rings holds the neighbours'. Three things about them are the contract:

• Rings are field-major. rings[j] is need j's value for every clipped neighbour, so a two-need request gives two rings. Slot i of every ring names the same neighbour, in the order neighbors states — counter-clockwise seen from outside the sphere, clipped to membership. A kernel that wants neighbour-major records writes zip(rings...). Each ring is a SmallVector where the system declares a maxneighbors and a Vector where it does not, which is the same rule the Values pass follows.

• Index(Local()) uses the caller's collection. On a complete grid it equals the global index; on a subset it uses the subset's numbering. A chunked sweep translates requests back to that caller axis, including stored Values. Automatic chunked execution applies to needs with StorageOrder(); a permutation order uses the whole-axis path.

• Centroid() supplies unit-sphere centroids. Each task computes a centroid on first use and keeps recent results in a bounded window. A locality-preserving order benefits from the window; a random permutation reduces its reuse. Use Value(table) when a complete precomputed field is faster for the workload. Distance and bearing remain the kernel's choice of radius.

Steepest descent, in one pass

A drainage kernel is the case the request form was built for: it reads one data array and the geometry, and it reads both for every neighbour. On the unit sphere the fall between two cells is an elevation difference over an angle, and the direction is a unit tangent vector — no radius enters, so nothing here commits to a datum.

using DiscreteGlobalGrids
using DiscreteGlobalGrids: Value, Centroid     # public, not exported
using LinearAlgebra: dot, normalize

sys   = IGeo7System()
cells = CellVector(subtree(sys, cellindex(levelgrid(sys, 1), 3), 5))

# A cone rising away from the first cell, so every cell drains towards it.
sink = cell_centroid(sys, cells[1])
dem  = [1000 * acos(clamp(dot(sink, cell_centroid(sys, c)), -1.0, 1.0))
        for c in cells]
length(cells), extrema(dem)

(2401, (2.1073424255447017e-5, 284.7536732137887))

function steepest(center, rings)
    z, p   = center                   # this cell's elevation and centroid
    zs, ps = rings                    # one ring per need, field-major
    best, dir = 0.0, zero(p)
    for (zn, pn) in zip(zs, ps)       # neighbour-major records, on our side
        c = clamp(dot(p, pn), -1.0, 1.0)
        drop = (z - zn) / acos(c)     # metres per radian: no radius involved
        drop > best && ((best, dir) = (drop, normalize(pn - c * p)))
    end
    return (best, dir)
end

fall, direction = mapneighbors(steepest, cells; needs = (Value(dem), Centroid()))

The kernel returned a concrete tuple, so the sweep split it the way it splits any other: one vector per component, both in collection index order.

typeof(fall), typeof(direction)

(Vector{Float64}, Vector{UnitSphericalPoint{Float64}})

Cell 1 sits at the bottom of the cone and has nowhere to go: no neighbour is lower, the loop never fires, and it keeps the zero fall and zero direction it started with. Every other cell drains, and its direction is a unit vector in the tangent plane at its own centroid.

(fall[1], direction[1]), (fall[end], direction[end])

((0.0, UnitSphericalPoint(0.0, 0.0, 0.0)), (978.7829236081853, UnitSphericalPoint(-0.5018483241213358, 0.8245920172468688, 0.2611441453859272)))

Any per-cell quantity

Centroid() is not a special case inside the sweep. It resolves to a cell field — a vector over the collection whose entries are computed by a per-cell function instead of stored — and the sweep then reads that field like any other Value. cellfield builds one, so these two requests are the same request:

using DiscreteGlobalGrids: cellfield          # public, not exported

byname  = mapneighbors(steepest, cells; needs = (Value(dem), Centroid()))
spelled = mapneighbors(steepest, cells;
    needs = (Value(dem), Value(cellfield(cell_centroid, cells))))
byname == spelled

true

Any function of (grid, cell) can be a field — cell_area as readily as cell_centroid. A field is pure: it never mutates and never remembers, so one field is safely read by every task of a threaded sweep, and what remembers is the bounded window the sweep gives each task.

known hands the field what you have already computed. A vector on the collection's cell axis is the complete case, and a complete field is read straight through with no window at all:

table = [cell_centroid(sys, c) for c in cells]
whole = cellfield(cell_centroid, cells; known = table)
mapneighbors(steepest, cells; needs = (Value(dem), Value(whole))) == byname

true

A one-dimensional cube on a Cells dimension is the partial case: the cells it carries are read from it, and every other cell is computed. Nothing requires the whole table, so precompute only the part that pays — the border, say, whose neighbours lie outside the collection's own index range and so miss the window most often:

using DimensionalData: DimArray

edge = cells[collect(border(cells))]
part = cellfield(cell_centroid, cells;
    known = DimArray([cell_centroid(sys, c) for c in edge],
                     (Cells(CellLookup(edge)),)))
length(edge), mapneighbors(steepest, cells;
    needs = (Value(dem), Value(part))) == byname

(240, true)

A field is read by local index, so it must be over the collection being swept. One built over other cells is an ArgumentError raised before the sweep starts, not a silently wrong answer.

The needs

DiscreteGlobalGrids.Engine.AbstractNeed Type

AbstractNeed

One per-neighbour quantity a neighbourhood sweep streams: Cell, Index, Value or Centroid.

A tuple of these is the needs keyword of mapneighbors and foreachneighbors. The callback is then f(center, rings), where center holds one entry per need for the visited cell and rings holds one ring per need — field-major, so rings[j] is need j's value for every clipped neighbour and slot i of every ring names the same neighbour.

source

DiscreteGlobalGrids.Engine.Cell Type

Cell()

Request each neighbour's cell identity, in the system's canonical id scheme. The ring's element type is the collection's own (eltype(cv)), and the center entry is the visited cell. Use Index to ask for the same cell in another scheme.

source

DiscreteGlobalGrids.Engine.Index Type

Index(Local())
Index(Global())
Index(T::Type{<:AbstractCellIndex})

Request each neighbour's index in one named space: Local for the collection's own 1:length(cv), Global for the complete grid at that level, or an id type listed in cellindextypes(system(cv)) for the same cell re-encoded, as reindex answers.

Index(Local()) names the collection the caller passed — the CellVector, or the cube's cell axis — not any chunk a sweep splits it into. The two integer spaces coincide on a complete grid and differ on every subset.

source

DiscreteGlobalGrids.Engine.Local Type

Local()

The index space of the collection the sweep was called on: 1:length(cv), the same numbers localindex answers with. An argument of Index.

source

DiscreteGlobalGrids.Engine.Global Type

Global()

The index space of the complete grid at the collection's level: the numbers globalindex answers with. An argument of Index.

source

DiscreteGlobalGrids.Engine.Value Type

Value(data::AbstractVector)

Request each neighbour's entry in data, which must be laid out against the collection: index k of the vector is index k of the collection, so axes(data) == (Base.OneTo(length(cv)),). The ring's element type is eltype(data).

Any number of Values may appear in one request, each with its own element type; the sweep reads them all from the one membership clip it already made for the cell.

source

DiscreteGlobalGrids.Engine.Centroid Type

Centroid()

Request each neighbour's cell_centroid on the unit sphere, as a GO.UnitSphericalPoint{Float64}. Distances and bearings need a radius and stay downstream of the sweep.

Centroid() is Value(cellfield(cell_centroid, cv)) asked for by name, and the two answer identically: the sweep gives each task a bounded window over the field, keyed by local index, and computes an entry on its first read inside that window. The values are exactly what cell_centroid answers; what the window changes is how often it is called — about once per cell wherever the visit order keeps neighbours close in local index, which storage order and a locality-preserving order do and a random permutation does not.

Spell the field out to hand the sweep what you have already computed: Value(cellfield(cell_centroid, cv; known = table)) for the whole collection, read straight through with no window at all, or a cube over a subset for the part of it you have.

source

DiscreteGlobalGrids.Engine.cellfield Function

cellfield(f, cv::CellVector; known = nothing) -> CellField

A vector over cv's cell axis whose kth entry is f(cv.grid, cv[k]) — cellfield(cell_centroid, cv) is the collection's centroid field — with known naming entries the caller has already computed, which are read instead of called for.

known accepts

• nothing, the default: every entry is computed on read;

• an AbstractVector with axis 1:length(cv): the field is complete and nothing is ever computed;

• a one-dimensional cube over a subset of cv — a DimensionalData array on a Cells dimension at the same system and level — whose cells are read from it and whose absent cells are computed.

The element type is f's return type for one cell of cv, inferred once when the field is built; an f whose return type cannot be inferred is an ArgumentError, as is a known whose elements are not of that type or whose shape or extent is none of the three above.

Reading is pure: the field never mutates and never remembers, so one field is safely shared by concurrent readers and reading the same index twice costs what reading it once did. Remembering belongs to the reader — a sweep handed Value(field) gives each of its tasks a bounded window over the field and computes an entry on its first read inside that window, which is why the field must be over the collection being swept and is rejected when it is not.

source

Index