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:
fasterize was fast — much faster than
raster::rasterize() — but it required a dense raster
template and could not express partial coverage at polygon
boundaries.
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.
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 tablesThis 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.
The next major step was replacing every GEOS dependency with native C++17 code. This happened in stages:
WKB parser (~230 lines): reads ISO WKB and EWKB,
both byte orders, silently skips Z/M ordinates. No external
dependencies. Replaces GEOSWKBReader_read_r.
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).
Shoelace orientation: determines ring winding
(CCW vs CW) via the signed-area formula. Replaces
geos_is_ccw.
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.
With the native core in place, the engine was extended beyond
polygons to handle all geometry types through a single
burn() entry point:
$runs (interior
cells, full coverage) and $edges (boundary cells, exact
coverage fractions).$lines with
(row, col, length, id) — the absolute length of line within
each cell, in CRS units.$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.
The mode parameter on burn() selects
between two code paths:
Uses the full exactextract walker to compute exact coverage fractions for every boundary cell. The walker traces each ring through the grid, building coordinate-level traversal records, then computes analytical cell-intersection areas. O(perimeter) time and memory.
This is the mode to use when exact area conservation matters — the total burned coverage equals the polygon’s true area to floating-point precision.
Reimplements the fasterize cell-centre rule with a dedicated lightweight sweep that bypasses the walker entirely. For each polygon edge, it computes x-intercepts at each row’s scanline y-midpoint, accumulates winding numbers per row, and sweeps left-to-right to emit runs. ~120 lines of C++.
No traversal coordinates, no CellRecord objects, no
BoundaryCellRecord bookkeeping — just edge-row
intersections and winding arithmetic. The output is runs only, no
$edges for polygons.
The boundary convention is left-inclusive (a cell whose centre falls exactly on a polygon edge is classified as “inside”), matching fasterize’s behaviour. On constructed geometries (aligned rectangles, offset rectangles, triangles, holes), approx mode produces cell-for-cell identical results to fasterize.
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.
After removing burn_sparse() and the GEOS
dependency:
cpp/): zero external
dependencies. Pure C++17.cpp11 (LinkingTo),
wk (Imports). No GEOS, no sf, no Armadillo.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.