Build management zones for a field from several years of Sentinel-2 imagery, optionally combined with yield, soil EC and terrain. Plan the soil sampling, bring the lab results back, and turn them into a variable-rate prescription. Built for Alberta and Saskatchewan.
Input a field boundary and your Copernicus credentials, and Waxwing pulls the growing-season NDVI history, works out where the field is consistently above or below its own average, and clusters that into as many zones as you ask for.
There is a full user guide in docs/Waxwing User Guide.docx —
sixteen sections covering everything from installing Python through to building a
prescription, with figures produced by the real pipeline. It is a Word document
so you can edit it for your own operation before handing it to someone.
It is generated, not hand-maintained. Every reference table in it — crop coefficients, sampling designs, colour ramps, soil zones, nutrients — is read out of the live code, so it cannot quietly drift from what the software does. To rebuild after changing a default:
python docs/build_figures.pypython docs/build_guide.pyThe first needs the usual project dependencies; the second additionally needs
python-docx.
Double-click run.bat. It finds Python, installs anything missing on first
run, and opens your browser once the server is ready.
To use a different port:
run.bat 8502Or start it by hand:
pip install -r requirements.txtstreamlit run app.pyEither way the app opens at http://localhost:8501. Leave the console window open while you work — closing it stops Waxwing.
Imagery comes from the Copernicus Data Space Ecosystem (CDSE). Registration is free, and the volume this app uses sits far inside the free monthly allowance — a season of history over a quarter section costs a fraction of a percent of it.
- Create an account at dataspace.copernicus.eu.
- Open the Sentinel Hub dashboard → User settings.
- Under OAuth clients, click Create new.
- Copy the Client ID and the Client secret.
The secret is displayed once — copy it before closing the dialog. People often call this pair an "API key"; it is two values and both are needed.
Paste them into the sidebar and hit Test connection. Optionally set them as environment variables instead:
setx CDSE_CLIENT_ID "your-client-id"
setx CDSE_CLIENT_SECRET "your-client-secret"The sidebar also offers Remember on this computer, which writes both values as
plain text to ~/.waxwing/credentials.json. Anything running under your Windows
account can read that file, so leave it off on a shared machine.
Draw the field on satellite imagery, or upload a boundary (GeoJSON, KML, GeoPackage, or a zipped shapefile). There is also an approximate Township-Range-Meridian jump to get the map near the right quarter — the Alberta Township System has correction jogs every fourth township, so it lands within a kilometre or two, then you draw against the imagery.
Two settings matter here:
- Headland trim (default 10 m) shrinks the boundary inward. A 10 m Sentinel-2 pixel sitting on the field edge is part crop and part road, ditch or shelterbelt; trimming keeps that mixture out of the clustering.
- Pixel size (default 10 m, Sentinel-2 native). Coarser smooths the result.
Pick years and a season window, then search, then download.
| Preset | Window | What it is good for |
|---|---|---|
| Peak biomass (default) | Jun 15 – Aug 15 | Closed canopy, differences show most clearly. Best single window for zoning. |
| Full growing season | May 1 – Sep 30 | Most images, but bare soil and stubble add noise. |
| Early season | May 1 – Jun 15 | Drainage and seedbed problems. Heavily soil-influenced. |
| Vegetative to early flower | Jun 1 – Jul 20 | Good compromise when peak has too few clear images. |
| Late season | Aug 15 – Sep 30 | Moisture-holding capacity, but uneven maturity. |
Winter and shoulder-season dates are never requested — there is no canopy signal under snow or on frozen ground.
Every image passes three gates before it is used:
- Catalog cloud cover — a coarse tile-level first pass.
- Clear pixels over your field (default ≥ 85 %) — a partly clouded image biases the field mean that every other image is normalised against, so partial scenes are dropped outright rather than half-used.
- Field mean NDVI (default ≥ 0.25) — rejects bare soil and harvested stubble.
Then Compute layers produces:
- Relative productivity — the multi-year average of normalised NDVI. The main zoning layer. Scenes are weighted by how well they agree with the consensus of all the others, so an undetected hazy or mis-registered date fades out instead of dragging the average. On clean imagery this changes nothing; where one date in three was spoiled it took zone agreement with the noise-free answer from 0.94 to 0.98.
- Variability — how much a spot's standing swings from year to year. Low is dependable ground; high changes rank between seasons.
- Mean NDVI and CV — for reference.
Upload a yield monitor CSV (or a point shapefile/GeoJSON). Columns are auto-detected and you can override them. Cleaning removes zeros, points outside the boundary, and extreme percentiles, then interpolates to the analysis grid by IDW.
The app reports the correlation between your yield map and the satellite productivity pattern — a weak correlation is a flag that the yield file still has pass artefacts in it.
Set weights on the layers, pick a zone count, and build. Score every zone count runs silhouette analysis across 2–8 zones. Patches below the minimum workable area are absorbed into their neighbour, because a sprayer cannot act on a two-pixel island.
Zones are numbered by productivity, so Zone 1 is always the best ground.
Method defaults to balanced breaks, which matters most for the common 3-zone case. A field with a slough, a saline patch or an eroded knoll has a tail of extreme pixels, and a plain variance-minimising split hands that tail its own zone — on benchmark quarter sections that came out at 12–16 ac of a 150 ac field, too little to justify its own rate. Balanced breaks solve the same split exactly (one dimension makes it a dynamic programme, not a heuristic) subject to a floor on zone size. A 20% floor took the smallest zone from 12 ac to 30 ac for one point of explained variance.
On a field that varies smoothly with no such patch, the floor changes nothing — the best split is already balanced.
| Method | Areas (ac, 6-field mean) | Balance | R² |
|---|---|---|---|
| k-means | 58 / 79 / 13 | 0.17 | 0.832 |
| Balanced, 20% floor — default | 54 / 67 / 29 | 0.43 | 0.821 |
| Equal-area quantiles | 50 / 50 / 50 | 0.99 | 0.812 |
Two things were tested and deliberately not changed: k-means is already finding the global optimum here (its gap to the exact 1-D dynamic programme is +0.0004%), and weighting each year equally instead of each scene made results worse, because a year with more images has a better-averaged estimate and deserves its weight.
Zone colours are chosen here, not at build time, so recolouring never re-runs the clustering. Ramps on offer:
| Ramp | Kind | Colourblind-safe |
|---|---|---|
| Green to red (agronomic) — default | diverging | no |
| Blue to red | diverging | yes |
| Single-hue green / blue | sequential | yes |
| Viridis | sequential | yes |
| Greyscale (for printing) | sequential | yes |
| Bauhaus primaries | categorical — labels zones, doesn't rank them | yes |
| Custom | one colour picker per zone | — |
Each ramp is sampled across its whole range for the zone count you asked for. That matters: taking the first three colours off an eight-step green-to-red ramp would hand back three greens and throw the ordering away. The picker warns you when a ramp isn't colourblind-safe, or when it doesn't encode order at all.
The chosen colours flow through everything — the map, the legend, the swatch
chips, and the color column in every vector export.
GeoPackage is the recommended download: one file, all layers, real CRS
metadata, and no 10-character cap on field names, so productivity_index stays
intact instead of becoming prod_idx. Layers are zones, field_boundary, and
yield_points when yield was used. Opens in QGIS, ArcGIS and most current farm
software. Shapefile is still there for older rate controllers.
One zip containing:
waxwing.gpkg All vector layers in one file (recommended)
sampling/ Field sheet, GPX waypoints, CSV and GeoJSON points
zones.shp (+.dbf/.shx/.prj) Same zones for older farm software
zones.geojson WGS84 lat/lon
zones.tif Zone raster, values 1..n, 0 = outside
productivity.tif Relative productivity
mean_ndvi.tif Mean raw NDVI
variability.tif Year-to-year swing (high = less dependable)
yield_interpolated.tif If yield was used
field_boundary.geojson The trimmed boundary
zone_statistics.csv Area and layer means per zone
scene_log.csv Every candidate image, and why it was kept or dropped
zone_map.png Printable map
productivity_map.png Printable map
metadata.json Full provenance
README.txt How this specific map was made
A zone map is a hypothesis about where the ground differs; soil sampling is how you test it. The design decides what the lab results can actually support, so all three families from the literature are here rather than just the convenient one.
| Design | Family | What it is for |
|---|---|---|
| Zone composite — default | directed | The standard precision-ag approach. Cheapest way to get a number you can vary rates on. |
| Stratified random | design-based | Unbiased per-zone means with honest confidence intervals. |
| Grid point | design-based | Best characterisation of variability; the only design that finds a pattern you didn't predict. |
| Grid cell | design-based | Fewer samples, each describing its cell rather than a point. |
| Systematic unaligned grid | design-based | Immune to aliasing against tramlines, drill passes and tile runs. |
| Spatial coverage | spatial | Most even geographic coverage for a given count (Brus, de Gruijter & Walvoort 2006). |
| Conditioned Latin hypercube | spatial | Reproduces the field's covariate distribution (Minasny & McBratney 2006). |
| Benchmark sites | directed | Fixed spots revisited yearly — the only way to separate a real trend from sampling noise. |
| Simple random | design-based | The baseline, and a fair independent check on whether the zones mean anything. |
| Transect | directed | Investigating a gradient: a salinity front, a slope sequence, a drainage run. |
Depths default to the prairie standard, 0–6″ plus 6–24″. Nitrate and sulphate move down the profile, so a surface test alone under-reads what the crop can reach. Three-way and deep splits are available for salinity or rooting work.
Cores per sample default to the published ranges — 8–10 for a grid point, 10–15 per management zone, 15–20 for a whole-field composite — and the app warns if you go below them.
How many samples? A built-in calculator solves n = (t·CV/D)² iteratively.
It also answers the inverse: at CV 40%, a 20-core composite gives about ±19%,
which is why the extension literature quotes ~±17% for nitrate. The formula
assumes independent samples; soil is spatially autocorrelated, so treat it as a
floor.
Placement rules. An edge buffer keeps samples off headlands, where compaction, overlapped fertiliser and road dust make the soil unrepresentative. If the buffer would erase a narrow field it is backed off in halves and the shortfall is reported — never silently dropped. Minimum spacing and a fixed random seed are available so a plan is reproducible.
Every plan is scored for how well it represents the field, with a two-sample Kolmogorov–Smirnov test per layer against every eligible pixel. This is what cLHS optimises directly, reported for all designs so they compare on equal terms.
Exports: a printable field sheet with tick boxes, GPX waypoints for a handheld
GPS, CSV, GeoJSON, and a sampling_points layer in the GeoPackage.
The step that closes the loop. Upload the lab report, matched to the sampling plan by sample ID (case and punctuation insensitive), and get three things.
Validation — a one-way ANOVA per nutrient answering whether the zones actually separated fertility. This is the honest test of the whole exercise, and a low number is reported as plainly as a high one. Nitrate usually scores badly: it is mobile and answers to last year's crop and this spring's weather more than to the stable pattern a satellite sees. Organic matter usually tracks the zones closely. Both come with a p-value, because with a handful of samples per zone noise alone can produce a respectable R².
Per-zone nutrient levels — what the lab found, summarised by zone.
A prescription — rates from a transparent balance:
N = target yield × lb N/bu + 30 lb residual target
− soil nitrate-N (0–24″) − pulse credit − manure credit
Phosphorus follows the Alberta response thresholds (build below 20 ppm Modified Kelowna, replace removal to 40, taper above). Potassium is skipped on soils already testing adequate. Every coefficient is editable and the formula is printed above the table. Exports a GeoJSON with N/P₂O₅/K₂O/S per zone for a rate controller.
Rates are starting points from published guidelines, not a recommendation for your field. Check with an agronomist before anything goes in the ground.
Almost nothing about the geometry differs across the prairies — the survey system, the imagery, the clustering and the sampling designs are shared. What differs is agronomy, so two things are chosen explicitly rather than inferred.
Province decides which meridians a land description can name. Alberta runs W4/W5/W6, Saskatchewan W1/W2/W3. There is no W4 land in Saskatchewan because the 4th Meridian is the border, and naming it there is refused with that explanation. It is defaulted from the field centroid but never silently assumed — the border runs through productive land.
Soil zone is the stronger axis. Nitrogen per bushel turns out to be flat across zones — Saskatchewan's own assumed rates work out to 1.87, 1.88 and 1.88 lb N/bu for canola in Brown, Dark Brown and Black — so the zone drives a prescription almost entirely through target yield. Five zones: Brown, Dark Brown, Black, Dark Gray, Gray.
Zone is deliberately not guessed from location. The obvious rule — dry southwest means Brown — is right for Saskatchewan and inverted for Alberta, whose Brown belt is the southeast Special Areas while the southwest is cool, wet foothills.
The survey grid is 486 × 483 chains, not 6 × 6 miles. Canadian road allowances are added to the 36 sections rather than carved out of them, so a township is 6.075 miles east–west and 6.0375 north–south. The old six-mile assumption under-ran by 121 m per range, compounding to kilometres. The grid also re-anchors at each base line rather than accumulating from township 1.
Measured against Red Deer, the one properly anchored landmark: 3.55 km → 1.28 km. Across six prairie landmarks the miss runs 1.3–6.4 km with the scatter dominating any trend, so the map jump quotes a flat range rather than a false formula — a city centre is not a township centre, and those checks cannot resolve the model any finer.
Versions before 1.4 prescribed roughly 60% too much nitrogen, in both provinces.
Crop.n_per_bu held total crop uptake including straw and the balance used it
as a fertiliser rate. Canola at 3.0 lb N/bu over a 50 bu target gave 150 lb N/ac,
against Saskatchewan's published assumed rate of 94 for the same target. Total
uptake counts nitrogen recycled from residue and mineralised in season, none of
which a spring nitrate test sees, so budgeting from uptake while crediting only
nitrate over-prescribes by half as much again.
The coefficients are now derived from the provincial guide's own rate and target pairs, and reproduce the published rates to 1.1% across fourteen crop-zone combinations. Uptake is retained as a display-only figure and never enters a rate.
Durum and mustard were also missing entirely, which matters in Saskatchewan.
Other indices. NDVI saturates once the canopy closes, exactly when you want to separate good ground from better. NDRE uses the red edge and keeps responding. GNDVI is more sensitive to chlorophyll and so to nitrogen. MSAVI2 suppresses the soil background early in the season.
Bare soil. Two indices read the ground instead of the crop — a Bare Soil Index and soil brightness, which is a usable organic-matter proxy on a clean bare date. Pick one and the season presets switch to the windows when a prairie field is uncovered (spring before emergence, after harvest) and the quality gate inverts: scenes with a crop growing on them are rejected. Soil properties do not change between a wet year and a dry one, which makes them a steadier zoning input than any single season of canopy.
Soil EC. Veris and EM38 exports ride the same cleaning and interpolation path as yield, with their own defaults — no zero-dropping, since a zero EC is just low conductivity rather than a combine that was not cutting.
Terrain. One DEM request gives slope, relative elevation and a topographic wetness index built on real D8 flow accumulation. Relative elevation is usually the most useful on prairie ground: it removes the regional tilt and leaves the knolls and hollows, which is what actually moves water. The Copernicus DEM is 30 m, so it finds a slough but not a dead furrow — Alberta's LiDAR is better if you have it.
Zone the odd years and the even years separately and see whether you get the same map. Agreement is reported as Cohen's kappa, which discounts the overlap you would get by chance from the zone sizes alone.
That distinction is not academic. On a test field of pure noise, raw pixel overlap reads 96.6% — which looks like a resounding success — while kappa reads −0.017, correctly calling it worthless. Raw overlap is the number that misleads, so kappa is the headline.
A leave-one-year-out pass then shows how far the map moves when each year is dropped. A year that changes everything is worth asking about: hail, a flooded spring, a different crop, or simply fewer clear images.
Both run on imagery already downloaded, so neither costs Copernicus quota.
Boundaries and the settings that produced a run are saved under data/fields/
as GeoPackages that open in QGIS, with a JSON sidecar for the settings. A farm is
not one quarter, and redrawing a boundary every time is the difference between a
demo and something used in August.
Field data stays on your machine. data/ is git-ignored in full.
A single NDVI image tells you about that day's weather, crop stage and sensor geometry as much as it tells you about the ground. What is stable across many dates and several years is the relative productivity pattern — which parts of the field consistently run above or below the field average.
So each scene is rescaled within the field before anything is averaged. Without that, a wet year with a heavy canopy dominates the average and a dry year barely registers, even though both describe the same underlying ground.
Z-score within field is the default. Ratio to field mean is more
interpretable (1.0 = field average), Percentile rank is the most robust to
outliers, and Raw NDVI is available for comparison.
Variation explained on the results page is the number to judge a zone map by. It is the share of the field's productivity variation the zones account for. Above 50 % is a well-separated field. Below 30 % means the ground varies as a smooth gradient rather than in blocks, and the zone boundaries are soft lines through that gradient rather than real agronomic edges. Most fields are somewhere in between.
If the silhouette score is low at every zone count, the field has no natural breaks. That is not a failure — it means the zone count should come from what your equipment and inputs can actually act on, not from the metric.
Productivity index in the statistics table is a zone's mean NDVI as a percentage of the field mean. 100 is average ground; 110 means the canopy there ran about 10 % denser than field average across the years analysed. It is a relative index, not a yield prediction.
These zones show where the field has been consistently different. They do not show why. A low zone could be thin topsoil, salinity, compaction, poor drainage, or the spot where the old yard site sat.
Soil sample by zone before varying rates. That is the step that turns a productivity map into an agronomic decision.
The interface is styled after the Bauhaus: warm paper, near-black rules, the three primaries as flat blocks, and one geometric typeface (Jost, a Futura descendant) carrying everything through weight and tracking rather than a second family. Step markers use the square, circle and triangle on Kandinsky's 1923 Bauhaus shape-to-colour assignment — triangle yellow, square red, circle blue.
Nothing is rounded except deliberate circles. baseRadius is zero.
The rendered maps are not themed. The NDVI, productivity, variability and yield layers keep their scientific colour ramps. Zone colours are yours to choose in step 5 — see the ramp table above — and default to the familiar agronomic green-to-red. The theme stops at the chrome.
Colour tokens live in waxwing/theme.py; the structural theme is
.streamlit/config.toml. Text contrast was measured, not eyeballed — body text
clears 4.5:1 on paper, and red as text uses a darker step (#C41E1E, 5.14:1)
because the fill red #D62828 lands at 4.36:1, just under the requirement.
To go back to stock Streamlit styling, delete .streamlit/config.toml and remove
the theme.inject() call near the top of app.py.
run.bat Windows launcher (deps, port check, browser)
app.py Streamlit UI — the seven steps
.streamlit/config.toml Bauhaus structural theme
docs/
Waxwing User Guide.docx The full user guide, editable
build_guide.py Rebuilds the guide from the live code
build_figures.py Rebuilds the guide's figures from the real pipeline
images/ Generated figures
waxwing/
config.py Endpoints, season presets, SCL classes, defaults
theme.py Bauhaus stylesheet, masthead, step headings
palettes.py Zone colour ramps and custom-colour handling
regional.py Provinces, soil zones, survey grid, target yields
sampling.py Soil sampling designs, sample-size stats, GPX/field sheet
lab.py Lab results, zone validation, VRA prescriptions
terrain.py DEM derivatives: slope, relative elevation, wetness
robustness.py Split-half and leave-one-year-out zone stability
library.py Saved field boundaries and their settings
boundary.py Boundary load/clean, UTM grid, headland trim, ATS helper
copernicus.py CDSE OAuth, Catalog search, Process API NDVI fetch
ndvi.py Scene selection, stack building, normalisation, temporal layers
yielddata.py Yield ingest, cleaning, IDW interpolation
zones.py Feature matrix, clustering, smoothing, polygonize, statistics
exports.py GeoTIFF/shapefile/GeoJSON/PNG/bundle
credentials.py Optional local credential storage
tests/ Four hand-run suites: maths, pipeline, UI, contrast
data/cache/ Downloaded NDVI, cached per field+date (never committed)
python tests/test_pipeline.pyRuns the whole analysis path — boundary → grid → synthetic NDVI stack → temporal layers → yield → clustering → every export format — with no API calls and no credentials. Uses synthetic imagery with a known productivity pattern, then checks the pipeline recovers it: acreages reconcile, Zone 1 is the most productive, the constructed low-productivity blob lands in the worst zone, and every file in the bundle is present and valid.
Downloaded NDVI is cached under data/cache keyed by field geometry, grid and date.
Re-running with different normalisation, weights or zone counts costs no additional
Copernicus quota. Delete the folder to force a fresh download.
- Sentinel-2 L2A on CDSE goes back to 2017. Four to six years is plenty; more years average over more crops and more weather, which is what makes the pattern dependable.
- Works anywhere Sentinel-2 covers. Outside Alberta, check that a season preset suits your growing season.
- The 10 m pixel is the real resolution limit. A 160-acre quarter is about 6,000 pixels, which is ample for clustering.
Imagery contains modified Copernicus Sentinel data, processed via the Copernicus Data Space Ecosystem.