ConservativeRegridding.jl

ConservativeRegridding.jl provides functionality to regrid between two arbitrary grids. A grid is a tessellation of a space into polygons (or grid cells), each with an associated value. Data on a grid is referred to as a field, whereas the grid itself defines the tessellation.

Regridding is performed conservatively, meaning the area-weighted mean is preserved. This is achieved by computing the intersection areas between all combinations of grid cells from the source and destination grids. These intersection areas provide the weights for averaging from neighboring cells.

Quick Start: SpeedyWeather Grid Transfer

Here's an example of regridding between two different geodesic grids from SpeedyWeather.jl:

using SpeedyWeather
using ConservativeRegridding

import GeoInterface as GI
import GeometryOps as GO

# Create random data on two different geodesic grid types
field1 = rand(OctaHEALPixGrid, 4)
field2 = rand(OctaminimalGaussianGrid, 4)

# Get the polygon vertices for each grid (each column is a polygon as (lon, lat) tuples)
# The Regridder constructor handles wrapping these as GeoInterface polygons and
# fixing antimeridian crossings internally
verts1 = RingGrids.get_gridcell_polygons(field1.grid)
verts2 = RingGrids.get_gridcell_polygons(field2.grid)

5×80 Matrix{Tuple{Float64, Float64}}:
 (90.0, 73.7992)  (180.0, 73.7992)  (270.0, 73.7992)  …  (0.0, -73.7992)
 (45.0, 52.8129)  (135.0, 52.8129)  (225.0, 52.8129)     (315.0, -90.0)
 (0.0, 73.7992)   (90.0, 73.7992)   (180.0, 73.7992)     (270.0, -73.7992)
 (45.0, 90.0)     (135.0, 90.0)     (225.0, 90.0)        (315.0, -52.8129)
 (90.0, 73.7992)  (180.0, 73.7992)  (270.0, 73.7992)     (0.0, -73.7992)

# Build the regridder (precomputes intersection areas)
R = ConservativeRegridding.Regridder(verts1, verts2)

# Regrid from field2 to field1
ConservativeRegridding.regrid!(field1, R, field2)

# Regrid in the reverse direction using the transpose
ConservativeRegridding.regrid!(field2, transpose(R), field1)

The key advantage is that the Regridder only needs to be constructed once. After that, both forward regridding (via R) and backward regridding (via transpose(R)) can be performed efficiently without recomputing intersection areas.

Mathematics of Regridding

The Intersection Area Matrix

Conservative regridding is built around a matrix $A$ of intersection areas between source and destination grid cells. For a source grid with $m$ cells and a destination grid with $n$ cells, the matrix $A$ is $n \times m$, where each entry $A_{ij}$ represents the area of intersection between destination cell $i$ and source cell $j$.

The algorithm uses an efficient spatial tree structure (STRtree) to compute only the non-zero intersections, avoiding the $O(nm)$ cost of checking all cell pairs.

Forward Regridding

Let $s$ be a vector of field values on the source grid and $d$ the destination field values. The forward regrid operation computes:

\[d_i = \frac{\sum_j A_{ij} s_j}{a^d_i}\]

or in matrix form:

\[d = \frac{A s}{a^d}\]

where $a^d_i$ is the area of destination cell $i$, and the division is element-wise. The matrix-vector product $As$ yields values weighted by intersection areas, and the division by $a^d$ normalizes to obtain the regridded field values.

Backward Regridding

The same intersection matrix $A$ can be reused for backward regridding by transposing it:

\[\tilde{s} = \frac{A^T d}{a^s}\]

where $a^s_j$ is the area of source cell $j$. The tilde on $\tilde{s}$ emphasizes that the round-trip operation $s \to d \to \tilde{s}$ does not recover the original field exactly—conservative regridding preserves the mean but generally reduces variance.

Conservation Property

Conservative regridding preserves the area-weighted mean:

\[\frac{\sum_i d_i a^d_i}{\sum_i a^d_i} = \frac{\sum_j s_j a^s_j}{\sum_j a^s_j}\]

This property holds when both grids cover the same total area. The areas can be computed from the regridder via row and column sums of $A$:

\[a^d_i = \sum_j A_{ij} \quad \text{and} \quad a^s_j = \sum_i A_{ij}\]

Implementation Details

In ConservativeRegridding.jl:

Sparse storage: The intersection matrix $A$ is stored as a sparse matrix since most grid cell pairs do not intersect.
STRtree acceleration: A Sort-Tile-Recursive tree enables efficient spatial queries to find intersecting cell pairs.
GeometryOps integration: Polygon intersections and area calculations are handled by GeometryOps.jl.
~Antimeridian handling: Polygons crossing the antimeridian are automatically fixed via GeometryOps.fix.~
Normalization: By default, the intersection areas are normalized to improve numerical conditioning.

API Reference

Missing docstring.

Missing docstring for ConservativeRegridding.Regridder. Check Documenter's build log for details.

Missing docstring.

Missing docstring for ConservativeRegridding.regrid!. Check Documenter's build log for details.

Intersection operator interface

These functions let you plug a custom weight kernel into the shared parallel assembly. See Customizing the weights: the intersection operator interface for a walkthrough.

ConservativeRegridding.intersection_areas — Function

intersection_areas(manifold, threaded, dst_tree, src_tree;
                   intersection_operator = DefaultIntersectionOperator(manifold),
                   npartitions = Threads.nthreads() * 4, progress = false,
                   cache = nothing)

Assemble the sparse intersection matrix between src_tree and dst_tree on manifold. Returns a SparseMatrixCSC of the output of the intersection operator.

This is more of a developer level function, which pulls together the intersection operator interface. Users should go through Regridder(…; intersection_operator = …).

This calls out to five functions, which dispatch on intersection_operator:

IntersectionReturnStyle(intersection_operator): return an IntersectionReturnStyle trait object, defining how the operator wants to receive and store the results of its computation. This is usually either OutOfPlaceSingleResult or InPlace.
work_items(intersection_operator, candidate_pairs): return a vector of "work items". For the regular area-of-intersection operator, this is a vector of (src_index, dst_index) pairs. But it can also be more complex, as in the spectral element regridder.
output_matrix_size(intersection_operator, src_tree, dst_tree): return the (nrows, ncols) shape of the sparse matrix. For the regular operator, this is (cell_index_count(dst_tree), cell_index_count(src_tree)). For the spectral element regridder, this is more complex.
output_eltype(intersection_operator, src_tree, dst_tree): return the element type of the sparse matrix. This is usually Float64, but may be different, especially if you wish to build up e.g. a matrix of intersection polygons, rather than just areas.
task_local_operator(intersection_operator): return a private operator for each assembly task, for operators that carry mutable state such as caches.

threaded is a GeometryOpsCore.BoolsAsTypes (True()/False(); convert via booltype(::Bool)). When threaded, work items are partitioned into npartitions chunks assembled on separate tasks via ChunkSplitters.jl.

Pass a SparseMatrixAssemblyCache with cache=... to reuse candidate and serial COO buffers across calls.

source

ConservativeRegridding.IntersectionReturnStyle — Type

abstract type IntersectionReturnStyle

Trait for how an intersection operator delivers its contribution per work item during sparse-matrix assembly. Resolved once via IntersectionReturnStyle(op) at the top of intersection_areas and threaded through the parallel assembly, so the lookup never happens in the hot loop.

Subtypes:

OutOfPlaceSingleResult: kernel op(src_cell, dst_cell) -> area; driver stores the COO triplet.
InPlace: kernel op(rows, cols, vals, item, src_tree, dst_tree) -> nothing pushes its own COO.

Defaults to OutOfPlaceSingleResult, matching DefaultIntersectionOperator.

source

ConservativeRegridding.OutOfPlaceSingleResult — Type

OutOfPlaceSingleResult <: IntersectionReturnStyle

In this IntersectionReturnStyle, the operator is a function op(src_cell, dst_cell) -> result.

The harness around it will store the result in a COO triplet (dst_index, src_index, result), so the operator can remain relatively pure.

The operator must also implement the should_store_result(op, result) -> Bool method, which determines whether the result should be stored in the sparse matrix.

source

ConservativeRegridding.InPlace — Type

InPlace <: IntersectionReturnStyle

In this IntersectionReturnStyle, the operator is a function op(rows, cols, vals, item, src_tree, dst_tree), pushing to the vectors rows, cols, and vals in place.

Here, item is a single entry from the list returned by work_items. Usually, that's a single candidate pair (src_index, dst_index).

This provides maximum flexibility, at the cost of a slightly more complex implementation being required.

source

ConservativeRegridding.work_items — Function

work_items(op, candidate_pairs) -> items

Map candidate (src_index, dst_index) pairs to the input you want to pass to the intersection operator. This can be used to e.g. run a grouping pass on the source index over the candidates, as is done for the spectral element regridder.

By default, this is a no-op and returns the candidate pairs as is.

source

ConservativeRegridding.output_matrix_size — Function

output_matrix_size(op, src_tree, dst_tree) -> (nrows, ncols)

Shape of the sparse matrix intersection_areas assembles for op.

Defaults to (cell_index_count(dst_tree), cell_index_count(src_tree)) — the dense global cell-index domains for dst rows and src columns. Operators whose counts differ from cell index counts (e.g. spectral-element node counts) may override this.

source

ConservativeRegridding.output_eltype — Function

output_eltype(op, [src_tree, dst_tree]) -> eltype

Element type of the sparse matrix intersection_areas assembles for op.

Defaults to Float64. Operators may override this to e.g. return a matrix of intersection polygons, rather than just areas.

source

ConservativeRegridding.should_store_result — Function

should_store_result(op, result) -> Bool

Determine whether the result should be stored in the sparse matrix, after it has been computed.

There is a default implementation for ::Number results across all operators, which is simply result > 0. All other combinations of operator and result type must have explicit dispatches implemented.

source

These are the various intersection operators already implemented:

ConservativeRegridding.DefaultIntersectionOperator — Type

DefaultIntersectionOperator(manifold::GeometryOps.Manifold)

Intersection operator that computes the area of intersection between a source and a destination cell. This is the operator that Regridder uses by default, and it produces the standard first-order conservative regridding weights.

Dispatches to the appropriate intersection algorithm based on the manifold: Foster-Hormann clipping on Planar, and convex-convex Sutherland-Hodgman on Spherical. It uses all the defaults of the intersection-operator interface (OutOfPlaceSingleResult return style, Float64 element type).

On Spherical, task_local_operator hands each assembly task a copy carrying a private clipping-buffer cache, reused by repeated builds on that task. Caches must not be shared across tasks.

source

ConservativeRegridding.IntersectionGridOperator — Type

IntersectionGridOperator(manifold::GeometryOps.Manifold)

Intersection operator that stores the raw polygons of intersection between source and destination cells, rather than their areas.

Pass it to Regridder via the intersection_operator keyword (with normalize = false), and regridder.intersections will then be a sparse matrix whose entry [i, j] holds the polygon of intersection between destination cell i and source cell j.

Since this preserves the exact geometry of intersection between the source and target grids, it's notably useful for determining flux transport during regridding in climate models.

Warning

In the spherical domain, this operator only supports convex cells.

source

and there are two more in the ClimaCore extension, for spectral element based regridding.