Rasterization, extraction, and zonal statistics

These ten verbs are owned by DiscreteGlobalGrids and share the names Rasters.jl uses, so call them qualified as DGG.rasterize, DGG.extract, DGG.zonal and so on when both packages are loaded. They work in memory on a DGGS cell axis.

Every verb accepts the same geometry inputs: GeoInterface geometries, features, feature collections, nested iterables of those, Tables tables, longitude/latitude Extents.Extent regions, and GO.UnitSpherical.SphericalCap caps. Table inputs read :geometry by default; geometrycolumn=:geom selects another column and geometrycolumn=(:longitude, :latitude) reads point coordinates.

Input coordinates follow the package geometry contract: longitude/latitude for ordinary GeoInterface points, or UnitSphericalPoint coordinates. Edges are great-circle arcs, and these functions do not transform projected coordinates. Cell intersection uses each system's published cell_boundary, including its approximation of curved edges.

Boundary rules

boundary	Polygon membership
:center (default)	The CentroidCovered rule: the polygon covers the cell's canonical interior representative, cell_centroid.
:intersects	Any intersection, including shared edges or vertices.
:touches	Alias for :intersects; different from the DE9IM Touches predicate.
:inside	The entire cell lies within the polygon. Coincident polygon boundaries are permitted.

Lines select intersected cells and points use deterministic containing-cell lookup. shape=:point, :line, or :polygon reinterprets a geometry as that kind — the rings of a polygon become lines, its vertices become points — and plural aliases are accepted. The canonical cell representative is not necessarily the mathematical area centroid.

Destinations

to accepts a grid, partial grid, cell vector or lookup, a cell dimension, dimension tuple, or an existing dimensional array or stack. A system requires a level: to=DGG.H3System(), level=5. A MultiOrderCellSet mixes levels and is not a destination; pass CellVector(set) or a level grid instead. Core allocation returns DimArray/DimStack, and Rasters templates keep their wrappers through the optional Rasters extension. The cell dimension can occur anywhere in a cube; spatial values broadcast over the other dimensions.

Rasterize values

import DiscreteGlobalGrids as DGG
import GeoInterface as GI
using Statistics

grid = DGG.levelgrid(DGG.HEALPixSystem(), 5)
points = [GI.Point((10.0, 25.0)), GI.Point((10.0, 25.0))]
counts = DGG.rasterize(count, points; to=grid)
sums = DGG.rasterize(sum, points; to=grid, fill=[2, 4])
DGG.rasterize!(sums, points; op=+, fill=1)

mean is sum ./ count, as in Rasters: init joins the sum, not the count.

Here is the vector-valued pattern from Rasters' crazy rasterization tutorial:

regions = [GI.Polygon([[(-5.0, 5.0), (15.0, 5.0), (15.0, 30.0),
                       (-5.0, 30.0), (-5.0, 5.0)]]),
           GI.Polygon([[(5.0, 15.0), (25.0, 15.0), (25.0, 40.0),
                       (5.0, 40.0), (5.0, 15.0)]])]
ids = DGG.rasterize(regions; to=grid, op=vcat,
    fill=[[i] for i in eachindex(regions)], eltype=Vector{Int},
    init=Int[], missingval=Int[], boundary=:intersects)

Every cell owns its mutable state, and overlapping cells hold region IDs in input order. Mutating binary operations such as append! also work.

Extract labelled slices

rows = DGG.extract(sums, points; id=true, index=true)

index=true reports the local cell-axis position, not a global cell ID. Points retain their input coordinates; polygon and line rows report cell representative longitude/latitude. Unsampled time and band dimensions stay labelled slices, and skipmissing=true drops a row when any element of its slice is missing.

Reduce zones without repeating selection per slice

stats = DGG.zonal(mean, sums; of=regions, emptyval=NaN)

The default spatialslices=true reduces the cell dimension independently for every nonspatial slice, so Ti × Cells × Band becomes Ti × Band × Zone. This extends Rasters' whole-zone behavior. spatialslices=false reduces the whole selected cube and answers an array with a plain Vector, one entry per zone, and a stack with a NamedTuple of those vectors. A single geometry omits the Zone dimension, and stack layers retain their own nonspatial dimensions.

Reference

DiscreteGlobalGrids.rasterize Function

rasterize([reducer,] data; to, fill, boundary=:center, kw...)

Burn spherical geometries into a Cells array.

• fill: a scalar, one value per feature, a property Symbol, a tuple of Symbols or NamedTuple for several layers, or a function updating each cell.

• Multiple features need a reducer, binary op, or function fill, applied in input order.

• count needs no fill; mean is sum ./ count, with init added to the sum.

• to: grid, cell axis, dimensional array, or stack; the spatial result broadcasts over the other dimensions.

• init, eltype, missingval set the cell state; mutable values are copied per cell.

• threaded parallelises cell writes; a custom op, reducer, or fill function also needs threadsafe=true.

• Pass eltype= for mutating custom reducers: an uninferred element type costs one extra fold per touched cell.

source

DiscreteGlobalGrids.rasterize! Function

rasterize!([reducer,] destination, data; fill, kw...)

Update selected cells in place; untouched cells keep their values.

• Binary op, streaming reducers, and function fills fold onto an existing nonmissing value; gathering reducers such as mean replace it.

source

DiscreteGlobalGrids.extract Function

DiscreteGlobalGrids.extract(A, geometries; boundary=:center, kw...)

Extract cell values as NamedTuple rows.

• Points use their containing cell; lines intersect cells; polygons follow boundary (:center, :intersects/:touches, :inside).

• Rows carry :geometry (the input point, else the cell's lon/lat), with id=true the feature number and index=true the cell-axis position.

• Unsampled dimensions stay labelled slices; name selects stack layers.

• skipmissing=true drops rows whose slice holds any missing element; flatten=false groups non-point rows by feature.

source

DiscreteGlobalGrids.zonal Function

DiscreteGlobalGrids.zonal(f, A; of, spatialslices=true, skipmissing=true, kw...)

Reduce the selected cells of each zone; selection runs once per zone.

• spatialslices: true reduces the cell dimension per other slice, false the whole selected cube, a tuple those dimensions (cell dimension required).

• spatialslices=false returns a Vector (one entry per zone) for an array and a NamedTuple of such vectors for a stack.

• Several zones append Dim{:Zone}; a single geometry returns its result.

• Zones outside the holding give missing; emptyval replaces empty slices, else f receives an empty iterator.

• Means are unweighted cell means.

source

DiscreteGlobalGrids.mask Function

mask(A; with, missingval=missing, invert=false, kw...) replaces uncovered cells.

source

DiscreteGlobalGrids.mask! Function

mask!(A; with, missingval, invert=false, kw...) replaces uncovered cells in place.

source

DiscreteGlobalGrids.boolmask Function

boolmask(data; to, invert=false, kw...) marks cells covered by any geometry.

source

DiscreteGlobalGrids.boolmask! Function

boolmask!(A, data; invert=false, kw...) overwrites A with that coverage.

source

DiscreteGlobalGrids.missingmask Function

missingmask(data; to, kw...) is true on covered cells and missing elsewhere.

source

DiscreteGlobalGrids.missingmask! Function

missingmask!(A, data; kw...) overwrites A with that coverage and missing.

source

Compatibility and execution

The initial compatibility reference is Rasters v0.15. Geometry selection uses a spherical grid/edge dual-tree traversal with prepared point location, conservative subtree acceptance, and exact leaf predicates. Common reducers stream accumulators; arbitrary iterable reducers gather ordered values. Threading uses disjoint output ownership: built-in reducers and operations run threaded, and a custom op, reducer, or fill function runs serially unless threadsafe=true.

progress and verbose are accepted compatibility controls; this version displays no progress bars. File output (filename, suffix, force), raster res/size, crs/mappedcrs, reprojection, fractional coverage, and chunk execution are outside the in-memory API, and unsupported options raise errors. Future chunk execution can reuse the partitioning API once selection plans are stable.