--- title: "Landscape composition in point buffers" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Landscape composition in point buffers} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") ``` ```{r setup} library(bufferscape) ``` ## 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: | field | points | classes that matter | |---|---|---| | air-quality exposure | monitors, home addresses | road surface, industry, tree cover | | environmental epidemiology | residential addresses in a cohort | greenspace, water, built surface | | food environment | schools, homes | outlet types within walking distance | | vector surveillance | ovitraps, light traps, tick drags | roofing, vegetation, standing water | | WASH | water points, households | sanitation infrastructure, drainage | | landscape ecology | camera traps, nest sites, quadrats | habitat classes, edge, canopy | 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 ```{r} kml <- system.file("extdata", "example_site.kml", package = "bufferscape") res <- buffer_composition(kml, radii = 50, grid_res = 5, verbose = FALSE) ``` Composition, largest first: ```{r} d <- res$long[res$long$area_m2 > 0, c("label_en", "area_m2", "area_w")] head(d[order(-d$area_m2), ], 6) ``` `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: ```{r} res$tanks ``` ## 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: ```{r} 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) ``` 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 ```{r} decay_kernel(c(0, 25, 50, 100), kernel = "exponential", lambda = 45) ``` 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: ```{r, eval = FALSE} 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: ```{r, eval = FALSE} # 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. ```{r} 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) ``` Pass it to `buffer_composition(categories = own)` and every output column, map colour and legend label follows it. Validate first if you like: ```{r, error = TRUE} validate_dictionary(data.frame(id = c(1, 1), category = "a", description = "b")) ``` 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. ```{r, eval = FALSE} 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 ``` ```{r} head(class_palette(scheme = "colorblind"), 4) ``` 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: ```{r, eval = FALSE} 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 ```{r, eval = FALSE} 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: ```{r, eval = FALSE} 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.