Port Operational Resilience to Coastal Hazards

Using ACES for large-scale Geospatial analysis in Disaster Management and Sustainability

Problem Statement

Ports are critical to trade and coastal communities, but how much storms and high water disrupt daily port activity is difficult to quantify.

Using AIS transit data, water-level, and storm-track datasets, this project asks How can we measure and predict port disruption from coastal hazards at scale?

Authors: Shoibolina Kaushik (Graduate Research Assistant), Zhe Zhang (Associate Professor)
CIDI-Spatial Lab
Department of Geography
Texas A&M University

---

# Setup
import os, time, warnings
from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf

from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_auc_score, roc_curve
from joblib import Parallel, delayed

warnings.filterwarnings("ignore")

FAST_DEMO = False # used later to speed up the parallel computing demo, if needed

try:
    import xgboost as xgb
    HAVE_XGB = True
except ImportError:
    from sklearn.ensemble import GradientBoostingClassifier
    HAVE_XGB = False

try:
    import shap
    HAVE_SHAP = True
except ImportError:
    HAVE_SHAP = False

WORK    = Path.cwd()
OUT_DIR = WORK / "data"
FIG_DIR = WORK / "figures"
FIG_DIR.mkdir(parents=True, exist_ok=True)


def allocated_cores():
    """How many cores this session actually has.

    os.cpu_count() reports the whole compute node, not the share the portal gave you.
    """
    try:
        n = len(os.sched_getaffinity(0))
    except AttributeError:
        n = os.cpu_count() or 1
    slurm = os.environ.get("SLURM_CPUS_PER_TASK")
    return min(n, int(slurm)) if slurm and slurm.isdigit() else n


N_CORES = allocated_cores()
print(f"xgboost   : {'yes' if HAVE_XGB else 'no -- using sklearn GradientBoosting'}")
print(f"shap      : {'yes' if HAVE_SHAP else 'no -- section 8 falls back to importances'}")
print(f"cores     : {N_CORES}   (os.cpu_count says {os.cpu_count()} -- that is the node)")
print(f"fast demo : {FAST_DEMO}")

xgboost   : yes
shap      : yes
cores     : 8   (os.cpu_count says 96 -- that is the node)
fast demo : False

1. Data

What are we doing?

We begin with an analysis-ready dataset that was created in separate preprocessing notebooks.

Those notebooks downloaded and processed:

• AIS ship-location data
• NOAA water-level data
• Tropical storm-track data

The large raw datasets (~871GB) were processed on ACES and saved as compact Parquet files. This notebook starts from those Parquet files.

---

Each row now represents one port on one day.

For each port-day, we have information about:

• how many ships were near the port,
• water-level conditions,
• whether a storm was nearby,
• storm distance and wind speed.

panel = pd.read_parquet(OUT_DIR / "panel_daily_6ports.parquet")
panel["date"] = pd.to_datetime(panel["date"]).dt.normalize()
panel = panel.sort_values(["port", "date"]).reset_index(drop=True)
panel

	date	port	n_points	unique_vessels	unique_commercial	n_moored	n_anchored	n_underway	n_stationary	mean_sog	...	storm_pts_200km	log_unique_vessels	lat	lon	frac_moored	frac_anchored	frac_underway	frac_stationary	commercial_share	activity_ratio
0	2019-01-04	Charleston	36466	101	21	12404	817	11284	30848	0.998316	...	0	4.624973	32.78	-79.92	0.340152	0.022404	0.309439	0.845939	0.207921	NaN
1	2019-01-05	Charleston	39449	100	19	12109	661	12956	34365	0.812723	...	0	4.615121	32.78	-79.92	0.306953	0.016756	0.328424	0.871125	0.190000	NaN
2	2019-01-06	Charleston	40287	101	21	12184	716	11211	34754	0.872778	...	0	4.624973	32.78	-79.92	0.302430	0.017772	0.278278	0.862660	0.207921	NaN
3	2019-01-07	Charleston	36480	106	22	11670	544	11670	30883	0.964361	...	0	4.672829	32.78	-79.92	0.319901	0.014912	0.319901	0.846573	0.207547	NaN
4	2019-01-08	Charleston	40467	106	23	11039	529	14345	33764	1.111713	...	0	4.672829	32.78	-79.92	0.272790	0.013072	0.354486	0.834359	0.216981	1.049505
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
6829	2024-10-12	Tampa	26929	60	4	3371	24	11185	24465	0.508066	...	0	4.110874	27.95	-82.45	0.125181	0.000891	0.415351	0.908500	0.066667	0.662983
6830	2024-10-13	Tampa	35447	72	8	5560	546	15412	31661	0.641442	...	0	4.290459	27.95	-82.45	0.156854	0.015403	0.434790	0.893193	0.111111	0.800000
6831	2024-10-14	Tampa	40738	88	13	6136	870	20557	36395	0.616508	...	0	4.488636	27.95	-82.45	0.150621	0.021356	0.504615	0.893392	0.147727	0.983240
6832	2024-10-15	Tampa	36899	91	14	7177	1163	15705	32668	0.653411	...	0	4.521789	27.95	-82.45	0.194504	0.031518	0.425621	0.885336	0.153846	1.011111
6833	2024-10-16	Tampa	39432	87	15	9482	924	15543	34427	0.737749	...	0	4.477337	27.95	-82.45	0.240465	0.023433	0.394172	0.873073	0.172414	0.977528

6834 rows × 27 columns

panel = pd.read_parquet(OUT_DIR / "panel_daily_6ports.parquet")
panel["date"] = pd.to_datetime(panel["date"]).dt.normalize()
panel = panel.sort_values(["port", "date"]).reset_index(drop=True)

print(f"{len(panel):,} port-days  |  {panel['port'].nunique()} ports  |  "
      f"{panel['date'].min():%Y-%m-%d} to {panel['date'].max():%Y-%m-%d}\n")

panel.groupby("port").agg(
    lat=("lat", "first"),
    days=("date", "nunique"),
    storm_days=("storm_day_200km", "sum"),
    ships_per_day=("unique_vessels", "mean"),
    level_1_in_20_m=("wl_mean", lambda s: s.quantile(0.95)),
).round(2).sort_values("lat")

6,834 port-days  |  6 ports  |  2019-01-04 to 2024-10-16

	lat	days	storm_days	ships_per_day	level_1_in_20_m
port
Tampa	27.95	1139	13	77.05	0.36
Houston	29.73	1139	10	419.50	0.48
NewOrleans	29.95	1139	17	389.34	0.53
Charleston	32.78	1139	26	132.68	0.46
Norfolk	36.85	1139	13	243.61	0.50
NewYorkNJ	40.67	1139	8	336.26	0.42

What should you notice?

Look at how different the six ports are.

Some ports normally see hundreds of ships per day, while others see far fewer. This is important later because the same change in ship count can mean very different things at different ports.

SHOW_PORT = "Charleston" # the port with the most storm events

g = panel[panel["port"] == SHOW_PORT]

fig, ax = plt.subplots(figsize=(14, 4))
ax.plot(g["date"], g["unique_vessels"], lw=0.8, color="tab:orange", label="ships per day")
ax.set_ylabel("distinct ships within 20 km")

ax2 = ax.twinx()
ax2.plot(g["date"], g["wl_anom"], lw=0.8, ls="--", color="tab:blue",
         label="water above normal (m)")
ax2.set_ylabel("metres above the 30-day typical level")

for d in g.loc[g["storm_day_200km"] == 1, "date"]:
    ax.axvspan(d, d + pd.Timedelta(days=1), color="red", alpha=0.15)

h1, l1 = ax.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax.legend(h1 + h2, l1 + l2, loc="upper left", fontsize=8)
ax.set_title(f"{SHOW_PORT}: ships vs water level   (red = storm within 200 km)")
plt.tight_layout()
plt.show()

print("The panel holds only days around hazard events, so the line sometimes joins days weeks apart.")

The panel holds only days around hazard events, so the line sometimes joins days weeks apart.

---

2. Effect of storm

What are we asking?

What happens to ship activity when a storm reaches a port?

For each storm event, we compare activity on the storm-arrival day with activity before the storm.

We calculate:

• how much vessel activity dropped, and
• how long it took activity to return to at least 90% of normal.

TAU_PRE, TAU_POST = 7, 7 # we take 7 days before the storm event, and 7 days after
RECOVERY_FRAC = 0.90 # 90% of normal operation is set as the port recovery fraction here


def storm_events(df):
    """Collapse runs of consecutive storm days into one event each."""
    out = []
    for port, g in df.groupby("port"):
        g = g.sort_values("date")
        flag = g["storm_day_200km"].to_numpy()
        dates = g["date"].to_numpy()
        prev = np.concatenate([[0], flag[:-1]])
        gap = np.concatenate([[True],
              (dates[1:] - dates[:-1]).astype("timedelta64[D]").astype(int) > 1])
        for ed in dates[(flag == 1) & ((prev == 0) | gap)]:
            out.append({"port": port, "event_date": pd.Timestamp(ed)})
    return pd.DataFrame(out)


def event_metrics(df, events):
    """Per event: the drop on day 0, and days until 90% of baseline is back."""
    rows = []
    for ev in events.itertuples(index=False):
        g = df[df["port"] == ev.port].set_index("date")
        win = g.reindex(pd.date_range(ev.event_date - pd.Timedelta(days=TAU_PRE),
                                      ev.event_date + pd.Timedelta(days=TAU_POST)))
        win["tau"] = (win.index - ev.event_date).days

        base = win.loc[win["tau"].between(-TAU_PRE, -1), "unique_vessels"].mean()
        day0 = win.loc[win["tau"] == 0, "unique_vessels"].iloc[0]
        if pd.isna(base) or base <= 0 or pd.isna(day0):
            continue

        post = win[win["tau"] >= 0]
        rec = post[post["unique_vessels"] >= RECOVERY_FRAC * base]
        at0 = win.loc[win["tau"] == 0].iloc[0]
        rows.append({
            "port": ev.port, "event_date": ev.event_date,
            "baseline_ships": base, "day0_ships": day0,
            "drop_frac": (base - day0) / base,
            "recovery_days": int(rec["tau"].iloc[0]) if len(rec) else np.nan,
            "wl_anom_day0": at0["wl_anom"],
            "wind_kt_day0": at0["storm_max_wind_kt_200km"],
            "dist_km_day0": at0["storm_min_dist_km"],
        })
    return pd.DataFrame(rows)


events = storm_events(panel)
ev = event_metrics(panel, events)

print(f"{len(events)} storm events across {events['port'].nunique()} ports, "
      f"{len(ev)} measurable\n")
print(ev.groupby("port").agg(
    events=("event_date", "count"),
    median_drop=("drop_frac", "median"),
    worst_drop=("drop_frac", "max"),
    median_recovery_days=("recovery_days", "median"),
).round(3).to_string())

print(f"\nall events: typical drop {ev['drop_frac'].median()*100:.0f}%, "
      f"worst {ev['drop_frac'].max()*100:.0f}%")

56 storm events across 6 ports, 56 measurable

            events  median_drop  worst_drop  median_recovery_days
port                                                             
Charleston      17        0.154       0.384                   1.0
Houston          6        0.193       0.297                   1.0
NewOrleans      10        0.151       0.282                   2.0
NewYorkNJ        5        0.106       0.263                   1.0
Norfolk          9        0.152       0.291                   1.0
Tampa            9        0.200       0.502                   2.0

all events: typical drop 16%, worst 50%

What should you notice?

Do not focus on every individual storm. Instead, look for the overall pattern:

Do ports usually become quieter during storms, and how large is the typical decline?

INK, ACCENT, STORM, MUTED = "#15242B", "#0B6E8C", "#B8433A", "#8FA3AB"

med   = ev.groupby("port")["drop_frac"].median().sort_values()
cnt   = ev.groupby("port").size()
worst = ev.groupby("port")["drop_frac"].max()
overall = ev["drop_frac"].median()

fig, ax = plt.subplots(figsize=(9, 4.6), facecolor="white")
bars = ax.barh(med.index, med.values * 100,
               color=[STORM if v >= overall else ACCENT for v in med.values],
               height=.62, zorder=3)

# each port's worst single event, as a faint tick behind the bar
ax.scatter(worst[med.index].values * 100, med.index,
           marker="|", s=260, color=MUTED, zorder=4, label="worst single event")

for y, (port, v) in enumerate(med.items()):
    ax.text(v * 100 + .6, y, f"{v*100:.1f}%", va="center", fontsize=11,
            color=INK, fontweight="bold")
    ax.text(-.6, y, f"{cnt[port]} events", va="center", ha="right",
            fontsize=9, color=MUTED)

ax.axvline(overall * 100, color=INK, ls="--", lw=1.2, zorder=5)
ax.annotate(f"all 56 events: {overall*100:.0f}%",
            xy=(overall * 100, len(med) - .35),
            xytext=(overall * 100 + 4, len(med) - .35),
            fontsize=10, color=INK, va="center",
            arrowprops=dict(arrowstyle="-", color=INK, lw=1))

ax.set_xlabel("ships lost on the day the storm arrived (%)", fontsize=11)
ax.set_xlim(-12, max(worst) * 100 + 6)
ax.set_title("Disruption of port traffic when a storm arrives", fontsize=14,
             fontweight="bold", color=INK, loc="left", pad=14)
ax.legend(loc="lower right", frameon=False, fontsize=9)
ax.grid(axis="x", alpha=.25, zorder=0)
for s in ("top", "right", "left"):
    ax.spines[s].set_visible(False)
ax.tick_params(left=False)

plt.tight_layout()
plt.savefig(FIG_DIR / "fig1_storm_decline.png", dpi=200, facecolor="white")
plt.show()

Across the events in this dataset, the typical decline is about 16%, although some storms produce much larger changes.

This gives us a descriptive answer. Next, we ask which part of the hazard is associated with that decline.

---

3. Which hazard actually stops a port?

What are we asking?

A port may become quieter for several related reasons.

For example:

• a storm may be nearby,
• water may be unusually high,
• or both may occur together.

We use regression to try to separate these effects.

The model also accounts for the fact that every port has a different normal level of activity.

First attempt

The first regression suggests that ship activity is lower when a storm is within 200 km.

# 3.1  First attempt: one storm effect, shared by all six ports.
reg = panel.dropna(subset=["wl_anom", "unique_vessels"]).copy()
reg["log_ships"] = np.log1p(reg["unique_vessels"])

daily = smf.ols("log_ships ~ wl_anom + storm_day_200km + C(port)", data=reg)\
           .fit(cov_type="HAC", cov_kwds={"maxlags": 7})
print(daily.summary().tables[1])


def plain(term, label):
    b, p = daily.params[term], daily.pvalues[term]
    pct = (np.exp(b) - 1) * 100 # log coefficient -> percent
    ptxt = "p < 0.001" if p < 0.001 else f"p = {p:.3f}"
    verdict = "convincing" if p < 0.05 else "NOT convincing"
    print(f"{label}: {pct:+.1f}% ships")
    print(f"    {ptxt}  ->  {verdict}")


print()
plain("storm_day_200km", "storm within 200 km    ")
plain("wl_anom",         "water 1 m above normal ")
print(f"\n{int(daily.nobs):,} port-days, {reg['port'].nunique()} ports. "
      "The usual bar is p < 0.05 -- how often chance alone would produce this.")

=========================================================================================
                            coef    std err          z      P>|z|      [0.025      0.975]
-----------------------------------------------------------------------------------------
Intercept                 4.8765      0.016    314.448      0.000       4.846       4.907
C(port)[T.Houston]        1.1593      0.017     67.892      0.000       1.126       1.193
C(port)[T.NewOrleans]     1.0871      0.017     65.398      0.000       1.055       1.120
C(port)[T.NewYorkNJ]      0.9195      0.023     39.418      0.000       0.874       0.965
C(port)[T.Norfolk]        0.6101      0.020     30.123      0.000       0.570       0.650
C(port)[T.Tampa]         -0.5324      0.020    -27.184      0.000      -0.571      -0.494
wl_anom                  -0.0351      0.031     -1.122      0.262      -0.097       0.026
storm_day_200km          -0.1188      0.042     -2.834      0.005      -0.201      -0.037
=========================================================================================

storm within 200 km    : -11.2% ships
    p = 0.005  ->  convincing
water 1 m above normal : -3.5% ships
    p = 0.262  ->  NOT convincing

6,791 port-days, 6 ports. The usual bar is p < 0.05 -- how often chance alone would produce this.

But should we automatically trust the first statistically significant result?

Check the result separately across ports... Notice something strange?

# 3.2  Let every port have its own storm effect, then ask whether they really differ.
inter = smf.ols("log_ships ~ wl_anom + storm_day_200km * C(port)", data=reg)\
           .fit(cov_type="HAC", cov_kwds={"maxlags": 7})


def storm_effect(port, model, frame):
    """One port's storm effect = the shared term + that port's own adjustment."""
    names = list(model.params.index)
    base = "storm_day_200km"
    vec = np.zeros(len(names))
    vec[names.index(base)] = 1.0
    hit = [n for n in names if n.startswith(base + ":") and f"T.{port}]" in n]
    if hit: # the reference port has no adjustment term
        vec[names.index(hit[0])] = 1.0
    t = model.t_test(vec) # gets the uncertainty of the sum right
    lo, hi = t.conf_int()[0]
    p = float(np.ravel(t.pvalue)[0])
    return {"port": port,
            "storm_days": int(frame.loc[frame["port"] == port, base].sum()),
            "pct": (np.exp(t.effect.item()) - 1) * 100,
            "lo": (np.exp(lo) - 1) * 100,
            "hi": (np.exp(hi) - 1) * 100,
            "p": p,
            "verdict": "convincing" if p < 0.05 else "not convincing"}


def ports_differ(model):
    """One test: do the six ports react differently at all?"""
    names = list(model.params.index)
    terms = [n for n in names if n.startswith("storm_day_200km:")]
    R = np.zeros((len(terms), len(names)))
    for i, n in enumerate(terms):
        R[i, names.index(n)] = 1.0
    return float(np.ravel(model.f_test(R).pvalue)[0])


eff = pd.DataFrame([storm_effect(p, inter, reg) for p in sorted(reg["port"].unique())])
print("storm effect per port -- % change in daily ships, with 95% range\n")
print(eff.round(2).to_string(index=False))
print(f"\ndo the ports react differently?  p = {ports_differ(inter):.3f}")
print("\n3.1 said about -11% for everyone. Compare that with the column of numbers.")

storm effect per port -- % change in daily ships, with 95% range

      port  storm_days    pct     lo    hi    p        verdict
Charleston          26  -8.44 -19.80  4.53 0.19 not convincing
   Houston          10 -12.80 -20.28 -4.62 0.00     convincing
NewOrleans          17 -27.74 -44.22 -6.39 0.01     convincing
 NewYorkNJ           8  12.92   1.95 25.08 0.02     convincing
   Norfolk          13   8.25  -3.86 21.88 0.19 not convincing
     Tampa          13 -22.45 -34.96 -7.54 0.00     convincing

do the ports react differently?  p = 0.000

3.1 said about -11% for everyone. Compare that with the column of numbers.

Two ports get busier during storms ( NY/NJ +12.9% and Norfolk +8.3% )

That tells us to investigate the comparison more carefully.

# 3.3  When do storm days happen, and what else is different about those days?
chk = reg.assign(month=reg["date"].dt.month)

nyc = chk[chk["port"] == "NewYorkNJ"].copy()
nyc["day type"] = np.where(nyc["storm_day_200km"] == 1, "storm day", "other day")

print("NewYorkNJ -- how many days of each kind fall in each month\n")
print(nyc.groupby("day type")["month"].value_counts().unstack(fill_value=0).to_string())

print("\nNewYorkNJ -- median ships by month, the port's own seasonal cycle\n")
print(nyc.groupby("month")["unique_vessels"].median().round(0).to_string())

NewYorkNJ -- how many days of each kind fall in each month

month       1    2    3   4   5   6   7   8    9   10  11  12
day type                                                     
other day  175  168  179  98  29  57  77  78  133  69  36  26
storm day    0    0    0   0   0   0   3   4    1   0   0   0

NewYorkNJ -- median ships by month, the port's own seasonal cycle

month
1     286.0
2     291.0
3     294.0
4     314.0
5     251.0
6     399.0
7     430.0
8     403.0
9     440.0
10    386.0
11    324.0
12    286.0

Finding the problem

Tropical storms mainly occur during hurricane season.

But our original comparison included quiet days from the whole year.

That means we may partly be comparing:

summer storm days vs. winter non-storm days

rather than simply:

storm days vs. similar non-storm days

We therefore restrict the comparison to the relevant season and control for month.

# 3.4  Compare like with like: hurricane season only (June-November).
season = chk[chk["month"].between(6, 11)]

inter3 = smf.ols("log_ships ~ wl_anom + storm_day_200km * C(port)", data=season)\
            .fit(cov_type="HAC", cov_kwds={"maxlags": 7})

eff3 = pd.DataFrame([storm_effect(p, inter3, season)
                     for p in sorted(season["port"].unique())])

print(f"hurricane season only -- {len(season):,} of {len(reg):,} port-days\n")
print(eff3.round(2).to_string(index=False))

p_diff = ports_differ(inter3)
n_neg = int((eff3["pct"] < 0).sum())
print(f"\n{n_neg} of {len(eff3)} ports now show a drop.")
print(f"do the ports react differently?  p = {p_diff:.3f}  --  "
      + ("yes, keep them separate." if p_diff < 0.05
         else "no, one shared number is defensible again."))

hurricane season only -- 2,741 of 6,791 port-days

      port  storm_days    pct     lo    hi    p        verdict
Charleston          25 -17.59 -27.23 -6.68 0.00     convincing
   Houston          10 -13.27 -20.41 -5.49 0.00     convincing
NewOrleans          17 -21.77 -37.74 -1.71 0.04     convincing
 NewYorkNJ           8  -6.13 -15.63  4.43 0.24 not convincing
   Norfolk          13  -2.99 -13.53  8.84 0.61 not convincing
     Tampa          13 -16.39 -31.76  2.44 0.08 not convincing

6 of 6 ports now show a drop.
do the ports react differently?  p = 0.248  --  no, one shared number is defensible again.

# 3.5  The final model: hurricane season, each port AND each month against its own normal.
final = smf.ols("log_ships ~ wl_anom + storm_day_200km + C(port) + C(month)",
                data=season).fit(cov_type="HAC", cov_kwds={"maxlags": 7})

b_s = final.params["storm_day_200km"]
lo_s, hi_s = final.conf_int().loc["storm_day_200km"]
p_s = final.pvalues["storm_day_200km"]
ptxt = "p < 0.001" if p_s < 0.001 else f"p = {p_s:.3f}"

print(f"{int(final.nobs):,} port-days. Change in daily ship count:\n")
print(f"storm within 200 km : {(np.exp(b_s)-1)*100:+.1f}%   "
      f"[{(np.exp(lo_s)-1)*100:+.1f}%, {(np.exp(hi_s)-1)*100:+.1f}%]   {ptxt}")

b_w = final.params["wl_anom"]
lo_w, hi_w = final.conf_int().loc["wl_anom"]
p_w = final.pvalues["wl_anom"]
q = season["wl_anom"].quantile([.5, .9, .99, 1.0])

print(f"\nwater level         : p = {p_w:.3f}. The coefficient is per metre, but half of all")
print(f"                      days sit within {q.iloc[0]:.2f} m of normal -- so read it at")
print( "                      levels the data actually contains:\n")
for m, tag in [(q.iloc[1], "1 day in 10"), (q.iloc[2], "1 day in 100"),
               (q.iloc[3], "worst in 6 years")]:
    print(f"    +{m:.2f} m  ({tag:16s}): {(np.exp(b_w*m)-1)*100:+5.1f}%   "
          f"[{(np.exp(lo_w*m)-1)*100:+.1f}%, {(np.exp(hi_w*m)-1)*100:+.1f}%]")

if hi_w > 0:
    print("\nEvery one of those ranges still touches zero. Not enough evidence is a")
    print("different statement from evidence of no effect.")

2,741 port-days. Change in daily ship count:

storm within 200 km : -14.8%   [-20.3%, -8.9%]   p < 0.001

water level         : p = 0.085. The coefficient is per metre, but half of all
                      days sit within 0.01 m of normal -- so read it at
                      levels the data actually contains:

    +0.20 m  (1 day in 10     ):  -2.1%   [-4.5%, +0.3%]
    +0.53 m  (1 day in 100    ):  -5.5%   [-11.3%, +0.8%]
    +1.30 m  (worst in 6 years): -12.9%   [-25.5%, +1.9%]

Every one of those ranges still touches zero. Not enough evidence is a
different statement from evidence of no effect.

Corrected result

After making the comparison more appropriate, a storm within 200 km is associated with roughly a 15% reduction in daily ship activity.

The water-level result is less certain in this dataset.

We also examine storm events themselves and ask whether stronger storms produce larger declines.

---

Visualizing the final result

The regression gives us an estimated effect and a range of plausible values.

In the figure below:

• the dot is the estimated change in daily ship activity;
• the horizontal line is the 95% confidence interval;
• the vertical line at 0% means “no change.”

If a confidence interval crosses 0, the data do not clearly distinguish that effect from no effect.

The top panel compares the two hazards in the final model.
The bottom panel checks whether the storm effect looks reasonably similar across individual ports.

Look for: Does the overall storm effect stay below zero, and do most ports point in the same direction?

# 3.6  Visualize the final hazard effects and check consistency across ports.

HAIR = "#DDE3E6"

# Show water level at a value the dataset actually experiences:
# the 99th percentile of hurricane-season water-level anomalies.
wl_step = season["wl_anom"].quantile(0.99)

# ---- Overall effects from the final model --------------------------------
top = []
for term, label, scale in [
    ("storm_day_200km", "Storm within 200 km", 1.0),
    ("wl_anom", f"Water +{wl_step:.2f} m", wl_step),
]:
    b = final.params[term]
    lo, hi = final.conf_int().loc[term]

    top.append({
        "label": label,
        "pct": (np.exp(b * scale) - 1) * 100,
        "lo":  (np.exp(lo * scale) - 1) * 100,
        "hi":  (np.exp(hi * scale) - 1) * 100,
        "sig": final.pvalues[term] < 0.05,
    })

# Per-port storm effects from the hurricane-season check
bot = eff3.sort_values("pct").reset_index(drop=True)
pooled_pct = top[0]["pct"]

# ---- Plot ---------------------------------------------------------------
fig, (axA, axB) = plt.subplots(
    2, 1,
    figsize=(9, 6.4),
    sharex=True,
    gridspec_kw={"height_ratios": [2, 5], "hspace": 0.30},
)
# Panel A: overall hazard effects
for i, r in enumerate(top):
    col = STORM if r["sig"] else MUTED

    axA.plot(
        [r["lo"], r["hi"]], [i, i],
        color=col, lw=4, solid_capstyle="round"
    )
    axA.scatter(
        r["pct"], i,
        s=110, color=col,
        edgecolor="white", linewidth=1.5, zorder=3
    )

    # Only label the estimated effect
    axA.text(
        r["pct"], i - 0.20,
        f"{r['pct']:+.1f}%",
        ha="center", va="bottom",
        fontsize=11, fontweight="bold", color=INK
    )
axA.set_yticks(range(len(top)))
axA.set_yticklabels([r["label"] for r in top], fontsize=11)
axA.set_ylim(len(top) - 0.35, -0.65)
axA.set_title(
    "Overall effects",
    loc="left", fontsize=11, fontweight="bold", color=INK
)
# Panel B: storm effect separately at each port
axB.axvline(
    pooled_pct,
    color=STORM, ls="--", lw=1.2, alpha=0.7
)

for i, r in bot.iterrows():
    col = STORM if r["p"] < 0.05 else MUTED

    axB.plot(
        [r["lo"], r["hi"]], [i, i],
        color=col, lw=3, solid_capstyle="round"
    )
    axB.scatter(
        r["pct"], i,
        s=75, color=col,
        edgecolor="white", linewidth=1.2, zorder=3
    )
axB.set_yticks(range(len(bot)))
axB.set_yticklabels(bot["port"], fontsize=10.5)
axB.set_ylim(len(bot) - 0.4, -0.65)
axB.set_title(
    "Storm effect estimated separately for each port",
    loc="left", fontsize=11, fontweight="bold", color=INK
)
axB.set_xlabel("Estimated change in daily ship activity (%)", fontsize=11)

# Shared plt formatting
span_lo = min([r["lo"] for r in top] + list(bot["lo"]))
span_hi = max([r["hi"] for r in top] + list(bot["hi"]))
width = span_hi - span_lo

for ax in (axA, axB):
    ax.axvline(0, color=INK, lw=1)
    ax.set_xlim(span_lo - width * 0.05,
                span_hi + width * 0.08)
    ax.grid(axis="x", color=HAIR, linewidth=0.8)
    ax.tick_params(left=False)

    for spine in ("top", "right", "left"):
        ax.spines[spine].set_visible(False)

axA.spines["bottom"].set_visible(False)
axB.spines["bottom"].set_color(HAIR)

fig.suptitle(
    "Estimated effects of coastal hazards on daily port activity",
    x=0.12, ha="left",
    fontsize=15, fontweight="bold", color=INK
)

fig.text(
    0.12, 0.92,
    "Dots are estimates; horizontal lines are 95% confidence intervals.",
    fontsize=9.5, color=MUTED
)

plt.tight_layout(rect=(0, 0, 1, 0.91))
plt.savefig(FIG_DIR / "fig2_hazards.png",
            dpi=200, facecolor="white", bbox_inches="tight")
plt.show()

What makes one storm more disruptive than another?

We now look only at storm events.

The previous analysis asked:

Does having a nearby storm reduce port activity?

This analysis asks a different question:

Once a storm is already nearby, which storm characteristics are associated with a larger drop?

For each storm event, we compare the observed activity drop with:

• storm wind speed,
• storm distance from the port, and
• water-level anomaly.

This is a simple linear regression, so it tests whether larger values are associated with consistently larger or smaller disruptions.

Look for: Which storm characteristic has the clearest relationship with the size of the activity drop?

Keep in mind that this is based on a relatively small number of storm events, so these results should be interpreted as associations rather than exact physical rules.

# 3.7  Given that a storm hit, what made the drop bigger? One row per storm event.
evr = ev.dropna(subset=["drop_frac", "wind_kt_day0", "dist_km_day0", "wl_anom_day0"]).copy()
storm_lin = smf.ols("drop_frac ~ wind_kt_day0 + dist_km_day0 + wl_anom_day0",
                    data=evr).fit(cov_type="HC3")

print(f"{len(evr)} storm events\n")
print(storm_lin.summary().tables[1])

p_dist = storm_lin.pvalues["dist_km_day0"]
p_wind = storm_lin.pvalues["wind_kt_day0"]
b_wind = storm_lin.params["wind_kt_day0"]

print(f"\nwind     : {b_wind*100:+.2f}% extra drop per knot"
      f"   -> a 100 kt storm costs about {b_wind*100*100:.0f}% more activity"
      f"   (p = {p_wind:.3f})")
print(f"distance : p = {p_dist:.3f}"
      + ("  -- not convincing" if p_dist > 0.05 else "  -- convincing"))
print("\nThis model fits a straight line: it asks whether every extra kilometre changes the")
print("drop by the same amount. Keep that in mind for section 7.")

56 storm events

================================================================================
                   coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------
Intercept        0.0216      0.055      0.396      0.692      -0.085       0.129
wind_kt_day0     0.0027      0.001      4.318      0.000       0.001       0.004
dist_km_day0    -0.0002      0.000     -0.753      0.451      -0.001       0.000
wl_anom_day0    -0.0383      0.046     -0.826      0.409      -0.129       0.053
================================================================================

wind     : +0.27% extra drop per knot   -> a 100 kt storm costs about 27% more activity   (p = 0.000)
distance : p = 0.451  -- not convincing

This model fits a straight line: it asks whether every extra kilometre changes the
drop by the same amount. Keep that in mind for section 7.

Circling back to the question: Which hazard actually stops a port?

After correcting for seasonality, the main result becomes much clearer.

• A storm within 200 km is associated with about 15% fewer ships per day.
• The estimated range is roughly 9% to 20% fewer ships, and it does not cross zero.
• After comparing storm days only with similar hurricane-season days, all six ports show a decline in activity.

Water level may also matter, but the evidence is weaker.

For a very high-water day of about +0.53 m above normal, the model estimates roughly 5.5% fewer ships, but the confidence interval still crosses zero. This means the data are not strong enough to clearly separate the effect from no effect.

Among storm events, stronger winds are associated with larger activity declines, while exact distance of storm within the 200 km window shows a less clear relationship.

Takeaway

Nearby storms are consistently associated with lower port activity, and stronger storms tend to produce larger declines. Water level may contribute as well, but this dataset does not provide equally strong evidence for its independent effect.

Main lesson

A regression can give a precise-looking answer even when the comparison behind it is unfair.

Always ask what observations are actually being compared.

---

4. Predicting Port Operation Disruption (given hazards data for a particular day)

What are we asking now?

So far we have described what happened during past storms.

Now we switch to a prediction question:

Can hazard information tell us whether a port is likely to have an unusually quiet day?

We define a disrupted day as a day when vessel activity falls:

below 90% of that port's recent normal activity

The model receives information about:

• water level,
• recent water-level conditions,
• storm presence,
• storm distance,
• storm wind,
• and season.

df = panel.copy()

# surge builds over days, so give the model recent history as well as today
g = df.groupby("port")["wl_anom"]
df["wl_anom_lag1"]  = g.shift(1)
df["wl_anom_lag2"]  = g.shift(2)
df["wl_anom_max3d"] = g.transform(lambda s: s.rolling(3, min_periods=1).max())

# Blank storm fields are informative, not missing. No storm nearby means no wind reading and no track points
df["storm_max_wind_kt_200km"] = df["storm_max_wind_kt_200km"].fillna(0.0)
df["storm_pts_200km"]         = df["storm_pts_200km"].fillna(0)
df["storm_min_dist_km"]       = df["storm_min_dist_km"].fillna(5000.0)

doy = df["date"].dt.dayofyear
df["doy_sin"] = np.sin(2 * np.pi * doy / 365.25)
df["doy_cos"] = np.cos(2 * np.pi * doy / 365.25)

FEATURES = [
    "wl_anom", "wl_anom_lag1", "wl_anom_lag2", "wl_anom_max3d",
    "storm_day_200km", "storm_min_dist_km", "storm_max_wind_kt_200km", "storm_pts_200km",
    "doy_sin", "doy_cos",
]
TARGET = "disrupted"
DISRUPTION_THRESHOLD = 0.90

df[TARGET] = (df["activity_ratio"] < DISRUPTION_THRESHOLD).astype(int)

# Called data_all rather than data: section 6 finds a reason to revise it.
data_all = df.dropna(subset=FEATURES + [TARGET, "activity_ratio"]).reset_index(drop=True)
data_all["year"] = data_all["date"].dt.year
ports = sorted(data_all["port"].unique())

print(f"rows usable : {len(data_all):,}   (dropped {len(df) - len(data_all):,})")
print(f"disrupted   : {data_all[TARGET].sum():,} days  ({data_all[TARGET].mean()*100:.1f}%)\n")
print(data_all.groupby("port")[TARGET].agg(days="size", disrupted="sum", rate="mean").round(3))

rows usable : 6,777   (dropped 57)
disrupted   : 909 days  (13.4%)

            days  disrupted   rate
port                              
Charleston  1131        216  0.191
Houston     1131         88  0.078
NewOrleans  1122         76  0.068
NewYorkNJ   1131        211  0.187
Norfolk     1131        150  0.133
Tampa       1131        168  0.149

Testing geographic transfer

We train the model on five ports and predict the sixth port.** The held-out port is completely unseen during training. We repeat this until every port has been held out once.

This is called leave-one-port-out cross-validation.

def make_model(scale_pos_weight=1.0, **kw):
    """Gradient-boosted trees; xgboost when available."""
    if HAVE_XGB:
        p = dict(n_estimators=300, max_depth=4, learning_rate=0.05, subsample=0.8,
                 colsample_bytree=0.8, scale_pos_weight=scale_pos_weight,
                 eval_metric="auc", random_state=42, n_jobs=1)
        p.update(kw)
        return xgb.XGBClassifier(**p)
    p = dict(n_estimators=300, max_depth=4, learning_rate=0.05, subsample=0.8, random_state=42)
    p.update({k: v for k, v in kw.items() if k in p})
    return GradientBoostingClassifier(**p)


def fit_one_port(held_out, frame, features=FEATURES, **kw):
    """Train on every other port, score the one held out."""
    train = frame[frame["port"] != held_out]
    test  = frame[frame["port"] == held_out]
    pos = max(int(train[TARGET].sum()), 1)
    model = make_model(scale_pos_weight=(len(train) - pos) / pos, **kw)   # rebalance classes
    model.fit(train[features], train[TARGET])
    prob = model.predict_proba(test[features])[:, 1]
    truth = test[TARGET].to_numpy()
    auc = roc_auc_score(truth, prob) if truth.min() != truth.max() else np.nan
    return {"port": held_out, "n_test": len(test), "disrupted": int(truth.sum()),
            "auc": auc, "prob": prob, "truth": truth}


def score_range(truth, prob, n_boot=300, seed=0):
    """Resample the held-out days to see how far the score could move by chance."""
    rng = np.random.default_rng(seed)
    vals = []
    for _ in range(n_boot):
        idx = rng.integers(0, len(truth), len(truth))
        if truth[idx].min() != truth[idx].max():
            vals.append(roc_auc_score(truth[idx], prob[idx]))
    return np.percentile(vals, [2.5, 97.5])


def leave_one_port_out(frame, features=FEATURES, with_range=True):
    fs = [fit_one_port(p, frame, features=features) for p in ports]
    if with_range:
        for f in fs:
            f["lo"], f["hi"] = score_range(f["truth"], f["prob"])
    cols = ["port", "n_test", "disrupted", "auc"] + (["lo", "hi"] if with_range else [])
    return fs, pd.DataFrame([{k: f[k] for k in cols} for f in fs])


folds_all, cv_all = leave_one_port_out(data_all)
pooled_all = roc_auc_score(np.concatenate([f["truth"] for f in folds_all]),
                           np.concatenate([f["prob"] for f in folds_all]))

print(cv_all.round(3).to_string(index=False))
print(f"\nmean of the six : {cv_all['auc'].mean():.3f}")
print(f"all days pooled : {pooled_all:.3f}")
print("\n0.5 = no better than guessing. lo/hi = where the score could plausibly land on a")
print("different sample of days. A range containing 0.5 means the port is unproven.")

      port  n_test  disrupted   auc    lo    hi
Charleston    1131        216 0.715 0.676 0.753
   Houston    1131         88 0.694 0.636 0.747
NewOrleans    1122         76 0.652 0.584 0.720
 NewYorkNJ    1131        211 0.621 0.577 0.670
   Norfolk    1131        150 0.806 0.765 0.848
     Tampa    1131        168 0.512 0.465 0.553

mean of the six : 0.667
all days pooled : 0.648

0.5 = no better than guessing. lo/hi = where the score could plausibly land on a
different sample of days. A range containing 0.5 means the port is unproven.

def plot_lopo(folds, cv, title):
    fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))

    ax = axes[0]
    for f in folds:
        if np.isnan(f["auc"]):
            continue
        fpr, tpr, _ = roc_curve(f["truth"], f["prob"])
        ax.plot(fpr, tpr, lw=1.5, label=f"{f['port']} ({f['auc']:.2f})")
    ax.plot([0, 1], [0, 1], "k--", lw=1, label="guessing")
    ax.set_xlabel("false alarms"); ax.set_ylabel("disruptions caught")
    ax.set_title(f"Predicting a port the model never saw\n{title}")
    ax.legend(fontsize=8); ax.grid(alpha=0.3)

    ax = axes[1]
    o = cv.sort_values("auc")
    ax.barh(o["port"], o["auc"], xerr=[o["auc"] - o["lo"], o["hi"] - o["auc"]],
            color="steelblue", capsize=4)
    ax.axvline(0.5, color="k", ls="--", lw=1)
    ax.set_xlim(0, 1); ax.set_xlabel("score on the held-out port   (0.5 = guessing)")
    ax.set_title(f"Score, with plausible range\n{title}")
    ax.grid(alpha=0.3, axis="x")

    plt.tight_layout()
    return fig


fig = plot_lopo(folds_all, cv_all, "first attempt, all days")
fig.savefig(FIG_DIR / "lopo_first_attempt.png", dpi=150)
plt.show()

Reading AUC

The main score is AUC:

• 0.5 = no better than guessing
• 1.0 = perfect ranking
• higher values indicate better discrimination

The first average score is about 0.67.

That looks promising.

But before celebrating, we need to ask:

What did the model actually learn?

---

5. Is that number real?

Why check the model?

A machine-learning model can predict the correct answer for the wrong reason.

So before trusting the AUC score, we run several checks.

We test whether the model still works when:

• trained on 2019–2022 and tested on 2023–2024,
• the pandemic years 2020–2021 are removed,
• only calendar variables are used,
• and only hazard variables are used.

CALENDAR_ONLY      = ["doy_sin", "doy_cos"]
FEATURES_NO_SEASON = [f for f in FEATURES if not f.startswith("doy_")]

# 1. hold out time rather than a port
tr = data_all[data_all["year"] <= 2022]
te = data_all[data_all["year"] >= 2023]
pos = max(int(tr[TARGET].sum()), 1)
m_time = make_model(scale_pos_weight=(len(tr) - pos) / pos)
m_time.fit(tr[FEATURES], tr[TARGET])
auc_time = roc_auc_score(te[TARGET], m_time.predict_proba(te[FEATURES])[:, 1])

# 2. without the pandemic years
no_covid = data_all[~data_all["year"].isin([2020, 2021])]
auc_nocovid = np.nanmean([fit_one_port(p, no_covid)["auc"] for p in ports])

# 3 and 4. calendar alone, hazards alone
auc_calendar = np.nanmean([fit_one_port(p, data_all, features=CALENDAR_ONLY)["auc"]
                           for p in ports])
auc_hazards  = np.nanmean([fit_one_port(p, data_all, features=FEATURES_NO_SEASON)["auc"]
                           for p in ports])

base = cv_all["auc"].mean()
print(f"everything, all years          : {base:.3f}   (section 5)")
print(f"train 2019-22, predict 2023-24 : {auc_time:.3f}   ({len(tr):,} / {len(te):,} days)")
print(f"everything, without 2020-21    : {auc_nocovid:.3f}")
print(f"calendar only                  : {auc_calendar:.3f}")
print(f"hazards only                   : {auc_hazards:.3f}")
print("\nmeasured against guessing (0.5):")
print(f"  calendar {auc_calendar - 0.5:+.3f}   hazards {auc_hazards - 0.5:+.3f}"
      f"   together {base - 0.5:+.3f}")

everything, all years          : 0.667   (section 5)
train 2019-22, predict 2023-24 : 0.655   (4,791 / 1,986 days)
everything, without 2020-21    : 0.655
calendar only                  : 0.617
hazards only                   : 0.566

measured against guessing (0.5):
  calendar +0.117   hazards +0.066   together +0.167

Something suspicious appears

The first two checks are reassuring: the score stays around 0.66.

But the surprising result is:

Calendar-only prediction performs better than hazards alone.

This means the model may be learning a seasonal pattern instead of only learning storm effects.

So next we investigate what part of the calendar is driving the prediction.

by_month = (data_all.assign(month=data_all["date"].dt.month)
            .groupby("month")[TARGET]
            .agg(days="size", disrupted="sum", rate="mean").round(3))
by_month["hurricane_season"] = np.where(by_month.index.isin(range(6, 12)), "yes", "")

print("disruption rate by month:")
print(by_month.to_string())

peak = by_month["rate"].idxmax()
print(f"\nworst month              : {peak}  ({by_month.loc[peak, 'rate']*100:.0f}% of days)")
print(f"hurricane season average : {by_month.loc[6:11, 'rate'].mean()*100:.0f}%")

disruption rate by month:
       days  disrupted   rate hurricane_season
month                                         
1      1038        236  0.227                 
2      1008        126  0.125                 
3      1074         91  0.085                 
4       588         30  0.051                 
5       174         12  0.069                 
6       342         30  0.088              yes
7       480         86  0.179              yes
8       492         72  0.146              yes
9       795         86  0.108              yes
10      414         40  0.097              yes
11      216         28  0.130              yes
12      156         72  0.462                 

worst month              : 12  (46% of days)
hurricane season average : 12%

MONTHS = ["Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"]
HURR   = set(range(6, 12))          # June-November

rate = by_month["rate"].reindex(range(1, 13))
days = by_month["days"].reindex(range(1, 13))
season_avg = rate.loc[6:11].mean()

colors = [STORM if m == 12 else ACCENT if m in HURR else MUTED for m in range(1, 13)]

fig, ax = plt.subplots(figsize=(10, 4.6), facecolor="white")
ax.bar(MONTHS, rate.values * 100, color=colors, zorder=3, width=.68)

# shade the hurricane season behind the bars
ax.axvspan(4.5, 10.5, color=ACCENT, alpha=.07, zorder=0)
ax.text(7.5, rate.max() * 100 * .96, "hurricane season",
        ha="center", fontsize=10, color=ACCENT, style="italic")

ax.axhline(season_avg * 100, color=ACCENT, ls="--", lw=1.4, zorder=4)
ax.text(-0.4, season_avg * 100 + 1.2,
        f"hurricane-season average {season_avg*100:.0f}%",
        fontsize=9.5, color=ACCENT)

for i, (v, d) in enumerate(zip(rate.values, days.values)):
    ax.text(i, v * 100 + .9, f"{v*100:.0f}%", ha="center", fontsize=10,
            color=INK, fontweight="bold" if i == 11 else "normal")
    ax.text(i, -3.2, f"{d:,}", ha="center", fontsize=8, color=MUTED)

ax.annotate("holiday shipping slowdown\n— nothing to do with weather",
            xy=(11, rate.iloc[11] * 100), xytext=(8.1, rate.max() * 100 * .80),
            fontsize=10.5, color=STORM, ha="center",
            arrowprops=dict(arrowstyle="->", color=STORM, lw=1.4))

ax.set_ylabel("share of days the port ran quiet (%)", fontsize=11)
ax.set_ylim(-5, rate.max() * 100 * 1.16)
ax.set_title("The model's strongest signal was December, not September",
             fontsize=14, fontweight="bold", color=INK, loc="left", pad=14)
ax.text(0, -5, "days in panel", fontsize=8, color=MUTED, ha="center", va="top")
ax.grid(axis="y", alpha=.25, zorder=0)
for s in ("top", "right"):
    ax.spines[s].set_visible(False)

plt.tight_layout()
plt.savefig(FIG_DIR / "fig3_monthly.png", dpi=200, facecolor="white")
plt.show()

The model learned Christmas!🎄🎁

December is disrupted on 46% of days - four times the hurricane-season rate. That pattern is largely a holiday shipping slowdown, not a coastal-hazard signal.

We therefore remove December 18 – January 8 and rerun the evaluation.

# Everything below this point uses `data`, not `data_all`.
HOLIDAY = (((data_all["date"].dt.month == 12) & (data_all["date"].dt.day >= 18)) |
           ((data_all["date"].dt.month == 1)  & (data_all["date"].dt.day <= 8)))

data = data_all[~HOLIDAY].reset_index(drop=True)

print(f"dropped {int(HOLIDAY.sum()):,} holiday days -- {len(data):,} remain "
      f"({int(HOLIDAY.sum()) / len(data_all) * 100:.0f}% of the data)")
print(f"disruption rate {data_all[TARGET].mean()*100:.1f}%  ->  {data[TARGET].mean()*100:.1f}%\n")

rows = []
for label, feats in [("calendar only", CALENDAR_ONLY),
                     ("hazards only",  FEATURES_NO_SEASON),
                     ("both",          FEATURES)]:
    before = np.nanmean([fit_one_port(p, data_all, features=feats)["auc"] for p in ports])
    after  = np.nanmean([fit_one_port(p, data,     features=feats)["auc"] for p in ports])
    rows.append({"inputs": label, "with holidays": before,
                 "without holidays": after, "change": after - before})
print(pd.DataFrame(rows).round(3).to_string(index=False))

# the corrected result -- this `cv` and `folds` are what sections 7-10 use
folds, cv = leave_one_port_out(data)
all_truth = np.concatenate([f["truth"] for f in folds])
all_prob  = np.concatenate([f["prob"]  for f in folds])

print("\ncorrected leave-one-port-out:")
print(cv.round(3).to_string(index=False))
print(f"\nmean of the six : {cv['auc'].mean():.3f}   (was {cv_all['auc'].mean():.3f})")
print(f"all days pooled : {roc_auc_score(all_truth, all_prob):.3f}")

dropped 396 holiday days -- 6,381 remain (6% of the data)
disruption rate 13.4%  ->  11.6%

       inputs  with holidays  without holidays  change
calendar only          0.617             0.571  -0.045
 hazards only          0.566             0.556  -0.010
         both          0.667             0.622  -0.045

corrected leave-one-port-out:
      port  n_test  disrupted   auc    lo    hi
Charleston    1065        173 0.684 0.645 0.730
   Houston    1065         81 0.675 0.618 0.741
NewOrleans    1056         71 0.646 0.579 0.732
 NewYorkNJ    1065        157 0.488 0.439 0.540
   Norfolk    1065         95 0.717 0.663 0.762
     Tampa    1065        162 0.520 0.472 0.567

mean of the six : 0.622   (was 0.667)
all days pooled : 0.591

The mean AUC decreases from about 0.67 → 0.62

Is that bad?

No.

The lower score is more trustworthy because we removed an easy but irrelevant shortcut.

fig = plot_lopo(folds, cv, "corrected, holiday window excluded")
fig.savefig(FIG_DIR / "lopo.png", dpi=150)
plt.show()

Where does the model work?

Performance is not equally strong at every port.

The table compares each port's:

• average ships per day,
• number of storm days,
• held-out AUC,
• and uncertainty range.

diag = (data.groupby("port")
        .agg(ships_per_day=("unique_vessels", "mean"),
             storm_days=("storm_day_200km", "sum"))
        .join(cv.set_index("port")[["auc", "lo", "hi"]]))
diag["ships_per_day"] = diag["ships_per_day"].round(0).astype(int)
diag["beats_guessing"] = np.where(diag["lo"] > 0.5, "yes", "no")

print(diag.round(2).sort_values("auc", ascending=False).to_string())

            ships_per_day  storm_days   auc    lo    hi beats_guessing
port                                                                  
Norfolk               246          13  0.72  0.66  0.76            yes
Charleston            134          26  0.68  0.65  0.73            yes
Houston               420          10  0.68  0.62  0.74            yes
NewOrleans            391          16  0.65  0.58  0.73            yes
Tampa                  77          13  0.52  0.47  0.57             no
NewYorkNJ             340           8  0.49  0.44  0.54             no

Four ports perform clearly above guessing:

Norfolk, Charleston, Houston, and New Orleans.

Two ports do not:

• Tampa: only about 77 ships per day, so normal day-to-day variation makes the disruption threshold noisy.
• New York/New Jersey: only 8 storm days, so there are too few hazard examples for the model to learn from.

This suggests that transfer works best when a port has both:

enough vessel traffic for stable daily counts + enough storm events to learn a hazard signal

What did we learn from this section?

The original AUC of 0.67 looked promising, but part of that performance came from learning the holiday shipping slowdown.

After removing that shortcut: mean AUC: 0.67 → 0.62

The lower score is more meaningful because it better reflects the hazard-related prediction problem.

Furthermore, the model works at 4 of 6 ports.

The two failures also tell us something useful about where the method can be applied:

• daily vessel counts must be stable enough to measure disruption,
• and the training data must contain enough storm events to learn from.

Main lesson

Model validation is not only about asking “How accurate is it?”, it is also about asking “Why does it work, and where does it fail?”

---

6. Scaling - one core vs many

Machine-learning models have settings called hyperparameters. To find good settings, we may need to train the model many times.

Here we test:

48 settings × 6 held-out ports = 288 model fits

Each fit is independent, which makes this a good task for parallel computing.

The comparison

We run exactly the same work using:

• one CPU core, and
• multiple CPU cores.

The result shows how independent jobs can be distributed across available computing resources.

if FAST_DEMO:
    GRID = [dict(max_depth=d, learning_rate=lr, n_estimators=n)
            for d in (3, 6) for lr in (0.03, 0.1) for n in (300, 600)]
else:
    GRID = [dict(max_depth=d, learning_rate=lr, n_estimators=n, subsample=s)
            for d in (3, 4, 6, 8) for lr in (0.03, 0.05, 0.1)
            for n in (400, 800) for s in (0.7, 0.9)]

#sized so the serial run takes minutes -- a few seconds of work would be swamped by process start-up and would understate the gain
jobs = [(p, prm) for prm in GRID for p in ports] 
print(f"{len(GRID)} settings x {len(ports)} ports = {len(jobs)} model fits, all independent")

def run_job(port, params):
    return {**params, "port": port, "auc": fit_one_port(port, data, **params)["auc"]}

t0 = time.perf_counter()
_ = Parallel(n_jobs=1)(delayed(run_job)(p, prm) for p, prm in jobs)
t_serial = time.perf_counter() - t0

nw = max(2, N_CORES)
t0 = time.perf_counter()
results = Parallel(n_jobs=nw, backend="loky")(delayed(run_job)(p, prm) for p, prm in jobs)
t_par = time.perf_counter() - t0

print(f" 1 core   : {t_serial:6.1f} s")
print(f"{nw:2d} cores  : {t_par:6.1f} s")
print(f" speed-up : {t_serial / t_par:.2f}x   ({t_serial / len(jobs) * 1000:.0f} ms per fit)")
print("\nThis is a SHARED node, so the figure moves with whoever else is running.")

keys = list(GRID[0])
best = (pd.DataFrame(results).groupby(keys)["auc"].mean()
        .sort_values(ascending=False).reset_index())
print("\nbest settings by mean held-out score:")
print(best.head(5).round(3).to_string(index=False))

BEST = {k: (int(v) if k in ("max_depth", "n_estimators") else float(v))
        for k, v in best.iloc[0][keys].items()}

48 settings x 6 ports = 288 model fits, all independent
 1 core   :  131.8 s
 8 cores  :   22.1 s
 speed-up : 5.97x   (458 ms per fit)

This is a SHARED node, so the figure moves with whoever else is running.

best settings by mean held-out score:
 max_depth  learning_rate  n_estimators  subsample   auc
         3           0.03           400        0.9 0.639
         3           0.03           400        0.7 0.637
         3           0.03           800        0.7 0.628
         4           0.03           400        0.9 0.627
         3           0.05           400        0.9 0.627

High-performance computing is useful not only because datasets are large, but also because scientific analyses often require repeating the same computation hundreds or thousands of times.

This is the same basic idea that made processing the original 871GB of AIS data practical.

---

7. What drives the prediction

What are we asking?

A model can make useful predictions without clearly explaining why. We use SHAP values to examine how different inputs influence the model's predictions. A larger SHAP magnitude means that a variable has a stronger influence on the model's output for that observation.

All days vs. storm days

There is an important complication: Storms occur on only a small fraction of days. If we average feature importance over every day, storm-related variables may look weak simply because most days have no storm.

So we examine feature importance twice:

1. across all days, and
2. across storm days only.

pos = max(int(data[TARGET].sum()), 1)
full_model = make_model(scale_pos_weight=(len(data) - pos) / pos, **BEST)
full_model.fit(data[FEATURES], data[TARGET])

if HAVE_SHAP:
    sv = shap.TreeExplainer(full_model).shap_values(data[FEATURES])
    rank = pd.Series(np.abs(sv).mean(axis=0), index=FEATURES).sort_values(ascending=False)
    plt.figure(figsize=(9, 6))
    shap.summary_plot(sv, data[FEATURES], show=False)
    plt.title("What pushes a day toward 'disrupted'")
    plt.tight_layout()
    plt.savefig(FIG_DIR / "shap.png", dpi=150)
    plt.show()
else:
    sv = None
    rank = pd.Series(full_model.feature_importances_, index=FEATURES).sort_values(ascending=False)
    rank.sort_values().plot.barh(figsize=(8, 5), color="steelblue")
    plt.tight_layout(); plt.show()

print("averaged over ALL days:")
print(rank.round(4).to_string())
print(f"\nstorm distance ranks #{list(rank.index).index('storm_min_dist_km') + 1} "
      f"of {len(rank)} -- mid-table, and the season leads. Read on.")

averaged over ALL days:
doy_sin                    0.3116
doy_cos                    0.2284
wl_anom                    0.1893
wl_anom_max3d              0.1500
storm_min_dist_km          0.1203
wl_anom_lag1               0.0869
wl_anom_lag2               0.0670
storm_max_wind_kt_200km    0.0236
storm_day_200km            0.0003
storm_pts_200km            0.0000

storm distance ranks #5 of 10 -- mid-table, and the season leads. Read on.

if sv is None:
    print("shap is not installed -- skipping. Install with: pip install shap")
else:
    storm_rows = data["storm_day_200km"].to_numpy() == 1

    cmp = pd.DataFrame({
        "all days":        pd.Series(np.abs(sv).mean(axis=0), index=FEATURES),
        "storm days only": pd.Series(np.abs(sv[storm_rows]).mean(axis=0), index=FEATURES),
    })
    cmp["x"] = (cmp["storm days only"] / cmp["all days"].replace(0, np.nan)).round(1)

    print(f"{storm_rows.sum()} storm days out of {len(data):,}")
    print()
    print(cmp.sort_values("storm days only", ascending=False).round(4).to_string())
    print()
    print("On the days that count, the storm inputs dominate and the season does not move")
    print("at all -- confirming the season is constant background, not a storm signal.")

86 storm days out of 6,381

                         all days  storm days only          x
storm_min_dist_km          0.1203           1.4017  11.700000
storm_max_wind_kt_200km    0.0236           0.5966  25.299999
doy_sin                    0.3116           0.3130   1.000000
wl_anom                    0.1893           0.2678   1.400000
wl_anom_max3d              0.1500           0.2116   1.400000
doy_cos                    0.2284           0.2042   0.900000
wl_anom_lag2               0.0670           0.1657   2.500000
wl_anom_lag1               0.0869           0.0691   0.800000
storm_day_200km            0.0003           0.0016   6.000000
storm_pts_200km            0.0000           0.0000        NaN

On the days that count, the storm inputs dominate and the season does not move
at all -- confirming the season is constant background, not a storm signal.

What should you notice?

On all days, seasonal and water-related variables can appear important. When we look specifically at storm days, variables such as storm distance, and storm wind become much more influential.

Therefore:

Feature importance depends on which observations you are asking about.

Does our 200 km cutoff make sense?

We originally defined a storm day using a 200 km radius. Now we can use the trained model to examine whether the data support that choice.

if sv is None:
    print("shap is not installed -- skipping. Install with: pip install shap")
else:
    i = FEATURES.index("storm_min_dist_km")

    sd = pd.DataFrame({"dist_km": data.loc[storm_rows, "storm_min_dist_km"].to_numpy(),
                       "shap": sv[storm_rows][:, i]})
    print("within the 200 km storm window:")
    print(sd.groupby(pd.cut(sd["dist_km"], [0, 50, 100, 150, 200]))["shap"]
            .agg(["count", "mean"]).round(3).to_string())
    print()
    print("Flat. A storm 20 km away and one 190 km away push the prediction by the same")
    print("amount -- so within the window, proximity does NOT matter. The regression was right.")

    # so where does distance actually start to matter?
    a = pd.DataFrame({"dist": data["storm_min_dist_km"].to_numpy(), "shap": sv[:, i]})
    a = a[a["dist"] < 4500]        # drop the 5000 "no Atlantic storm at all" sentinel

    plt.figure(figsize=(7.5, 4.5))
    plt.scatter(a["dist"], a["shap"], alpha=0.3, s=12, color="steelblue")
    plt.axvline(200, color="r", ls="--", lw=1, label="our 200 km cutoff")
    plt.axhline(0, color="k", lw=1)
    plt.xscale("log")
    plt.xlabel("distance to nearest storm (km, log scale)")
    plt.ylabel("push toward 'disrupted'")
    plt.title("Where the model actually draws the line")
    plt.legend(); plt.grid(alpha=0.3); plt.tight_layout()
    plt.savefig(FIG_DIR / "distance_effect.png", dpi=150)
    plt.show()

    for lo, hi in [(0, 200), (200, 500), (500, 1000), (1000, 4500)]:
        m = a["dist"].between(lo, hi)
        print(f"{lo:>5}-{hi:<5} km : {m.sum():>5} days, mean push {a.loc[m, 'shap'].mean():+.3f}")

within the 200 km storm window:
            count   mean
dist_km                 
(0, 50]        14  1.387
(50, 100]      24  1.347
(100, 150]     18  1.449
(150, 200]     30  1.424

Flat. A storm 20 km away and one 190 km away push the prediction by the same
amount -- so within the window, proximity does NOT matter. The regression was right.

    0-200   km :    86 days, mean push +1.402
  200-500   km :   183 days, mean push +0.708
  500-1000  km :   420 days, mean push -0.078
 1000-4500  km :  1446 days, mean push -0.161

Within 200 km, the effect of storm distance is surprisingly flat. A storm 20 km away and a storm near the edge of the 200 km window can influence the model similarly.

When we look farther away, the storm-related signal declines gradually rather than disappearing exactly at 200 km.

Therefore:

The model can help us question assumptions that we originally put into the analysis.

And the model finds signal out to about 1,000 km; our cutoff is probably too tight, a finding only the model could surface.

---

8. What sea-level rise does

In this section, we ask:

How often would today's unusually high water levels occur if average sea level were higher?

For each port, we define its high-water mark as the water level reached on roughly the highest 5% of observed days. So under present conditions, the high-water mark is reached about 5% of the time.

A. Sea-level-rise scenario

We then add hypothetical sea-level rise to the historical observations. For example: What if every observed water level were 0.25 m higher? We then count how often the port would cross its current high-water mark.

What should you notice?

The same amount of sea-level rise does not affect every port equally. For example, under the +0.25 m scenario, Tampa reaches its current high-water mark far more frequently than it does today. The differences depend on each port's local water-level distribution, not simply its latitude.

B. Stronger-storm scenario

The notebook also asks how additional storm surge changes the model's disruption-risk score. This part does use the trained model.

So keep the two ideas separate:

• Higher seas: change how often today's high-water level is reached.
• Stronger surge: change the model's predicted disruption risk.

# --- A. higher seas: no model, absolute water level -----------------------
# Each port's HIGH-WATER MARK: the level it reaches on its worst 5% of days.
SLR_M = [0.0, 0.25, 0.50, 1.00]
high_water_mark = data.groupby("port")["wl_mean"].quantile(0.95)

emp = (pd.DataFrame([{"slr_m": s, "port": port,
                      "share": float(((g["wl_mean"] + s) >= high_water_mark[port]).mean())}
                     for s in SLR_M for port, g in data.groupby("port")])
       .pivot(index="port", columns="slr_m", values="share"))

print("A. HIGHER SEAS -- share of days reaching each port's high-water mark")
print("   (columns are metres of sea-level rise added)")
print(emp.round(3).to_string())

print("\neach port's high-water mark, m above mean sea level:")
print(high_water_mark.round(3).to_string())


for port in emp.index:
    print(f"  {port:<11} 5% of days  ->  {emp.loc[port, 0.25]*100:.0f}% at +0.25 m"
          f"   ({emp.loc[port, 0.25] / emp.loc[port, 0.0]:.0f}x more often)")

# --- B. stronger storms: model-based, inside the observed range ------------
print("\n\nobserved surge (m) -- the model knows nothing beyond this:")
print(data["wl_anom"].quantile([0.5, 0.9, 0.99, 1.0]).round(3).to_string())

SURGE_M = [0.0, 0.05, 0.10, 0.20]
WL_FEATURES = ["wl_anom", "wl_anom_lag1", "wl_anom_lag2", "wl_anom_max3d"]

rows = []
for s in SURGE_M:
    scen = data.copy()
    scen[WL_FEATURES] = scen[WL_FEATURES] + s
    scen["prob"] = full_model.predict_proba(scen[FEATURES])[:, 1]
    for port, g in scen.groupby("port"):
        rows.append({"surge_m": s, "port": port, "mean_prob": g["prob"].mean()})

scen_df = pd.DataFrame(rows)
pivot = scen_df.pivot(index="port", columns="surge_m", values="mean_prob")

print("\nB. STRONGER STORMS -- change in disruption risk if surges run +N m larger (%)")
print(((pivot.div(pivot[0.0], axis=0) - 1) * 100).round(1).to_string())

A. HIGHER SEAS -- share of days reaching each port's high-water mark
   (columns are metres of sea-level rise added)
slr_m        0.00   0.25   0.50   1.00
port                                  
Charleston  0.051  0.380  0.920  1.000
Houston     0.051  0.403  0.914  0.998
NewOrleans  0.050  0.270  0.792  1.000
NewYorkNJ   0.051  0.480  0.900  0.998
Norfolk     0.051  0.329  0.880  1.000
Tampa       0.051  0.581  0.963  1.000

each port's high-water mark, m above mean sea level:
port
Charleston    0.466
Houston       0.482
NewOrleans    0.533
NewYorkNJ     0.413
Norfolk       0.503
Tampa         0.366
  Charleston  5% of days  ->  38% at +0.25 m   (8x more often)
  Houston     5% of days  ->  40% at +0.25 m   (8x more often)
  NewOrleans  5% of days  ->  27% at +0.25 m   (5x more often)
  NewYorkNJ   5% of days  ->  48% at +0.25 m   (9x more often)
  Norfolk     5% of days  ->  33% at +0.25 m   (6x more often)
  Tampa       5% of days  ->  58% at +0.25 m   (11x more often)


observed surge (m) -- the model knows nothing beyond this:
0.50    0.008
0.90    0.187
0.99    0.422
1.00    1.299

B. STRONGER STORMS -- change in disruption risk if surges run +N m larger (%)
surge_m     0.00  0.05  0.10       0.20
port                                   
Charleston   0.0   2.8   7.8  13.600000
Houston      0.0   5.3  10.2  17.200001
NewOrleans   0.0   5.2   9.0  16.100000
NewYorkNJ    0.0   3.7   8.4  14.500000
Norfolk      0.0   4.1   9.1  15.500000
Tampa        0.0   4.3  11.9  20.400000

Main lesson

The same environmental change can produce very different operational exposure at different locations.

SHOW_SLR = 0.25       # +0.50 saturates: every port lands 0.74-0.91 and the map goes flat

fig, axes = plt.subplots(1, 3, figsize=(17, 4.5))

ax = axes[0]
for port in emp.index:
    ax.plot(emp.columns, emp.loc[port], marker="o", label=port)
ax.set_xlabel("sea-level rise (m)")
ax.set_ylabel("share of days reaching the high-water mark")
ax.set_title("A. Higher seas (no model)"); ax.legend(fontsize=8); ax.grid(alpha=0.3)

ax = axes[1]
rel = (pivot.div(pivot[0.0], axis=0) - 1) * 100
for port in rel.index:
    ax.plot(rel.columns, rel.loc[port], marker="o", label=port)
ax.set_xlabel("added surge (m)"); ax.set_ylabel("increase in disruption risk (%)")
ax.set_title("B. Stronger storms (model)"); ax.legend(fontsize=8); ax.grid(alpha=0.3)

# lat/lon are used here for DRAWING only -- they were kept out of the model
ax = axes[2]
latlon = data.groupby("port")[["lat", "lon"]].first()
delta = emp[SHOW_SLR] - emp[0.0]
sc = ax.scatter(latlon["lon"], latlon["lat"], c=delta[latlon.index],
                s=300, cmap="OrRd", edgecolor="k", zorder=3)
# print the value beside each port: colour scales rescale to whatever range they
# are given, so two very different maps can look identical
for p, r in latlon.iterrows():
    ax.annotate(f"{p}\n{delta[p]:.2f}", (r["lon"], r["lat"]),
                textcoords="offset points", xytext=(9, -6), fontsize=8)
plt.colorbar(sc, ax=ax, label=f"extra share of days at +{SHOW_SLR:.2f} m")
ax.set_xlabel("longitude"); ax.set_ylabel("latitude")
ax.set_xlim(latlon["lon"].min() - 4, latlon["lon"].max() + 6)
ax.set_title(f"C. Where +{SHOW_SLR:.2f} m affects the most"); ax.grid(alpha=0.3)

plt.tight_layout()
plt.savefig(FIG_DIR / "scenarios.png", dpi=150)
plt.show()

print(f"extra share of days at the high-water mark, +{SHOW_SLR:.2f} m, most to least:")
print(delta.sort_values(ascending=False).round(3).to_string())
print("\nNo north-south pattern. Exposure follows each port's own high-water mark -- Tampa 0.37 m, New Orleans 0.53 m -- which reflects local tidal range, not latitude.")

extra share of days at the high-water mark, +0.25 m, most to least:
port
Tampa         0.531
NewYorkNJ     0.429
Houston       0.352
Charleston    0.330
Norfolk       0.278
NewOrleans    0.220

No north-south pattern. Exposure follows each port's own high-water mark -- Tampa 0.37 m, New Orleans 0.53 m -- which reflects local tidal range, not latitude.

# ============================================================
# Interactive map: which ports sea-level rise reaches first
# ============================================================

import folium
from IPython.display import display

# Esri basemap
TILES = (
    "https://server.arcgisonline.com/ArcGIS/rest/services/"
    "World_Street_Map/MapServer/tile/{z}/{y}/{x}"
)

ATTR = (
    "Tiles © Esri — Sources: Esri, DeLorme, NAVTEQ, "
    "USGS, Intermap, iPC, NRCAN, Esri Japan, METI, "
    "Esri China (Hong Kong), Esri Thailand, TomTom"
)

# Years represented in the dataset
SPAN = f"{data['date'].min():%Y}\u2013{data['date'].max():%Y}"

# Difference between exposure under selected sea-level rise and the present-day baseline
exposure = emp[SHOW_SLR] - emp[0.0]

if "port" in cv.columns:
    cvx = cv.set_index("port")
else:
    cvx = cv.copy()

m2 = folium.Map(
    location=[34, -84],
    zoom_start=5,
    tiles=None
)
folium.TileLayer(
    tiles=TILES,
    attr=ATTR,
    name="Esri World Street Map",
    overlay=False,
    control=True
).add_to(m2)
for port, r in data.groupby("port")[["lat", "lon"]].first().iterrows():

    e = exposure[port]
    s = cvx.loc[port]
    usable = s["lo"] > 0.5

    folium.CircleMarker(
        location=[r["lat"], r["lon"]],
        radius=10 + 30 * e,
        color="black",
        weight=1,
        fill=True,
        fill_opacity=0.80,
        fill_color=(
            "#b30000" if e > 0.45
            else "#fc8d59" if e > 0.30
            else "#fee8c8"
        ),
        tooltip=(
            f"{port}: "
            f"{emp.loc[port, SHOW_SLR] * 100:.0f}% of days "
            f"at +{SHOW_SLR:.2f} m"
        ),
        popup=folium.Popup(
            f"<b>{port}</b><br>"
            f"High-water mark <b>{high_water_mark[port]:.2f} m</b><br>"
            f"Reached on 5% of days in {SPAN}<br>"
            f"\u2192 <b>{emp.loc[port, SHOW_SLR] * 100:.0f}% of days</b> "
            f"at +{SHOW_SLR:.2f} m<br><hr>"
            f"Model score <b>{s['auc']:.2f}</b> "
            f"({s['lo']:.2f}\u2013{s['hi']:.2f})<br>"
            f"<i>{'usable here' if usable else 'unproven here'}</i>",
            max_width=250
        )
    ).add_to(m2)
folium.LayerControl().add_to(m2)

print(
    f"Circle size and colour = how much more often the high-water "
    f"mark is reached at +{SHOW_SLR:.2f} m of sea-level rise."
)

print(
    f"Baseline is 5% of days across {SPAN}. "
    "Click a port for its numbers."
)

display(m2)

Circle size and colour = how much more often the high-water mark is reached at +0.25 m of sea-level rise.
Baseline is 5% of days across 2019–2024. Click a port for its numbers.

Make this Notebook Trusted to load map: File -> Trust Notebook

---

9. Saving the model

A model file on its own is not reusable. The next person needs to know which columns to feed it and in what order, what "disrupted" was defined as, which ports it can be trusted at, and where it must not be used at all. All of that gets saved alongside the model.

# 9.1  Save the model together with everything needed to use it correctly.
import joblib

MODEL_DIR = WORK / "models"
MODEL_DIR.mkdir(parents=True, exist_ok=True)

cv_idx = cv.set_index("port") # leave-one-port-out scores from section 4

artifact = {
    "model":                full_model,
    "features":             FEATURES,  # order matters at predict time
    "target":               TARGET,
    "disruption_threshold": DISRUPTION_THRESHOLD,
    "trained_on_ports":     ports,
    "excluded_window":      "18 Dec - 8 Jan (holiday slowdown, not hazard-driven)",
    "storm_radius_km":      200,
    "cv_scores":            cv_idx[["auc", "lo", "hi"]].round(3).to_dict("index"),
    "trustworthy_at":       [p for p in ports if cv_idx.loc[p, "lo"] > 0.5],
    "high_water_mark_m":    high_water_mark.round(3).to_dict(),   # from section 8
    "trained_utc":          pd.Timestamp.utcnow().isoformat(),
    "predicts":             ("P(port below 90% of its own 30-day median ship count) "
                             "from hazard + season inputs only"),
    "scope":                ("North Atlantic storms only -- not valid for Pacific ports. "
                             "Unproven at ports under ~100 ships/day."),
}

path = MODEL_DIR / "port_disruption_model.joblib"
joblib.dump(artifact, path)
print(f"saved {path.name}  ({path.stat().st_size / 1024:.0f} KB)\n")

saved port_disruption_model.joblib  (458 KB)

(Optional) Reload model and make example predictions

Look for: the difference between an in-sample score and a held-out score from a model that had never seen that port.

# Reload the model to view the file descriptions
bundle = joblib.load(path)
unproven = [p for p in ports if p not in bundle["trustworthy_at"]]

print("predicts        :", bundle["predicts"])
print("inputs          :", len(bundle["features"]), "columns, order fixed")
print("disrupted means : below "
      f"{bundle['disruption_threshold']:.0%} of the port's own normal")
print("trustworthy at  :", ", ".join(bundle["trustworthy_at"]))
print("unproven at     :", ", ".join(unproven) if unproven else "-")
print("scope           :", bundle["scope"])

predicts        : P(port below 90% of its own 30-day median ship count) from hazard + season inputs only
inputs          : 10 columns, order fixed
disrupted means : below 90% of the port's own normal
trustworthy at  : Charleston, Houston, NewOrleans, Norfolk
unproven at     : NewYorkNJ, Tampa
scope           : North Atlantic storms only -- not valid for Pacific ports. Unproven at ports under ~100 ships/day.

# 9.2  The honest (held-out) score for specific port-days.
# full_model was trained on all six ports, so its scores are in-sample. The leave-one-port-out folds hold predictions from a model that had NOT seen that port
fold_by_port = {f["port"]: f for f in folds}

def held_out_score(port, date):
    f = fold_by_port[port]
    rows = data[data["port"] == port].reset_index(drop=True)
    hit = rows.index[rows["date"] == pd.Timestamp(date)]
    return float(f["prob"][hit[0]]) if len(hit) else np.nan


# three storm days, plus one ordinary day at the same port for contrast
calm = (data[(data["port"] == "NewOrleans") & (data["storm_day_200km"] == 0)]
        .sample(1, random_state=0).iloc[0])

CASES = [("NewOrleans", "2020-10-28"),
         ("NewOrleans", "2019-10-26"),
         ("Houston",    "2019-09-18"),
         ("NewOrleans", calm["date"].strftime("%Y-%m-%d"))]

for port, date in CASES:
    sel = data[(data["port"] == port) & (data["date"] == pd.Timestamp(date))]
    if sel.empty:
        print(f"{port} {date}: not in the panel\n")
        continue

    r = sel.iloc[0]
    # predict from the DataFrame, not from a row Series -- a Series taken from a
    # mixed-type frame is object dtype, and xgboost rejects it
    in_sample = full_model.predict_proba(sel[FEATURES])[0, 1]

    print(f"{port} {date}")
    print(f"   wind {r['storm_max_wind_kt_200km']:.0f} kt, "
          f"{r['storm_min_dist_km']:.0f} km away, surge {r['wl_anom']:+.2f} m")
    print(f"   {r['unique_vessels']:.0f} ships = {r['activity_ratio']:.2f} of normal"
          f"  ->  quiet_day = {int(r[TARGET])}")
    print(f"   in-sample {in_sample:.3f}   |   held-out {held_out_score(port, date):.3f}")
    print()

NewOrleans 2020-10-28
   wind 100 kt, 98 km away, surge +0.14 m
   316 ships = 0.86 of normal  ->  quiet_day = 1
   in-sample 0.951   |   held-out 0.924

NewOrleans 2019-10-26
   wind 45 kt, 51 km away, surge +0.06 m
   373 ships = 0.99 of normal  ->  quiet_day = 0
   in-sample 0.485   |   held-out 0.572

Houston 2019-09-18
   wind 30 kt, 55 km away, surge +0.16 m
   411 ships = 0.96 of normal  ->  quiet_day = 0
   in-sample 0.280   |   held-out 0.100

NewOrleans 2022-02-27
   wind 0 kt, 5000 km away, surge +0.03 m
   399 ships = 1.00 of normal  ->  quiet_day = 0
   in-sample 0.196   |   held-out 0.149

A high score means the day looks more like previously observed disruption conditions, not that the port must shut down.

The New Orleans storm on October 28, 2020 had 100 kt winds, was 98 km away, and ship activity fell to 86% of normal, producing a high held-out score.

Similar storms can still produce different operational responses.

---

(Optional) Putting the prediction on the map

The interactive map shows:

• the storm track,
• the port,
• the 200 km storm radius,
• and the 20 km ship-count radius.

Click the track and port markers to inspect the event and model output.

The map connects the prediction back to the actual spatial event.

# ============================================================================
# Interactive map: the storm, the port, and the two radii that define the data
# ============================================================================

try:
    import folium
    HAVE_FOLIUM = True
except ImportError:
    HAVE_FOLIUM = False
    print("folium not installed: pip install folium")

TILES = (
    "https://server.arcgisonline.com/ArcGIS/rest/services/"
    "World_Street_Map/MapServer/tile/{z}/{y}/{x}"
)

ATTR = (
    "Tiles © Esri — Sources: Esri, DeLorme, NAVTEQ, "
    "USGS, Intermap, iPC, NRCAN, Esri Japan, METI, "
    "Esri China (Hong Kong), Esri Thailand, TomTom"
)

SPAN = f"{data['date'].min():%Y}\u2013{data['date'].max():%Y}"

EVENT_PORT = "NewOrleans"
EVENT_DATE = "2020-10-28"

ib_path = OUT_DIR / "ibtracs.NA.list.v04r01.csv"


def _hav_km(lon, lat, lon0, lat0):
    """Calculate great-circle distance in kilometres."""

    R = 6371.0088

    lon = np.deg2rad(np.asarray(lon, dtype=float))
    lat = np.deg2rad(np.asarray(lat, dtype=float))

    lon0 = np.deg2rad(float(lon0))
    lat0 = np.deg2rad(float(lat0))

    a = (
        np.sin((lat - lat0) / 2) ** 2
        + np.cos(lat)
        * np.cos(lat0)
        * np.sin((lon - lon0) / 2) ** 2
    )

    return 2 * R * np.arcsin(np.sqrt(a))

# Store each port's held-out validation fold
fold_by_port = {f["port"]: f for f in folds}


def held_out_score(port, date):
    """Return the prediction from the fold where the port was not in training."""

    rows = data[data["port"] == port].reset_index(drop=True)
    hit = rows.index[rows["date"] == pd.Timestamp(date)]

    if len(hit):
        return float(fold_by_port[port]["prob"][hit[0]])

    return np.nan


if not HAVE_FOLIUM:
    pass

elif not ib_path.exists():
    print(f"Storm tracks not found at {ib_path}")

else:
    # Get the port coordinates
    plat, plon = data.loc[
        data["port"] == EVENT_PORT,
        ["lat", "lon"]
    ].iloc[0]

    avail = pd.read_csv(
        ib_path,
        nrows=0,
        skiprows=[1]
    ).columns

    want = [
        c for c in
        ["SID", "NAME", "ISO_TIME", "LAT", "LON", "USA_WIND"]
        if c in avail
    ]

    ib = pd.read_csv(
        ib_path,
        skiprows=[1],
        low_memory=False,
        usecols=want
    )

    ib["ISO_TIME"] = pd.to_datetime(
        ib["ISO_TIME"],
        errors="coerce"
    )

    for c in ["LAT", "LON", "USA_WIND"]:
        ib[c] = pd.to_numeric(
            ib[c],
            errors="coerce"
        )

    ib = ib.dropna(
        subset=["ISO_TIME", "LAT", "LON"]
    )

    # Find the storm that came nearest the port during the event window
    d0 = pd.Timestamp(EVENT_DATE)

    win = ib[
        ib["ISO_TIME"].between(
            d0 - pd.Timedelta(days=6),
            d0 + pd.Timedelta(days=3)
        )
    ].copy()

    win["dist"] = _hav_km(
        win["LON"],
        win["LAT"],
        plon,
        plat
    )

    sid = win.loc[
        win["dist"].idxmin(),
        "SID"
    ]

    # Extract the complete track for the selected storm
    track = (
        ib[ib["SID"] == sid]
        .sort_values("ISO_TIME")
    )

    if "NAME" in track.columns:
        name = str(track["NAME"].iloc[0]).title()
    else:
        name = str(sid)

    m = folium.Map(
        location=[plat - 2, plon],
        zoom_start=6,
        tiles=None
    )

    # Add the Esri basemap
    folium.TileLayer(
        tiles=TILES,
        attr=ATTR,
        name="Esri World Street Map",
        overlay=False,
        control=True
    ).add_to(m)

    # Add the two radii that define the dataset
    folium.Circle(
        location=[plat, plon],
        radius=200_000,
        color="#c00",
        weight=2,
        dash_array="6",
        fill=False,
        tooltip=(
            "200 km \u2014 inside this, the day counts "
            "as a storm day"
        )
    ).add_to(m)

    folium.Circle(
        location=[plat, plon],
        radius=20_000,
        color="#06c",
        weight=1,
        fill=True,
        fill_opacity=0.30,
        tooltip="20 km \u2014 where we count ships"
    ).add_to(m)

    # Add the storm track
    folium.PolyLine(
        track[["LAT", "LON"]].values,
        color="#444",
        weight=2,
        tooltip=f"{name} storm track"
    ).add_to(m)

    # Add storm observations along the track
    for _, t in track.iterrows():

        if pd.notna(t.get("USA_WIND")):
            w = t["USA_WIND"]
        else:
            w = 0

        if w >= 64:
            point_color = "#d7301f"
        elif w >= 34:
            point_color = "#fdae61"
        else:
            point_color = "#9ecae1"

        folium.CircleMarker(
            location=[t["LAT"], t["LON"]],
            radius=3 + w / 18,
            color=None,
            fill=True,
            fill_opacity=0.85,
            fill_color=point_color,
            popup=(
                f"{t['ISO_TIME']:%Y-%m-%d %H:%M}<br>"
                f"{w:.0f} kt"
            )
        ).add_to(m)

    # Get port information for the selected event date
    r = data[
        (data["port"] == EVENT_PORT)
        & (data["date"] == d0)
    ].iloc[0]

    # Add the port marker and event details
    folium.Marker(
        location=[plat, plon],
        icon=folium.Icon(
            color="darkblue",
            icon="ship",
            prefix="fa"
        ),
        popup=folium.Popup(
            f"<b>Port of {EVENT_PORT}</b><br>"
            f"{EVENT_DATE}<br><hr>"
            f"{r['unique_vessels']:.0f} ships = "
            f"<b>{r['activity_ratio']:.2f}</b> of normal<br>"
            f"Storm {r['storm_min_dist_km']:.0f} km away at "
            f"{r['storm_max_wind_kt_200km']:.0f} kt<br>"
            f"Water {r['wl_anom']:+.2f} m<br><hr>"
            f"Held-out score "
            f"<b>{held_out_score(EVENT_PORT, EVENT_DATE):.3f}</b><br>"
            f"<i>From a model never trained on this port</i>",
            max_width=280
        )
    ).add_to(m)

    # Add a layer control
    folium.LayerControl().add_to(m)

    print(
        f"{name} \u2014 red = hurricane force, "
        "orange = tropical storm, blue = weaker"
    )

    print(
        "Click the track points for time and wind; "
        "click the ship for the port's day."
    )

    display(m)

Zeta — red = hurricane force, orange = tropical storm, blue = weaker
Click the track points for time and wind; click the ship for the port's day.

Make this Notebook Trusted to load map: File -> Trust Notebook

---

10. Optional - does another method agree?

Gradient boosting builds trees one after another; a Random Forest builds them independently and votes.Main idea is that if different algorithms reach similar conclusions, the pattern is more likely to come from the data rather than one specific model.

# ANSWER
def fit_one_port_rf(held_out, frame, features=FEATURES):
    train = frame[frame["port"] != held_out]
    test  = frame[frame["port"] == held_out]
    rf = RandomForestClassifier(n_estimators=300, max_depth=8, class_weight="balanced",
                                random_state=42, n_jobs=N_CORES)
    rf.fit(train[features], train[TARGET])
    prob = rf.predict_proba(test[features])[:, 1]
    truth = test[TARGET].to_numpy()
    auc = roc_auc_score(truth, prob) if truth.min() != truth.max() else np.nan
    return {"port": held_out, "auc": auc, "model": rf}


rf_folds = [fit_one_port_rf(p, data) for p in ports]
rf_cv = pd.DataFrame([{"port": f["port"], "auc_rf": f["auc"]} for f in rf_folds])

compare = cv[["port", "auc"]].rename(columns={"auc": "auc_boosted"}).merge(rf_cv, on="port")
compare["gap"] = (compare["auc_rf"] - compare["auc_boosted"]).round(3)
print(compare.round(3).to_string(index=False))
print(f"\nmean   boosted {compare['auc_boosted'].mean():.3f}   "
      f"random forest {compare['auc_rf'].mean():.3f}")

rf_imp = pd.Series(rf_folds[0]["model"].feature_importances_, index=FEATURES)\
           .sort_values(ascending=False)
print("\ntop inputs")
print(pd.DataFrame({"random_forest": rf_imp.head(5).index,
                    "boosted":       rank.head(5).index}).to_string(index=False))
print("\nSame ports succeed, same ports fail, same inputs lead. The result is a property of the data, not of the algorithm.")

      port  auc_boosted  auc_rf   gap
Charleston        0.684   0.695 0.011
   Houston        0.675   0.680 0.005
NewOrleans        0.646   0.666 0.020
 NewYorkNJ        0.488   0.524 0.036
   Norfolk        0.717   0.743 0.027
     Tampa        0.520   0.527 0.007

mean   boosted 0.622   random forest 0.639

top inputs
    random_forest           boosted
          doy_cos           doy_sin
    wl_anom_max3d           doy_cos
          doy_sin           wl_anom
storm_min_dist_km     wl_anom_max3d
     wl_anom_lag1 storm_min_dist_km

Same ports succeed, same ports fail, same inputs lead. The result is a property of the data, not of the algorithm.

---

Summary

• Nearby storms matter: when a storm is within 200 km, ports handle about 15% fewer ships on average. Stronger winds are associated with larger declines.

• High water may matter too: the estimated effect is negative, but the evidence is not strong enough to clearly separate it from no effect.

• The model transfers to 4 of 6 ports: it works best where there is enough daily vessel traffic and enough storm examples to learn from.

• Sea-level rise changes exposure quickly: with +0.25 m, Tampa would reach today's high-water mark on about 58% of days instead of 5%.

Main lesson

Do not trust the first good-looking result.
We found two misleading patterns — seasonality in the regression and the Christmas slowdown in the prediction model — and corrected both before interpreting the results.

Who this is for

• Port authorities — understand which hazards disrupt operations
• Coastal planners — assess changing high-water exposure
• Researchers — reuse the workflow for AIS-based port resilience studies

Next steps

• include winter storms,
• test storm distances beyond 200 km,
• and add more ports.