Skip to contents

Weighted spatial sampling

Source: R/wt_sample.R
wt_sample.Rd

Take a spatial sample based on a weights dataset.

Usage

wt_sample(x, size, weights, ...)

# S3 method for class 'sfc'
wt_sample(x, size, weights, ...)

# S3 method for class 'SpatVector'
wt_sample(x, size, weights = NULL, ...)

# S3 method for class 'sf'
wt_sample(
  x,
  size,
  weights = NULL,
  power = 1,
  method = "random",
  area_weight = 0,
  replace = TRUE,
  of_largest_polygon = TRUE,
  ...
)

# S3 method for class 'SpatRaster'
wt_sample(
  x,
  size,
  weights = NULL,
  power = 1,
  method = "random",
  replace = TRUE,
  ...
)

Arguments

x

A spatial dataset containing either vector data (sf, sfc, SpatVector) or raster data (SpatRaster). This dataset both describes the extent of the sampling area and contains the feature used to weight the sampling.

size

Number of points to sample.

weights

If x is of class sf, SpatVector or SpatRaster, specifies the column by which to weight the sampling. Can also be NULL in which case sampling will be based on the first column or layer in the dataset. If x is of class sfc, is used to specify the weights vector directly, as sfc objects contain no information about features.

...

If x is of class sf, further arguments passed to the custom function passed to method if a function was provided. Otherwise, arguments passed to st_sample.sf.

power

Numeric value that scales the impact of weighting.

  • power = 0 draws an unweighted sample

  • 0 < power < 1 reduces the impact of weighting

  • power = 1 draws normal weighted sample

  • power > 1 increases the impact of weighting

method

How to generate points from sampled polygons/cells.

  • "centroid" computes the geographic center.

  • "random" performs random polygon sampling. If x is a vector dataset, uses st_sample(..., type = "random"). Can be very slow for larger polygons. If x is a raster dataset, generates uniformly distributed raster values within the cell extent using runif.

  • "random_parallel" is a parallelized alternative to "random". Only available if x consists of vector data and the future.apply package is installed. This can be very useful if traditional random spatial sampling is too slow. Only available for polygon sampling.

  • "random_cpp" is a C++ alternative to "random". It is much faster than both other random methods but I wrote it myself so there's a reasonable chance it produces junk. Only available for polygon sampling.

  • Can also be a function or a purrr-style lambda taking at least one argument (the sampled points/polygons) and returning an object that can be coerced to sfc. If x is a SpatRaster, the function can take an additional argument (the weights raster passed to weights). Further optional arguments to these custom functions can be provided using ....

Defaults to "centroid" for raster data and "random" for vector data. Generally, "centroid" makes more sense for regular high-resolution data (e.g. grids) while "random" is more natural for irregular low-resolution data (e.g. administrative divisions).

area_weight

Weight scaling applied to compensate for polygon size. If not 0, computes an additional weight based on the polygon area. Smaller polygons are weighted less than larger polygons to avoid a sort of "double weighting" where busy (i.e. economically stronger or more populated) areas are divided into more and smaller divisions. Defaults to 0, i.e. no area weighting. Ignored if x is a SpatRaster.

replace

Passed to the replace argument in sample. If TRUE, allows the sampling to select the same polygons/cells multiple times. In this case, method should be set to "random" or a custom function if sampling identical locations is not desired. Defaults to FALSE for raster data and TRUE for vector data. See also the method argument.

of_largest_polygon

Passed to st_centroid. Ignored if method is not "centroid". If TRUE and x consists of multipolygons, computes the centroid only of the largest polygon. Generally, this should be TRUE to prevent sampling outside of the polygon boundaries.

Value

An object of class sfc containing the sampled point geometries.

Examples

library(ggplot2)
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE
library(terra)

# Retrieve weights data
bel_pop <- wp_data("BEL")

# Demonstrate different weighting powers
samp_0 <- wt_sample(bel_pop, size = 1000, power = 0)
samp_5 <- wt_sample(bel_pop, size = 1000, power = 0.5)
samp_10 <- wt_sample(bel_pop, size = 1000, power = 1)
samp_15 <- wt_sample(bel_pop, size = 1000, power = 1.5)
groups <- paste("Weight:", rep(seq(0, 1.5, 0.5), each = length(samp_0)))
geom <- c(samp_0, samp_5, samp_10, samp_15)

gr_samp <- st_sf(group = groups, geometry = geom)
ggplot(gr_samp) +
  geom_sf() +
  facet_wrap(~group) +
  theme_void()


# Demonstrate impact of method
bel_pop_agg <- aggregate(bel_pop, fact = 30)
n <- nrow(as.data.frame(bel_pop_agg))
samp_cent <- wt_sample(bel_pop_agg, size = n, method = "centroid")
samp_rand <- wt_sample(bel_pop_agg, size = n, method = "random")
plot(bel_pop_agg, main = "Red dots = centroid, green dots = random")
plot(samp_cent, add = TRUE, pch = 16, col = "red")
plot(samp_rand, add = TRUE, pch = 16, col = "green")


bel_pop_agg <- project(as.polygons(bel_pop_agg), "EPSG:3035")
# Random sampling can be sped up by using parallelization
# This only works for polygons
if (requireNamespace("future.apply")) {
  wt_sample(bel_pop_agg, size = n, method = "random_parallel")
}
#> Geometry set for 36 features 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 3850320 ymin: 3030500 xmax: 4000022 ymax: 3141552
#> Projected CRS: ETRS89-extended / LAEA Europe
#> First 5 geometries:
#> POINT (3938234 3110710)
#> POINT (3936725 3109539)
#> POINT (3928676 3099564)
#> POINT (3932532 3106660)
#> POINT (3923066 3106341)

# ... or by using C++
wt_sample(bel_pop_agg, size = n, method = "random_cpp")
#> Geometry set for 36 features 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 3853743 ymin: 3039158 xmax: 4017302 ymax: 3138570
#> Projected CRS: ETRS89-extended / LAEA Europe
#> First 5 geometries:
#> POINT (3900729 3086620)
#> POINT (3893619 3077160)
#> POINT (3976659 3082146)
#> POINT (3930415 3123362)
#> POINT (3929091 3123687)

if (FALSE) { # \dontrun{
library(tidycensus)

co_pop <- get_decennial(
  geography = "tract",
  variables = "P1_001N",
  state = "CO",
  year = 2020,
  geometry = TRUE
)

# The `area_weight` argument can be useful if polygon sizes are of
# drastically different sizes
co_unweighted <- wt_sample(co_pop["value"], size = 100, power = 0)
co_area_weighted <- wt_sample(co_pop["value"], size = 100, power = 0, area_weight = 1)
plot(st_geometry(co_pop))
plot(co_unweighted, pch = 16)
plot(co_area_weighted, pch = 16)} # }