The coon-way: Fragmented raster computations

This document showcases and explains how raster computations are performed the coon-way, providing all the information necessary to use GeoRacoon’s parallelization tools efficiently.

Loosely put, the coon-way consists in decomposing raster-based operations into sets of independent block-wise computations and providing a high-level interface to keep the user out of the weeds that such a decomposition brings about. Instead of handling the complete raster in one operation, the coon-way divides it into spatial blocks that can then be processed sequentially or in parallel, depending on the available hardware.

There is no single, universal approach for the decomposition of raster-based operations. For each application, the feasibility and benefits of a decomposition depend entirely on the internal logic and the underlying tools implementing it. In GeoRacoon, we exploit GDAL’s capacity to read and write only blocks in a raster file and combine this with specific formulations of the logic in two particular applications:

1. Applying a filter in the form of a convolution with a distribution (Gaussian in most cases) on a finite-sized kernel.

2. Performing multiple linear regression over the pixels in a raster.

Trading time for memory

The benefits of such decompositions fall along two axes: 1. the amount of memory (RAM) required, and 2. the duration (“wall time”) needed to carry out the computation.

As we will see, the decompositions can be configured in a way to maximise the benefits along one axis at the expense of the other. In a way, the coon-way allows the user to trade RAM for runtime and vice-versa.

This makes it possible to adapt the computation to hardware restrictions: Smaller blocks and fewer workers reduce the amount of memory required at any one time but generally increase the number of tasks and their scheduling overhead. Larger blocks and more workers can reduce runtime at the cost of a larger peak memory footprint.

Before we present the fragmented raster computations for a finite-kernel convolution and multiple linear regression, along with benchmark results showcasing some real-world benefits of the coon-way, note that the benchmark results are taken from GeoRacoon’s own benchmark suite. The benchmark results are machine-dependent and should therefore be read as an illustration of the available trade-offs, rather than as universal performance constants.

The benchmark results shown below are based on a square synthetic raster of size RASTER_SIZE x RASTER_SIZE with:

RASTER_SIZE = 20000

Applying a Gaussian filter block-wise

A Gaussian filter requires neighboring pixels, so processing independent blocks requires overlap or halo regions at their boundaries. convster.parallel.apply_filter() computes a halo from the filter parameters and calls riogrande.prepare.create_views() to construct two corresponding view lists. Each worker reads an expanded view including the halo, applies the filter, and returns only its inner view for writing. The expanded views overlap, but the inner views tile the output without overlap. At the outer edge of the raster, the halo is reduced to fit.

filter_border = (2, 2)
views, inner_views = create_views(
    view_size=(6, 6),
    border=filter_border,
    size=(24, 24),
)
inner_xs = sorted({view[0] for view in inner_views})
inner_ys = sorted({view[1] for view in inner_views})
center_xs = set(inner_xs[1:3])
center_ys = set(inner_ys[1:3])
center = [
    (view, inner_view)
    for view, inner_view in zip(views, inner_views)
    if inner_view[0] in center_xs and inner_view[1] in center_ys
    and inner_view[:2] != (12, 12)
]
corner = [
    (view, inner_view)
    for view, inner_view in zip(views, inner_views)
    if inner_view[:2] == (18, 18)
]
selected = [*center, *corner]

fig, ax = plt.subplots(figsize=(7, 7))
ax.set_xlim(0, 24)
ax.set_ylim(24, 0)
ax.set_aspect("equal")
ax.tick_params(which="both", bottom=False, left=False,
               labelbottom=False, labelleft=False)

for view, _ in selected:
    x, y, width, height = view
    ax.add_patch(Rectangle(
        (x, y), width, height, facecolor="#f4a261", alpha=0.4,
        edgecolor="#e76f51", linewidth=1.5,
    ))

for _, inner_view in selected:
    x, y, width, height = inner_view
    ax.add_patch(Rectangle(
        (x, y), width, height, facecolor="#8ecae6",
        edgecolor="black", linewidth=1.5,
    ))

for position in range(1, 24):
    ax.axvline(position, color="0.75", linewidth=0.6, zorder=10)
    ax.axhline(position, color="0.75", linewidth=0.6, zorder=10)

ax.text(12, -1.0, f"border: {filter_border[0]} pixel",
        ha="center", va="bottom")
ax.legend(handles=[
    Patch(facecolor="#8ecae6", edgecolor="black", label="output block"),
    Patch(facecolor="#f4a261", alpha=0.4, edgecolor="#e76f51",
          label="worker view with border"),
], loc="upper center", bbox_to_anchor=(0.5, -0.05), frameon=False)
ax.set_title("Block-wise Gaussian filtering")
fig.tight_layout()

The illustration shows example blocks that can be processed independently. An orange region represents the complete view sent to a worker, including the border pixels needed by the filter. The inner blue region is the corresponding inner block: it is the part of the worker’s result that contributes to the final raster. Because only the blue regions are written, these contributions tile the output without overlap even though the orange worker views overlap.

Benchmark results

These heatmaps show the runtime and peak process-tree memory of the block-wise filter relative to the native implementation. The native value is repeated as the n_jobs = 1 reference column.

filter_fig = plot_routine(
    FILTER_TIME,
    FILTER_MEMORY,
    FILTER_NATIVE_TIME,
    FILTER_NATIVE_MEMORY,
    RASTER_SIZE,
    "Gaussian filtering (apply_filter + bpgaussian)",
)

The figure shows benchmark results for the synthetic filter benchmarks. The two plots show runtime and peak process-tree memory relative to the native, full-raster implementation The native result is the reference case and is shown as n_jobs = 1 (see the apply_filter API documentation). Consequently, a value of 0.5 means half the native runtime or memory, while 1.5 means 150%.

Multiple linear regression as additive block contributions

To decompose multiple linear regression, we use the normal equations. It does not construct one design matrix for the complete raster. Instead, it accumulates the two sufficient statistics needed by ordinary least squares in two passes over the blocks:

Let y be the response and let X be the design matrix. With an intercept, X has one column for every predictor plus a final column of ones. The ordinary least-squares estimate is obtained from the normal equations:

\[\hat{\beta} = (X^T X)^{-1} X^T y.\]

For a large raster, each usable pixel is one observation. Loading all observations into X would make the largest array proportional to the number of pixels. Instead, split the observations into spatial blocks X_1, X_2, ..., X_B and corresponding response vectors y_b. Matrix multiplication is additive over rows, so:

\[X^T X = \sum_{b=1}^{B} X_b^T X_b, \qquad X^T y = \sum_{b=1}^{B} X_b^T y_b.\]

Each block can therefore be processed independently. The order in which the partial matrices are added does not change the mathematical result, although floating-point addition can introduce tiny rounding differences.

The current implementation makes this decomposition explicit in two passes:

1. coonfit.parallel.get_XT_X() computes X_b.T @ X_b for every block and adds the results. The result is the global X.T @ X matrix.

2. coonfit.parallel.compute_weights() inverts this accumulated matrix, then coonfit.parallel.get_optimal_betas() computes and adds X_b.T @ y_b for every block.

3. The coefficient vector is calculated as inv(X.T @ X) @ (X.T @ y).

The regression state that must be retained between blocks is thus only the (p + 1) x (p + 1) cross-product matrix, where p is the number of predictors, plus two vectors of length p + 1. This is the key memory benefit of the formulation: the accumulated mathematical state depends on the number of predictors, not on the number of raster pixels.

This does not mean that every byte of peak memory is independent of raster size. A worker temporarily holds the data for its current block, and the corresponding design matrix has approximately block_pixels x (p + 1) entries. Peak memory also includes the selector, multiprocessing workers, and queued partial results. Block size and the number of workers therefore remain important memory parameters.

Two independent tuning dimensions

block_size controls how much data one worker handles at a time. Smaller blocks reduce the temporary working set, but increase the number of blocks and consequently the scheduling and I/O overhead. Larger blocks usually reduce that overhead, while requiring more memory per worker.

n_jobs controls how many blocks may be processed concurrently. More workers can reduce elapsed time when CPU and I/O resources are available, but each worker can have its own block-sized working set. Increasing n_jobs can therefore increase peak memory substantially.

These controls allow the computation to be adapted to the machine:

• with limited RAM, use fewer workers and smaller blocks, effectively processing a longer sequence of smaller tasks;

• with more RAM, use larger blocks to reduce overhead; and

• with more available CPU resources, increase n_jobs to process more blocks concurrently.

Thus the implementation is not merely a parallel version of the same computation. It exposes a practical runtime–memory trade-off in two dimensions while keeping the regression’s accumulated state small.

Relative runtime and memory

plot_routine(
    PARALLEL_TIME,
    PARALLEL_MEMORY,
    NATIVE_TIME,
    NATIVE_MEMORY,
    RASTER_SIZE,
    "Multiple linear regression (compute_weights)",
)

<Figure size 1300x500 with 4 Axes>

Reading the plots

The left heatmap shows whether a block-size/worker combination reduces elapsed time compared with the native implementation. The right heatmap shows the memory cost of achieving that runtime. A cell with value 0.75 on the runtime plot completes in 75 percent of the native time; a cell with value 1.5 on the memory plot uses 150 percent of native peak memory.

The useful configuration is consequently not necessarily the fastest cell. On a memory-constrained machine, a slightly slower cell below the available memory limit is preferable to a faster cell that cannot run. Conversely, when sufficient RAM and CPU capacity are available, larger blocks and more workers may reduce runtime at the cost of a larger peak working set.

Note

The top row (block_fraction = 1) is not a meaningful parallel configuration: a block as large as the raster leaves a single task, yet GeoRacoon’s parallel routines create and tear down a pool of n_jobs workers, leading to fixed costs independent of the actual work carried out. When the block spans the full raster, the native single-process implementation is the one to choose.

The fixed cost of that worker pool grows with n_jobs in every configuration and shows up in both the wall time and peak memory figures. It is an inherent multiprocessing overhead, not something specific to GeoRacoon. Extra workers only pay off when there are enough blocks to run in parallel; once the pool outgrows the available work, the fixed cost is paid with no offsetting gain. The block_fraction = 1 row makes this plain, since a single block leaves nothing to parallelize.