C12 · Crosswinds at Santa Barbara — Bayesian Modelling in Practice

The brief

A flight school at Santa Barbara (SBA) sends student pilots on solo flights. Its safety rules for a student solo, as far as wind goes:

"The crosswind component on the runway in use must not exceed 8 knots. Calm is fine. Students only use runway 7/25 (the long runway); the tower gives them the end that faces into the wind, so tailwind is not an issue. First-solo students have a stricter 6-knot limit."

The chief instructor asks:

"By month and hour of the day, what is the probability that a solo slot is flyable? We want a default two-hour solo window between 08:00 and 18:00 local time. Which one should it be, how sure are you, and how many bookings a month should we expect to cancel because of the wind? Does the answer change for first-solo students?"

You have three years of hourly wind reports from the airport's own weather station. The obvious approach - fit a model for the wind speed, fit a model for the wind direction, and combine them - gives a confident, plausible, wrong answer. Find out why, fix it, and then check your probabilities on a year the model has not seen.

This is the applied sequel to example E81 (circular data). Do E81 first: this challenge assumes you know the mean resultant vector, the canonical (natural-parameter) von Mises link, circular predictors as harmonics, and what the Santa Barbara wind does over a day. E81 modelled the direction; here the question is about a vector - speed and direction together - and about a small tail probability.

How this notebook works

Each task states what to deliver, not how. Write your code in the YOUR CODE HERE cells.
Stuck? h.hint("task2") reveals hints one level at a time: nudge → approach → code skeleton. Try to get by on nudges.
h.check("task2", ...) compares your numbers with the reference solution.
A full worked solution lives in notebooks/solutions/. Open it only when you are done (or truly stuck).

import logging
import time

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import pytensor.tensor as pt
from scipy import special, stats

from pymc_challenges import Hints, data

RANDOM_SEED = 1932
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")
logging.getLogger("pymc").setLevel(logging.WARNING)  # several fits: no sampler banner per fit
pd.set_option("display.width", 160)
pd.set_option("display.max_columns", 30)
T_START = time.time()

h = Hints("C12")
h.tasks()

C12 - Crosswinds at Santa Barbara: when can a student fly solo?

task0 - Look at the data (2 hints)
task1 - The crosswind, and what the raw data can say (3 hints, 2 check(s))
task2 - The obvious model, speed and direction separately (3 hints, 2 check(s))
task3 - Let the direction depend on the speed (3 hints, 2 check(s))
task4 - Is the model honest about 2023? (3 hints, 3 check(s))
task5 - The schedule (3 hints, 4 check(s))

The data

Routine hourly reports (METARs), issued at minute 53 of each hour. drct is the direction the wind blows from, in degrees clockwise from true north, reported in steps of 10; sknt is the speed in knots. A calm is reported as direction 0 and speed 0. An empty direction with a non-zero speed is a variable wind (no steady direction).

The school books slots in local clock time, so we convert the UTC timestamps to America/Los_Angeles (daylight-saving time included). We fit on 2021-2022 and keep 2023 aside: nothing you do before Task 4 may look at it.

data.describe("ksba_wind")
raw = data.load("ksba_wind")
local = pd.to_datetime(raw.valid).dt.tz_localize("UTC").dt.tz_convert("America/Los_Angeles")
wind = raw.assign(hour=local.dt.hour, month=local.dt.month, year=local.dt.year, date=local.dt.date)
wind = wind[(wind.year >= 2021) & wind.sknt.notna()].reset_index(drop=True)   # 8 reports fall on 31 Dec 2020 local
train = wind[wind.year < 2023].reset_index(drop=True)
test = wind[wind.year == 2023].reset_index(drop=True)
print(f"training reports (2021-2022): {len(train):,}; held-out 2023: {len(test):,}")
train.head()

ksba_wind: Routine hourly METAR reports, 2021-01-01 to 2023-12-31 (26,101 rows, at minute 53 of each hour): valid (UTC timestamp), drct (direction the wind blows FROM, degrees clockwise from true north, reported in steps of 10; 0 together with sknt = 0 means calm), sknt (speed, knots), tmpf and dwpf (air and dew-point temperature, Fahrenheit). Empty = missing (e.g. variable direction).
  source: Iowa Environmental Mesonet (IEM, Iowa State University) ASOS/METAR archive, station SBA (Santa Barbara Municipal Airport, California; an FAA/NWS Automated Surface Observing System). Public US government observations redistributed by IEM, which asks to be acknowledged. BUILT by tools/build_e81_circular.py (the URL returns an uncompressed CSV) - keep the cached file.
  url:    https://mesonet.agron.iastate.edu/cgi-bin/request/asos.py?station=SBA&data=drct&data=sknt&data=tmpf&data=dwpf&year1=2021&month1=1&day1=1&year2=2024&month2=1&day2=1&tz=Etc/UTC&format=onlycomma&latlon=no&missing=empty&trace=empty&direct=no&report_type=3
training reports (2021-2022): 17,354; held-out 2023: 8,677

	valid	drct	sknt	tmpf	dwpf	hour	month	year	date
0	2021-01-01 08:53	0.0	0.0	47.0	32.0	0	1	2021	2021-01-01
1	2021-01-01 09:53	0.0	0.0	45.0	33.0	1	1	2021	2021-01-01
2	2021-01-01 10:53	0.0	0.0	41.0	32.0	2	1	2021	2021-01-01
3	2021-01-01 11:53	70.0	3.0	42.0	33.0	3	1	2021	2021-01-01
4	2021-01-01 12:53	0.0	0.0	39.0	32.0	4	1	2021	2021-01-01

Runway 7/25, from the FAA airport record (form 5010, as listed e.g. on AirNav): 6,052 ft long, magnetic headings 075° / 255°, true headings 089° / 269°. (Runways are named after their magnetic heading rounded to 10°; the local magnetic variation is about 12-14° east.)

RWY_MAGNETIC = 75.0
RWY_TRUE = 89.0
LIMIT = 8.0          # knots of crosswind, solo students
LIMIT_FIRST = 6.0    # knots of crosswind, first-solo students
DAY_HOURS = range(8, 18)   # the school is open 08:00-18:00

Task 0 · Look at the data

Deliver

How many reports are calm, how many have no direction, and how fast is the wind when the direction is missing? What are the possible values of drct and sknt?
Rose diagrams of the afternoon wind (12:00-16:00), split by speed class (light, moderate, strong). Does the direction depend on the speed?

Your turn Work this task in the notebook.

Task 1 · The crosswind, and what the raw data can say

Deliver

For every report, the crosswind component on runway 7/25 and whether the slot is flyable under the 8-knot rule. Decide what calms and missing directions mean. Why does the choice of runway end not matter for the crosswind?
The share of 2021-2022 reports that exceed the limit, with the true heading - and with the magnetic heading 075° used by mistake. Which error does the mistake make, and how big is it?
A month x hour table of the observed exceedance rate. Is it good enough to schedule from?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Split the wind vector into a component along the runway and one across it. Which reference direction do METAR reports use, and which does the runway's name use? For calms and for variable winds, ask what the largest possible crosswind is.

2 Approach

Crosswind = speed x |sin(direction - runway heading)|, with both angles TRUE. The other runway end adds 180 degrees, which only flips the sign of the sine. Calm: 0. Variable (missing direction, speed > 0): use the speed itself as a worst case - and notice that it never matters at these speeds. Then pivot_table(index="month", columns="hour", values="exceed") and print the number of reports and of exceedances per cell.

3 Skeleton

def crosswind(df, heading):
    xw = df.sknt * np.abs(np.sin(np.deg2rad(df.drct - ___)))
    return xw.where(df.drct.notna(), ___)        # calm gives 0 automatically

train["xw"] = crosswind(train, 89.0)
train["exceed"] = train.xw > ___
exceed_true = train.exceed.mean()
exceed_magnetic = (crosswind(train, ___) > 8).mean()

Self-checks in the notebook (h.check("task1", …)):

exceed_true — share of ALL 2021-2022 reports with more than 8 kt of crosswind on the TRUE runway heading 089, about 0.024
exceed_magnetic — the same share computed (wrongly) with the magnetic heading 075, about 0.034

Task 2 · The obvious model: speed and direction separately

A crosswind exceedance needs a strong wind and the wrong direction. The natural plan is to model each and combine them.

Deliver

A model for the wind speed as a smooth function of hour and month that treats the calms honestly (they are speeds below a threshold, not a separate kind of wind) and respects that speeds are reported in whole knots. Check it on the share of calms and of strong winds by hour.
A model for the reported direction as a smooth function of hour and month. Make it flexible enough to have more than one preferred direction (E81 showed that one von Mises is not).
Combine them into the probability of exceeding 8 knots of crosswind for every month x hour, and compare it, hour by hour, with the observed rate. What goes wrong, and why?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

A calm is "speed below about 3 knots", and a report of 7 knots is "between 6.5 and 7.5". What kind of likelihood describes an observation that is only known to lie in an interval? For the direction: 36 reported values, more than one lobe at some hours. And when every observation in a group shares the same predictors, what do you actually need to keep from the group?

2 Approach

Speed: log S ~ Normal(m, sigma) with m and log sigma linear in harmonic features of hour (several daily harmonics) and month; the likelihood of a report s >= 3 is Phi((log(s + 0.5) - m) / sigma) - Phi((log(s - 0.5) - m) / sigma), of a calm Phi((log 2.5 - m) / sigma). Use pm.logcdf(pm.Normal.dist(), z) and pm.math.logdiffexp, weighted by the count of reports per (month, hour, speed).

Direction: a categorical over the 36 bins with logits sum_k a_k cos(k theta) + b_k sin(k theta), each a_k, b_k linear in the features (K = 1 is E81's canonical von Mises; take K = 3). The log-likelihood of a cell is sum_j count_j * logit_j - N * logsumexp(logits): concave, cheap, and it handles the 10-degree rounding exactly.

Combine per (month, hour) and posterior draw: P(exceed) = sum over speeds and bins of P(speed) P(bin) 1[speed |sin(bin - 89)| > 8]. Compare with the observed rate per hour.

3 Skeleton

with pm.Model() as speed_model:
    a = pm.Normal("a", 0, ___, shape=P)
    b = pm.Normal("b", 0, ___, shape=P)
    m = 1.5 + pt.dot(F_sp, a)
    sigma = pt.exp(pt.dot(F_sp, b))
    lcdf_hi = pm.logcdf(pm.Normal.dist(), (log_hi - m) / sigma)
    lcdf_lo = pm.logcdf(pm.Normal.dist(), (log_lo - m) / sigma)
    ll = pt.switch(s_rep == 0, ___, pm.math.logdiffexp(lcdf_hi, lcdf_lo))
    pm.Potential("ll", pt.sum(counts * ll))

basis = np.concatenate([np.stack([np.cos(k * TH), np.sin(k * TH)]) for k in range(1, K + 1)])
with pm.Model() as direction_model:
    W = pm.Normal("W", 0, ___, shape=(P, 2 * K))
    logits = pt.dot(pt.dot(F_cells, W), basis)                    # (cells, 36)
    pm.Potential("ll", pt.sum(Y * logits) - pt.sum(___ * pt.logsumexp(logits, axis=1)))

exc = np.einsum("csv,csd,vd->cs", p_speed, p_dir, (XW_GRID > 8).astype(float))

Self-checks in the notebook (h.check("task2", …)):

calm_0700 — predicted share of calm reports at 07:53, averaged over months, about 0.49 (observed 0.57)
exceed_afternoon_indep — exceedance rate the INDEPENDENT model predicts for 11:00-15:59, averaged over months, about 0.10 (observed about 0.03)

Task 3 · Let the direction depend on the speed

Deliver

A direction model that conditions on the reported speed as well as on hour and month. Try the canonical von Mises of E81 first, then a more flexible circular family. Combined with your speed model, this is a model of the full wind vector.
The exceedance check of Task 2, by hour and by month, for both. Which one can you use?
Show where the inadequate one goes wrong: compare the predicted and observed directions of strong afternoon winds.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Write the joint distribution of speed and direction as p(speed) x p(direction | speed). Which conditional did Task 0 show you? Then ask what a crosswind limit is a question about: the bulk of the direction distribution, or its shoulders at high speed?

2 Approach

Add the standardised log speed, its square and their products with the leading hour and month harmonics to the direction model's features, and refit on counts per (month, hour, speed, bin). Fit it with K = 1 (a von Mises) and with K = 3. For prediction, compute per cell and speed q = P(crosswind > 8 | speed) from the direction probabilities, then sum P(speed) q over speeds. For the diagnosis, overlay the predicted rose for strong (13 kt and more) afternoon winds on the observed one.

3 Skeleton

def features_speed(hour, month, speed):
    F = harmonics(hour, month)
    zs = (np.log(speed) - ___) / ___
    return np.column_stack([F, zs[:, None] * F[:, :7], zs ** 2, (zs ** 2)[:, None] * F[:, 1:3]])

q = np.zeros((n_cells, n_draws, 41))
for v in range(9, 41):                                     # 8 kt or less can never exceed
    p = dir_probs(features_speed(cells.hour, cells.month, np.full(n_cells, v)), W, K)
    q[:, :, v] = (p * (XW_GRID[v] > 8)).sum(___)
exc = (p_speed * q).sum(-1)

Self-checks in the notebook (h.check("task3", …)):

exceed_afternoon_vm — 11:00-15:59 exceedance rate with a single von Mises conditional on speed, about 0.10 (still far too high)
exceed_afternoon_flex — the same with the flexible (K = 3 log-Fourier) direction model conditional on speed, about 0.03

Task 4 · Is the model honest about 2023?

The posterior bands in Task 3 are narrow. But the model treats 17,000 hourly reports as independent pieces of evidence, and wind is not like that.

Deliver

For each month of 2023, the observed share of daytime (08:00-17:59) reports above the limit, and your model's 90% predictive interval for it. How many of the 12 months fall inside?
Evidence from 2021-2022 that bad hours are clustered (compare what the model expects with what happened at the level of days).
A model extension that represents the clustering, refitted on 2021-2022 only, and the same 2023 check.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

A predictive interval for a month's share has to include the noise of the reports, not just the uncertainty in the probabilities. Then compare, in 2021-2022, how many DAYS have at least one bad hour with what independent hours would give.

2 Approach

Simulate each 2023 daytime report as a Bernoulli draw with its cell probability, per posterior draw, and take the monthly share; count months inside the 5-95% range. For the clustering, expected share of days with at least one exceedance = mean over days of 1 - prod(1 - p_h). The extension: a Bernoulli regression on the daytime reports with the climatological logit as an offset, an intercept, and a random effect per day (non-centred). Simulate a new day effect per day of 2023 for the check.

3 Skeleton

with pm.Model() as day_layer:
    c = pm.Normal("c", 0, 1.5)
    sigma_day = pm.HalfNormal("sigma_day", ___)
    u = pm.Normal("u", 0, 1, shape=n_days)
    pm.Bernoulli("y", logit_p=offset + c + ___ * u[day_code], observed=exceed)
    idata_day = pm.sample(var_names=["c", "sigma_day"])

# 2023, one month: lc (draws, reports) = climatological logits + c
u = rng.normal(size=(n_draws, n_days_in_month)) * sigma_day[:, None]
share = (rng.random(lc.shape) < special.expit(lc + u[:, ___])).mean(1)

Self-checks in the notebook (h.check("task4", …)):

covered_iid — number of 2023 months (of 12) whose daytime exceedance share lies inside the i.i.d. model's 90% predictive interval, about 6
covered_day — the same count with day-level random effects, 11 or 12
sigma_day — posterior mean of the day-effect sd on the logit scale, about 2.2

Task 5 · The schedule

A solo slot is a two-hour window $[h, h+2)$; call it flyable if both reports made inside it (at $h$:53 and $h+1$:53) are within the limit.

Deliver

For every month and every window starting between 08:00 and 16:00, the probability that the window is flyable under the 8-knot rule, from your final model (with its uncertainty).
The recommended default window for the whole year: its average probability, the posterior probability that it is the best window, and how much better it is than the runner-up and the worst. Would a month-specific window be worth the complication?
For a student with 20 solo bookings in March in the default window, the expected number cancelled because of the wind, with a 90% predictive interval.
The same window ranking for first-solo students (6-knot limit). Does the advice change?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Two hours of the same day share their weather. With the day effect, what is the probability that both reports of a window are within the limit? Then think about what the instructor can act on: a ranking, the probability that the ranking is right, and the size of the difference.

2 Approach

Given the day effect u, P(window) = (1 - p1(u)) (1 - p2(u)); average over u ~ N(0, 1) with Gauss-Hermite nodes (np.polynomial.hermite_e.hermegauss, weights normalised), for each posterior draw, month and window start 8..16. Average over months for the yearly figure; P(best) = share of draws in which a window has the highest average. Cancellations of 20 bookings: Binomial(20, 1 - P(window)) per draw. For 6 knots, recompute q with the new limit from the same posterior and refit the small day layer.

3 Skeleton

x, w = np.polynomial.hermite_e.hermegauss(30)
w = w / w.sum()
for mo in range(1, 13):
    for j, h0 in enumerate(range(8, 17)):
        l1 = logit_clim[cell(mo, h0)][:, None] + c[:, None] + sigma_day[:, None] * x
        l2 = ___
        pw[mo - 1, j] = ((1 - special.expit(l1)) * (1 - special.expit(l2)) * w).sum(___)
year = pw.mean(0)                                  # (windows, draws)
p_best = np.bincount(year.argmax(0), minlength=9) / year.shape[1]

Self-checks in the notebook (h.check("task5", …)):

p_default — year-average probability that the recommended window (08:00-10:00) is flyable under the 8-kt rule, about 0.965
p_best_default — posterior probability that 08:00-10:00 is the best of the nine windows, about 1
cancel_march — expected cancellations of 20 March bookings in that window, about 1.1
p_default_first — year-average probability that 08:00-10:00 is flyable under the 6-kt first-solo rule, about 0.84

Going further

Gusts. The METAR also reports gusts (not in this extract). Download them from the IEM archive, add a gust model conditional on the mean speed and direction, and redo Task 5 with a gust crosswind limit. How much of the morning advantage survives?
Persistence. Replace the day effect with a latent AR(1) process over hours (or an HMM over weather regimes) and check whether "the next two hours" can be forecast from the last report - a much more useful product for a dispatcher than climatology.
Other runway. Runway 15/33 (true headings 166/346) is shorter but nearly across 7/25. Under a rule "either runway", how much flyable time does the school gain, and which wind regimes does the second runway rescue?
Continuous direction. Replace the 36-bin categorical with a continuous log-Fourier density plus a rounding model (E81 handled rounding by jittering). Do the tail probabilities change?

h.progress()

C12: 0/17 hints used, 13/13 checks passed.

Reference solution

A complete, executed solution with commentary on every task — including the ways the obvious model fails. Try the hint ladder first.

Open the worked solution →

#The brief

#How this notebook works

#The data

#Task 0 · Look at the data

#Task 1 · The crosswind, and what the raw data can say

#Task 2 · The obvious model: speed and direction separately

#Task 3 · Let the direction depend on the speed

#Task 4 · Is the model honest about 2023?

#Task 5 · The schedule

#Going further

Reference solution

The brief

How this notebook works

The data

Task 0 · Look at the data

Task 1 · The crosswind, and what the raw data can say

Task 2 · The obvious model: speed and direction separately

Task 3 · Let the direction depend on the speed

Task 4 · Is the model honest about 2023?

Task 5 · The schedule

Going further