controlledburn

Rasterize polygons without materializing any pixel values. controlledburn computes exact coverage fractions for polygon-grid intersections and returns results as sparse tables: run-length encoded interior cells and individually weighted boundary cells.

The scanline algorithm is O(perimeter) in both time and memory.

Installation

remotes::install_github("hypertidy/controlledburn")

Usage

library(controlledburn)
library(geos)

poly <- as_geos_geometry("POLYGON ((1 1, 9 1, 9 9, 1 9, 1 1))")

# Sparse output — no dense matrix allocated
r <- burn_scanline(poly, extent = c(0, 10, 0, 10), dimension = c(20L, 20L))
r
#> <controlledburn> 20 x 20 grid, 1 geometry
#>   runs:  74 (256 interior cells)
#>   edges: 0 boundary cells
#>   sparsity: 36.0% empty
#> <controlledburn> 20 x 20 grid, 1 geometry
#>   runs:  N (interior cells)
#>   edges: N boundary cells

# Materialise only when you need it
mat <- materialise_chunk(r)

Default grid parameters

With no extent or dimension, controlledburn derives both from the geometry:

r <- burn_scanline(poly)
# extent from wk::wk_bbox(), 256 cells on the long axis

Or specify resolution:

r <- burn_scanline(poly, resolution = 0.5)

Geometry input

Accepts geos_geometry, sfc, wk::wkb(), blob, or raw WKB list (compatible with vapour/gdalraster):

Note that geos::as_geos_geometry() provides interchange for terra::vect().

Shared boundary complementarity

Adjacent polygons with shared edges produce complementary coverage fractions that sum to exactly 1.0 in every boundary cell — no gaps, no overlaps:

left  <- as_geos_geometry("POLYGON ((0 0, 5 0, 5 10, 0 10, 0 0))")
right <- as_geos_geometry("POLYGON ((5 0, 10 0, 10 10, 5 10, 5 0))")

r <- burn_scanline(c(left, right), extent = c(0,10,0,10), dimension = c(20L,20L))

# Coverage sums to 1.0 in every touched cell
mat1 <- materialise_chunk(r, id = 1)
mat2 <- materialise_chunk(r, id = 2)
max(mat1 + mat2)
#> [1] 1
#> [1] 1

Output format

burn_scanline() and burn_sparse() return a list with class "controlledburn" containing:

• runs: data.frame(row, col_start, col_end, id) — interior cells with coverage = 1.0, run-length encoded by row.
• edges: data.frame(row, col, weight, id) — boundary cells with exact partial coverage in (0, 1).
• extent: c(xmin, xmax, ymin, ymax)
• dimension: c(ncol, nrow)

This is the natural output of scanline rasterization — no dense matrix is allocated until materialise_chunk() is called.

Performance

Scanline algorithm scales with perimeter, not area:

Shape	Resolution	Scanline	Dense (burn_sparse)	Speedup
Star	3200×3200	13 ms	225 ms	17×
Jagged	3200×3200	15 ms	136 ms	9×
NC counties	2000×800	29 ms	61 ms	2×

Memory for real-world grids (CGAZ at 32K×16K, ~500M cells): ~50 MB sparse vs ~2 GB dense.

Fast

v <- new(gdalraster::GDALVector, "/vsicurl/https://github.com/mdsumner/geoboundaries/releases/download/latest/geoBoundariesCGAZ_ADM0.parquet")
v$returnGeomAs ## WKB
#> [1] "WKB"
gcol <- v$getGeometryColumn()
v$setIgnoredFields( setdiff(v$getFieldNames(), gcol))
wkbgeom <- wk::wkb(v$fetch(-1)[[gcol]])
v$close()

system.time(burn_scanline(wkbgeom))
#>    user  system elapsed 
#>   0.524   0.015   0.538

system.time(r1 <- burn_scanline(wkbgeom, dimension = c(8192, 4096)))
#>    user  system elapsed 
#>   0.953   0.000   0.953
str(r1)
#> List of 4
#>  $ runs     :'data.frame':   81149 obs. of  4 variables:
#>   ..$ row      : int [1:81149] 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 ...
#>   ..$ col_start: int [1:81149] 5712 5706 5706 5704 5704 5703 5702 5702 5701 5699 ...
#>   ..$ col_end  : int [1:81149] 5712 5712 5715 5717 5719 5719 5719 5719 5718 5718 ...
#>   ..$ id       : int [1:81149] 1 1 1 1 1 1 1 1 1 1 ...
#>  $ edges    :'data.frame':   469461 obs. of  4 variables:
#>   ..$ row   : int [1:469461] 1065 1066 1066 1066 1066 1066 1066 1066 1067 1067 ...
#>   ..$ col   : int [1:469461] 5712 5707 5708 5709 5710 5711 5712 5713 5705 5706 ...
#>   ..$ weight: num [1:469461] 0.0126 0.1343 0.0792 0.296 0.1872 ...
#>   ..$ id    : int [1:469461] 1 1 1 1 1 1 1 1 1 1 ...
#>  $ extent   : num [1:4] -180 180 -90 83.6
#>  $ dimension: int [1:2] 8192 4096
#>  - attr(*, "class")= chr "controlledburn"

Really fast.

system.time(r1 <- burn_scanline(wkbgeom, dimension = c(8192, 4096) * 20))
#>    user  system elapsed 
#>  17.005   0.625  17.629
pryr::object_size(r1)
#> 278.56 MB
tibble::as_tibble(r1$runs)
#> # A tibble: 2,402,331 × 4
#>      row col_start col_end    id
#>    <int>     <int>   <int> <int>
#>  1 21300    114229  114229     1
#>  2 21301    114228  114231     1
#>  3 21302    114227  114232     1
#>  4 21303    114226  114232     1
#>  5 21304    114225  114232     1
#>  6 21305    114224  114232     1
#>  7 21306    114210  114215     1
#>  8 21306    114221  114234     1
#>  9 21307    114209  114235     1
#> 10 21308    114209  114236     1
#> # ℹ 2,402,321 more rows
tibble::as_tibble(r1$edges)
#> # A tibble: 12,006,186 × 4
#>      row    col   weight    id
#>    <int>  <int>    <dbl> <int>
#>  1 21299 114228 0.49016      1
#>  2 21299 114229 0.52762      1
#>  3 21299 114230 0.029348     1
#>  4 21300 114227 0.14377      1
#>  5 21300 114228 0.91647      1
#>  6 21300 114230 0.95528      1
#>  7 21300 114231 0.61698      1
#>  8 21300 114232 0.36007      1
#>  9 21301 114226 0.44189      1
#> 10 21301 114227 0.95691      1
#> # ℹ 12,006,176 more rows
r1[c("extent", "dimension")]
#> $extent
#> [1] -180.00000  180.00000  -90.00000   83.63339
#> 
#> $dimension
#> [1] 163840  81920

History

controlledburn was derived from fasterize by Noam Ross (EcoHealth Alliance), removing Armadillo and raster package dependencies to return bare scanline indexes.

Version 0.1.0 is a complete rewrite using the exactextract algorithm (Daniel Baston, vendored from exactextractr) for exact coverage fractions via a new O(perimeter) scanline sweep.

See also: vaster for grid cell indexing, exactextractr for raster extraction with polygon weights.

Code of Conduct

Please note that the controlledburn project is released with a Contributor Code of Conduct. By contributing to this project, you agree to abide by its terms.