Landscape composition in point buffers

The problem

Given a set of point locations and a set of land-cover polygons, how much of each class lies near each point?

This is the shape of a great many study designs:

In land-use regression the buffer-and-weight step is the model; in the others it produces covariates. Either way the computation is identical.

“Near” is the hard part. A plain area-within-buffer treats a feature touching the point and one at the buffer edge as equivalent. Weighting by distance fixes that, but only if the weighting is done properly, and the usual shortcuts are not.

bufferscape computes, for every point and every class, the exact area inside a buffer and a distance-decay weighted effective area in which the kernel is integrated over polygon geometry.

It works from vector polygons. That matters at fine scales: rasterising a 3 m wide alley or a 2 m water tank onto a 10 m grid destroys it, and the global land-cover products that operate at that resolution cannot represent the features these designs depend on.

A worked example

kml <- system.file("extdata", "example_site.kml", package = "bufferscape")
res <- buffer_composition(kml, radii = 50, grid_res = 5, verbose = FALSE)

Composition, largest first:

d <- res$long[res$long$area_m2 > 0, c("label_en", "area_m2", "area_w")]
head(d[order(-d$area_m2), ], 6)
#>                     label_en  area_m2    area_w
#> 3               Mature trees 755.8604 323.71040
#> 17      Internal favela road 696.6825 376.50506
#> 28     Narrow internal alley 384.5187 211.43289
#> 12                Vacant lot 383.8599 151.55558
#> 13 Waste accumulation / dump 248.6434  94.77581
#> 4             Clay tile roof 207.0477 118.70774

area_m2 is the exact surface inside the buffer; area_w is that surface after distance weighting. Their ratio is how close to the trap a category sits.

Point features are counted separately:

res$tanks
#> # A tibble: 1 x 6
#>   Ovitrap_ID radius_m tank_open tank_sealed tank_total pool_count
#>   <chr>         <dbl>     <int>       <int>      <int>      <int>
#> 1 SITE_1           50         1          20         21          1

Why not weight at the centroid

The cheap way to weight a polygon is to measure the distance to its centroid. For compact features this is fine. For an elongated feature that passes close to the sampling point it is not: the centroid can sit almost on the point while most of the polygon’s area is far away, so the whole polygon is weighted as though it were adjacent.

Both are reported, so the difference can be inspected directly:

w <- res$long[res$long$area_m2 > 0, ]
w$bias_pct <- 100 * (w$area_w_centroid - w$area_w) / w$area_w
head(w[order(-abs(w$bias_pct)), c("label_en", "area_m2", "bias_pct")], 5)
#>                     label_en  area_m2  bias_pct
#> 17      Internal favela road 696.6825 39.687922
#> 28     Narrow internal alley 384.5187 34.165319
#> 12                Vacant lot 383.8599  1.999895
#> 3               Mature trees 755.8604  1.751168
#> 10 Concrete slab with debris 134.6372 -1.527390

In this buffer the internal road is the worst case. Roads, drainage channels and alleys have exactly the geometry that breaks the centroid approximation, and they are usually the features of epidemiological interest.

Choosing the kernel

decay_kernel(c(0, 25, 50, 100), kernel = "exponential", lambda = 45)
#> [1] 1.0000000 0.5737534 0.3291930 0.1083680

Weights are bounded on [0, 1] and equal 1 at the sampling point, so a weighted area is always between zero and the true area.

The default lambda = 45 m comes from close-kin genetic estimates of mean Aedes aegypti dispersal. For another taxon, set it from the relevant dispersal literature. Because the choice is an assumption, refit across a grid and compare by AIC rather than trusting a single value:

fits <- lapply(c(15, 30, 45, 60, 100), function(L) {
  r <- buffer_composition(kml, radii = 50, lambda = L, verbose = FALSE)
  # ... build the model matrix from r$wide and fit ...
})

Choosing the buffer size

There is no cap on radii. The worked example uses 50 m because that is the scale at which the vector-surveillance question was posed, but the same call takes any radius:

# fine-scale: a few tens of metres
buffer_composition(kml, radii = c(20, 30, 40, 50))

# neighbourhood scale, e.g. walkability or food environment
buffer_composition(kml, radii = c(100, 250, 500))

# land-use regression around an air-quality monitor
buffer_composition(kml, radii = c(50, 100, 300, 1000), lambda = 300)

Two things to keep in mind as the radius grows. Set lambda to the process you are modelling, not to the buffer: a kernel with lambda = 45 inside a 1 km buffer gives almost all the weight to the innermost tenth, which may be right or may be a mistake. And grid_res is the integration step, so leaving it at 1 m over a 1 km buffer is roughly 3 million points per site; scale it with the radius, for instance grid_res = radius / 50.

Scale bars on the figures adapt automatically.

Using your own classification

The 29-class schema shipped with the package is specific to informal settlements in Rio de Janeiro. Nothing else in the package depends on it.

own <- data.frame(
  id          = 1:4,
  category    = c("water", "water", "built", "vegetation"),
  description = c("pond", "channel", "roof", "canopy"),
  label_en    = c("Standing water", "Drainage channel", "Roof", "Tree canopy"),
  fill        = c("#2C7FB8", "#41B6C4", "#BDBDBD", "#31A354"),
  pattern     = c("none", "none", "none", "dots")
)
class_dictionary(own)
#>   id   category description         label_en    fill pattern
#> 1  1      water        pond   Standing water #2C7FB8    none
#> 2  2      water     channel Drainage channel #41B6C4    none
#> 3  3      built        roof             Roof #BDBDBD    none
#> 4  4 vegetation      canopy      Tree canopy #31A354    dots
#>             full_name
#> 1        1_water_pond
#> 2     2_water_channel
#> 3        3_built_roof
#> 4 4_vegetation_canopy

Pass it to buffer_composition(categories = own) and every output column, map colour and legend label follows it. Validate first if you like:

validate_dictionary(data.frame(id = c(1, 1), category = "a", description = "b"))
#> Error: `id` must be unique; duplicated: 1

Only id, category and description are required. Supplying fill is strongly preferred for land cover: without it, colours are generated on a ramp ordered by id and carry no meaning.

Palettes

Four schemes ship with the package, and the same palette argument is taken by the map and the composition chart, so a figure pair can be made to match.

map_composition(res, "SITE_1", palette = "aerial")      # appearance-matched
map_composition(res, "SITE_1", palette = "colorblind")  # colour-vision-safe
map_composition(res, "SITE_1", palette = "greyscale")   # print

head(class_palette(scheme = "colorblind"), 4)
#>   id    fill pattern
#> 1  1 #008B65    none
#> 2  2 #28AD89    diag
#> 3  3 #6EC8AF   cross
#> 4  4 #BB5300    none

A palette of many nominal colours cannot be made safe for colour-vision deficiency – the perceptual space is not large enough, and any such palette contains pairs that converge. The "colorblind" scheme therefore uses colour for the coarse group only, and separates members within a group by a lightness step and a texture, so no class depends on hue as its only cue. That redundancy is what makes a figure accessible, and it survives greyscale printing.

Everything not exposed as an argument can be added afterwards, because the return value is an ordinary ggplot object:

map_composition(res, "SITE_1", legend = FALSE) +
  ggplot2::labs(title = "My title")

One exception: the legend is drawn as text inside the panel rather than as a ggplot guide, so theme(legend.*) has no effect on it. Use the legend, top_n and title arguments instead.

Overlapping polygons

Polygons may overlap, because real surfaces do – a tree crown over a roof is both. Category areas are therefore not constrained to sum to the buffer area, and a buffer can exceed 100% classified. Proportions are deliberately not computed; if you need them, decide your own denominator.

A whole folder

out <- batch_composition("path/to/kml/folder", radii = c(20, 30, 40, 50))
out$summary

This writes a workbook with one modelling-ready sheet per radius, a map per site, a composition chart per site, and a point-feature chart per group. One malformed file is logged and skipped rather than stopping the run.

The workbook can also be written on its own, with control over what it carries:

write_composition_report(res, "composition.xlsx",
                         radii   = c(30, 50),
                         metrics = c("exact", "weighted"),
                         digits  = 2)

Carrying all four metrics for a large dictionary makes for a very wide sheet; dropping the centroid comparison when the methods analysis is not needed roughly halves it.

Point features and reference distances

Alongside the areas, the package counts point features falling inside each buffer and measures straight-line distances from the site to reference features digitised outside it. In the worked example those are water containers and landmarks such as an expressway and a drainage channel; the patterns that recognise them are arguments, so the same machinery serves any point layer.

Distances that were never digitised stay NA. They are not coerced to zero, which would place the site at the feature and invert the sign of its coefficient.