⚙️ How-To Agriculture & vegetation Interpolation & surfaces Agriculture & Crops · step 13 Notebook By · Updated

How to lay out a soil sampling grid in Python

On this page

Problem statement

Soil sampling is expensive per point and the layout decides what the samples can tell you. A regular grid estimates a mean well and a spatial pattern poorly; a zone-based design characterises known zones and cannot discover new ones; a design with no short-distance pairs cannot fit a variogram, so nothing can be interpolated from it afterwards.

The number of samples is usually fixed by budget. The question this guide answers is where to put them, given what the samples are for — and the answer differs for a mean, a map, and a set of management zones.

Quick answer

A regular grid with a random offset, clipped to the field and pulled in from the boundary:

import numpy as np, geopandas as gpd
from shapely.geometry import Point

def grid_samples(field, spacing_m, crs, inset_m=15, seed=0):
    """field: a shapely polygon in a projected CRS, passed separately — shapely has no .crs."""
    rng = np.random.default_rng(seed)
    poly = field.buffer(-inset_m)
    minx, miny, maxx, maxy = poly.bounds
    ox, oy = rng.uniform(0, spacing_m, 2)               # random origin, not a fixed one
    xs = np.arange(minx + ox, maxx, spacing_m)
    ys = np.arange(miny + oy, maxy, spacing_m)
    pts = [Point(x, y) for x in xs for y in ys]
    g = gpd.GeoDataFrame(geometry=[p for p in pts if poly.contains(p)], crs=crs)
    g["sample_id"] = [f"S{i + 1:03d}" for i in range(len(g))]
    return g

The inset matters: headlands are compacted, differently fertilised and unrepresentative, and a sample taken fifteen metres from a hedge is a sample of the hedge.

A regular grid, a zone-stratified design, a nested design and a space-filling design. regular gridmean yes, variogram nozone-stratifiedprecise zone meansnestedfits a variogramspace-fillingfinds new patternsFour designs; each answers a different question spatialworkflow.io
Four designs; each answers a different question and none answers all of them.

What you get — from The data used on this page

field 19.5 ha; zone areas, ha: {'high': 6.7, 'low': 6.7, 'mid': 6.7}

Purpose first, spacing from budget, randomised origin

1. Decide what the samples are for

  • A field mean for a lime or phosphate recommendation: a regular or random grid, enough points to hit the precision you need.
  • A map to interpolate: a grid plus clustered pairs at short lags, or the variogram cannot be fitted.
  • Zone characterisation: stratified by zone, with enough points per zone to estimate each mean.
  • Discovering new patterns: a space-filling design that does not assume the zones.

2. Choose the spacing from the budget and the field

A 100 m grid on a 25 ha field is about 25 points and a 75 m grid about 44 — before the edges take their share. On the 19.5 ha Ferrara field, after a 15 m inset, 100 m gave 17 points and 75 m gave 29 (Example 1).

3. Randomise the grid origin

A fixed origin aligned with the field's corner risks aligning with tramlines, drainage runs or previous applications. A random offset costs nothing and removes the risk.

4. Inset from the boundary

Fifteen to twenty metres. Headlands are consistently different and they are usually managed separately anyway.

5. Add short-lag pairs if a map is wanted

A regular grid has a minimum lag equal to its spacing, so the variogram has no points below it and the nugget cannot be separated from the short-range structure. Adding a handful of close pairs — two points 10 and 25 m apart at a subset of grid nodes — fixes it cheaply.

6. Stratify by zone if the zones already exist

Allocate points proportionally to area, with a minimum per zone, and randomise within each zone. That gives a defensible mean per zone and keeps the coverage.

7. Produce a navigable file with stable identifiers

Points with an identifier, in a format the field computer reads, with the identifier printed on the bag label. The commonest cause of a ruined sampling campaign is the labels not matching the points.

A grid-only lag distribution against one with added short-lag pairs. the lags a design produces75 m grid alone✕no pairs below 75 m✕nugget cannot be estimated✕the variogram is extrapolated✕interpolation smoothing is a guessgrid plus close pairs✓pairs at 10 and 25 m✓nugget constrained directly✓40 extra points on a 29-point grid✓to reach 24 pairs under 20 mMeasured on a 19.5 ha field: short lags are the expensive part of the budget spatialworkflow.io
A regular grid has nothing below its spacing, which is exactly where the variogram is decided.

Grid, stratified and space-filling designs

The data used on this page

A real field: the 22 ha corn field at Porto Felloni near Ferrara whose three yield-monitor harvests are published on Zenodo (CC BY 4.0), used on yield monitor data explained. The field outline is the concave hull of the 2017 yield points — 19.5 ha of it, because the authors had trimmed the borders — and the zones are the mean of the three kriged yield maps in the archive, split into thirds. The files' declared CRS is wrong; they are in Gauss-Boaga (EPSG:3003).

import os
import zipfile

import geopandas as gpd
import numpy as np
import rasterio
import requests
import shapely
from rasterio.features import shapes
from shapely.geometry import shape

if not os.path.exists("FIELD_YIELD"):
    reply = requests.get("https://zenodo.org/records/14216165/files/FIELD_YIELD.zip?download=1", timeout=300)
    reply.raise_for_status()
    open("FIELD_YIELD.zip", "wb").write(reply.content)
    zipfile.ZipFile("FIELD_YIELD.zip").extractall(".")

CRS = "EPSG:3003"
points = gpd.read_file("FIELD_YIELD/yield2017/yield2017.shp").set_crs(CRS, allow_override=True)
field = shapely.concave_hull(points.union_all(), ratio=0.05)

maps = []
for year in (2016, 2017, 2018):
    with rasterio.open(f"FIELD_YIELD/yield{year}krig/{year}yieldkrig.tif") as src:
        maps.append(src.read(1, masked=True))
        transform = src.transform
mean = np.ma.mean(maps, axis=0)
breaks = np.quantile(mean.compressed(), [1 / 3, 2 / 3])
classes = np.digitize(mean.filled(np.nan), breaks).astype("int16")
classes[mean.mask] = -1
zones = gpd.GeoDataFrame(
    [{"zone": ("low", "mid", "high")[int(v)], "geometry": shape(g)}
     for g, v in shapes(classes, mask=classes >= 0, transform=transform)], crs=CRS).dissolve("zone").reset_index()
print(f"field {field.area / 1e4:.1f} ha; zone areas, ha: "
      f"{dict(zip(zones['zone'], (zones.area / 1e4).round(1)))}")
field 19.5 ha; zone areas, ha: {'high': 6.7, 'low': 6.7, 'mid': 6.7}

Example 1 — grid, stratified and space-filling designs

import numpy as np, geopandas as gpd
from shapely.geometry import Point

def random_points_in(poly, n, rng, max_tries=200):
    minx, miny, maxx, maxy = poly.bounds
    out = []
    for _ in range(max_tries * n):
        p = Point(rng.uniform(minx, maxx), rng.uniform(miny, maxy))
        if poly.contains(p):
            out.append(p)
            if len(out) == n:
                break
    return out

def stratified_samples(field, zones, n_total, min_per_zone=4, inset_m=15, seed=0):
    rng = np.random.default_rng(seed)
    z = gpd.overlay(zones, gpd.GeoDataFrame(geometry=[field.buffer(-inset_m)],
                                            crs=zones.crs), how="intersection")
    z["area"] = z.area
    share = z["area"] / z["area"].sum()
    n = np.maximum((share * n_total).round().astype(int), min_per_zone)

    rows = []
    for (_, row), k in zip(z.iterrows(), n):
        for p in random_points_in(row.geometry, int(k), rng):
            rows.append({"zone": row.get("zone"), "geometry": p})
    g = gpd.GeoDataFrame(rows, crs=zones.crs)
    g["sample_id"] = [f"S{i + 1:03d}" for i in range(len(g))]
    print(g.groupby("zone").size().to_string())
    return g

def space_filling(field, n, crs, inset_m=15, seed=0, candidates=5000):
    """Maximin: greedily pick points that are as far as possible from those chosen."""
    rng = np.random.default_rng(seed)
    poly = field.buffer(-inset_m)
    cand = np.array([[p.x, p.y] for p in random_points_in(poly, candidates, rng)])
    chosen = [int(rng.integers(len(cand)))]
    d = np.hypot(*(cand - cand[chosen[0]]).T)
    for _ in range(n - 1):
        nxt = int(d.argmax())
        chosen.append(nxt)
        d = np.minimum(d, np.hypot(*(cand - cand[nxt]).T))
    return gpd.GeoDataFrame(
        {"sample_id": [f"S{i + 1:03d}" for i in range(n)]},
        geometry=[Point(*cand[i]) for i in chosen], crs=crs)

The field and its zones are real: the outline of the Ferrara corn field from its yield points, and three zones from the mean of its three kriged yield maps, split at the tertiles:

for spacing in (100, 75):
    grid = grid_samples(field, spacing, CRS)
    print(f"{spacing} m grid on {field.area / 1e4:.1f} ha: {len(grid)} points")
strat = stratified_samples(field, zones, n_total=len(grid))
spread = space_filling(field, len(grid), CRS)
print(f"stratified {len(strat)}, space-filling {len(spread)}")
100 m grid on 19.5 ha: 17 points
75 m grid on 19.5 ha: 29 points
zone
high    10
low      9
mid     10
stratified 29, space-filling 29

Example 2 — check the lag distribution before sampling

import numpy as np
from scipy.spatial.distance import pdist

def lag_profile(points, bins=(0, 20, 50, 100, 200, 400, 800)):
    xy = np.c_[points.geometry.x, points.geometry.y]
    d = pdist(xy)
    counts, edges = np.histogram(d, bins=bins)
    print(f"{len(points)} points, {len(d):,} pairs, "
          f"minimum separation {d.min():.1f} m")
    for lo, hi, c in zip(edges[:-1], edges[1:], counts):
        print(f"  {lo:>4.0f}–{hi:<4.0f} m: {c:5,} pairs")
    if counts[0] < 20:
        print("  ! too few short-lag pairs to fit a nugget")
    return counts, edges
for name, design in (("75 m grid", grid), ("space-filling", spread)):
    print(name)
    lag_profile(design)
75 m grid
29 points, 406 pairs, minimum separation 75.0 m
     0–20   m:     0 pairs
    20–50   m:     0 pairs
    50–100  m:    43 pairs
   100–200  m:   120 pairs
   200–400  m:   165 pairs
   400–800  m:    78 pairs
  ! too few short-lag pairs to fit a nugget
space-filling
29 points, 406 pairs, minimum separation 65.7 m
     0–20   m:     0 pairs
    20–50   m:     0 pairs
    50–100  m:    40 pairs
   100–200  m:    96 pairs
   200–400  m:   170 pairs
   400–800  m:   100 pairs
  ! too few short-lag pairs to fit a nugget

Twenty to thirty pairs in the shortest bin is a reasonable minimum for fitting a nugget. A pure 75 m grid has none: its closest points are 75 m apart.

Example 3 — add short-lag pairs cheaply

import numpy as np, geopandas as gpd
from shapely.geometry import Point

def add_close_pairs(points, n_pairs=6, distances=(10, 25), seed=0):
    rng = np.random.default_rng(seed)
    pick = rng.choice(len(points), size=min(n_pairs, len(points)), replace=False)
    extra = []
    for i in pick:
        p = points.geometry.iloc[i]
        for dist in distances:
            angle = rng.uniform(0, 2 * np.pi)
            extra.append(Point(p.x + dist * np.cos(angle), p.y + dist * np.sin(angle)))
    g = gpd.GeoDataFrame(geometry=extra, crs=points.crs)
    g["sample_id"] = [f"P{i + 1:03d}" for i in range(len(g))]
    out = gpd.pd.concat([points, g], ignore_index=True)
    print(f"{len(points)} grid points + {len(g)} close points = {len(out)}")
    return out
nested = add_close_pairs(grid, n_pairs=6)
lag_profile(nested)
from scipy.spatial.distance import pdist

for n in (12, 20, 29):
    more = add_close_pairs(grid, n_pairs=n)
    short = int((pdist(np.c_[more.geometry.x, more.geometry.y]) < 20).sum())
    print(f"  close pairs at {n} grid nodes: {short} pairs under 20 m")
29 grid points + 12 close points = 41
41 points, 820 pairs, minimum separation 10.0 m
     0–20   m:     7 pairs
    20–50   m:    11 pairs
    50–100  m:    90 pairs
   100–200  m:   223 pairs
   200–400  m:   339 pairs
   400–800  m:   150 pairs
  ! too few short-lag pairs to fit a nugget
29 grid points + 24 close points = 53
  close pairs at 12 grid nodes: 17 pairs under 20 m
29 grid points + 40 close points = 69
  close pairs at 20 grid nodes: 24 pairs under 20 m
29 grid points + 58 close points = 87
  close pairs at 29 grid nodes: 31 pairs under 20 m

On this field the cheap fix was not enough. Pairs at six grid nodes — twelve extra points, 41% more cost — gave only 7 pairs under 20 m, below the check's own threshold of 20; it took pairs at 20 nodes, 40 extra points, to reach 24. Short lags are expensive on a small field, so budget for them from the start rather than bolting a handful on, or accept a nugget estimated from very few pairs and say so.

Why a regular grid cannot fit a variogram

Why a regular grid cannot fit a variogram

A variogram is estimated from pairs binned by separation. A grid with 75 m spacing has no pairs closer than 75 m, so the model must be extrapolated to the origin to estimate the nugget — and the nugget is what determines how much the interpolation smooths. Adding pairs at 10 and 25 m constrains it directly, which is why nested and clustered designs exist.

Why the grid origin should be random

A field's features are often periodic: tramlines at the boom width, drainage runs at a fixed spacing, previous applications aligned with the same tramlines. A grid whose origin is the field corner and whose spacing happens to be a multiple of the tramline spacing samples the same position in every pass, which is the worst possible alignment. A random offset makes that a coincidence rather than a certainty.

Why to inset from the boundary

Headlands are turned on, compacted, double-fertilised and shaded by hedges. Their soil is genuinely different, and including them in a field mean biases it. If the headland is of interest, sample it as its own stratum rather than letting it contaminate the field.

Why stratification and discovery pull in opposite directions

Stratifying by existing zones gives precise estimates of those zones' means and, by construction, no information about whether the zones are right. A space-filling design spreads points to cover the field's variation without assuming a structure, which is what you want when the zones are the question. Doing both means splitting the budget.

The design decisions a soil sampling campaign has to make explicitly. decide these before samplingwhat the samples are fora mean, a map or zone meansrandomised grid originnot the field cornerinset from the boundary15–20 mshort-lag pairs if a map is wantedstable sample identifiersmatching the bag labelsdepth and compositing recordedLabels not matching the points is the commonest campaign failure spatialworkflow.io
Six decisions, all cheap to make and expensive to discover afterwards.

Narrow fields, composites and depth

  • Buffering inwards can empty a narrow field. Check for an empty polygon.
  • Composite samples of several cores per point are standard; record the pattern.
  • Depth matters. A 0–15 cm and a 15–30 cm sample are different variables.
  • Sample identifiers must match the labels. This is the commonest campaign failure.
  • Record the sampling date. Nutrient levels move through a season.
  • Keep the design file. Resampling the same points later is far more informative.
  • Access constraints — wet spots, standing crop — will move points; record where they went.
  • Report the achieved spacing, not the designed one.

FAQ

What spacing should I use for soil sampling?

Sixty to a hundred metres is common for phosphorus and potassium mapping, but the right answer depends on the variogram range, which you do not know until you have sampled.

Should the grid origin be fixed?

No. Randomise it, so the grid cannot align with tramlines, drainage runs or previous applications.

Why inset from the field boundary?

Headlands are compacted, differently fertilised and shaded. Sample them as their own stratum if they matter, not as part of the field.

Why can I not interpolate my grid samples?

Because a regular grid has no pairs closer than its spacing, so the variogram nugget cannot be estimated. Add a few pairs at 10 and 25 m.

Should I stratify by management zone?

If the zones are established and you want their means, yes. If you want to discover whether the zones are right, use a space-filling design instead.

How many samples per zone?

Proportional to area, with a minimum of about four so each zone mean has some precision.