95 KiB
95 KiB
In [9]:
import os
from pathlib import Path
import pandas as pd
from dotenv import load_dotenv
from sqlalchemy import create_engine, text
load_dotenv(Path.cwd() / ".env")
engine = create_engine(
"postgresql+psycopg2://{user}:{pw}@{host}:{port}/{db}".format(
user=os.environ["PG_USER"],
pw=os.environ["PG_PASSWORD"],
host=os.environ["PG_HOST"],
port=os.environ["PG_PORT"],
db=os.environ["PG_DB"],
)
)
# LZ_attractions_io carries one row per (timestamp, ride) and uses the
# Item _id that matches both ride_master and LZ_attractions_io_poi.
# (LZ_attractions_io_queuetimes is a same-data duplicate keyed by QueueLine
# id, which doesn't join to anything — don't use it.)
QUERY = text(
"""
SELECT a.time_stamp,
a._id AS ride_id,
r.ride AS ride_name,
r.capacity,
r.duration,
a.queue_time AS queue_time_sec
FROM knoebels."LZ_attractions_io" a
JOIN knoebels.ride_master r
ON r.attractions_dot_io_id = a._id
WHERE a.queue_time IS NOT NULL
ORDER BY a.time_stamp
"""
)
df = pd.read_sql(QUERY, engine)
df["queue_time_min"] = df["queue_time_sec"] / 60.0
# Scraper writes naive timestamps that represent America/New_York wall clock.
# Localize so downstream time-of-day / day-of-week features are unambiguous.
df["time_stamp"] = (
pd.to_datetime(df["time_stamp"])
.dt.tz_localize("America/New_York", ambiguous="NaT", nonexistent="NaT")
)
df = df.dropna(subset=["time_stamp"]).reset_index(drop=True)
def get_ride_time_matrix(df: pd.DataFrame, freq: str = "5min") -> pd.DataFrame:
"""Pivot the long df into a (timestamp x ride_name) matrix of queue times.
Rows are snapped to `freq` bins so all rides share a common index.
Values are mean queue time in minutes within each bin.
"""
snapped = df.assign(ts=df["time_stamp"].dt.floor(freq))
return (
snapped.pivot_table(
index="ts",
columns="ride_name",
values="queue_time_min",
aggfunc="mean",
)
.sort_index().ffill().fillna(0) # Impute missing values: forward fill first (assume line stayed the same), then fill remaining NaNs (like morning startup) with 0
)
print(f"Loaded {len(df):,} rows across {df['ride_name'].nunique()} rides")
print(f"Time range: {df['time_stamp'].min()} → {df['time_stamp'].max()}")
Loaded 635,281 rows across 57 rides Time range: 2024-05-15 12:23:34.954315-04:00 → 2026-05-26 16:35:16.363800-04:00
In [ ]:
'''
if we were to run SVD on raw minutes, the algorithm will just chase the rides
with the biggest raw numbers.
Frist, need to standardize the data so every ride has a mean of 0 and a variance of 1.
scikit-learn makes this trivial
'''
from sklearn.preprocessing import StandardScaler
# Rides to exclude from the PCA:
# powersurge — torn down years ago, ~90% of its rows are imputed zeros and
# its loadings are pure artifact.
EXCLUDED_RIDES = ["powersurge"]
df_matrix = get_ride_time_matrix(df).drop(columns=EXCLUDED_RIDES, errors="ignore")
# Standardize the columns (rides)
scaler = StandardScaler()
scaled_matrix = scaler.fit_transform(df_matrix)
pd.DataFrame(scaled_matrix, index=df_matrix.index, columns=df_matrix.columns).head()
In [12]:
from sklearn.decomposition import PCA
# Let's extract the top 3 components ("Park Mood", "Capacity", etc.)
pca = PCA(n_components=3)
principal_components = pca.fit_transform(scaled_matrix)
# Look at how much variance (signal) each component explains
print(f"Variance explained by PC1 (Park Mood): {pca.explained_variance_ratio_[0]:.1%}")
print(f"Variance explained by PC2: {pca.explained_variance_ratio_[1]:.1%}")
print(f"Variance explained by PC3: {pca.explained_variance_ratio_[2]:.1%}")Variance explained by PC1 (Park Mood): 40.0% Variance explained by PC2: 12.9% Variance explained by PC3: 7.0%
In [14]:
# Interpret each PC by looking at (1) which rides load onto it and (2) when its score is high/low.
loadings = pd.DataFrame(
pca.components_.T,
index=df_matrix.columns,
columns=["PC1", "PC2", "PC3"],
)
scores = pd.DataFrame(
principal_components,
index=df_matrix.index,
columns=["PC1", "PC2", "PC3"],
)
# Top-magnitude rides per PC — these are the ones "shaping" the axis.
for pc in loadings.columns:
print(f"\n=== {pc} — top 8 rides by |loading| ===")
top = loadings[pc].reindex(loadings[pc].abs().sort_values(ascending=False).index).head(8)
print(top.to_string())
# Daily-mean PC scores over time — reveals seasonality each PC captures.
scores.resample("D").mean().plot(subplots=True, figsize=(12, 6), title=["PC1", "PC2", "PC3"])
Out [14]:
=== PC1 — top 8 rides by |loading| === ride_name stratosfear 0.179701 merry_mixer 0.173788 flyer 0.168306 s&g_carousel 0.168079 tornado 0.166412 kiddie_whip 0.166132 spanish_bambini 0.165741 red_baron 0.165697 === PC2 — top 8 rides by |loading| === ride_name bumper_cars 0.230985 giant_wheel 0.209181 kiddie_boats -0.196950 pete's_fleet -0.195781 spanish_bambini -0.195189 s&g_carousel -0.195019 red_baron -0.190284 helicopters -0.189507 === PC3 — top 8 rides by |loading| === ride_name paratrooper 0.335514 kiddie_bumper_cars -0.327763 motor_boats -0.244509 powersurge 0.232955 tumbling_timbers -0.207141 roto_jets 0.199833 stratosfear 0.184504 tornado 0.180223
array([<Axes: title={'center': 'PC1'}, xlabel='ts'>,
<Axes: title={'center': 'PC2'}, xlabel='ts'>,
<Axes: title={'center': 'PC3'}, xlabel='ts'>], dtype=object)
In [ ]:
# Join ride metadata (area + category) to the loadings so the labels write themselves.
# NOTE: attractions.io category names carry a trailing space ("Kiddie ", "Family ",
# "Thrill "), so trim() here — otherwise downstream == / .isin() / dict lookups on
# the bare names silently miss (e.g. the PC2×PC3 scatter colors all fall through).
metadata = pd.read_sql(
text(
"""
SELECT r.ride AS ride_name,
trim(c.name) AS category,
p.parent AS area
FROM knoebels.ride_master r
JOIN knoebels."LZ_attractions_io_poi" p
ON p._id = r.attractions_dot_io_id
LEFT JOIN knoebels."LZ_attractions_io_categories" c
ON c._id = p.category
"""
),
engine,
).set_index("ride_name")
annotated = loadings.join(metadata)
print("=== Mean loading by area ===")
print(annotated.groupby("area")[["PC1", "PC2", "PC3"]].mean().round(3))
print("\n=== Mean loading by category ===")
print(annotated.groupby("category")[["PC1", "PC2", "PC3"]].mean().round(3))
In [ ]:
# Validate PC interpretations from the time side.
# PC1 (busyness) should peak weekend afternoons.
# PC2 (Kiddieland vs main midway) should fade in the evening as the kids leave.
# PC3 (daypart) should run positive midday — school/daycare group rates flood
# the daytime family rides — and go negative after dark, when the evening
# coaster/bumper-car/Ferris-wheel crowd takes over.
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
scores.groupby(scores.index.hour).mean().plot(ax=axes[0])
axes[0].set_title("Mean PC score by hour of day")
axes[0].set_xlabel("Hour (local)")
axes[0].axhline(0, color="grey", linewidth=0.5)
scores.groupby(scores.index.dayofweek).mean().plot(ax=axes[1])
axes[1].set_title("Mean PC score by day of week (0=Mon, 6=Sun)")
axes[1].set_xlabel("Day of week")
axes[1].axhline(0, color="grey", linewidth=0.5)
plt.tight_layout()
In [ ]:
# Low-rank reconstruction: how much of any one ride's queue is explained by the top 3 PCs?
# pca.inverse_transform takes us back to standardized space; scaler.inverse_transform back to minutes.
reconstructed = pd.DataFrame(
scaler.inverse_transform(pca.inverse_transform(principal_components)),
index=df_matrix.index,
columns=df_matrix.columns,
)
# Per-ride R^2 of the rank-3 model: 1 - SS_residual / SS_total.
ss_resid = ((df_matrix - reconstructed) ** 2).sum()
ss_total = ((df_matrix - df_matrix.mean()) ** 2).sum()
ride_r2 = (1 - ss_resid / ss_total).rename("rank3_R2")
print("=== Best explained by top 3 PCs (move with the park) ===")
print(ride_r2.sort_values(ascending=False).head(10).round(3).to_string())
print("\n=== Worst explained (ride-specific dynamics) ===")
print(ride_r2.sort_values().head(10).round(3).to_string())
# Visualize one well-explained ride vs the reconstruction.
RIDE = ride_r2.sort_values(ascending=False).index[0]
fig, ax = plt.subplots(figsize=(14, 4))
df_matrix[RIDE].resample("D").mean().plot(ax=ax, label="actual", alpha=0.6)
reconstructed[RIDE].resample("D").mean().plot(ax=ax, label="top-3 PC reconstruction", alpha=0.9)
ax.set_title(f"{RIDE}: actual vs 3-PC reconstruction (daily-mean queue, min)")
ax.set_ylabel("Queue (min)")
ax.legend()
In [ ]:
# Blog-ready scatter: rides positioned in PC2 x PC3 space, colored by category.
# PC2 = Kiddieland corner vs main midway (spatial; see GPS-cluster validation)
# PC3 = evening crowd (−) vs daytime / school-group daypart (+); NOT throughput.
fig, ax = plt.subplots(figsize=(10, 8))
category_colors = {"Family": "tab:blue", "Kiddie": "tab:orange", "Thrill": "tab:red"}
for category, group in annotated.dropna(subset=["category"]).groupby("category"):
ax.scatter(
group["PC2"], group["PC3"],
label=category,
color=category_colors.get(category, "grey"),
s=80, alpha=0.7, edgecolor="white",
)
# Label only the rides farthest from the origin so the plot stays readable.
distance = (annotated["PC2"] ** 2 + annotated["PC3"] ** 2) ** 0.5
to_label = distance.sort_values(ascending=False).head(15).index
for ride in to_label:
row = annotated.loc[ride]
ax.annotate(ride, (row["PC2"], row["PC3"]), fontsize=8, alpha=0.85,
xytext=(4, 2), textcoords="offset points")
ax.axhline(0, color="grey", linewidth=0.5)
ax.axvline(0, color="grey", linewidth=0.5)
ax.set_xlabel("PC2 Kiddieland corner ← → main midway")
ax.set_ylabel("PC3 evening crowd ← → daytime / school-group daypart")
ax.set_title("Knoebels rides in PC2 × PC3 loading space")
ax.legend(title="Category")
plt.tight_layout()
In [ ]:
# Derive park areas from GPS instead of a missing `parent` column.
# Each POI has a (lat, lon) point; kiddie rides are physically clustered in
# Kiddieland. KMeans on coords gives us an area label derived from geometry
# alone, which we can then cross-check against the PCA loadings.
#
# k=5 carves out a tight 19-ride Kiddieland. k=3 was too coarse — it swept
# main-midway rides near the park entrance (giant_wheel, grand_carousel,
# bumper_cars) into the kiddie blob. IMPORTANT: lat and lon must be standardized
# before KMeans — on raw degrees the clustering is seed-unstable (the contrast r
# wobbles 0.60-0.71 between random_states); standardizing makes it identical
# across seeds at r=0.78. The kiddie cluster is identified by lowest mean PC2
# loading, so this stays self-contained (no hand-typed ride list).
from sklearn.cluster import KMeans
poi = pd.read_sql(
text(
"""
SELECT r.ride AS ride_name,
(p.location)[0]::float AS lat,
(p.location)[1]::float AS lon
FROM knoebels.ride_master r
JOIN knoebels."LZ_attractions_io_poi" p
ON p._id = r.attractions_dot_io_id
WHERE p.location IS NOT NULL
"""
),
engine,
).set_index("ride_name")
located = annotated.join(poi, how="inner")
# KMeans on the geometry only — no PCA information leaks in. Standardize the
# coordinates first so the clustering is reproducible (see note above).
N_AREAS = 5
coords = StandardScaler().fit_transform(located[["lat", "lon"]])
kmeans = KMeans(n_clusters=N_AREAS, n_init=10, random_state=0)
located["area_cluster"] = kmeans.fit_predict(coords)
# Identify the Kiddieland cluster as the one with the lowest mean PC2 loading,
# then export ride-name lists for the validation cell to build a contrast from
# geometry-derived groups (rather than the misfiling `category` tag).
kiddie_cluster = located.groupby("area_cluster")["PC2"].mean().idxmin()
kiddie_cluster_rides = located.index[located["area_cluster"] == kiddie_cluster]
main_cluster_rides = located.index[located["area_cluster"] != kiddie_cluster]
fig, axes = plt.subplots(1, 2, figsize=(15, 7))
# Left: physical map colored by PC2 loading. If the kiddie hypothesis is right,
# one corner of the park should be uniformly blue.
sc = axes[0].scatter(
located["lon"], located["lat"],
c=located["PC2"], cmap="RdBu_r", s=140,
edgecolor="black", linewidth=0.4,
)
axes[0].set_title("Rides colored by PC2 loading\n(blue = kiddie side, red = main park)")
axes[0].set_xlabel("Longitude")
axes[0].set_ylabel("Latitude")
plt.colorbar(sc, ax=axes[0], label="PC2 loading")
for ride, row in located.iterrows():
axes[0].annotate(ride, (row["lon"], row["lat"]),
fontsize=6, alpha=0.65,
xytext=(3, 3), textcoords="offset points")
# Right: same map, colored by KMeans cluster on coords alone.
for cluster, group in located.groupby("area_cluster"):
label = f"Area {cluster}" + (" (Kiddieland)" if cluster == kiddie_cluster else "")
axes[1].scatter(group["lon"], group["lat"], label=label, s=140, alpha=0.85)
axes[1].set_title(f"Spatial clusters (KMeans k={N_AREAS} on standardized lat/lon)")
axes[1].set_xlabel("Longitude")
axes[1].set_ylabel("Latitude")
axes[1].legend()
plt.tight_layout()
print(f"Kiddieland = cluster {kiddie_cluster} ({len(kiddie_cluster_rides)} rides):")
print(" " + ", ".join(sorted(kiddie_cluster_rides)))
print("\n=== Mean PC loading per spatial cluster ===")
print(located.groupby("area_cluster")[["PC1", "PC2", "PC3"]].mean().round(3))
print("\n=== Ride counts per cluster ===")
print(located.groupby("area_cluster").size())
In [ ]:
# Diagnostics: refine the PC labels, and check whether ffill+0 imputation is
# corrupting any rides into a "broken-ride" axis.
for pc in ["PC1", "PC2", "PC3"]:
print(f"\n=== {pc} — full sorted loadings (most negative → most positive) ===")
print(loadings[pc].sort_values().to_string())
silent_rate = (df_matrix == 0).mean().sort_values(ascending=False).rename("silent_rate")
print("\n=== Per-ride silent rate (fraction of rows = 0; high = closed often) ===")
print(silent_rate.head(20).round(3).to_string())
In [ ]:
# Validate PC3 as a *time-of-day* (daypart) axis.
#
# The original "queue turnover / throughput" hypothesis is untestable here:
# ride_master.capacity and .duration are 100% NULL (0 of 60 rides), so there's
# no throughput to correlate against. And splitting rides by PC3 sign shows no
# queue-magnitude difference (neg-PC3 mean queue 9.8 min / busy-rate 0.37 vs
# pos-PC3 8.8 / 0.35) — PC3 is NOT "builds a line vs walk-on."
#
# What PC3 actually tracks is WHEN a ride is busy. The PC3 score runs strongly
# positive midday and negative after dark. Mechanism (park ops): discounted
# school/daycare group rates run from open to ~mid-afternoon and flood the
# daytime family flat-rides (tea_cups, tilt-a-whirl, roto_jets, italian_trapeze);
# the evening crowd (teens / date-night) shifts to coasters, bumper cars, motor
# boats and the Ferris wheel after dark (impulse, bumper_cars, giant_wheel).
afternoon_rides = loadings.index[loadings["PC3"] > 0.1]
evening_rides = loadings.index[loadings["PC3"] < -0.1]
print(f"daytime rides (PC3 > +0.1): {len(afternoon_rides)} -> {sorted(afternoon_rides)}")
print(f"evening rides (PC3 < -0.1): {len(evening_rides)} -> {sorted(evening_rides)}")
hour = df_matrix.index.hour
# De-trend by the park-wide hourly mean so we see *shape* (when each group is
# busy relative to the whole park), not PC1's overall nightly ramp.
park_by_hour = df_matrix.mean(axis=1).groupby(hour).mean()
fig, axes = plt.subplots(1, 2, figsize=(15, 5))
# Left: the PC3 score's own daypart profile — this IS the axis.
scores.groupby(scores.index.hour)["PC3"].mean().plot(ax=axes[0], marker="o")
axes[0].axhline(0, color="grey", linewidth=0.5)
axes[0].set_title("PC3 score by hour of day\n(+ = daytime pattern, − = evening pattern)")
axes[0].set_xlabel("Hour (local)")
axes[0].set_ylabel("Mean PC3 score")
# Right: each loading group's queue as a share of the park-wide level, by hour.
for name, grp, color in [("daytime rides (PC3 +)", afternoon_rides, "tab:orange"),
("evening rides (PC3 −)", evening_rides, "tab:blue")]:
share = df_matrix[grp].mean(axis=1).groupby(hour).mean() / park_by_hour
share.plot(ax=axes[1], marker="o", label=name, color=color)
axes[1].axhline(1.0, color="grey", linewidth=0.5)
axes[1].set_title("Group queue ÷ park-wide queue, by hour\n(de-trended — isolates daypart shape from PC1)")
axes[1].set_xlabel("Hour (local)")
axes[1].set_ylabel("Share of park-wide mean queue")
axes[1].legend()
plt.tight_layout()
In [ ]:
# Validate the PC labels against the most direct possible ground truth:
# the average queue time of the rides each PC is supposed to represent.
# PC1 should track park-wide mean queue.
# PC2 should track the Kiddieland-vs-main-midway contrast (NOT raw kiddie mean —
# on busy days both go up together, so raw kiddie mean is muddied by PC1).
#
# The kiddie/main split here comes from the GPS-cluster cell above
# (`kiddie_cluster_rides` / `main_cluster_rides`), NOT the `category` tag. The
# taxonomy misfiles ~3 rides across the divide (the big giant_wheel and
# grand_carousel are tagged Family but sit on the main midway; the small
# s&g_carousel sits in Kiddieland), which capped the category-based contrast at
# r=0.73. The geometry-derived grouping (k=5) lifts it to ~0.78 — and uses
# spatial information the PCA never saw, so it's an independent validation.
#
# Run the GPS-cluster cell first; this depends on its exported ride lists.
kiddie_cols = df_matrix.columns.intersection(kiddie_cluster_rides)
main_cols = df_matrix.columns.intersection(main_cluster_rides)
print(f"Kiddieland rides in df_matrix: {len(kiddie_cols)} | main-park rides: {len(main_cols)}")
park_mean = df_matrix.mean(axis=1)
kiddie_mean = df_matrix[kiddie_cols].mean(axis=1)
main_mean = df_matrix[main_cols].mean(axis=1)
contrast = main_mean - kiddie_mean
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
axes[0].scatter(park_mean, scores["PC1"], s=4, alpha=0.25)
axes[0].set_xlabel("Park-wide mean queue (min)")
axes[0].set_ylabel("PC1 score")
axes[0].set_title(f"PC1 vs park-wide mean\nPearson r = {park_mean.corr(scores['PC1']):.3f}")
axes[1].scatter(kiddie_mean, scores["PC2"], s=4, alpha=0.25)
axes[1].set_xlabel("Kiddieland mean queue (min)")
axes[1].set_ylabel("PC2 score")
axes[1].set_title(f"PC2 vs raw Kiddieland mean\nPearson r = {kiddie_mean.corr(scores['PC2']):.3f}")
axes[2].scatter(contrast, scores["PC2"], s=4, alpha=0.25)
axes[2].set_xlabel("Main-midway mean − Kiddieland mean (min)")
axes[2].set_ylabel("PC2 score")
axes[2].set_title(f"PC2 vs main-vs-Kiddieland contrast\nPearson r = {contrast.corr(scores['PC2']):.3f}")
plt.tight_layout()