INLASPDERasters.jl

INLASPDERasters

Raster glue for INLASPDE: extract covariate values from a Rasters.Raster at mesh vertices. Prediction-to-raster and uncertainty surfaces (M2, M3) will land here as separate milestones.

See plans/plan.md for the package roadmap and CLAUDE.md for scope and style rules. This is not a standalone package: it depends on INLASPDE and is only meaningful in conjunction with a fitted SPDE model.

Exceedance(c)

Wrapper requesting per-cell posterior exceedance probabilities P(u(s) > c | y) from the sample-based predict_raster overload.

Probability functionals do not compose linearly under barycentric averaging, so exceedance rasters cannot be derived from Gaussian- approximation summaries — the honest path is sample-based: draw joint posterior samples of the latent SPDE field, project each draw through the same mesh→raster projector, and reduce with mean(η_samples .> c, dims = 2).

Example

rng = Xoshiro(1234)
ex = predict_raster(rng, model, res, template;
    component = "SPDE2[3]", quantity = Exceedance(log(500.0)),
    n_samples = 2000)

_check_crs(rs_crs, mesh_crs) -> Nothing

Assert mesh_crs == rs_crs when mesh_crs !== nothing. The keyword defaults to nothing everywhere it appears — back-compat with v0.2.x "trust the caller" CRS handling. When supplied, CRS mismatches at the API boundary raise an informative error rather than producing silent geographic foot-guns (ADR-041).

_resolve_spde_component(model::LatentGaussianModel, id; time_index = nothing)
    -> (i::Int, mesh::INLAMesh, slice::AbstractRange{Int})

Locate a raster-projectable spatial component inside model.components and return its 1-based index, the retained INLAMesh, and the component-local index range that picks out the spatial slice the projector consumes. Used by the predict_raster(model, res, template) overloads to bridge a random_effects block (SPDE-flavored vertex vector or Kronecker-structured space-time block) to a barycentric mesh→raster projector.

Accepted spatial components

• SPDE2 and SPDE2NonStationary constructed via the mesh-aware constructor (so the mesh field is populated per ADR-036). Slice is 1:n_v (identity).
• KroneckerComponent(spatial, temporal) where spatial is one of the two SPDE flavours. Requires time_index::Integer ∈ 1:length(temporal); slice is time_index:n_t:n_v · n_t, picking the spatial vector at the requested time slot under KroneckerMapping's time-inner-index flattening (x[(s - 1) · n_t + t]).

Accepted id shapes

• Int — 1-based component index.
• String — matches LatentGaussianModels._component_name(c, i) (e.g. "SPDE2[3]", "SPDE2NonStationary[2]", "KroneckerComponent[1]").
• Type{<:SPDE2}, Type{<:SPDE2NonStationary}, or Type{<:KroneckerComponent} — auto-locate the unique component of that type; throws if zero or multiple match.

Errors

Raises ArgumentError with a user-readable enumeration of the model's components when the id cannot be resolved or the resolved component is not raster-projectable. KroneckerComponent paths additionally require time_index; non-Kronecker components reject time_index ≠ nothing.

extract_at_mesh(raster::Raster, mesh::INLAMesh;
                method::Symbol = :bilinear,
                outside::Symbol = :error,
                missingval = NaN,
                mesh_crs = nothing) -> Vector{Float64}

Sample raster at each vertex of mesh and return a length-num_vertices vector of values.

Arguments

• raster::Raster — a 2D raster with X and Y dimensions. The raster must be defined on regular, monotonically-ordered coordinates along both axes (ascending or descending — both are supported).
• mesh::INLAMesh — the SPDE mesh. Mesh vertex coordinates are assumed to be in the same CRS as raster. INLAMesh does not currently carry CRS metadata; pass mesh_crs to assert against Rasters.crs(raster) at the API boundary.

Keywords

• method = :bilinear — one of :bilinear or :nearest. Bilinear interpolation reproduces affine raster fields exactly at mesh vertices. Nearest-neighbour is useful for categorical covariates.
• outside = :error — policy for vertices outside the raster extent: :error throws, :missing substitutes missingval.
• missingval = NaN — the sentinel inserted for outside-domain vertices when outside = :missing.
• mesh_crs = nothing — optional CRS assertion (ADR-041). When supplied, must equal Rasters.crs(raster) or an ArgumentError is raised. Default nothing preserves the v0.2.x "trust the caller" behaviour. Reprojection is out of scope; if the CRSs differ, pre-project one side before calling.

Returns

Vector{Float64} of length num_vertices(mesh) with one extracted value per vertex, ordered by vertex index.

Notes

• For bilinear sampling, a vertex at the raster edge still has a well-defined value (the bracketing cell collapses to the edge).
• If the raster itself contains missing values and a bracketing cell has any missing corner, the returned value is NaN for that vertex under :bilinear, and missingval under :nearest if the nearest cell itself is missing.

predict_raster(values::AbstractVector{<:Real},
               mesh::INLAMesh,
               template::Raster;
               outside::Symbol = :missing,
               missingval::Real = NaN) -> Raster

Project a vertex-valued field onto a raster grid matching template.

Builds a barycentric INLASPDE.MeshProjector from mesh to the cell centres of template and applies it to values. The returned raster has the same dimensions, extent, resolution, and dim order as template; only the underlying array is replaced.

Arguments

• values — length-num_vertices(mesh) vector: the per-vertex field to project. Typically the posterior mean (or any linear functional of it) of an SPDE component.
• mesh — the SPDE mesh values live on.
• template — a 2D Raster with X and Y dimensions that defines the target grid.

Keywords

• outside = :missing — policy for raster cells that fall outside the mesh domain: :error throws, :missing substitutes missingval.
• missingval = NaN — sentinel used when outside = :missing.

Returns

A Raster with the same dims as template whose data is the projected field. Cells outside the mesh domain carry missingval.

predict_raster(model::LatentGaussianModel, res::INLAResult,
               template::Raster;
               component,
               quantity::Symbol = :mean,
               level::Real = 0.95,
               outside::Symbol = :missing,
               missingval::Real = NaN,
               mesh_crs = nothing,
               time_index = nothing) -> Raster

Project a posterior summary of a raster-projectable component onto a raster grid matching template. Wraps the vertex-vector predict_raster primitive — the user-facing entry point lifts the per-component slice out of LatentGaussianModels.random_effects and forwards the corresponding vector through the barycentric mesh→raster projector.

Arguments

• model — a fitted LatentGaussianModel containing at least one raster-projectable spatial component: SPDE2, SPDE2NonStationary, or KroneckerComponent(spatial::SPDE2-flavored, temporal). The spatial child must be constructed via …(mesh::INLAMesh; …) so the mesh is retained per ADR-036.
• res — the INLAResult returned by inla(model, y; …).
• template — a 2D Raster defining the target grid; the returned raster shares its dims, extent, resolution, and dim order.

Keywords

• component — selector identifying the spatial component to project. Accepts Int (1-based index), String (matching the random_effects key, e.g. "SPDE2[3]", "SPDE2NonStationary[2]", "KroneckerComponent[1]"), or Type{<:SPDE2} / Type{<:SPDE2NonStationary} / Type{<:KroneckerComponent} (auto-locates the unique component of that type, errors if there are zero or multiple).
• quantity = :mean — which posterior summary to project: one of :mean, :sd, :lower, :upper.
• level = 0.95 — credible level forwarded to random_effects for :lower / :upper.
• outside, missingval — forwarded to the vertex-vector predict_raster form.
• mesh_crs = nothing — optional CRS assertion (ADR-041). When supplied, must equal Rasters.crs(template) or an ArgumentError is raised. Default nothing preserves the v0.2.x "trust the caller" behaviour.
• time_index = nothing — required when component resolves to a KroneckerComponent; specifies which 1-based time slot to project. Must be nothing for stationary SPDE2 / non-stationary SPDE2NonStationary components.

Returns

A Raster matching template's dims with the projected per-pixel posterior summary.

predict_raster(rng::AbstractRNG, model::LatentGaussianModel,
               res::INLAResult, template::Raster;
               component, quantity,
               n_samples::Integer = 1000,
               outside::Symbol = :missing,
               missingval::Real = NaN,
               mesh_crs = nothing) -> Raster

Sample-based projection of an SPDE component onto a raster grid. Draws n_samples joint posterior samples via LatentGaussianModels.posterior_sample, slices the SPDE block out of each draw, projects all draws through the barycentric mesh→raster projector in a single sparse-dense GEMM, then reduces the resulting n_cells × n_samples matrix column-wise per quantity.

Quantity dispatch

• quantity = :mean — per-cell posterior mean. Asymptotically agrees with the Gaussian-approximation predict_raster(model, res, template; quantity = :mean) overload up to MC error.
• quantity isa Real with 0 ≤ quantity ≤ 1 — per-cell empirical quantile (Statistics.quantile). E.g. quantity = 0.5 is the posterior median; quantity = 0.025 / 0.975 are the 95% credible bounds.
• quantity isa Exceedance(c) — per-cell exceedance probability P(u(s) > c | y) ≈ mean(η_samples .> c, dims = 2). The honest path for tail probabilities; cannot be derived from Gaussian-approximation summaries.

Performance

P.A * x_samples[spde_range, :] is one sparse-dense GEMM yielding the full n_cells × n_samples projected block. All reductions are column-wise scans of this block. Memory cost is n_cells × n_samples Float64s; pick n_samples accordingly for very large rasters.

Reproducibility

Pass a seeded rng (e.g. Xoshiro(1234)) to obtain bit-wise identical rasters across calls.

KroneckerComponent space-time

Like the Gaussian-approximation overload, this entry point accepts a KroneckerComponent(spatial::SPDE2-flavored, temporal) via component plus a 1-based time_index keyword. The slice picks the spatial vector at the requested time slot from each posterior draw before the sparse-dense GEMM, so the per-cell statistics (mean / quantile / exceedance) reflect the conditional posterior at that single time slot.

See also: the predict_raster(model, res, template; …) overload for the Gaussian-approximation overload (cheaper, no sampling, but cannot produce exceedance rasters).

quantile_rasters(mean::AbstractVector{<:Real},
                 sd::AbstractVector{<:Real},
                 mesh::INLAMesh,
                 template::Raster;
                 z::Real = 1.959963984540054,
                 outside::Symbol = :missing,
                 missingval::Real = NaN) -> NamedTuple

Project per-vertex posterior mean and standard deviation onto a raster grid, together with symmetric Gaussian credible-interval rasters.

Semantics

Each output raster is a linear projection of the corresponding vertex quantity via the barycentric mesh-to-pixel projector P:

• mean = P * mean_v
• sd = P * sd_v
• lower = P * (mean_v - z * sd_v)
• upper = P * (mean_v + z * sd_v)

Projecting sd linearly is a convenient but approximate "vertex-level-quantile" view: the SPDE posterior at vertices is not diagonal, so a cell-level exact standard deviation would require diag(P Σ P'). The per-vertex interval is sharp at the vertex and interpolated linearly inside each triangle; this matches the Gaussian summary fields callers plot on R-INLA outputs.

Keywords

• z = 1.96 — the Gaussian quantile half-width (default = 97.5th percentile, giving a 95% interval). Pass z = 1.645 for 90%, etc.
• outside, missingval — forwarded to predict_raster.

Returns

A NamedTuple (; mean, sd, lower, upper) of four Rasters, all sharing the dims of template.