Metadata-Version: 2.4
Name: spatialsample
Version: 0.1.1
Summary: Design-based spatial sampling in Python: GRTS, BAS, LPM, stratified and systematic designs, Horvitz-Thompson estimation, spatial balance and Olofsson map accuracy assessment.
Author-email: Yusuf Eminoğlu <yusufeminoglu@gmail.com>
License-Expression: MIT
Project-URL: Homepage, https://github.com/YusufEminoglu/spatialsample_sdk
Project-URL: Repository, https://github.com/YusufEminoglu/spatialsample_sdk
Project-URL: Issues, https://github.com/YusufEminoglu/spatialsample_sdk/issues
Project-URL: Changelog, https://github.com/YusufEminoglu/spatialsample_sdk/blob/main/CHANGELOG.md
Keywords: spatial-sampling,survey-design,grts,balanced-acceptance-sampling,local-pivotal-method,spatially-balanced,horvitz-thompson,stratified-sampling,accuracy-assessment,land-cover,ecology,field-survey,gis
Classifier: Development Status :: 4 - Beta
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Operating System :: OS Independent
Classifier: Topic :: Scientific/Engineering :: GIS
Classifier: Topic :: Scientific/Engineering :: Mathematics
Classifier: Topic :: Scientific/Engineering :: Information Analysis
Classifier: Intended Audience :: Science/Research
Classifier: Intended Audience :: Developers
Classifier: Typing :: Typed
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.20.0
Provides-Extra: scipy
Requires-Dist: scipy>=1.7.0; extra == "scipy"
Provides-Extra: shapely
Requires-Dist: shapely>=2.0.0; extra == "shapely"
Provides-Extra: geo
Requires-Dist: shapely>=2.0.0; extra == "geo"
Requires-Dist: geopandas>=0.12.0; extra == "geo"
Provides-Extra: all
Requires-Dist: scipy>=1.7.0; extra == "all"
Requires-Dist: shapely>=2.0.0; extra == "all"
Requires-Dist: geopandas>=0.12.0; extra == "all"
Provides-Extra: dev
Requires-Dist: pytest>=7.0.0; extra == "dev"
Requires-Dist: pytest-cov>=4.0.0; extra == "dev"
Requires-Dist: ruff>=0.1.0; extra == "dev"
Requires-Dist: mypy>=1.0.0; extra == "dev"
Dynamic: license-file

# spatialsample

[![PyPI version](https://img.shields.io/pypi/v/spatialsample.svg?color=10b981)](https://pypi.org/project/spatialsample/)
[![Python version support](https://img.shields.io/pypi/pyversions/spatialsample.svg?color=3b82f6)](https://pypi.org/project/spatialsample/)
[![License: MIT](https://img.shields.io/badge/License-MIT-blue.svg)](https://github.com/YusufEminoglu/spatialsample_sdk/blob/main/LICENSE)
[![Code style: Ruff](https://img.shields.io/badge/code%20style-ruff-D7FF64.svg)](https://docs.astral.sh/ruff/)

**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

```bash
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

```python
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

| Function | Design | Inclusion probabilities | Reference |
|---|---|---|---|
| `srs` | Simple random sampling of locations | exact, `n / area` | |
| `systematic` | Square or hexagonal grid, uniformly random start, optional random rotation | exact, `1 / cell area` (random `n`) | |
| `stratified` | Stratified by polygons or raster classes; proportional, equal or Neyman allocation; SRS, BAS, GRTS or LPM within strata | exact, `n_h / area_h` | Cochran (1977) |
| `stratified_pixels` | Stratified random sampling of raster pixels without replacement | exact, `n_h / N_h` | Olofsson et al. (2014) |
| `grts` | Generalized Random Tessellation Stratified: randomised quadrant-recursive addresses, systematic selection on the resulting line, reverse hierarchical ordering, oversample sites | exact (finite frames), equal or unequal | Stevens & Olsen (2004) |
| `bas` | Balanced Acceptance Sampling from a random-start Halton sequence (bases 2, 3); oversample continues the sequence | exact, `n / area` | Robertson et al. (2013) |
| `lpm` | Local pivotal method (LPM2), equal or unequal probabilities | exact | Grafström et al. (2012) |
| `poisson_disk` | Minimum-distance (simple sequential inhibition) sampling | approximate, `n / area` | |
| `cluster` | Cluster plots (FIA-style centre + 3 at 120 degrees, 2 x 2, line, cross or custom) around SRS/BAS/GRTS centres | `n_clusters / area` per subplot | |
| `two_stage` | PPS-by-area (or equal) selection of primary units, then SRS inside each | exact; self-weighting with PPS | |
| `points_along_lines` | Systematic points along roads/rivers/hedgerows with random start | exact, `1 / spacing` per unit length | |
| `random_transects` | Straight transects of fixed length and random (or fixed) bearing | `n / area` | |

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:

| Design | Voronoi balance, mean | sd |
|---|---:|---:|
| systematic, square (random rotation) | 0.041 | 0.011 |
| systematic, hexagonal (random rotation) | 0.044 | 0.008 |
| LPM | 0.074 | 0.017 |
| BAS | 0.087 | 0.021 |
| GRTS | 0.105 | 0.024 |
| Poisson-disk (default distance) | 0.118 | 0.033 |
| simple random | 0.352 | 0.106 |

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):

| Design | bias | true SD of the estimate | mean estimated SE |
|---|---:|---:|---:|
| simple random | +0.026 | 0.905 | 0.880 (SRS variance) |
| GRTS | -0.015 | 0.331 | 0.399 (local variance) |
| LPM | -0.026 | 0.282 | 0.423 (local variance) |
| BAS | +0.011 | 0.180 | 0.426 (local variance) |
| systematic, hexagonal | -0.003 | 0.185 | 0.456 (local variance) |

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

```python
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)

```python
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:

```python
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

```bash
# 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

| Area | Names |
|---|---|
| Domains | `Region`, `Polygon`, `RasterMask` |
| Designs | `srs`, `systematic`, `stratified`, `stratified_pixels`, `grts`, `bas`, `lpm`, `poisson_disk`, `cluster`, `two_stage`, `points_along_lines`, `random_transects` |
| Building blocks | `inclusion_probabilities`, `allocate`, `grts_order`, `reverse_hierarchical_order`, `halton`, `lpm2_select` |
| Estimation | `ht_total`, `ht_mean`, `srs_variance`, `local_variance`, `local_weights`, `estimate`, `Estimate` |
| Balance | `voronoi_balance`, `spatial_balance`, `compare_designs` |
| Accuracy | `olofsson`, `AccuracyAssessment`, `sample_size`, `allocate_accuracy_sample` |
| I/O | `Sample.to_geojson`, `Sample.to_csv`, `Sample.to_geodataframe`, `read_region`, `read_strata`, `read_points`, `read_lines`, `read_matrix_csv` |

## Development and testing

```bash
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](https://github.com/YusufEminoglu/spatialsample_sdk/blob/main/LICENSE).
If you use spatialsample in research, please cite it using
[CITATION.cff](https://github.com/YusufEminoglu/spatialsample_sdk/blob/main/CITATION.cff).
