How to Build a Flight Delay Predictor in Python with Real Flight and Weather Data

Most flight delay prediction tutorials report accuracy in the 80s or 90s. Most of them are measuring something other than what they claim — usually a train/test split that lets the model see the answer, or input columns that only exist after the flight has already left.

This tutorial builds a delay predictor on real data and reports honest numbers, which are more modest than you might expect. That is the useful result: it shows how hard the problem is, which features actually carry signal, and which single addition moves the model more than everything else combined.

Everything below was run end to end on real data before publishing. The numbers are what the code printed.

What we're building

A classifier that predicts, two hours before scheduled departure, whether a flight will leave the gate 15 or more minutes late. The two-hour lead matters. A predictor is only useful if it runs before the event, so every feature has to be something you'd know at that moment.

Flight data: the US Bureau of Transportation Statistics Reporting Carrier On-Time Performance dataset — every domestic flight by a reporting carrier, with scheduled and actual times. We use January 2024: 547,271 flights.

Weather data: hourly METARs from the Iowa Environmental Mesonet ASOS archive, which is free and scriptable.

Scope: departures from four big hubs — Chicago O'Hare, Atlanta, Dallas/Fort Worth and Denver. That's 89,896 flown departures, 26.5% of which left 15+ minutes late.

You need Python 3.10+, pandas, numpy and scikit-learn.

Step 0: Get the data

The BTS publishes monthly zips at a stable URL:

curl -O "https://transtats.bts.gov/PREZIP/On_Time_Reporting_Carrier_On_Time_Performance_1987_present_2024_1.zip"
unzip On_Time_Reporting_Carrier_On_Time_Performance_1987_present_2024_1.zip -d bts/

The unzipped CSV is about 250 MB and has 109 columns. We'll load ten of them.

For weather, the IEM service returns one CSV per station. The report_type=3 flag keeps routine hourly observations and tz=Etc/UTC returns timestamps in UTC, which matters in a moment:

mkdir -p wx
for st in ORD ATL DFW DEN; do
  curl -o "wx/$st.csv" "https://mesonet.agron.iastate.edu/cgi-bin/request/asos.py?\
station=$st&data=tmpf&data=sknt&data=gust&data=vsby&data=p01i&data=wxcodes\
&year1=2024&month1=1&day1=1&year2=2024&month2=2&day2=1\
&tz=Etc%2FUTC&format=onlycomma&latlon=no&missing=M&trace=T&report_type=3"
  sleep 3   # be polite to a free service
done

Step 1: Load flights

import glob
import numpy as np
import pandas as pd

HUBS = {"ORD": "America/Chicago", "ATL": "America/New_York",
        "DFW": "America/Chicago", "DEN": "America/Denver"}

cols = ["FlightDate", "DayOfWeek", "Reporting_Airline", "Origin", "Dest",
        "CRSDepTime", "Distance", "DepDel15", "Cancelled", "Diverted"]
flights = pd.read_csv(glob.glob("bts/*.csv")[0], usecols=cols)

flights = flights[flights["Origin"].isin(HUBS)]
flights = flights[(flights["Cancelled"] == 0) & (flights["Diverted"] == 0)]
flights = flights.dropna(subset=["DepDel15"]).copy()
flights["DepDel15"] = flights["DepDel15"].astype(int)

DepDel15 is the target — BTS sets it to 1 when the gate departure was 15 or more minutes late. We drop cancelled and diverted flights because "was it late" isn't a meaningful question for a flight that didn't depart normally. That's a real modelling choice, though. A production system would probably want to predict cancellations too.

Notice what we did not load: DepTime, DepDelay, TaxiOut, ArrDelay, and the delay-cause columns like LateAircraftDelay and WeatherDelay. Every one of those is only known after the flight has departed. Feeding them to a model is the most common way delay tutorials end up with suspiciously high accuracy, because the model learns to read the answer.

Step 2: Convert scheduled departures to UTC

This step is where tutorials most often go quietly wrong.

CRSDepTime is the scheduled departure in local time, stored as an integer like 1425. The weather timestamps are UTC. If you join them without converting, every Chicago flight gets matched to weather from six hours away. Nothing errors, the join succeeds, the model trains, and it learns from misaligned data.

hhmm = flights["CRSDepTime"].astype(int)
hours = (hhmm // 100) % 24          # BTS writes midnight as 2400
minutes = hhmm % 100
local_naive = (pd.to_datetime(flights["FlightDate"])
               + pd.to_timedelta(hours, unit="h")
               + pd.to_timedelta(minutes, unit="m"))

parts = []
for origin, tz in HUBS.items():
    mask = flights["Origin"] == origin
    parts.append(local_naive[mask]
                 .dt.tz_localize(tz, ambiguous="NaT", nonexistent="NaT")
                 .dt.tz_convert("UTC"))
flights["dep_utc"] = pd.concat(parts)
flights["dep_hour_local"] = hours
flights = flights.dropna(subset=["dep_utc"])

The ambiguous and nonexistent arguments handle the daylight-saving transitions, which don't occur in January but will bite the moment you extend this to a full year. Keep the local hour as its own feature too. Delays build through the day, and that pattern follows local clocks, not UTC.

Step 3: Join weather known two hours ahead

The tempting approach is to join each flight to the weather observed at its departure time. Don't. Two hours before departure you don't know that yet, so a model trained on it is being handed information it will never have in production.

Instead, join each flight to the most recent observation at or before its prediction time. merge_asof does exactly this:

def load_metars(path):
    wx = pd.read_csv(path, na_values=["M"])
    wx["valid"] = pd.to_datetime(wx["valid"], utc=True)
    wx["p01i"] = pd.to_numeric(wx["p01i"].replace("T", 0.001), errors="coerce")
    codes = wx["wxcodes"].fillna("")
    wx["snow"] = codes.str.contains("SN").astype(int)
    wx["thunder"] = codes.str.contains("TS").astype(int)
    wx["fog"] = codes.str.contains("FG").astype(int)
    wx["gust"] = wx["gust"].fillna(wx["sknt"])     # no gust reported = steady wind
    return wx[["station", "valid", "tmpf", "sknt", "gust", "vsby",
               "p01i", "snow", "thunder", "fog"]]

wx = pd.concat([load_metars(p) for p in glob.glob("wx/*.csv")]).sort_values("valid")

LEAD = pd.Timedelta(hours=2)
flights["wx_asof"] = flights["dep_utc"] - LEAD
flights = flights.sort_values("wx_asof")
flights = pd.merge_asof(flights, wx, left_on="wx_asof", right_on="valid",
                        left_by="Origin", right_by="station", direction="backward",
                        tolerance=pd.Timedelta(hours=3))
flights = flights.dropna(subset=["tmpf", "sknt", "vsby"])

direction="backward" means "the latest observation at or before this moment". The tolerance stops a flight from being matched to a stale observation if the station had a reporting gap. T in the precipitation column means a trace amount, and we treat it as a small positive number rather than zero.

The honest version of this, in production, would use the forecast for departure time rather than the last observation, because a TAF is what you actually know about conditions two hours out. The distinction is covered in METAR vs TAF. For a tutorial, the last observation is a reasonable and leak-free stand-in.

Step 3b: The feature that matters most

Here's the addition that turns out to beat everything else: how the airport is running right now.

For each flight, compute the share of departures from the same airport that left late, among flights scheduled between five and two-and-a-half hours earlier. Those flights have already gone, so their outcome is known by prediction time. The feature uses no future information.

def recent_delay_rate(group, start=pd.Timedelta(hours=5), end=pd.Timedelta(hours=2.5)):
    group = group.sort_values("dep_utc")
    t = group["dep_utc"].to_numpy()
    cum = np.concatenate([[0], np.cumsum(group["DepDel15"].to_numpy())])
    lo = np.searchsorted(t, (group["dep_utc"] - start).to_numpy(), side="left")
    hi = np.searchsorted(t, (group["dep_utc"] - end).to_numpy(), side="right")
    n = hi - lo
    rate = np.where(n > 0, (cum[hi] - cum[lo]) / np.maximum(n, 1), np.nan)
    return pd.Series(rate, index=group.index)

flights["origin_recent_delay_rate"] = (flights.groupby("Origin", group_keys=False)
                                       .apply(recent_delay_rate))

The cumulative-sum-and-searchsorted pattern computes a sliding window over ninety thousand rows in well under a second, where a row-by-row loop would take minutes. Flights early in the morning have no earlier flights in the window and get NaN, which the gradient boosting model handles natively.

Step 4: Split by time, not at random

flights["FlightDate"] = pd.to_datetime(flights["FlightDate"])
cutoff = pd.Timestamp("2024-01-24")
train = flights[flights["FlightDate"] < cutoff]
test = flights[flights["FlightDate"] >= cutoff]

This gives 66,360 training flights from January 1–23 and 23,491 test flights from January 24–31.

Why not train_test_split? Because flights on the same day share weather, congestion and knock-on delays. A random split puts some of Tuesday's flights in training and the rest in test, and the model effectively memorises that Tuesday was bad. We measured it on this data: the same model scores 0.743 ROC AUC on a random split and 0.661 on the time-based one. The random-split number is not a better model. It's leakage, and it would disappear the moment you deployed.

There's a second lesson hiding in the split. The training weeks include the mid-January Arctic outbreak and run at 29.7% delayed, while the test week is calmer at 17.3%. The model is being evaluated on conditions different from those it learned on, which is exactly what happens in production and is the right way to test.

Step 5: Train two models

from sklearn.compose import ColumnTransformer
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler

CATEGORICAL = ["Reporting_Airline", "Origin", "Dest"]
NUMERIC = ["dep_hour_local", "DayOfWeek", "Distance",
           "tmpf", "sknt", "gust", "vsby", "p01i", "snow", "thunder", "fog",
           "origin_recent_delay_rate"]

X_train, y_train = train[CATEGORICAL + NUMERIC], train["DepDel15"]
X_test, y_test = test[CATEGORICAL + NUMERIC], test["DepDel15"]

# Baseline: logistic regression
logit = Pipeline([
    ("prep", ColumnTransformer([
        ("cat", OneHotEncoder(handle_unknown="ignore", min_frequency=20), CATEGORICAL),
        ("num", StandardScaler(), NUMERIC),
    ])),
    ("model", LogisticRegression(max_iter=2000)),
])
logit.fit(X_train.fillna(0), y_train)
print("logistic regression AUC:", roc_auc_score(y_test, logit.predict_proba(X_test.fillna(0))[:, 1]))

# Gradient boosting, with native categorical support
Xg_train, Xg_test = X_train.copy(), X_test.copy()
for c in CATEGORICAL:
    cats = pd.CategoricalDtype(sorted(Xg_train[c].unique()))
    Xg_train[c] = Xg_train[c].astype(cats)
    Xg_test[c] = Xg_test[c].astype(cats)

gbm = HistGradientBoostingClassifier(categorical_features="from_dtype",
                                     max_iter=300, learning_rate=0.05, random_state=42)
gbm.fit(Xg_train, y_train)
p_gbm = gbm.predict_proba(Xg_test)[:, 1]
print("gradient boosting AUC:", roc_auc_score(y_test, p_gbm))

Always fit a simple baseline first. If a gradient-boosted model can't beat logistic regression by much, that tells you the signal is limited, which is worth knowing before you spend a week tuning. Fixing the categories from the training set (CategoricalDtype) stops an unseen destination in the test week from silently shifting the encoding.

The gradient boosting model trains in about a second on a laptop.

Results

On the held-out final week of January:

ROC AUC of 0.5 is a coin toss and 1.0 is perfect. So schedule and weather alone get you a modest step above chance, and one live operational feature adds as much as everything else combined.

Permutation importance tells the same story. Without the live feature, the top signals were local departure hour, visibility and airline. With it, origin_recent_delay_rate goes straight to first place, ahead of airline and hour.

What the numbers mean in practice

AUC is abstract. A threshold makes it concrete:

from sklearn.metrics import precision_score, recall_score

for thr in (0.3, 0.5):
    pred = (p_gbm >= thr).astype(int)
    print(f"threshold {thr}: flags {pred.mean():.1%}, "
          f"precision {precision_score(y_test, pred):.1%}, "
          f"recall {recall_score(y_test, pred):.1%}")

The base rate in the test week is 17.3%. At a threshold of 0.5 the model flags one flight in twenty, and 40% of those are actually delayed, which is more than double the base rate. That's useful for ranking which flights to watch. It is nowhere near good enough to tell a passenger "your flight will be late".

Precision falls when you lower the threshold to catch more delays, and there's no setting that gets both high. That trade-off is the honest shape of this problem with this information.

Why it isn't better, and what would help

The model is missing the largest single cause of delay. In BTS's own delay-cause columns, a late-arriving aircraft is consistently one of the biggest contributors. If the aircraft scheduled to operate your 14:00 departure is still in the air 400 miles away at 13:30, your flight will be late. Nothing in the schedule or the weather knows that.

Three things would help, in roughly descending order:

1. The inbound aircraft's live status. Where is the airframe that will operate this flight, and is it running late? That's the most predictive single input available, and it's live data, not historical.
2. Air traffic control programmes. A ground delay programme or ground stop at the destination holds departures regardless of how everything else looks. These are published in real time.
3. Forecast weather at departure time rather than the last observation, particularly for thunderstorms, which build fast and cause outsized disruption.

All three are live operational data. The pattern from Step 3b generalises: the best predictor of delay is delay that's already happening somewhere in the system.

Taking it to production

If you want to run this on live flights rather than historical files, the inputs above map onto real endpoints. To be clear about what exists: SkyLink doesn't offer a single "will this flight be delayed" classifier. What it provides are the live inputs this tutorial shows you need.

• Flight status gives current departure and arrival delay for a flight, which is the live version of the inbound-aircraft signal.
• FAA delays gives active ground stops and ground delay programmes at US airports.
• METAR and TAF give the observation and the forecast.
• ML flight time gives a gate-to-gate block-time estimate with min and max bounds, a statistical baseline for how long the flight should take.

Flight delay prediction with an API shows how those combine into a working delay estimate without training anything. All of them are available on the free trial. A sensible path is to train on the historical BTS data, as in this tutorial, and then feed the model live features from the API at prediction time.

One caution if you do: the BTS dataset covers US domestic flights by reporting carriers only. A model trained on it has never seen an international flight or a small carrier, and won't transfer to them without new training data.