spatialsample

Design-based spatial sampling in Python. Decide where to put field plots, household survey points or validation samples for a land-cover map, with spatially balanced designs (GRTS, BAS, LPM), classical designs, exact inclusion probabilities and design weights, Horvitz-Thompson estimation, a spatial-balance metric to compare designs, and the Olofsson et al. (2014) good-practice workflow for map accuracy and area estimation.

The core needs only NumPy. SciPy (faster nearest-neighbour search), Shapely and GeoPandas are optional.

Installation

pip install spatialsample              # NumPy only
pip install "spatialsample[scipy]"     # + SciPy KD-tree for large balance/frame computations
pip install "spatialsample[geo]"       # + Shapely and GeoPandas interoperability
pip install "spatialsample[all]"       # everything

Quick start

import spatialsample as ss

# Study area in a projected CRS (metres): an L-shaped block with a lake (hole).
area = ss.Region.from_coords(
    [(0, 0), (1000, 0), (1000, 600), (600, 600), (600, 1000), (0, 1000)],
    holes=[[(200, 200), (400, 200), (400, 400), (200, 400)]],
)

sites = ss.grts(area, n=50, n_over=10, seed=1)   # 50 base + 10 replacement sites
sites.to_geojson("sites.geojson")                # or .to_csv(), .to_geodataframe()

print(sites.base().weight.sum())                 # 800000.0  (= area; weights in m^2)
print(ss.spatial_balance(sites, domain=area).value)  # Voronoi balance, lower is better

# After the field work: estimate the mean of a variable y observed at the base sites
base = sites.base()
y = observe(base.xy)                             # your measurements
est = ss.estimate(base, y)                       # HT total, ratio mean, local variance
print(est.mean, est.mean_se)

Inputs can be a Region, a RasterMask, any GeoJSON mapping (Polygon, MultiPolygon, Feature, FeatureCollection) or any object with __geo_interface__ (Shapely, GeoPandas). GRTS and LPM also accept an (N, 2) array of candidate sites (a finite frame, such as households, lakes or pixels), optionally with a size measure prob= for unequal-probability sampling.

Coordinates are treated as planar. Reproject longitude/latitude data to a projected or equal-area CRS before sampling; distance-based designs and the balance metric assume linear units.

Designs

Every design returns a Sample with xy, pi (inclusion probability for finite frames, inclusion density per unit area for continuous designs, following Cordy 1993), weight = 1/pi, site_type (base / over), visiting order, optional stratum, index and extra columns (cluster, psu, chainage, pixel row/col, ...).

GRTS and LPM on a continuous area use a dense, randomly offset grid frame (default max(4000, 100 (n + n_over)) cells for GRTS) and move each selected cell centre to a uniform location within its cell, as recent versions of the R package spsurvey do.

Which design is most spread out?

Voronoi spatial balance (Stevens & Olsen 2004): the variance of the summed inclusion probability inside each site's Voronoi polygon, mean((v_i - 1)^2). Zero is perfectly balanced. Measured with ss.compare_designs on the L-shaped area above (800,000 m^2, with a hole), n = 50, 100 independent draws per design, 20,000-point evaluation grid:

Systematic grids are the most even but have a random sample size, no valid design-based variance estimator and are vulnerable to periodic patterns; GRTS, BAS and LPM give fixed sample sizes, exact inclusion probabilities and three to five times better balance than simple random sampling. Reproduce with spatialsample compare area.geojson --n 50 --reps 100.

Does balance pay off?

Estimating the mean of the smooth surface 10 + 0.01 x + 5 sin(y / 150) over the same area with n = 30 (400 draws per design, ratio estimator, ss.estimate defaults):

All biases are within Monte Carlo error. Spatially balanced designs cut the standard error by a factor of 2.7 (GRTS) to 5 (BAS) here; the local neighbourhood variance estimator captures much of that gain while staying on the conservative side for smooth surfaces.

Estimation

ss.ht_total(y, pi)                    # Horvitz-Thompson total (integral for continuous designs)
ss.ht_mean(y, pi, size=area)          # HT mean; without size: Hajek ratio mean
ss.srs_variance(y, pi)                # with-replacement variance of the HT total
ss.local_variance(y, xy, pi, k=4)     # Stevens & Olsen (2003) neighbourhood variance
ss.estimate(sample, y)                # all of the above with a CI, design-aware defaults

local_variance uses neighbourhoods of each site and its three nearest sample neighbours, initial weights proportional to 1/pi balanced to a doubly stochastic matrix by iterative proportional fitting. By default each neighbourhood term is divided by 1 - sum_j w_ij^2, the usual small-sample correction for a weighted variance, which makes the estimator unbiased (relative to the with-replacement variance) when the variable has no spatial structure; pass correction=False for the uncorrected form, which is about 25 % smaller with k = 4. Tests check both properties by simulation.

Map accuracy and area estimation (Olofsson et al. 2014)

import spatialsample as ss

# Section 5 worked example of Olofsson et al. (2014): sample counts, rows = map class,
# columns = reference class; mapped areas in 30 m pixels (0.09 ha each).
matrix = [[66, 0, 5, 4],
          [0, 55, 8, 12],
          [1, 0, 153, 11],
          [2, 1, 9, 313]]
pixels = [200_000, 150_000, 3_200_000, 6_450_000]
labels = ["Deforestation", "Forest gain", "Stable forest", "Stable non-forest"]

res = ss.olofsson(matrix, [p * 0.09 for p in pixels], labels=labels)
print(res.summary())

class                  mapped area     est. area    +/- 95% CI     user     +/-  producer     +/-
-------------------------------------------------------------------------------------------------
Deforestation               18,000        21,158         6,158    0.880   0.074     0.749   0.213
Forest gain                 13,500        11,686         3,756    0.733   0.101     0.847   0.254
Stable forest              288,000       285,770        15,510    0.927   0.040     0.935   0.034
Stable non-forest          580,500       581,386        16,281    0.963   0.021     0.962   0.018
overall accuracy 0.947 +/- 0.018 (95% CI), n = 640

These reproduce the published estimates (for example a deforestation area of 21,158 ha +/- 6,158 ha against 18,000 ha mapped; user's accuracy 0.88, producer's accuracy 0.75 +/- 0.21; overall accuracy 0.95). The paper's sample-size example (conjectured user's accuracies 0.7, 0.6, 0.9, 0.95 and a target standard error of 0.01 for overall accuracy) is ss.sample_size(pixels, [0.7, 0.6, 0.9, 0.95], 0.01) == 641.

To draw that sample from a classified raster:

mask = ss.RasterMask(landcover_array, origin=(x_min, y_max), cell_size=30, nodata=0)
n_h = ss.allocate_accuracy_sample(641, pixels, rare={0: 75, 1: 75})   # fix rare classes
validation = ss.stratified_pixels(mask, dict(zip(mask.class_counts(), n_h)), seed=7)
validation.to_csv("validation_points.csv")      # x, y, row, col, stratum, pi, weight

Command line

# 50 GRTS sites plus 10 replacements
spatialsample grts area.geojson --n 50 --over 10 --seed 1 -o sites.geojson

# other designs: srs, systematic (--shape hex --random-rotation), bas, lpm, poisson,
# stratified (--field zone --allocation proportional --method grts)
spatialsample bas area.geojson --n 40 --over 5 --seed 3 -o sites.csv

# spatial balance of any set of sites (GeoJSON points or CSV with x,y)
spatialsample balance sites.geojson --frame area.geojson

# compare designs on your own study area
spatialsample compare area.geojson --n 50 --reps 50

# Olofsson accuracy assessment from an error matrix CSV (header row and label column optional)
spatialsample olofsson matrix.csv --areas 200000,150000,3200000,6450000 --pixel-area 0.09
spatialsample samplesize --areas 200000,150000,3200000,6450000 --users 0.7,0.6,0.9,0.95 --target-se 0.01

Output is GeoJSON (default, also to stdout when -o is omitted) or CSV by file extension. Each site carries site_id, design, site_type, order, pi and weight.

API overview

Development and testing

pip install -e ".[dev,all]"
pytest --cov=spatialsample
ruff check . && ruff format --check src tests && mypy src

The test suite includes Monte Carlo checks: empirical inclusion frequencies of GRTS and LPM match unequal target probabilities, Horvitz-Thompson estimates are unbiased for finite and continuous designs, GRTS/BAS/LPM are better balanced than simple random sampling, holes are never sampled, and every design is reproducible with a seed.

References

• Cochran, W. G. (1977). Sampling Techniques (3rd ed.). Wiley.
• Cordy, C. B. (1993). An extension of the Horvitz-Thompson theorem to point sampling from a continuous universe. Statistics & Probability Letters, 18(5), 353-362.
• Dumelle, M., Kincaid, T., Olsen, A. R., & Weber, M. (2023). spsurvey: Spatial sampling design and analysis in R. Journal of Statistical Software, 105(3), 1-29.
• Grafström, A., Lundström, N. L. P., & Schelin, L. (2012). Spatially balanced sampling through the pivotal method. Biometrics, 68(2), 514-520.
• Horvitz, D. G., & Thompson, D. J. (1952). A generalization of sampling without replacement from a finite universe. Journal of the American Statistical Association, 47(260), 663-685.
• Olofsson, P., Foody, G. M., Herold, M., Stehman, S. V., Woodcock, C. E., & Wulder, M. A. (2014). Good practices for estimating area and assessing accuracy of land change. Remote Sensing of Environment, 148, 42-57.
• Robertson, B. L., Brown, J. A., McDonald, T., & Jaksons, P. (2013). BAS: Balanced acceptance sampling of natural resources. Biometrics, 69(3), 776-784.
• Stevens, D. L., & Olsen, A. R. (2003). Variance estimation for spatially balanced samples of environmental resources. Environmetrics, 14(6), 593-610.
• Stevens, D. L., & Olsen, A. R. (2004). Spatially balanced sampling of natural resources. Journal of the American Statistical Association, 99(465), 262-278.
• Tillé, Y. (2006). Sampling Algorithms. Springer.

License

MIT, see LICENSE. If you use spatialsample in research, please cite it using CITATION.cff.