Architecture: how controlledburn developed

Origins: fasterize

controlledburn began as a fork of fasterize by Noam Ross and EcoHealth Alliance. fasterize was a fast polygon rasterizer for R: given a set of polygons and a raster template, it produced a materialised raster using a scanline algorithm based on Wylie et al. (1967). Its key properties:

• Cell-centre rule: a cell is “inside” the polygon if the cell centre falls inside the polygon boundary. Binary in/out, no coverage fractions.
• Pixel functions: Sum, First, Last, Min, Max, Count, Any — applied per-cell across overlapping polygons.
• Dense output: always materialised a full raster object in memory.
• Dependencies: Rcpp, sf (for geometry input), Armadillo (for the raster matrix).

fasterize was fast — much faster than raster::rasterize() — but it required a dense raster template and could not express partial coverage at polygon boundaries.

The sparse insight

controlledburn’s founding insight was that the scanline sweep could emit sparse run-length encoded output instead of filling a dense matrix. For a polygon on a grid, most cells are either fully inside (runs of consecutive columns) or fully outside. Only the boundary cells need individual attention. The output:

• $runs: (row, col_start, col_end, id) — interior cells, full coverage. One record per horizontal span.
• $edges: (row, col, fraction, id) — boundary cells with a coverage fraction in (0, 1).

This is O(perimeter) in both time and memory, rather than O(area). A small polygon on a 100,000 × 100,000 grid produces a handful of run records, not 10 billion pixels.

Integrating exactextract

The exact coverage fractions came from Daniel Baston’s exactextract C++ library. The algorithm traces polygon rings through grid cells, computing analytically exact sub-pixel coverage fractions — not sampling, not approximation. It was the gold standard for zonal statistics.

The first integration vendored exactextract’s C++ code into controlledburn. This brought exact fractions but also brought the GEOS C API as a hard dependency. exactextract’s raster_cell_intersection() function took GEOSGeometry* pointers and called GEOS functions throughout: GEOSGetExteriorRing_r, GEOSCoordSeq_getXY_r, GEOSisEmpty_r, etc.

The pipeline at this stage:

WKB bytes (from R)
  → GEOS WKBReader → GEOSGeometry*
    → exactextract::raster_cell_intersection(grid, ctx, geom)
      → dense Raster<float> (bounding-box-sized matrix)
    → dense_to_sparse() → runs + edges tables

This was burn_sparse(). It worked, but it allocated a dense intermediate for each geometry’s bounding box, required GEOS at runtime, and only handled polygons.

Replacing GEOS: the native core

The next major step was replacing every GEOS dependency with native C++17 code. This happened in stages:

1. WKB parser (~230 lines): reads ISO WKB and EWKB, both byte orders, silently skips Z/M ordinates. No external dependencies. Replaces GEOSWKBReader_read_r.

2. Ring walker: walks polygon rings through the grid cell-by-cell, tracking entry/exit sides and building traversal coordinate lists. Replaces exactextract’s GEOS-dependent ring iteration (GEOSGetExteriorRing_r, GEOSGetInteriorRingN_r, GEOSCoordSeq_getXY_r).

3. Shoelace orientation: determines ring winding (CCW vs CW) via the signed-area formula. Replaces geos_is_ccw.

4. Coordinate min/max envelopes: simple loop over coordinates. Replaces GEOSEnvelope_r and the component-box queries.

The analytical coverage math — the actual geometric computation of how much of a cell is covered by a polygon edge — was kept from exactextract. These files are pure computational geometry with no GEOS dependency: cell.cpp, traversal.cpp, traversal_areas.cpp, box.cpp, grid.cpp, coordinate.cpp, side.cpp, measures.cpp, perimeter_distance.cpp. Nine files, vendored in src/exactextract/.

The result was burn() (originally burn_scanline()), a GEOS-free scanline engine producing the same sparse output.

Unified geometry rasterization

With the native core in place, the engine was extended beyond polygons to handle all geometry types through a single burn() entry point:

• Polygons produce $runs (interior cells, full coverage) and $edges (boundary cells, exact coverage fractions).
• Lines produce $lines with (row, col, length, id) — the absolute length of line within each cell, in CRS units.
• Points produce $points with (row, col, id) — no measure column (a point either falls in a cell or it doesn’t).

Each table’s measure means exactly one thing. No semantic confusion from mixing polygon fractions with line lengths in a shared column.

Two burn modes

The mode parameter on burn() selects between two code paths:

Performance

Benchmarked on CGAZ (218 country polygons, 10.1M vertices):

Grid	Cells	cb approx	fasterize	Winner
256 × 128	33K	1.0s	0.5s	fasterize
4096 × 2048	8.4M	1.1s	0.4s	fasterize
16384 × 8192	134M	1.2s	2.8s	cb
32768 × 16384	537M	1.4s	9.7s	cb
65536 × 32768	2.1B	1.7s	OOM	cb only
131072 × 65536	8.6B	2.4s	OOM	cb only

The crossover is around 134 million cells. Below that, fasterize’s simpler per-cell cost wins. Above it, controlledburn’s sparse output and O(perimeter) memory win decisively — fasterize runs out of memory above ~2 billion cells, while controlledburn reaches 8.6 billion in 2.4 seconds.

Current dependency tree

After removing burn_sparse() and the GEOS dependency:

• C++ core (cpp/): zero external dependencies. Pure C++17.
• R package: cpp11 (LinkingTo), wk (Imports). No GEOS, no sf, no Armadillo.
• Python package: pybind11 + scikit-build-core. No GEOS.
• Vendored exactextract: 9 GEOS-free analytical geometry files (coverage fractions, box intersections, traversal areas).

The three surfaces

The C++ core is canonical. R and Python are thin bindings:

                     ┌──────────────────────┐
                     │    C++ core (cpp/)    │
                     │  scanline engine      │
                     │  WKB parser           │
                     │  exactextract subset  │
                     │  materialize.hpp      │
                     └─────────┬────────────┘
                               │
               ┌───────────────┼───────────────┐
               │               │               │
        ┌──────▼──────┐ ┌─────▼──────┐ ┌──────▼──────┐
        │  R package  │ │   Python   │ │   C++ lib   │
        │  cpp11 shim │ │  pybind11  │ │  CMake      │
        │  burn()     │ │  burn()    │ │  burn()     │
        └─────────────┘ └────────────┘ └─────────────┘

tools/sync-core.sh copies the canonical C++ source from cpp/ into the R package’s src/ and inst/include/ directories. The Python package’s CMakeLists.txt references cpp/ directly.

Shared parity fixtures in fixtures/ (CSV with WKT, WKB hex, and expected results) are read by test suites in all three languages, ensuring cross-language consistency.

Acknowledgements

• Noam Ross and EcoHealth Alliance: the original fasterize scanline algorithm.
• Daniel Baston and ISciences, LLC: the exactextract C++ analytical coverage algorithm, vendored under Apache-2.0.
• The boundary handling work was informed by discussion in hypertidy/fasterize#6 and OSGeo/gdal#14615. Dan Baston has since folded the exactextract logic into GEOS itself.