Methods & code
How a water year becomes a piece of music, and how the synced visual is built
The idea in one paragraph
Sonification means turning data into sound. Instead of mapping one number straight to pitch, which sounds like an alarm, these pieces separate timescales and signal according to music theory. Slow signals (the season, lake level, how far the year sits from its 15-year normal) set the harmony. Medium signals (streamflow, plant productivity) set density and presence. Fast or eventful signals (a rainy day, ice forming, fish passing a dam) trigger short motifs. A repeating piano figure holds it all together: a minimalist form, where one pattern stays constant and everything around it changes. This means the sound isn’t exactly the data, or even exploratory data analysis, but is accompaniment to data visualization designed for deeper appreciation and exploration.
Design principles
- Separate the timescales, as above.
- Minimalism as the form. A piano ostinato (a short repeated figure) is the constant. The data changes its density, register, harmony and surroundings.
- Consonance is enforced. Every pitched note is drawn from the current chord, and chords pass two guards against clashing intervals (see Consonance guards below).
- Every layer must be audible when its data is active, and nothing may drown it out. This is checked by measurement (see Quality checks).
- Provenance is visible. Every layer is tagged R retrieved measurement, D derived, or S simulated timing.
Data
| Record | Source | Used for | |
|---|---|---|---|
| Granby River at Grand Forks, daily discharge + data flags | Water Survey of Canada 08NN002 | flow, anomaly, hydrologic phase, ice (the “Ice Conditions” flag) | R |
| Burrell Creek above Gloucester Creek | WSC 08NN023 | hand-drum drive | R |
| Billings climate station | ECCC daily climate, station 1100 | temperature, rain, snowfall, snow on ground | R |
| Okanagan Lake at Kelowna, level | WSC 08NM083 | lake drone, anomaly, lake phase | R |
| Okanagan River at Penticton (dam outflow) | WSC 08NM050 | piano density | R |
| Okanagan River near Oliver | WSC 08NM085 | soft percussion | R |
| Summerland climate station | ECCC daily climate, station 979 | temperature, precipitation, humidity | R |
| Gross primary productivity (photosynthesis) | NASA MODIS MOD17A2HGF, 8-day, mean of a 5×5 block of 500 m pixels (about 2.3 km on a side), via ORNL DAAC | string pad | R |
| Adult sockeye count and water temperature | Wells Dam, Columbia Basin Research DART | marimba runs, thermal stress | R |
| Day length | solar geometry at each site’s latitude | mode brightness | D |
Hydrometric and climate data come from ECCC’s GeoMet OGC API (daily means, 2010–2024, so every day of 2024 can be compared with its 15-year normal). Wells Dam sits on the Columbia just below the Okanagan confluence and above the Wenatchee, so its sockeye are overwhelmingly Okanagan-bound, but it is not an Okanagan gauge.
Data dictionaries
Every table the pieces read, and every table of what they decided, is described in a data dictionary in the data-dict format: what each value is, its units and range, how it was measured and its quirks, how the tables relate, and a glossary. The site’s build checks the data against the dictionaries on every push. There are three layers:
- The records. The water-year records live as tidy tables in
data/parquet/, beside the original API responses (browse). The Climate Pair’s monthly indices, storms, volcanoes and Kettle River peaks are indata/climate/parquet/(browse). The values are copied from what the agencies published; the pieces transform them in code. Quirks the dictionaries record include Summerland’s separate rain and snow values, missing on more than half its days (so the Okanagan piece splits total precipitation by temperature), lake levels measured above a local gauge datum rather than sea level, the AMO record ending in January 2023, and two definitions of storm landfall. - The mappings. Each track has a mapping dictionary in
tracks/mapping/<track>/: what the piece read for each day (or year), and what it decided for each bar (mode, chord, phase, density). Every transform between the two is named and bounded in the column that uses it, as a step in a modelling recipe would be, and simple rules are written as assertions the build checks: a storm’s category follows from its peak wind, a density stays between 0 and 3. The tables come from the score of the locked render, and only if its fingerprint matches the render lock, so a mapping always describes the audio of the version it names. Browse them for Granby, Okanagan, the Pacific, the Atlantic, Meridian Chorus (light, air, sea and ice; daylight only) and Superior Ice Year (main; shorter). - The words. One shared glossary holds every term the dictionaries define, grouped by field, and a test checks that each dictionary uses its definitions word for word.
Some records are played but not republished, and their values appear in no table here: the Wells Dam sockeye counts and water temperatures (Columbia River DART data; the listening page shows them as “heard only”), SILSO sunspot numbers (CC BY-NC) and IAP ocean heat content (terms unconfirmed). Their dictionaries stay with the local copies, because a dictionary’s ranges are values too; the fetch scripts download the records themselves.
Versions, variants and remixes
Every track has a version. A version changes when the music changes. Each track’s card and changelog are on the Tracks page. renders.lock.json records, for each track: its version, the script and arguments that render it, the random seed, content hashes of the engine code, of every data file and instrument sample the render read, and of every file it wrote, plus the package and ffmpeg versions. python scripts/renders.py verify re-renders every track into a scratch folder and compares; verify granby-wy2024 checks one. It says whether the audio still matches and, if not, which code or inputs changed since the lock. The lock is written only when a scratch render matches the published files, so it always describes the audio that exists. Code is pinned by content hash rather than by commit, so the check holds outside the repository too.
Variants are named changes to one track, kept for listening comparisons and version history: python pilot_granby_wy2024.py --variant frozen renders Granby with its mood held in Dorian all year. They live in variants/<piece>.yaml, and each piece lists the parameters it accepts (DEFAULTS at the top of its script: for Granby the year, tempo, key, a frozen mode, the mode range, muted layers, layer gains and the seed). A name the piece doesn’t declare is refused, so a typo can’t silently render the default.
Remixes are new tracks built from a published one. A short YAML file names the parent track and version, what it keeps, what it changes and who made it; the remix renders under its own name, and the lock records its lineage. See remixes/README.md and CONTRIBUTING.md.
The listening companion
The Listen page uses no framework, only features built into every browser:
| Job | How |
|---|---|
| Playback and timing | an HTML <audio> element; its current playback time is the single clock everything reads from |
| Animation | requestAnimationFrame, redrawing while the music plays |
| The rings | Canvas 2D. The full year is drawn once dim and once lit, and each frame shows the lit copy only up to the playhead |
| Linking sound to picture | how loud each instrument group is, measured 8 times a second from its own audio track ahead of time, then looked up by playback time. That drives the glows, the meters and the drifting particles |
| Data | one JSON file per piece: daily records, per-bar harmony, events and loudness |
The loudness is measured ahead of time rather than live because a browser can only analyse the finished mix, not each instrument separately.
Reproduce it
Everything runs in Python 3.13 with the pinned packages in requirements.txt (numpy, scipy, soundfile, matplotlib), plus ffmpeg. Run the scripts from scripts/.
To hear the engine without downloading the sample library, python demo.py renders a 45-second synthesizer sketch of the Granby year (no samples needed; about 10 seconds to render).
pip install -r ../requirements.txt
python fetch_granby.py # Granby, Burrell Creek and Billings pulls (cached)
python fetch_okanagan.py # Okanagan, MODIS and DART pulls (cached)
python fetch_samples.py # VSCO subset + bird recordings + licence manifest
python fetch_records.py wsc 08NN002 granby --years 2010-2024 # any WSC station (or: eccc STN_ID name)
python pilot_granby_wy2024.py # renders about 7 minutes of audio in about 30 s
python pilot_granby_wy2024.py wy2019 # the same piece, another year
python pilot_granby_wy2024.py --variant frozen
python renders.py verify # re-render every locked track and compare (needs the samples)
python check_sample_pitch.py # every sample's pitch against its label
python pilot_okanagan_wy2024.py
python pilot_check_plot.py granby_wy2024
python export_companion.py # data + audio for the listening page
python build_parquet.py # tidy tables from the cached pulls (data/parquet, data/climate/parquet)
python export_mapping.py # each locked track's mapping tables (tracks/mapping)By default the sample library goes to samples/ and the rendered audio to rendered/, both inside the repository and both git-ignored. Set SONIFICATION_SAMPLES and SONIFICATION_HEAVY to put them elsewhere. Renders are seeded: with the pinned packages, the same data gives the same audio, byte for byte. The raw Wells Dam file is not in the repository; fetch_okanagan.py downloads it from DART.
Limitations
- Simulated timing. Fish-spawning windows are modelled from air temperature and season, not observed. Bird windows are typical regional timing, not 2023–24 sightings.
- Air temperature stands in for water temperature on the Granby, where no stream-temperature record was available. The Wells temperature is Columbia River water, measured in the fishway’s attraction flow, which is pumped from the dam tailrace.
- The satellite productivity values average 25 pixels, a block about 2.3 km on a side, so they mix land covers. The Summerland point includes irrigated orchard.
Source code
Note
timegrid.py: the year, its framing, bars and seconds
"""
Time for the pieces: which days a track covers (Year), how days become bars,
and how bars become seconds (Grid). Split out of seasonal.py on 2026-09-27.
A Year is a run of consecutive days with a start set by its FRAMING:
water Oct 1 by default (northern temperate hydrology); named by the
calendar year it ends in, so WY2024 = Oct 1 2023 - Sep 30 2024.
Pass start=(month, day) where the water-supply regime differs.
calendar Jan 1 - Dec 31.
solstice one solstice to the next (December by default, or June).
lunar 12 lunations from the first new moon on or after an anchor
date (Jan 1 by default): 354 or 355 days.
equatorial a start the track must supply, set by its water-supply regime
(e.g. the onset of the rains); named like a water year.
Length follows the framing, so leap years, 365-day years and 354-day lunar
years all come out right. One BAR = 3 DAYS; a year whose length isn't a
multiple of 3 ends on a short bar (per_bar averages what that bar has).
Solstice and new-moon times use Meeus, Astronomical Algorithms (2nd ed.),
ch. 27 (mean solstice) and ch. 49 (new moon, main periodic terms); dates are
UTC and good to well under a day (tests/test_timegrid.py checks them against
published times).
"""
import datetime as dt
import math
import numpy as np
from dsp_core import SR, day_length_hours
DAYS_PER_BAR = 3
# ------------------------------------------------------------- astronomy ---
_J2000_ORDINAL_OFFSET = 1721424.5 # JD of 0001-01-01 00:00 is ordinal 1 + this
def _jd_to_datetime(jd):
"""Julian day -> UTC datetime (ignores the ~70 s TT-UT difference)."""
days = jd - _J2000_ORDINAL_OFFSET
ordinal = int(math.floor(days))
return dt.datetime.fromordinal(ordinal) + dt.timedelta(days=days - ordinal)
def solstice(year, which="december"):
"""UTC datetime of the June or December solstice in `year` (Meeus table 27.B,
valid 2000-3000; the periodic terms are left out, which costs minutes)."""
Y = (year - 2000) / 1000
if which == "december":
jde = 2451900.05952 + 365242.74049 * Y - 0.06223 * Y**2 - 0.00823 * Y**3 + 0.00032 * Y**4
elif which == "june":
jde = 2451716.56767 + 365241.62603 * Y + 0.00325 * Y**2 + 0.00888 * Y**3 - 0.00030 * Y**4
else:
raise ValueError(f"solstice: which must be 'december' or 'june', not {which!r}")
return _jd_to_datetime(jde)
def _new_moon_jde(k):
"""Meeus ch. 49: JDE of new moon number k (k = 0 near 2000-01-06)."""
T = k / 1236.85
jde = (2451550.09766 + 29.530588861 * k + 0.00015437 * T**2
- 0.000000150 * T**3 + 0.00000000073 * T**4)
E = 1 - 0.002516 * T - 0.0000074 * T**2
r = math.radians
M = r(2.5534 + 29.10535670 * k - 0.0000014 * T**2 - 0.00000011 * T**3)
Mp = r(201.5643 + 385.81693528 * k + 0.0107582 * T**2 + 0.00001238 * T**3 - 0.000000058 * T**4)
F = r(160.7108 + 390.67050284 * k - 0.0016118 * T**2 - 0.00000227 * T**3 + 0.000000011 * T**4)
Om = r(124.7746 - 1.56375588 * k + 0.0020672 * T**2 + 0.00000215 * T**3)
s = math.sin
jde += (-0.40720 * s(Mp) + 0.17241 * E * s(M) + 0.01608 * s(2 * Mp) + 0.01039 * s(2 * F)
+ 0.00739 * E * s(Mp - M) - 0.00514 * E * s(Mp + M) + 0.00208 * E**2 * s(2 * M)
- 0.00111 * s(Mp - 2 * F) - 0.00057 * s(Mp + 2 * F) + 0.00056 * E * s(2 * Mp + M)
- 0.00042 * s(3 * Mp) + 0.00042 * E * s(M + 2 * F) + 0.00038 * E * s(M - 2 * F)
- 0.00024 * E * s(2 * Mp - M) - 0.00017 * s(Om))
return jde
def new_moons_from(anchor, n):
"""UTC datetimes of the first `n` new moons whose UTC date is on or after `anchor`."""
k = math.floor((anchor.year + (anchor.timetuple().tm_yday - 1) / 365.25 - 2000) * 12.3685) - 1
out = []
while len(out) < n:
t = _jd_to_datetime(_new_moon_jde(k))
if t.date() >= anchor:
out.append(t)
k += 1
return out
# ------------------------------------------------------------------ year ---
FRAMINGS = ("water", "calendar", "solstice", "lunar", "equatorial")
class Year:
"""The days one track covers: `start` (inclusive) to `end` (exclusive)."""
def __init__(self, start, end, framing, label):
if framing not in FRAMINGS:
raise ValueError(f"unknown framing {framing!r}; one of {FRAMINGS}")
self.start, self.end, self.framing, self.label = start, end, framing, label
self.n_days = (end - start).days
self.days = [start + dt.timedelta(i) for i in range(self.n_days)]
self.n_bars = -(-self.n_days // DAYS_PER_BAR)
def __repr__(self):
return f"Year({self.label}: {self.start} to {self.days[-1]}, {self.n_days} days, {self.n_bars} bars)"
# --- framings -----------------------------------------------------------
@classmethod
def water(cls, year, start=(10, 1)):
"""Water year `year`, ending in calendar year `year` (WY2024 = Oct 2023 - Sep 2024)."""
s = dt.date(year if start == (1, 1) else year - 1, *start)
tag = "" if start == (10, 1) else f"@{start[0]:02d}-{start[1]:02d}"
return cls(s, s.replace(year=s.year + 1), "water", f"WY{year}{tag}")
@classmethod
def calendar(cls, year):
return cls(dt.date(year, 1, 1), dt.date(year + 1, 1, 1), "calendar", f"CY{year}")
@classmethod
def solstice(cls, year, which="december"):
"""From the `which` solstice of `year` to the same solstice a year later."""
s, e = solstice(year, which).date(), solstice(year + 1, which).date()
return cls(s, e, "solstice", f"{'DS' if which == 'december' else 'JS'}{year}")
@classmethod
def lunar(cls, year, anchor=None):
"""12 lunations from the first new moon on or after `anchor` (default Jan 1 of `year`)."""
moons = new_moons_from(anchor or dt.date(year, 1, 1), 13)
return cls(moons[0].date(), moons[12].date(), "lunar", f"LY{moons[0].date():%Y%m%d}")
@classmethod
def equatorial(cls, year, start):
"""A year set by the local water-supply regime; `start` = (month, day) is required."""
s = dt.date(year if start == (1, 1) else year - 1, *start)
return cls(s, s.replace(year=s.year + 1), "equatorial", f"EY{year}@{start[0]:02d}-{start[1]:02d}")
@classmethod
def parse(cls, spec):
"""'wy2024', 'wy2024@07-01', 'cy2024', 'ds2023', 'js2024', 'ly2024',
'ly2024@2024-02-10', 'ey2024@11-01' -> Year. Used by the pieces' --year flag."""
s = spec.strip().lower()
head, _, at = s.partition("@")
kind, num = head[:2], int(head[2:])
md = tuple(int(p) for p in at.split("-")) if at else None
if kind == "wy":
return cls.water(num, md or (10, 1))
if kind == "cy":
return cls.calendar(num)
if kind in ("ds", "js"):
return cls.solstice(num, "december" if kind == "ds" else "june")
if kind == "ly":
return cls.lunar(num, dt.date(*md) if md else None)
if kind == "ey":
if not md:
raise ValueError("an equatorial year needs its start: ey2024@MM-DD")
return cls.equatorial(num, md)
raise ValueError(f"can't read year spec {spec!r}")
# --- day helpers --------------------------------------------------------
def doy_index(self, doy):
"""Index into `days` of calendar day-of-year `doy` (1-366); the first
occurrence if the year spans two of them. ValueError if it isn't in the year."""
for i, d in enumerate(self.days):
if d.timetuple().tm_yday == doy:
return i
raise ValueError(f"day-of-year {doy} isn't in {self.label}")
def daylength(self, lat):
"""Astronomical day length (hours) at latitude `lat` for each day."""
return np.array([day_length_hours(d.timetuple().tm_yday, lat) for d in self.days])
WY2024 = Year.water(2024) # the default: every piece before 2026-09-27 was WY2024
# ------------------------------------------------------------------ grid ---
class Grid:
def __init__(self, bpm, n_bars=WY2024.n_bars, beats_per_bar=4, tail_s=10.0):
self.bpm = bpm
self.beat = 60.0 / bpm
self.bar = self.beat * beats_per_bar
self.n_bars = n_bars
self.dur = n_bars * self.bar + tail_s
self.n = int(self.dur * SR)
def t(self, bar, beat=0.0):
return (bar + beat / 4.0) * self.bar
def s(self, bar, beat=0.0, jitter_ms=0.0, rng=None):
j = (rng.normal(0, jitter_ms) / 1000.0) if (rng is not None and jitter_ms) else 0.0
return max(0, int((self.t(bar, beat) + j) * SR))
def day_to_time(self, day_float):
return day_float / DAYS_PER_BAR * self.bar
def per_bar(x, how="mean"):
"""Daily array -> one value per bar (mean, sum, max...). A short last bar
(years that aren't a multiple of 3 days) reduces over the days it has."""
x = np.asarray(x)
n_bars = -(-len(x) // DAYS_PER_BAR)
pad = n_bars * DAYS_PER_BAR - len(x)
if pad == 0:
return getattr(x.reshape(n_bars, DAYS_PER_BAR), how)(axis=1)
x = np.concatenate([x.astype(float), np.full(pad, np.nan)]).reshape(n_bars, DAYS_PER_BAR)
return getattr(np, "nan" + how)(x, axis=1)
# --------------------------------------------------------- daily helpers ---
def norm01(x, lo_pct=2, hi_pct=98):
lo, hi = np.percentile(x, lo_pct), np.percentile(x, hi_pct)
return np.clip((x - lo) / (hi - lo + 1e-12), 0, 1)
def smooth(x, win):
k = np.ones(win) / win
xp = np.concatenate([np.full(win, x[0]), x, np.full(win, x[-1])])
return np.convolve(xp, k, mode="same")[win:-win]
def events_from_flag(flag, min_gap=3):
"""Onset/offset day indices of runs in a boolean daily flag, merging gaps
shorter than `min_gap` days (so one mild afternoon doesn't 'break up' the
ice and re-freeze it the next day)."""
f = flag.astype(bool).copy()
idx = np.flatnonzero(f)
for a, b in zip(idx[:-1], idx[1:]):
if 1 < b - a <= min_gap:
f[a:b] = True
on = np.flatnonzero(f & ~np.r_[False, f[:-1]])
off = np.flatnonzero(f & ~np.r_[f[1:], False])
return list(zip(on, off)), f
def fmt_time(sec):
return f"{int(sec // 60)}:{int(sec % 60):02d}"
Note
harmony.py: modes, chords, voice leading, progressions
"""
The harmonic grammar shared by the pieces: pure functions, no I/O. Split out
of seasonal.py on 2026-09-27.
Music-theory model (written up for GW in the listening guide):
MODES ordered dark -> bright. Neighbours differ by exactly ONE note, so
sliding one step along this ladder is the smallest possible change of mood.
The season sets a base brightness (day length); the anomaly vs. that day's
15-year normal nudges it one step darker (drought/low) or brighter (high).
"""
import numpy as np
from timegrid import DAYS_PER_BAR
MODES = [ # dark -> bright; each differs from its neighbour by one note
("Phrygian", [0, 1, 3, 5, 7, 8, 10]),
("Aeolian", [0, 2, 3, 5, 7, 8, 10]),
("Dorian", [0, 2, 3, 5, 7, 9, 10]),
("Mixolydian", [0, 2, 4, 5, 7, 9, 10]),
("Ionian", [0, 2, 4, 5, 7, 9, 11]),
("Lydian", [0, 2, 4, 6, 7, 9, 11]),
]
def degree_to_midi(deg, mode_idx, tonic):
steps = MODES[mode_idx][1]
o, d = divmod(deg, 7)
return tonic + 12 * o + steps[d]
def chord(deg, mode_idx, tonic, extensions=(0, 2, 4, 8)):
"""Stack thirds on scale degree `deg`: 0,2,4 = triad; 6 = 7th; 8 = 9th.
Default is an add9 (no 7th) - the open, 'cinematic' colour.
Two consonance guards (the chime clashes in track 04 were this class of
problem):
- a DIMINISHED triad (root-to-5th = 6 semitones) is swapped for the
chord a third below, which contains the same upper notes plus a
stable root (Bdim -> G, the textbook substitution);
- a 9th that lands a half step above the root (b9) is swapped for the
7th, which is colour rather than clash."""
if (degree_to_midi(deg + 4, mode_idx, tonic) - degree_to_midi(deg, mode_idx, tonic)) % 12 == 6:
deg -= 2
root = degree_to_midi(deg, mode_idx, tonic)
out = []
for e in extensions:
m = degree_to_midi(deg + e, mode_idx, tonic)
if e == 8 and (m - root) % 12 == 1:
m = degree_to_midi(deg + 6, mode_idx, tonic)
out.append(m)
return out
def voice_lead(prev, target_pcs, lo, hi):
"""Place each target pitch-class in [lo,hi] near the previous voicing.
NOTE (2026-09-26 tests): tone i is placed near previous voice i, but the
previous voicing is sorted and the new chord isn't, so this is not the
minimal motion its name promises (tests/test_harmony.py, strict xfail).
The fix changes the music, so it waits for GW's A/B (sprint plan §3)."""
out = []
for i, pc in enumerate(target_pcs):
ref = prev[i] if prev and i < len(prev) else (lo + hi) // 2
cands = [m for m in range(lo, hi + 1) if m % 12 == pc % 12]
out.append(min(cands, key=lambda m: abs(m - ref)))
return sorted(out)
def mode_track(season_bright, anomaly, lo=1, hi=5, phrase=8, max_step=1, floor=1, ceil=5):
"""Per-bar mode index. season_bright/anomaly are per-bar 0..1 / -0.5..0.5.
Changes only at phrase boundaries (every `phrase` bars) and by at most one
step, so the listener hears a slow drift, never a lurch. `floor`=1 keeps
Phrygian (whose b2 sits a half step above the tonic) out by default."""
target = lo + season_bright * (hi - lo) + np.clip(anomaly * 2.4, -1.2, 1.2)
out = np.zeros(len(target), dtype=int)
cur = int(np.clip(round(target[0]), floor, ceil))
for b in range(len(target)):
if b % phrase == 0:
t = int(round(np.mean(target[b:b + phrase])))
cur += int(np.clip(t - cur, -max_step, max_step))
cur = int(np.clip(cur, floor, ceil))
out[b] = cur
return out
def phase_bars(phase_day, min_run=6):
"""Daily phase labels -> one label per bar (majority vote), then any run
shorter than `min_run` bars is absorbed into the run before it, so the
progression doesn't flip back and forth on a two-day wobble."""
n_bars = -(-len(phase_day) // DAYS_PER_BAR)
raw = []
for b in range(n_bars):
v, c = np.unique(phase_day[b * DAYS_PER_BAR:(b + 1) * DAYS_PER_BAR], return_counts=True)
raw.append(v[np.argmax(c)])
runs, i = [], 0
while i < n_bars:
j = i
while j < n_bars and raw[j] == raw[i]:
j += 1
runs.append([raw[i], i, j])
i = j
merged = []
for r in runs:
if merged and (r[2] - r[1] < min_run or merged[-1][0] == r[0]):
merged[-1][2] = r[2]
else:
merged.append(r)
return [ph for ph, a, z in merged for _ in range(a, z)]
def chord_track(bar_phase, prog, bars_per_chord):
"""Walk each phase's chord loop; a new phase restarts its loop at the top."""
chords, pos, cur = [], 0, None
for b, ph in enumerate(bar_phase):
if b % bars_per_chord == 0:
if ph != cur:
cur, pos = ph, 0
else:
pos = (pos + 1) % len(prog[cur])
chords.append(prog[cur][pos])
return chords
def fire_bars(intensity_bar, rate, phase=0.5):
"""Integrate-and-fire: accumulate intensity*rate per bar and fire a motif
each time the total crosses 1. Deterministic, so a sparse layer always
speaks in proportion to its data (random per-bar gating left the Granby
whitefish layer completely silent on one 2026-09-24 render: ~5% odds
over its short window)."""
acc, out = phase, []
for b, v in enumerate(intensity_bar):
acc += v * rate
if acc >= 1:
out.append(b)
acc -= 1
return set(out)
Note
records.py: loading the data records for a year
"""
Record loaders: cached pulls (data/*.json) -> date->value dicts -> daily
arrays for one Year. Split out of seasonal.py on 2026-09-27.
Every function that returns a daily array takes `year` (a timegrid.Year) and
defaults to WY2024, the year every piece before 2026-09-27 was built on.
Since 2026-09-27 the loaders read the tidy parquet tables in data/parquet/
(built by build_parquet.py and described by data/parquet/data-dict.yaml).
A station that isn't in the tables yet (fetched since the last build) is
read from its raw JSON pull instead, so fetch_records.py keeps working.
"""
import csv, io, json, re, datetime as dt
from functools import lru_cache
from pathlib import Path
import numpy as np
from timegrid import WY2024
HERE = Path(__file__).resolve().parent.parent
DATA = HERE / "data"
PARQUET = DATA / "parquet"
@lru_cache(maxsize=None)
def _rows(table, station):
"""One station's rows from a parquet table, in stored order, as dicts (None if absent)."""
p = PARQUET / f"{table}.parquet"
if not p.exists():
return None
import pyarrow.compute as pc
import pyarrow.parquet as pq
t = pq.read_table(p)
t = t.filter(pc.equal(t["station_id"], station))
return t.to_pylist() if t.num_rows else None
@lru_cache(maxsize=None)
def _load(fname):
"""Parse a cached pull once per run (the climate file is read for four variables)."""
with open(DATA / fname) as fh:
return json.load(fh)
def wsc_series(fname, key="DISCHARGE"):
"""dict date->value (and date->symbol) for a WSC station, named by its pull's
file name (e.g. granby_08NN002_2010_2024.json)."""
out, sym = {}, {}
m = re.search(r"_(\d\d[A-Z]{2}\d{3})_", fname)
rows = _rows("hydrometric_daily", m.group(1)) if m else None
if rows is not None:
col = key.lower()
for r in rows:
d = str(r["date"])
if r[col] is not None:
out[d] = float(r[col])
sym[d] = r[col + "_symbol"]
return out, sym
for f in _load(fname):
p = f["properties"]
if p.get(key) is not None:
out[p["DATE"][:10]] = float(p[key])
sym[p["DATE"][:10]] = p.get(key + "_SYMBOL_EN")
return out, sym
def climate_series(fname, key):
"""dict date->value of one daily climate variable (e.g. "MEAN_TEMPERATURE"),
for a station named by its pull's file name (billings_1100_climate_2010_2024.json)."""
out = {}
m = re.search(r"_(\d+)_climate_", fname)
rows = _rows("climate_daily", m.group(1)) if m else None
if rows is not None and key.lower() in rows[0]:
for r in rows:
if r[key.lower()] is not None:
out[str(r["date"])] = float(r[key.lower()])
return out
for f in _load(fname):
p = f["properties"]
if p.get(key) is not None:
out[p["LOCAL_DATE"][:10]] = float(p[key])
return out
def in_year(series, fill="interp", year=WY2024):
"""Pull `year`'s days out of a date->value dict as an array. Gaps are
interpolated (or zero-filled with fill='zero'; left NaN with fill=None)."""
x = np.array([series.get(str(d), np.nan) for d in year.days], dtype=float)
if fill == "interp" and np.isnan(x).any():
ok = ~np.isnan(x)
if ok.sum() == 0:
return np.zeros(year.n_days)
x = np.interp(np.arange(year.n_days), np.flatnonzero(ok), x[ok])
elif fill == "zero":
x = np.nan_to_num(x)
return x
def coverage(series, year=WY2024):
"""Fraction of `year`'s days with a value: check before rendering a new year."""
return float(np.mean([str(d) in series for d in year.days]))
def doy_percentile(series, halfwin=7, year=WY2024):
"""For each day of `year`: where does its value sit among ALL years'
values within +-halfwin days of the same calendar day? 0..1. This is the
'is it normal for the time of year' signal - the anomaly, not the level."""
by_doy = {}
for k, v in series.items():
d = dt.date.fromisoformat(k)
by_doy.setdefault(d.timetuple().tm_yday, []).append(v)
out = np.full(year.n_days, 0.5)
for i, d in enumerate(year.days):
v = series.get(str(d))
if v is None:
continue
doy = d.timetuple().tm_yday
pool = []
for o in range(-halfwin, halfwin + 1):
pool += by_doy.get((doy - 1 + o) % 366 + 1, [])
pool = np.array(pool)
out[i] = (pool < v).mean() + 0.5 * (pool == v).mean()
return out
def modis_gpp_daily(name, year=WY2024):
"""MODIS MOD17A2HGF 8-day GPP at a point -> daily array for `year`, plus a
0..1 version scaled by the 2015-2024 95th percentile.
Returns (None, None) if the pull isn't cached (caller falls back to a
growing-degree-day model and says so)."""
rows = None
tbl = _rows("modis_gpp", name)
if tbl is not None: # regroup the tidy rows into composites, pixels in stored order
by = {}
for r in tbl:
by.setdefault(str(r["date"]), []).append(r["gpp_raw"])
rows = [{"date": d, "data": v} for d, v in by.items()]
else:
p = DATA / f"modis_gpp_{name}_2015_2024.json"
if not p.exists():
return None, None
with open(p) as fh:
rows = json.load(fh)["rows"]
dates, vals = [], []
for r in rows:
v = [x for x in r["data"] if x < 32761]
if v:
dates.append(dt.date.fromisoformat(r["date"]) + dt.timedelta(4)) # composite centre
vals.append(np.mean(v) * 0.0001 / 8 * 1000) # g C m-2 day-1
if dates[0] > year.start + dt.timedelta(8) or dates[-1] < year.days[-1] - dt.timedelta(8):
return None, None # the cache doesn't cover this year: don't extrapolate a flat season
t = np.array([(d - year.start).days for d in dates])
v = np.array(vals)
daily = np.interp(np.arange(year.n_days), t, v)
return daily, np.clip(daily / np.percentile(v, 95), 0, 1)
def growing_degree_days(tmean, base=5.0):
return np.cumsum(np.maximum(0, tmean - base))
def dart_sockeye(year=WY2024):
"""Wells Dam daily adult sockeye + Wells water temperature (C) for `year`.
Counting season is May-Nov; off-season days = 0 / NaN."""
raw = _load("wells_dam_sockeye_2015_2024.json")
cnt, temp = {}, {}
for y in range(year.start.year, year.days[-1].year + 1):
if str(y) not in raw: # a silent zero would play as "no fish came"
raise ValueError(f"no Wells Dam counts cached for {y} ({year.label})")
for r in csv.DictReader(io.StringIO(raw[str(y)])):
d = r.get("Date") or ""
if not d.startswith(str(y)):
continue
cnt[d] = max(0.0, float(r["Sock"] or 0))
if r["TempC"] and float(r["TempC"]) != 0.0: # 0.0 is a missing code (2015-07-23): no 0 °C Columbia in season
temp[d] = float(r["TempC"])
return in_year(cnt, fill="zero", year=year), np.array([temp.get(str(d), np.nan) for d in year.days])
def dart_run_totals(year=WY2024):
"""Adult sockeye counted at Wells Dam in each calendar year `year` touches."""
raw = _load("wells_dam_sockeye_2015_2024.json")
out = {}
for y in range(year.start.year, year.days[-1].year + 1):
if str(y) not in raw:
raise ValueError(f"no Wells Dam counts cached for {y} ({year.label})")
out[y] = sum(max(0.0, float(r["Sock"] or 0)) for r in csv.DictReader(io.StringIO(raw[str(y)]))
if (r.get("Date") or "").startswith(str(y)))
return out
def check_coverage(named_series, year, need=0.9):
"""Refuse to render a year the records don't cover: {name: series} -> ValueError
listing every series below `need`, so a gap is never heard as a quiet year."""
short = {k: coverage(v, year) for k, v in named_series.items()}
short = {k: c for k, c in short.items() if c < need}
if short:
raise ValueError(f"{year.label}: records too thin to render: "
+ ", ".join(f"{k} {c:.0%}" for k, c in short.items()))
# ---------------------------------------------------------- Climate Pair ---
CLIMATE = DATA / "climate"
def _climate_table(name):
"""A Climate Pair parquet table as {column: list}, or None if it hasn't been built."""
p = CLIMATE / "parquet" / f"{name}.parquet"
if not p.exists():
return None
import pyarrow.parquet as pq
return pq.read_table(p).to_pydict()
def climate_monthly(name="indices_monthly"):
"""{'months': ['1958-01', ...], 'series': {name: [value or None]}} from the parquet table
(or the parsed JSON it was built from). 'fetched_only_monthly' is the local-only file."""
t = _climate_table(name)
if t is None:
return json.loads((CLIMATE / f"{name}.json").read_text())
months = [f"{d:%Y-%m}" for d in t.pop("month")]
return {"months": months, "series": t}
def volcanoes():
t = _climate_table("volcanoes")
if t is None:
return json.loads((CLIMATE / "indices_monthly.json").read_text())["volcanoes"]
return [{"date": str(d), "name": n, "vei": v, "cooling": c, **({"note": x} if x else {})}
for d, n, v, c, x in zip(t["date"], t["name"], t["vei"], t["cooling"], t["note"])]
def storms(basin):
"""One basin's tropical cyclones (atl, nepac, wpac), in stored order."""
t = _climate_table("storms")
if t is None:
return json.loads((CLIMATE / f"storms_{basin}.json").read_text())
keys = [k for k in t if k != "basin"]
return [{k: (t[k][i].isoformat() if k == "peak_time" else t[k][i]) for k in keys}
for i, b in enumerate(t["basin"]) if b == basin]
def kettle_annual_peak():
"""{'1929': {'max_daily_m3s': ..., 'date': 'YYYY-MM-DD'}, ...}"""
t = _climate_table("kettle_annual_peak")
if t is None:
return json.loads((CLIMATE / "kettle_annual_max.json").read_text())
return {str(y): {"max_daily_m3s": q, "date": str(d)} for y, q, d in zip(t["year"], t["max_daily_m3s"], t["date"])}
Note
renders.py and variants.py: the render lock, variants and remixes
"""
Render manifest and verification (2026-09-27): the backbone of track versions.
`renders.lock.json` (repo root) holds one entry per track (and variant): its
version, the script and arguments that render it, the seed, SHA-256 hashes of
the engine code, of every data file and sample the render opened, and of every
file it wrote. A version changes when the music changes; this file proves when
it didn't.
python renders.py lock TRACK VERSION -- SCRIPT [ARGS...] render, check, record
python renders.py verify [TRACK ...] re-render, compare (all if none named)
python renders.py list
Both `lock` and `verify` render into a scratch folder (SONIFICATION_HEAVY,
_LIGHT and _SITE_AUDIO point there), so the published files are never
overwritten. `lock` also requires the scratch render to match the files in
place, so a lock always describes the audio that exists. Code is pinned by
content hash rather than a git commit: the lock is written before the commit,
and the hashes also hold outside the repository.
"""
import argparse
import datetime as dt
import hashlib
import json
import os
import platform
import re
import subprocess
import sys
import tempfile
from pathlib import Path
HERE = Path(__file__).resolve().parent
ROOT = HERE.parent
LOCK = ROOT / "renders.lock.json"
ENGINE = ["dsp_core.py", "timegrid.py", "harmony.py", "records.py", "seasonal.py", "sampler.py", "synth.py", "solar.py", "meteogram_kit.py"]
PY = sys.executable
# Run SCRIPT with an audit hook that logs every file it opens for reading.
_RUNNER = r"""
import sys, os, runpy
log = open(os.environ["RENDERS_OPEN_LOG"], "w")
def hook(event, args):
if event == "open" and isinstance(args[0], (str, bytes, os.PathLike)):
mode = args[1] if len(args) > 1 and isinstance(args[1], str) else "r"
if "r" in mode and "+" not in mode:
log.write(os.fspath(args[0]) if not isinstance(args[0], bytes) else args[0].decode() )
log.write("\n"); log.flush()
sys.addaudithook(hook)
import soundfile # samples open through libsndfile, which the audit hook doesn't see
_init = soundfile.SoundFile.__init__
def _logged(self, file, mode="r", *a, **k):
if isinstance(file, (str, os.PathLike)) and "r" in mode and "+" not in mode:
log.write(os.fspath(file) + "\n"); log.flush()
_init(self, file, mode, *a, **k)
soundfile.SoundFile.__init__ = _logged
try: # parquet tables open through pyarrow's C++ readers, also unseen by the hook
import pyarrow.parquet as _pq
_rt = _pq.read_table
def _read_table(source, *a, **k):
if isinstance(source, (str, os.PathLike)):
log.write(os.fspath(source) + "\n"); log.flush()
return _rt(source, *a, **k)
_pq.read_table = _read_table
except ImportError:
pass
script, *rest = sys.argv[1:]
sys.argv = [script] + rest
sys.path.insert(0, os.path.dirname(os.path.abspath(script)))
runpy.run_path(script, run_name="__main__")
"""
def sha256(path):
h = hashlib.sha256()
with open(path, "rb") as fh:
for chunk in iter(lambda: fh.read(1 << 20), b""):
h.update(chunk)
return h.hexdigest()
def load_lock():
return json.loads(LOCK.read_text()) if LOCK.exists() else {"tracks": {}}
def save_lock(lock):
LOCK.write_text(json.dumps(lock, indent=1, sort_keys=True) + "\n")
def render(script, args, scratch):
"""Render into `scratch`; return ({relpath: sha} of outputs, sorted input paths, log text)."""
env = dict(os.environ)
env.update(SONIFICATION_HEAVY=str(scratch / "heavy"), SONIFICATION_LIGHT=str(scratch / "light"),
SONIFICATION_SITE_AUDIO=str(scratch / "site_audio"), RENDERS_OPEN_LOG=str(scratch / "opened.txt"))
r = subprocess.run([PY, "-c", _RUNNER, str(HERE / script), *args], cwd=ROOT, env=env,
capture_output=True, text=True)
if r.returncode:
raise SystemExit(f"render failed: {script} {' '.join(args)}\n{r.stdout[-2000:]}\n{r.stderr[-3000:]}")
outputs = {}
for kind in ("heavy", "light", "site_audio"):
base = scratch / kind
for p in sorted(base.rglob("*")) if base.exists() else []:
if p.is_file():
outputs[f"{kind}/{p.relative_to(base)}"] = sha256(p)
opened = sorted({Path(line).resolve() for line in (scratch / "opened.txt").read_text().splitlines() if line})
return outputs, opened, r.stdout
def in_place(rel):
"""Where an output's published copy lives (the same env overrides the pieces read)."""
sys.path.insert(0, str(HERE))
import dsp_core
kind, _, sub = rel.partition("/")
return {"heavy": dsp_core.HEAVY_OUT, "light": dsp_core.LIGHT_OUT, "site_audio": dsp_core.SITE_AUDIO}[kind] / sub
def inputs_record(opened):
"""Hash the data files and samples a render opened (not Python's own files)."""
import dsp_core
data, samples = {}, hashlib.sha256()
n_samples = 0
for p in opened:
if not p.is_file():
continue
if ROOT / "data" in p.parents:
data[str(p.relative_to(ROOT))] = sha256(p)
elif Path(dsp_core.SAMPLES).resolve() in p.parents:
samples.update(f"{p.relative_to(Path(dsp_core.SAMPLES).resolve())}:{sha256(p)}\n".encode())
n_samples += 1
return data, {"files": n_samples, "sha256": samples.hexdigest()}
def environment():
import numpy, scipy, soundfile
ff = subprocess.run(["ffmpeg", "-version"], capture_output=True, text=True).stdout.split("\n")[0]
return {"python": platform.python_version(), "numpy": numpy.__version__, "scipy": scipy.__version__,
"soundfile": soundfile.__version__, "ffmpeg": ff}
def cmd_lock(a):
script_text = (HERE / a.script).read_text()
seeds = sorted({int(a or b) for a, b in re.findall(r"(?:default_rng|random\.seed)\((\d+)\)|\"seed\":\s*(\d+)", script_text)})
with tempfile.TemporaryDirectory(prefix="renders-") as tmp:
outputs, opened, log = render(a.script, a.args, Path(tmp))
data, samples = inputs_record(opened)
said = re.search(r"^seeds: ([\d, ]+)$", log, re.M) # a piece whose seed depends on its arguments says which
if said:
seeds = [int(s) for s in said.group(1).split(",")]
bad = [rel for rel, h in outputs.items() if not in_place(rel).exists() or sha256(in_place(rel)) != h]
if bad and not a.new:
raise SystemExit(f"{a.track}: the scratch render differs from the files in place ({len(bad)} of "
f"{len(outputs)}), e.g. {bad[:3]}. Re-render in place first, or pass --new.")
lock = load_lock()
prev = lock["tracks"].get(a.track, {})
lock["tracks"][a.track] = {
"version": a.version, "locked": dt.date.today().isoformat(),
"script": a.script, "args": a.args, "seeds": seeds,
"code": {f: sha256(HERE / f) for f in [a.script] + ENGINE},
"data": data, "samples": samples, "outputs": outputs, "environment": environment(),
}
if "--remix" in a.args: # lineage: the parent it came from and who made it
sys.path.insert(0, str(HERE))
import variants
r = variants.load_remix(ROOT / a.args[a.args.index("--remix") + 1])
lock["tracks"][a.track].update(parent=r["parent"], credits=r["credits"], keeps=r.get("keeps"))
if "--variant" in a.args:
lock["tracks"][a.track]["variant"] = a.args[a.args.index("--variant") + 1]
if prev and prev.get("version") != a.version:
lock["tracks"][a.track]["previous"] = {"version": prev["version"], "locked": prev.get("locked")}
elif prev.get("previous"): # re-locking the same version keeps its history
lock["tracks"][a.track]["previous"] = prev["previous"]
save_lock(lock)
print(f"{a.track} v{a.version}: {len(outputs)} outputs, {len(data)} data files, {samples['files']} samples locked")
def cmd_verify(a):
lock = load_lock()["tracks"]
names = a.tracks or sorted(lock)
failed = []
for name in names:
e = lock[name]
with tempfile.TemporaryDirectory(prefix="renders-") as tmp:
outputs, opened, _ = render(e["script"], e["args"], Path(tmp))
data, samples = inputs_record(opened)
diff = sorted(set(outputs) ^ set(e["outputs"])) + [k for k in outputs if e["outputs"].get(k, outputs[k]) != outputs[k]]
code = [f for f, h in e["code"].items() if (HERE / f).exists() and sha256(HERE / f) != h]
inputs = [f"{f} (" + ("now read" if f not in e["data"] else "no longer read" if f not in data else "changed") + ")"
for f in set(data) | set(e["data"]) if data.get(f) != e["data"].get(f)]
if samples != e["samples"]:
inputs.append("samples (changed)")
ok = not diff
print(f"{'OK ' if ok else 'FAIL'} {name} v{e['version']}: {len(outputs) - len(diff)}/{len(e['outputs'])} outputs identical"
+ (f"; code changed since lock: {', '.join(code)}" if code else "")
+ (f"; inputs changed: {', '.join(sorted(inputs))}" if inputs else ""))
if not ok:
failed.append(name)
for k in diff[:10]:
print(f" differs: {k}")
if failed:
raise SystemExit(f"{len(failed)} of {len(names)} tracks did not reproduce")
def cmd_list(a):
for name, e in sorted(load_lock()["tracks"].items()):
print(f"{name:40s} v{e['version']:6s} {e['locked']} {e['script']} {' '.join(e['args'])}")
def main():
ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
sub = ap.add_subparsers(dest="cmd", required=True)
p = sub.add_parser("lock")
p.add_argument("track"), p.add_argument("version"), p.add_argument("script"), p.add_argument("args", nargs="*")
p.add_argument("--new", action="store_true", help="lock a render that has no files in place yet")
p.set_defaults(func=cmd_lock)
p = sub.add_parser("verify")
p.add_argument("tracks", nargs="*")
p.set_defaults(func=cmd_verify)
sub.add_parser("list").set_defaults(func=cmd_list)
a = ap.parse_args()
a.func(a)
if __name__ == "__main__":
main()"""
Variants and remixes (2026-09-27): named parameter overrides of a piece.
A piece declares its parameters and their defaults, then calls setup():
P, Y, NAME, META = V.setup("granby", DEFAULTS)
and reads P["bpm"], P["mute"] and so on. With no flags, P is the defaults
and NAME is the published name, so the published render is unchanged.
piece.py wy2019 another year (see timegrid.Year.parse)
piece.py --variant frozen a named variant from variants/<piece>.yaml
piece.py --remix remixes/dusk.yaml a remix: its own track, with lineage and credits
piece.py --set bpm=60 --set mute=[perc] ad hoc overrides (named __set-<hash>)
VARIANTS are for listening A/Bs and version history: variants/<piece>.yaml maps
a name to overrides plus an optional `note`. A variant renders as
<piece>_<year>__<variant>.
A REMIX is a new track built from a parent (remixes/*.yaml):
remix: granby-wy2024-dusk # the new track's id
parent: {piece: granby, track: granby-wy2024, version: "1.0"}
keeps: [the Granby and Burrell flow records, the ice record, the mapping]
changes: {bpm: 60, mute: [perc], mode: Dorian}
credits: [{name: ..., role: remix}]
notes: why this remix exists
It renders under its own id, and renders.py records its parent in the lock.
Every override key must be one the piece declares, so a typo fails loudly
instead of quietly rendering the default.
"""
import argparse
import hashlib
import json
from pathlib import Path
import yaml
import timegrid
ROOT = Path(__file__).resolve().parent.parent
VARIANTS = ROOT / "variants"
REMIX_KEYS = {"remix", "parent", "keeps", "changes", "credits", "notes"}
def load_remix(path):
r = yaml.safe_load(Path(path).read_text())
missing = {"remix", "parent", "changes", "credits"} - set(r)
extra = set(r) - REMIX_KEYS
if missing or extra:
raise SystemExit(f"{path}: missing {sorted(missing)} / unknown {sorted(extra)} (keys: {sorted(REMIX_KEYS)})")
return r
def _check(piece, defaults, over, where):
bad = set(over) - set(defaults)
if bad:
raise SystemExit(f"{where}: {piece} has no parameter(s) {sorted(bad)}; it has {sorted(defaults)}")
def setup(piece, defaults, argv=None):
"""-> (params, Year, output name, meta). meta is {} for the published default."""
ap = argparse.ArgumentParser(description=f"render {piece}")
ap.add_argument("year", nargs="?", help="e.g. wy2024 (default), wy2019, cy2023")
g = ap.add_mutually_exclusive_group()
g.add_argument("--variant", help=f"a name from variants/{piece}.yaml")
g.add_argument("--remix", help="a remix spec (remixes/*.yaml)")
ap.add_argument("--set", action="append", default=[], metavar="KEY=VALUE", help="override one parameter")
a = ap.parse_args(argv)
P, over, meta, tag, name = dict(defaults), {}, {}, "", None
if a.variant:
spec = yaml.safe_load((VARIANTS / f"{piece}.yaml").read_text())
if a.variant not in spec:
raise SystemExit(f"no variant {a.variant!r} in variants/{piece}.yaml; there are {sorted(spec)}")
over = {k: v for k, v in spec[a.variant].items() if k != "note"}
_check(piece, defaults, over, f"variants/{piece}.yaml:{a.variant}")
meta = {"variant": a.variant, "overrides": over, "note": spec[a.variant].get("note")}
tag = f"__{a.variant}"
if a.remix:
r = load_remix(a.remix)
if r["parent"].get("piece") != piece:
raise SystemExit(f"{a.remix} is a remix of {r['parent'].get('piece')!r}, not {piece!r}")
over = dict(r["changes"])
_check(piece, defaults, over, a.remix)
meta = {"remix": r["remix"], "parent": r["parent"], "overrides": over, "credits": r["credits"]}
name = r["remix"].replace("-", "_")
extra = {}
for s in a.set:
k, _, v = s.partition("=")
extra[k] = yaml.safe_load(v)
if extra:
_check(piece, defaults, extra, "--set")
over.update(extra)
meta.setdefault("overrides", {}).update(extra)
tag += "__set-" + hashlib.sha256(json.dumps(extra, sort_keys=True).encode()).hexdigest()[:6]
if a.year:
over["year"] = P["year"] = a.year
if meta:
meta["overrides"]["year"] = a.year
P.update(over)
Y = timegrid.Year.parse(P["year"])
if name is None:
name = f"{piece}_{Y.label.lower()}{tag}"
elif tag:
name += tag
return P, Y, name, meta
def bus_settings(base, P):
"""Apply P['gains'] to a {bus: settings} table, checking every name in
P['gains'] and P['mute'] is a layer. Muted layers are still built (the
score code writes to them); the piece leaves them out of the master."""
unknown = (set(P.get("mute", [])) | set(P.get("gains", {}))) - set(base)
if unknown:
raise SystemExit(f"no layer(s) {sorted(unknown)}; the layers are {sorted(base)}")
out = {}
for k, kw in base.items():
kw = dict(kw)
if k in P.get("gains", {}):
kw["gain"] = float(P["gains"][k])
out[k] = kw
return out
Note
seasonal.py: the older single module, now re-exporting the three above
"""
Compatibility shim (2026-09-27). The engine that lived here is now three
modules, and this file re-exports them so older scripts and the published
repo keep working:
timegrid.py Year (the year and its framing), Grid, per_bar, daily helpers
harmony.py the mode ladder, chords, voice leading, progressions (pure)
records.py cached-pull loaders and the per-year daily arrays
plus sampler.sustain (sample playback). The module-level DAYS, N_DAYS,
N_BARS and WY_START are the WY2024 defaults; new code should build a
timegrid.Year and pass it (the pilots take --year, e.g. wy2023, cy2024).
"""
import datetime as dt # noqa: F401 (pieces use W.dt)
from dsp_core import SR, HEAVY_OUT, LIGHT_OUT, SITE_AUDIO, day_length_hours # noqa: F401
from timegrid import ( # noqa: F401
DAYS_PER_BAR, FRAMINGS, WY2024, Grid, Year, events_from_flag, fmt_time, new_moons_from,
norm01, per_bar, smooth, solstice,
)
from harmony import ( # noqa: F401
MODES, chord, chord_track, degree_to_midi, fire_bars, mode_track, phase_bars, voice_lead,
)
from records import ( # noqa: F401
DATA, HERE, check_coverage, climate_series, coverage, dart_run_totals, dart_sockeye, doy_percentile, growing_degree_days,
in_year, modis_gpp_daily, wsc_series,
)
from sampler import sustain # noqa: F401
WY_START = WY2024.start
DAYS = WY2024.days
N_DAYS = WY2024.n_days
N_BARS = WY2024.n_bars
wy = in_year # the old name, from when every piece was a water year
def doy_index(doy):
return WY2024.doy_index(doy)
def daylength(lat):
return WY2024.daylength(lat)
Note
sampler.py: instruments, buses, reverb, master
"""
Sample-based instruments, mixing buses, and convolution reverb for the
seasonal pilots (2026-09-24). Replaces "synthesize every sound from sine
waves" with recorded acoustic instruments (VSCO 2 CE, CC0) + real bird
recordings (Commons, CC BY-SA). See samples_manifest.csv.
Pieces:
Instrument - folder of pitched samples -> note(midi, vel, dur)
OneShots - folder of unpitched hits -> hit(name_prefix, vel)
Bus - stereo track with gain/EQ/pan and a reverb send
make_ir - synthetic stereo impulse response (room/hall), bright/dark
master - sum buses, apply reverb returns, gentle glue compression
"""
import re
import subprocess
from functools import lru_cache
from pathlib import Path
import numpy as np
import soundfile as sf
from scipy.signal import butter, sosfilt, fftconvolve
from dsp_core import SR, SAMPLES # SAMPLES: repo samples/ unless SONIFICATION_SAMPLES is set
VSCO = SAMPLES / "vsco2ce"
NOTE_IDX = {"C": 0, "C#": 1, "Db": 1, "D": 2, "D#": 3, "Eb": 3, "E": 4, "F": 5, "F#": 6,
"Gb": 6, "G": 7, "G#": 8, "Ab": 8, "A": 9, "A#": 10, "Bb": 10, "B": 11}
def name_to_midi(name, octave_offset=0):
m = re.search(r"(?<![A-Za-z])([A-G][#b]?)(-?\d)(?!\d)", name)
if not m:
return None
return 12 * (int(m.group(2)) + 1 + octave_offset) + NOTE_IDX[m.group(1)]
@lru_cache(maxsize=None)
def load(path):
x, sr = sf.read(str(path), always_2d=True, dtype="float32")
if sr != SR:
raise ValueError(f"{path}: sr {sr}")
if x.shape[1] == 1:
x = np.repeat(x, 2, axis=1)
return x[:, :2]
def estimate_f0(x, sr=SR, fmin=25, fmax=2000):
"""Autocorrelation f0 over the loudest 0.4 s of the sample (after attack)."""
m = x.mean(axis=1)
start = int(0.15 * sr)
seg = m[start:start + int(0.4 * sr)]
seg = seg - seg.mean()
ac = fftconvolve(seg, seg[::-1])[len(seg) - 1:]
lo, hi = int(sr / fmax), int(sr / fmin)
neg = np.flatnonzero(ac[:hi] < 0) # skip the zero-lag lobe: for low notes it otherwise
lo = max(lo, int(neg[0])) if len(neg) else lo # wins (a 55 Hz tone read as 2 kHz; test_sound.py)
lag = lo + np.argmax(ac[lo:hi])
return sr / lag
class Instrument:
"""Pitched sampled instrument. Picks the nearest recorded note (and the
dynamic layer nearest `vel`), then repitches by resampling. Resampling
also shortens/lengthens the note - fine for shifts of +-2-3 semitones,
which is all a well-sampled set needs."""
def __init__(self, folder, glob="*.wav", octave_offset=0, mapping=None, dyn_regex=r"_v(\d)",
gain=1.0, attack_trim=0.0):
self.folder = VSCO / folder if not Path(folder).is_absolute() else Path(folder)
self.gain = gain
self.attack_trim = attack_trim
self.zones = [] # (midi, dyn, path)
for p in sorted(self.folder.glob(glob)):
if mapping:
midi, dyn = mapping(p.name)
else:
midi = name_to_midi(p.stem, octave_offset)
d = re.search(dyn_regex, p.stem)
dyn = int(d.group(1)) if d else 1
if midi is not None:
self.zones.append((midi, dyn, p))
if not self.zones:
raise FileNotFoundError(self.folder)
self.dyns = sorted({z[1] for z in self.zones})
def _pick(self, midi, vel, rng):
want_dyn = self.dyns[min(len(self.dyns) - 1, int(vel * len(self.dyns)))]
cands = [z for z in self.zones if z[1] == want_dyn] or self.zones
best = min(abs(z[0] - midi) for z in cands)
pool = [z for z in cands if abs(z[0] - midi) == best]
return pool[rng.integers(len(pool))] if rng is not None else pool[0]
def note(self, midi, vel=0.7, dur=None, release=0.6, rng=None):
src_midi, _, path = self._pick(midi, vel, rng)
x = load(path)
if self.attack_trim:
x = x[int(self.attack_trim * SR):]
ratio = 2 ** ((midi - src_midi) / 12.0)
if dur is not None: # don't resample audio we're about to cut off
x = x[: int((dur + release) * SR * ratio) + 2]
if abs(ratio - 1) > 1e-4:
n_out = int(len(x) / ratio)
idx = np.arange(n_out) * ratio
x = np.stack([np.interp(idx, np.arange(len(x)), x[:, c]) for c in range(2)], axis=1)
if dur is not None:
n = int((dur + release) * SR)
x = x[:n].copy()
r = int(release * SR)
if len(x) > r:
x[-r:] *= np.linspace(1, 0, r)[:, None] ** 2
# velocity -> level (dynamic layers already carry the timbre change)
return x * (0.35 + 0.65 * vel) * self.gain
class OneShots:
def __init__(self, folder=VSCO / "Percussion", glob="*.wav"):
self.files = sorted(Path(folder).glob(glob))
def hit(self, prefix, vel=0.7, rng=None):
pool = [f for f in self.files if f.name.startswith(prefix)]
if not pool:
raise FileNotFoundError(prefix)
p = pool[rng.integers(len(pool))] if rng is not None else pool[0]
return load(p) * vel
# ---------------------------------------------------------------- birds ----
@lru_cache(maxsize=None)
def bird_phrases(name, min_len=0.6, max_len=4.0, top=12):
"""Decode a Commons bird recording, high-pass it (drop wind/traffic rumble),
and cut it into its loudest song phrases. Returns a list of stereo arrays."""
src = next((SAMPLES / "birds").glob(name + ".*"))
wav = SAMPLES / "birds" / "_decoded" / (name + ".wav")
wav.parent.mkdir(exist_ok=True)
if not wav.exists():
subprocess.run(["ffmpeg", "-loglevel", "error", "-y", "-i", str(src), "-ar", str(SR), "-ac", "2", str(wav)],
check=True)
x = load(wav)
x = sosfilt(butter(4, 900, "highpass", fs=SR, output="sos"), x, axis=0)
env = np.convolve(np.abs(x.mean(axis=1)), np.ones(2205) / 2205, mode="same")
thr = np.percentile(env, 80)
on = env > thr
# merge gaps < 0.25 s
phrases, i, n = [], 0, len(on)
while i < n:
if on[i]:
j = i
gap = 0
while j < n and gap < int(0.25 * SR):
gap = 0 if on[j] else gap + 1
j += 1
a, b = max(0, i - int(0.08 * SR)), min(n, j)
if (b - a) / SR >= min_len:
seg = x[a:a + int(max_len * SR)].copy()
f = min(len(seg) // 4, int(0.05 * SR))
seg[:f] *= np.linspace(0, 1, f)[:, None]
seg[-f:] *= np.linspace(1, 0, f)[:, None]
phrases.append((env[a:b].mean(), seg))
i = j
else:
i += 1
phrases.sort(key=lambda t: -t[0])
out = [p / (np.abs(p).max() + 1e-9) for _, p in phrases[:top]]
if not out:
raise ValueError(f"no phrases found in {name}")
return out
# ---------------------------------------------------------------- buses ----
def eq(x, hp=None, lp=None):
if hp:
x = sosfilt(butter(2, hp, "highpass", fs=SR, output="sos"), x, axis=0)
if lp:
x = sosfilt(butter(2, lp, "lowpass", fs=SR, output="sos"), x, axis=0)
return x
def pan_stereo(x, pan):
"""pan -1..1, equal-power. x is (n,2)."""
a = (pan + 1) * np.pi / 4
return x * np.array([np.cos(a), np.sin(a)]) * np.sqrt(2)
class Bus:
def __init__(self, name, n, gain=1.0, hp=None, lp=None, send=0.0, ir="hall"):
self.name, self.gain, self.hp, self.lp, self.send, self.ir = name, gain, hp, lp, send, ir
self.buf = np.zeros((n, 2), dtype=np.float32)
def add(self, x, at, gain=1.0, pan=0.0):
at = int(at)
if at >= len(self.buf) or at < 0:
return
x = x[: len(self.buf) - at]
if pan:
x = pan_stereo(x, pan)
self.buf[at:at + len(x)] += x * gain
def processed(self):
return eq(self.buf, self.hp, self.lp) * self.gain
def make_ir(seconds=3.2, bright=0.5, predelay=0.02, seed=7):
"""Synthetic stereo hall IR: decorrelated noise per channel with a
frequency-dependent decay (highs die ~3x faster than lows, as in a real
room), plus a handful of early reflections. `bright` 0..1 sets the
high-band share - 'snow on the ground' uses a dark, short one."""
rng = np.random.default_rng(seed)
n = int(seconds * SR)
t = np.arange(n) / SR
ir = np.zeros((n, 2))
bands = [(None, 400, 1.0), (400, 2500, 0.7), (2500, None, 0.35)]
for lo, hi, life in bands:
for c in range(2):
b = rng.standard_normal(n)
b = eq(b[:, None], hp=lo, lp=hi)[:, 0]
t60 = seconds * life * (1.0 if lo is None else (0.5 + 0.5 * bright))
w = 1.0 if lo is None else (0.4 + 0.8 * bright) * (0.7 if hi is None else 1.0)
ir[:, c] += b * np.exp(-6.9 * t / t60) * w
for k in range(8):
d = int((predelay + rng.uniform(0.005, 0.07)) * SR)
ir[d, rng.integers(2)] += rng.uniform(0.3, 0.8) * (-1) ** k
pd = int(predelay * SR)
ir = np.concatenate([np.zeros((pd, 2)), ir])[:n]
fade = np.linspace(1, 0, int(0.2 * SR)) ** 2
ir[-len(fade):] *= fade[:, None]
return ir / np.sqrt((ir ** 2).sum(axis=0).mean())
def convolve_reverb(x, ir):
return np.stack([fftconvolve(x[:, c], ir[:, c])[: len(x)] for c in range(2)], axis=1)
def glue(x, thresh_db=-14, ratio=2.0, attack=0.02, release=0.3):
"""Very gentle RMS bus compressor, block-based, so the whole mix breathes
together without any one layer punching through."""
blk = 1024
n = len(x)
rms = np.sqrt(np.convolve((x ** 2).mean(axis=1), np.ones(blk) / blk, mode="same") + 1e-12)
db = 20 * np.log10(rms)
over = np.maximum(0, db - thresh_db)
gr_db = -over * (1 - 1 / ratio)
# smooth gain reduction (attack/release one-pole)
g = np.empty(n)
a_c, r_c = np.exp(-1 / (attack * SR)), np.exp(-1 / (release * SR))
cur = 0.0
for i in range(0, n, 64): # control-rate smoothing, 64-sample hops
tgt = gr_db[i]
c = a_c if tgt < cur else r_c
cur = c * cur + (1 - c) * tgt
g[i:i + 64] = cur
return x * (10 ** (g / 20))[:, None]
def master(buses, irs, out_wav, stems_dir=None, ir_env=None, wet_gain=0.35):
"""ir_env: optional {ir_name: per-sample 0..1 envelope}. Every bus send
goes to EVERY reverb, and the envelopes crossfade between them - e.g.
'open' hall vs 'snow' (dark, short) following snow on the ground."""
n = len(buses[0].buf)
dry = np.zeros((n, 2))
send_sum = np.zeros((n, 2))
stems = {}
for b in buses:
y = b.processed()
dry += y
send_sum += y * b.send
stems[b.name] = y
if stems_dir:
Path(stems_dir).mkdir(parents=True, exist_ok=True)
sf.write(Path(stems_dir) / f"{b.name}.flac", (y / 4).astype(np.float32), SR)
send_sum = eq(send_sum, hp=180)
wet = np.zeros((n, 2))
for k, ir in irs.items():
env = ir_env[k][:, None] if ir_env else 1.0 / len(irs)
wet += convolve_reverb(send_sum, ir) * env * wet_gain
mix = glue(dry + wet)
mix = mix / (np.abs(mix).max() + 1e-9) * 0.89
f = int(3 * SR)
mix[:f] *= np.linspace(0, 1, f)[:, None]
mix[-int(6 * SR):] *= np.linspace(1, 0, int(6 * SR))[:, None] ** 1.5
sf.write(out_wav, mix.astype(np.float32), SR)
return mix, stems
def to_mp3(wav, mp3):
subprocess.run(["ffmpeg", "-loglevel", "error", "-y", "-i", str(wav), "-af", "loudnorm=I=-16:TP=-1.5:LRA=11",
"-ar", "44100", "-b:a", "192k", str(mp3)], check=True)
def vsco_kit():
"""The pilot instrument set. octave_offset=1 where VSCO names notes with
C3 = middle C (verified 2026-09-24 by autocorrelation f0 against the
filename: cello, violins, contrabass, flute, marimba, glock read exactly
12 semitones high vs. their labels; piano + harp matched as-is)."""
import re as _re
piano_map = lambda n: (21 + 2 * int(_re.search(r"_(\d+)\.wav$", n).group(1)),
int(_re.search(r"dyn(\d)", n).group(1)))
return {
"piano": Instrument("Keys/Upright Piano", gain=0.9, mapping=piano_map),
"bass": Instrument("Strings/Solo Contrabass/SusNV", octave_offset=1),
"bpizz": Instrument("Strings/Solo Contrabass/Pizz", octave_offset=1),
"cello": Instrument("Strings/Cello Section/susvib", octave_offset=1),
"cpizz": Instrument("Strings/Cello Section/pizzT", octave_offset=1),
"vln": Instrument("Strings/Violin Section/susVib", octave_offset=1),
"vpizz": Instrument("Strings/Violin Section/Pizz", octave_offset=1),
"flute": Instrument("Woodwinds/Flute/susNV", octave_offset=1),
"harp": Instrument("Strings/Harp"),
"glock": Instrument("Percussion/Glock", octave_offset=1),
"marimba": Instrument("Percussion/Marimba", octave_offset=1),
"perc": OneShots(),
}
def sustain(ins, midi, vel, dur, rng, seg=4.5, xf=1.2, release=1.5):
"""Hold a sampled note longer than its recording by re-articulating
overlapping takes with crossfades (string sections do exactly this in
real life - players change bow at different times). Moved here from
seasonal.py on 2026-09-27."""
n = int((dur + release) * SR)
out = np.zeros((n, 2), dtype=np.float32)
t = 0.0
first = True
while t < dur:
piece = ins.note(midi, vel, dur=min(seg, dur - t) + xf, release=xf, rng=rng)
if not first:
f = min(len(piece), int(xf * SR))
piece = piece.copy()
# skip the new take's bow attack and fade it in under the old one
piece[:f] *= np.linspace(0, 1, f)[:, None]
a = int(t * SR)
piece = piece[: n - a]
out[a:a + len(piece)] += piece
t += seg
first = False
r = int(release * SR)
out[-r:] *= np.linspace(1, 0, r)[:, None] ** 2
return out
Note
pilot_granby_wy2024.py: the Granby score
"""
PILOT A - Granby River, water year 2023-24 (Oct 1 2023 - Sep 30 2024).
Any other year the records cover: `pilot_granby_wy2024.py wy2019` (WY2011-WY2024;
the arrangement was tuned on WY2024, and the phase thresholds assume a water year).
Key of D, 72 bpm, 1 bar = 3 days, ~7 minutes.
A warm El Nino winter: the Granby gauge flagged only 28 ice days, in broken
episodes, and the freshet came early and small (peak 147 m3/s vs 465 in the
2010-24 record). The piece should sound like a year that never fully
committed to winter.
LAYER DATA (R = retrieved, D = derived, S = simulated) SOUND
piano ostinato Granby flow, in-year level (R) density: half notes -> 8ths
harmony / mode day length (D) + Granby flow vs 15-yr normal (R) mode brightness, 1 step per 8 bars
progression hydrologic phase from Granby flow (D) which 4-chord loop plays
bass / cello Granby flow (R) contrabass root; cello 5th above at high flow
string pad MODIS GPP, Granby valley (R) [GDD fallback (D)] violins swell with plant productivity
ice WSC "Ice Conditions" flag (R) + Billings temp (R) freeze: glock + bowed cymbal; crackle + singing-ice
chirps while frozen; breakup: cymbal swell + harp
rain / snow Billings TOTAL_RAIN / TOTAL_SNOW (R) rain = harp drops, snow = glockenspiel flakes
reverb space Billings snow on ground (R) snow cover = dark, short room (snow absorbs sound)
hand percussion Burrell Creek flow + flashiness (R) tumba/conga/shaker, tresillo 3-3-2 at peak
redband trout spawning window from temp proxy + season (S) flute: short rising figure
whitefish fall spawning window from temp proxy (S) cello pizzicato: falling figure
birds typical phenology windows (S timing, R recordings) dipper (winter), thrush + bank swallow (summer)
"""
import json
import sys
from pathlib import Path
import numpy as np
sys.path.insert(0, str(Path(__file__).resolve().parent))
import sampler as S
from sampler import Bus, bird_phrases, make_ir, SR
import seasonal as W
import variants as V
from timegrid import per_bar, smooth, norm01
from dsp_core import ice_crinkle_grain
# Named parameters: variants/granby.yaml and remixes override these (variants.py)
DEFAULTS = {
"year": "wy2024",
"bpm": 72,
"tonic": 50, # MIDI; D3
"mode": None, # a mode name (e.g. "Dorian") freezes the harmony's mood
"mode_range": [1, 4], # Aeolian..Ionian, the ladder's span for season + anomaly
"mute": [], # layers left out of the mix
"gains": {}, # layer -> gain, replacing the defaults below
"seed": 2024,
}
def main():
OUT = W.LIGHT_OUT / "pilots"
OUT.mkdir(parents=True, exist_ok=True)
HEAVY = W.HEAVY_OUT / "pilots" # WAV + stems (large) live in runtime
HEAVY.mkdir(parents=True, exist_ok=True)
P, Y, NAME, META = V.setup("granby", DEFAULTS)
DAYS, N_DAYS, N_BARS = Y.days, Y.n_days, Y.n_bars
rng = np.random.default_rng(P["seed"])
np.random.seed(P["seed"]) # the ice grain helper uses the global RNG; seed it so renders repeat exactly
TONIC = P["tonic"]
BPM = P["bpm"]
g = W.Grid(BPM, Y.n_bars)
print(f"{N_BARS} bars, {g.dur / 60:.1f} min")
# --------------------------------------------------------------- features --
q_all, q_sym = W.wsc_series("granby_08NN002_2010_2024.json")
b_all, _ = W.wsc_series("burrell_08NN023_2010_2024.json")
clim = "billings_1100_climate_2010_2024.json"
clim_T = W.climate_series(clim, "MEAN_TEMPERATURE")
W.check_coverage({"Granby flow": q_all, "Burrell flow": b_all, "Billings temperature": clim_T}, Y)
q = W.in_year(q_all, year=Y)
qn = norm01(np.log(q)) # in-year level 0..1
qpct = W.doy_percentile(q_all, year=Y) # vs. 15-yr normal
ice = np.array([q_sym.get(str(d)) == "Ice Conditions" for d in DAYS])
bq = W.in_year(b_all, year=Y)
bn = norm01(np.log(bq))
bflash = norm01(np.abs(np.diff(np.log(bq), prepend=np.log(bq[0]))))
T = W.in_year(clim_T, year=Y)
rain = W.in_year(W.climate_series(clim, "TOTAL_RAIN"), fill="zero", year=Y)
snow = W.in_year(W.climate_series(clim, "TOTAL_SNOW"), fill="zero", year=Y)
sog = W.in_year(W.climate_series(clim, "SNOW_ON_GROUND"), fill="zero", year=Y)
DL = Y.daylength(49.03)
gpp, gppn = W.modis_gpp_daily("granby_valley", Y)
GPP_SOURCE = "MODIS MOD17A2HGF (retrieved)"
if gpp is None:
gdd = W.growing_degree_days(T)
gppn = np.clip(smooth(np.maximum(0, T - 5), 21) / 12, 0, 1) * (gdd > 60)
GPP_SOURCE = "growing-degree proxy from Billings temperature (derived; MODIS pull unavailable)"
doy = np.array([d.timetuple().tm_yday for d in DAYS])
Tproxy = smooth(T, 14) # stand-in for water temperature (no Granby stream temp on record)
# ---------------------------------------------------------------- harmony --
bright = per_bar(0.7 * norm01(DL, 0, 100) + 0.3 * norm01(smooth(T, 15)))
anom = per_bar(smooth(qpct, 15)) - 0.5
modes = W.mode_track(bright, anom, lo=P["mode_range"][0], hi=P["mode_range"][1])
if P["mode"]:
modes = np.full(len(modes), [m[0] for m in W.MODES].index(P["mode"]))
# phase uses LINEAR flow share of the year's peak (log flow, used for the
# piano, would call a 20 m3/s February bump a "peak")
qlin = smooth(q / q.max(), 9)
dq = np.gradient(qlin)
phase_day = np.where(ice | ((qlin < 0.1) & ((doy > 300) | (doy < 60))), "dormant",
np.where(qlin > 0.5, "peak",
np.where((dq > 0.004) & (qlin > 0.1), "rising",
np.where((dq < -0.003) & (qlin > 0.15), "recession", "low"))))
PROG = { # scale degrees (0 = tonic). Each loop is 4 chords x 2 bars.
"dormant": [0, 5, 0, 3], # i VI i iv - still, circling home
"rising": [0, 3, 4, 5], # i iv v VI - climbs, doesn't resolve
"peak": [5, 3, 0, 4], # VI iv i v - wide, 'anthemic'
"recession": [0, 6, 5, 3], # i VII VI iv - bass steps down: letting go
"low": [0, 3, 0, 5], # i iv i VI - sparse, late-summer
}
bar_phase = W.phase_bars(phase_day, min_run=6)
chords = W.chord_track(bar_phase, PROG, bars_per_chord=2)
# ------------------------------------------------------------ instruments --
K = S.vsco_kit()
piano, bass, cello, cpizz, vln = K["piano"], K["bass"], K["cello"], K["cpizz"], K["vln"]
flute, harp, glock, perc = K["flute"], K["harp"], K["glock"], K["perc"]
n = g.n
B = {k: Bus(k, n, **kw) for k, kw in V.bus_settings({
"piano": dict(gain=1.0, hp=70, send=0.35),
"bass": dict(gain=0.55, lp=900, send=0.12),
"pad": dict(gain=0.28, hp=220, send=0.6),
"ice": dict(gain=1.5, hp=500, send=0.5),
"precip": dict(gain=1.05, hp=250, send=0.55),
"perc": dict(gain=2.0, hp=60, send=0.15),
"fish": dict(gain=0.24, hp=150, send=0.5),
"birds": dict(gain=0.2, hp=900, send=0.35),
}, P).items()}
events = [] # (seconds, label) -> listening guide
def mark(sec, text):
events.append((round(sec, 1), text))
def tones(b, lo, hi, prev=None, ext=(0, 2, 4, 8)):
return W.voice_lead(prev, W.chord(chords[b], modes[b], TONIC, ext), lo, hi)
# ---------------------------------------------------------- piano + bass --
PATTERN = [0, 2, 3, 2, 1, 2, 3, 2] # indexes into the 4-note voicing, low->high
DENSITY = {0: [0, 4], 1: [0, 2, 4, 6], 2: list(range(8)), 3: list(range(8))}
qn_bar, dl_bar = per_bar(qn), per_bar(norm01(DL, 0, 100))
prev = None
for b in range(N_BARS):
v = tones(b, 55, 77, prev)
prev = v
lvl = int(np.sqrt(qn_bar[b]) * 3.99)
base_vel = 0.28 + 0.3 * dl_bar[b] + 0.15 * qn_bar[b]
for step in DENSITY[lvl]:
vel = np.clip(base_vel + (0.1 if step in (0, 4) else 0) + rng.normal(0, 0.03), 0.1, 0.95)
B["piano"].add(piano.note(v[PATTERN[step]], vel, dur=g.beat * 1.3, release=0.9, rng=rng),
g.s(b, step * 0.5, jitter_ms=7, rng=rng), pan=-0.1 + 0.05 * (PATTERN[step] - 1.5))
if lvl == 3 and step in (3, 7): # a soft upper-octave echo at full flow
B["piano"].add(piano.note(v[PATTERN[step]] + 12, vel * 0.5, dur=g.beat, release=0.8, rng=rng),
g.s(b, step * 0.5 + 0.25, jitter_ms=5, rng=rng), pan=0.25)
# left hand: low root each bar (the 'ground' under the ostinato)
root_lo = W.voice_lead(None, [W.degree_to_midi(chords[b], modes[b], TONIC)], 38, 50)[0]
B["piano"].add(piano.note(root_lo, 0.35 + 0.2 * qn_bar[b], dur=g.bar * 0.9, release=1.2, rng=rng),
g.s(b, 0, jitter_ms=4, rng=rng), pan=-0.2)
if b % 2 == 0: # chord change -> sustained strings
vel = 0.3 + 0.5 * qn_bar[b]
B["bass"].add(W.sustain(bass, root_lo - 12 if root_lo > 45 else root_lo, vel, g.bar * 2, rng),
g.s(b), gain=0.9)
if qn_bar[b] > 0.5:
B["bass"].add(W.sustain(cello, root_lo + 7, vel * 0.8, g.bar * 2, rng), g.s(b), gain=0.5, pan=0.2)
# ------------------------------------------------------ plant productivity --
gpp_bar = per_bar(gppn)
prev = None
for b in range(0, N_BARS, 2):
if gpp_bar[b] < 0.08:
prev = None
continue
v = tones(b, 64, 86, prev, ext=(2, 4, 8))
prev = v
for i, m in enumerate(v):
B["pad"].add(W.sustain(vln, m, 0.2 + 0.6 * gpp_bar[b], g.bar * 2, rng), g.s(b), gain=gpp_bar[b] ** 1.2,
pan=(-0.5, 0.1, 0.5)[i])
first_green = next((b for b in range(N_BARS) if gpp_bar[b] >= 0.08 and b > 40), None) # spring, not the October tail
# ------------------------------------------------------------------- ice --
runs, _ = W.events_from_flag(ice, min_gap=3)
for on, off in runs:
b = on // 3
top = tones(b, 79, 96)[::-1]
for i, m in enumerate(top[:3]): # descending glock: the surface closing
B["ice"].add(glock.note(m, 0.45, dur=1.5, release=2.0, rng=rng), g.s(b, i * 0.5), gain=0.7, pan=0.3 - 0.3 * i)
B["ice"].add(perc.hit("susCymb1-bow", 0.5, rng), g.s(b), gain=0.6)
mark(g.day_to_time(on), f"Ice forms on the Granby ({DAYS[on]:%b %d}) - glockenspiel + bowed cymbal")
for d in range(on, off + 1):
cold = np.clip(-T[d] / 15, 0.1, 1)
t0 = g.day_to_time(d)
for _ in range(int(3 + 10 * cold)):
grain = ice_crinkle_grain(SR, dur=0.03, intensity=cold)[:, None].repeat(2, 1)
B["ice"].add(grain * 0.25, (t0 + rng.uniform(0, g.day_to_time(1))) * SR, pan=rng.uniform(-0.8, 0.8))
if T[d] < -4 and rng.random() < 0.45: # 'singing ice': dispersive descending chirp
L = 1.6
tt = np.arange(int(L * SR)) / SR
f = 2400 * np.exp(-tt * 1.6) + 350
chirp = np.sin(2 * np.pi * np.cumsum(f) / SR) * np.exp(-tt * 2.2) * (1 - np.exp(-tt * 60))
B["ice"].add(chirp[:, None].repeat(2, 1) * 0.18, (t0 + rng.uniform(0, 1)) * SR, pan=rng.uniform(-0.6, 0.6))
bo = min(off + 1, N_DAYS - 1)
swell = perc.hit("susCymb1-cresc-Median", 0.6, rng)
t_break = g.day_to_time(bo)
t_freeze = g.day_to_time(on) + 1.0
B["ice"].add(swell, max(t_freeze * SR, t_break * SR - len(swell) * 0.8), gain=0.6 if off - on >= 3 else 0.3)
arp = tones(bo // 3, 57, 91)
arp = sorted(set(arp + [m + 12 for m in arp]))
for i, m in enumerate(arp): # rising harp: the river opening
B["ice"].add(harp.note(m, 0.5, rng=rng), t_break * SR + i * 0.09 * SR, gain=0.6, pan=-0.4 + 0.1 * i)
mark(t_break, f"Ice breaks up ({DAYS[bo]:%b %d}) - cymbal swell + rising harp")
# -------------------------------------------------------------- rain/snow --
for d in range(N_DAYS):
b = min(d // 3, N_BARS - 1)
t0 = g.day_to_time(d)
if rain[d] > 0.2:
pool = tones(b, 67, 91)
for _ in range(min(9, int(round(1.6 * np.sqrt(rain[d]))))):
B["precip"].add(harp.note(int(rng.choice(pool)), rng.uniform(0.25, 0.55), rng=rng),
(t0 + rng.uniform(0, g.day_to_time(1))) * SR, pan=rng.uniform(-0.7, 0.7))
if snow[d] > 0.2:
pool = tones(b, 81, 98)
for _ in range(min(7, int(round(1.3 * np.sqrt(snow[d]))))):
B["precip"].add(glock.note(int(rng.choice(pool)), 0.22, dur=1.0, release=2.5, rng=rng),
(t0 + rng.uniform(0, g.day_to_time(1))) * SR, gain=0.6, pan=rng.uniform(-0.8, 0.8))
big_rain = int(np.argmax(rain))
mark(g.day_to_time(big_rain), f"Wettest day, {rain[big_rain]:.0f} mm rain ({DAYS[big_rain]:%b %d}) - harp shower")
big_snow = int(np.argmax(snow))
mark(g.day_to_time(big_snow), f"Biggest snowfall, {snow[big_snow]:.0f} cm ({DAYS[big_snow]:%b %d}) - glockenspiel flurry")
# ------------------------------------------------------- hand percussion --
active = np.clip((smooth(qn, 7) - 0.25) / 0.2, 0, 1)
drive = np.sqrt(np.clip(0.6 * bn + 0.4 * bflash, 0, 1)) * active
drive_bar = per_bar(drive)
SW = 0.12 # swing: offbeat 8ths land a little late
for b in range(N_BARS):
lv = drive_bar[b]
if lv < 0.12:
continue
for step in range(8):
beat = step * 0.5 + (SW * 0.5 if step % 2 else 0)
at = lambda: g.s(b, beat, jitter_ms=6, rng=rng)
if step in (0, 5):
B["perc"].add(perc.hit("Tumba-HitN", 0.5 + 0.4 * lv, rng), at(), pan=-0.3)
if lv > 0.3 and (step in (2, 6) or (lv > 0.6 and step == 3)):
B["perc"].add(perc.hit("Conga-HitN", 0.4 + 0.4 * lv, rng), at(), pan=0.25)
if lv > 0.45 and step % 2:
B["perc"].add(perc.hit("Tamb1-Shake", 0.15 + 0.25 * lv, rng), at(), pan=0.45)
if lv > 0.5 and step == 7:
B["perc"].add(perc.hit("Quinto-Tap1", 0.3, rng), at(), pan=0.5)
if lv > 0.78 and step in (0, 3, 6): # tresillo 3-3-2 against the 4/4
B["perc"].add(perc.hit("LogDrumLo", 0.55, rng), at(), pan=-0.5)
peak_d = int(np.argmax(q))
mark(g.day_to_time(peak_d), f"Freshet peak, {q[peak_d]:.0f} m3/s ({DAYS[peak_d]:%b %d}) - full ostinato + tresillo drums")
first_perc = next((b for b in range(N_BARS) if drive_bar[b] >= 0.12), None)
if first_perc is not None:
mark(g.t(first_perc), f"Burrell Creek wakes ({DAYS[first_perc * 3]:%b %d}) - hand drums enter")
# ------------------------------------------------------------------- fish --
in_window = lambda a, b_: ((doy >= a) & (doy <= b_)).astype(float)
trout = np.exp(-((Tproxy - 8) / 2.5) ** 2) * in_window(91, 181) * (qn > 0.3)
whitefish = np.exp(-((Tproxy - 4.5) / 2.0) ** 2) * in_window(288, 344)
tr_bar, wf_bar = per_bar(trout), per_bar(whitefish)
first_t = first_w = None
tr_fire, wf_fire = W.fire_bars(tr_bar, 0.6), W.fire_bars(wf_bar, 0.75)
for b in range(N_BARS):
if b in tr_fire:
v = tones(b, 72, 90, ext=(4, 7, 8)) # 5th -> octave -> 9th, rising
for i, m in enumerate(sorted(v)):
B["fish"].add(flute.note(m, 0.5, dur=(0.45 if i < 2 else 1.6) * g.beat * 2, release=0.6, rng=rng),
g.s(b, 1 + i * 0.5, jitter_ms=5, rng=rng), pan=0.35)
first_t = first_t if first_t is not None else b
if b in wf_fire:
v = sorted(tones(b, 50, 64, ext=(0, 2, 4)), reverse=True)
for i, m in enumerate(v):
B["fish"].add(cpizz.note(m, 0.55, rng=rng), g.s(b, 2 + i * 0.5, jitter_ms=5, rng=rng), pan=-0.35)
first_w = first_w if first_w is not None else b
if first_t is not None:
mark(g.t(first_t), f"Redband trout spawning window opens ({DAYS[first_t * 3]:%b %d}, simulated) - flute figure")
if first_w is not None:
mark(g.t(first_w), f"Mountain whitefish spawning ({DAYS[first_w * 3]:%b %d}, simulated) - cello pizzicato")
# ------------------------------------------------------------------ birds --
BIRDS = [ # (recording, day-of-year window, max prob per bar, gain) - typical timing, not observations
("american_dipper", (335, 366), 0.2, 0.7), ("american_dipper", (1, 75), 0.2, 0.7),
("swainsons_thrush_a", (140, 205), 0.35, 0.9), ("swainsons_thrush_b", (140, 205), 0.2, 0.9),
("bank_swallow", (125, 215), 0.18, 0.5),
]
for name, (a, z), p, gain in BIRDS:
phr = bird_phrases(name)
fire = W.fire_bars(per_bar(in_window(a, z)), p, phase=rng.uniform(0.3, 0.9))
for b in range(N_BARS):
if b in fire:
x = phr[rng.integers(len(phr))]
B["birds"].add(x * gain, g.s(b, rng.uniform(0, 3)), pan=rng.uniform(-0.7, 0.7))
mark(g.day_to_time(Y.doy_index(140)), "Swainson's thrush arrives (typical ~May 20)")
# ---------------------------------------------------------------- master --
snowcov = np.clip(smooth(sog, 5) / 25, 0, 1)
t_day = np.arange(n) / SR / g.day_to_time(1)
snow_env = np.interp(t_day, np.arange(N_DAYS), snowcov)
irs = {"open": make_ir(3.4, bright=0.65), "snow": make_ir(1.4, bright=0.1, seed=11)}
mix, stems = S.master([b for k, b in B.items() if k not in P["mute"]], irs, HEAVY / f"{NAME}.wav", stems_dir=HEAVY / f"{NAME}_stems",
ir_env={"open": 1 - snow_env, "snow": snow_env})
S.to_mp3(HEAVY / f"{NAME}.wav", OUT / f"{NAME}.mp3")
# ------------------------------------------------------ score + guide ----
for b in range(0, N_BARS, 8):
if b == 0 or modes[b] != modes[b - 8]:
mark(g.t(b), f"Mode -> {W.MODES[modes[b]][0]} ({DAYS[b * 3]:%b %d})")
for b in range(1, N_BARS):
if bar_phase[b] != bar_phase[b - 1]:
mark(g.t(b), f"Phase -> {bar_phase[b]} ({DAYS[b * 3]:%b %d})")
if first_green is not None:
mark(g.t(first_green), f"Green-up: string pad enters ({DAYS[first_green * 3]:%b %d}; {GPP_SOURCE.split(' (')[0]})")
events.sort()
level = {k: float(20 * np.log10(np.sqrt((v ** 2).mean()) + 1e-12)) for k, v in stems.items()}
score = {
"name": NAME, "bpm": BPM, "tonic": "D", "bar_seconds": g.bar, "days_per_bar": 3,
"duration_s": g.dur, "gpp_source": GPP_SOURCE,
"days": [{"date": str(d), "t": round(g.day_to_time(i), 2), "flow": round(float(q[i]), 2),
"flow_pct_normal": round(float(qpct[i]), 3), "burrell": round(float(bq[i]), 3),
"ice": bool(ice[i]), "temp": round(float(T[i]), 1), "rain": float(rain[i]), "snow": float(snow[i]),
"snow_on_ground": float(sog[i]), "daylength": round(float(DL[i]), 2), "gpp_norm": round(float(gppn[i]), 3),
"trout_sim": round(float(trout[i]), 3), "whitefish_sim": round(float(whitefish[i]), 3)}
for i, d in enumerate(DAYS)],
"bars": [{"bar": b, "t": round(g.t(b), 2), "mode": W.MODES[modes[b]][0], "degree": int(chords[b]),
"phase": bar_phase[b], "piano_density": int(np.sqrt(qn_bar[b]) * 3.99),
"perc_drive": round(float(drive_bar[b]), 3)} for b in range(N_BARS)],
"events": events, "stem_rms_db": level,
}
if META: # variants and remixes say what they changed; the published render has no such key
score["variant"] = META
with open(OUT / f"{NAME}_score.json", "w") as fh:
json.dump(score, fh, indent=1)
print("\n".join(f"{W.fmt_time(t)} {e}" for t, e in events))
print("stem RMS dBFS:", {k: round(v, 1) for k, v in level.items()})
if __name__ == "__main__":
main()
Note
pilot_okanagan_wy2024.py: the Okanagan score
"""
PILOT B - Okanagan, water year 2023-24 (Oct 1 2023 - Sep 30 2024).
Any other year the records cover: `pilot_okanagan_wy2024.py wy2019` (WY2016-WY2024,
bounded by the Wells Dam counts; the arrangement was tuned on WY2024).
Key of F, 64 bpm, 1 bar = 3 days, ~7.7 minutes. Slower tempo and HALF the
harmonic rhythm of the Granby pilot (a chord every 4 bars, not 2), because a
big regulated lake changes slowly.
A drought year on the lake (it rose only ~0.7 m and peaked ~1 m below 2017)
that was ALSO a very large sockeye return at Wells Dam (491k, the most in the
2015-24 counts). Low water and many fish in the same summer.
LAYER DATA (R = retrieved, D = derived, S = simulated) SOUND
lake drone Okanagan Lake level at Kelowna 08NM083 (R) F pedal (contrabass), cello 5th while filling
harmony / mode day length (D) + lake level vs 15-yr normal (R) mode brightness
progression lake phase: winter / filling / full / drawdown (D) which 4-chord loop plays
piano ostinato Okanagan River at Penticton 08NM050 = dam outflow (R) density follows the release; when the
outflow is held flat the pattern is played
strictly, unvaried
soft percussion Okanagan River near Oliver 08NM085 (R) log drum / conga / shaker
sockeye Wells Dam daily adult sockeye count, DART (R) marimba runs climbing the scale
thermal stress Wells Dam water temperature, DART (R) a held high violin note above 18 C
sockeye spawning Oct window near Oliver, scaled by run size (S) cello pizzicato, low and falling
string pad MODIS GPP, Summerland bench (R) violins swell with plant productivity
rain / snow Summerland CS precip + temperature split (R/D) harp drops / glockenspiel flakes
frost Summerland CS mean temp < -8 C (R) ice crackle (Jan 2024 arctic outbreak)
air / space Summerland CS min relative humidity (R) humid = long hall, dry = short bright room
birds typical phenology windows (S timing, R recordings) meadowlark (spring), sandhill cranes (passage)
"""
import json
import sys
from pathlib import Path
import numpy as np
sys.path.insert(0, str(Path(__file__).resolve().parent))
import sampler as S
from sampler import Bus, bird_phrases, make_ir, SR
import seasonal as W
from timegrid import per_bar, smooth, norm01
from dsp_core import ice_crinkle_grain
def main():
OUT = W.LIGHT_OUT / "pilots"
OUT.mkdir(parents=True, exist_ok=True)
HEAVY = W.HEAVY_OUT / "pilots" # WAV + stems (large) live in runtime
HEAVY.mkdir(parents=True, exist_ok=True)
Y = W.Year.parse(sys.argv[1] if len(sys.argv) > 1 else "wy2024")
DAYS, N_DAYS, N_BARS = Y.days, Y.n_days, Y.n_bars
y0, y1 = DAYS[0].year, DAYS[-1].year
NAME = f"okanagan_{Y.label.lower()}"
rng = np.random.default_rng(1997)
np.random.seed(1997) # the ice grain helper uses the global RNG; seed it so renders repeat exactly
TONIC = 53 # F3
BPM = 64
g = W.Grid(BPM, Y.n_bars, tail_s=12.0)
print(f"{N_BARS} bars, {g.dur / 60:.1f} min")
# --------------------------------------------------------------- features --
lake_all, _ = W.wsc_series("okanagan_lake_kelowna_08NM083_2010_2024.json", key="LEVEL")
pen_all, _ = W.wsc_series("okanagan_river_penticton_08NM050_2010_2024.json")
oli_all, _ = W.wsc_series("okanagan_river_oliver_08NM085_2010_2024.json")
clim = "summerland_cs_979_climate_2010_2024.json"
clim_T = W.climate_series(clim, "MEAN_TEMPERATURE")
W.check_coverage({"Okanagan Lake level": lake_all, "Penticton flow": pen_all,
"Oliver flow": oli_all, "Summerland temperature": clim_T}, Y)
lake = W.in_year(lake_all, year=Y)
lake_n = norm01(lake, 0, 100)
lake_pct = W.doy_percentile(lake_all, year=Y)
pen = W.in_year(pen_all, year=Y)
pen_n = norm01(pen, 0, 100)
oli = W.in_year(oli_all, year=Y)
oli_n = norm01(np.log(oli))
T = W.in_year(clim_T, year=Y)
precip = W.in_year(W.climate_series(clim, "TOTAL_PRECIPITATION"), fill="zero", year=Y)
rain = np.where(T > 1.0, precip, 0.0)
snow = np.where(T <= 1.0, precip, 0.0) # mm water equivalent
hum = W.in_year(W.climate_series(clim, "MIN_REL_HUMIDITY"), year=Y)
DL = Y.daylength(49.6)
gpp, gppn = W.modis_gpp_daily("summerland", Y)
GPP_SOURCE = "MODIS MOD17A2HGF (retrieved)"
if gpp is None:
gppn = np.clip(smooth(np.maximum(0, T - 5), 21) / 12, 0, 1)
GPP_SOURCE = "temperature proxy (derived; MODIS pull unavailable)"
sock, wells_T = W.dart_sockeye(Y)
doy = np.array([d.timetuple().tm_yday for d in DAYS])
# ---------------------------------------------------------------- harmony --
bright = per_bar(0.75 * norm01(DL, 0, 100) + 0.25 * norm01(smooth(T, 15)))
anom = per_bar(smooth(lake_pct, 15)) - 0.5
modes = W.mode_track(bright, anom, lo=1, hi=4, phrase=8)
ls_ = smooth(lake, 7)
dl_ = np.gradient(ls_) * 100 # cm/day
phase_day = np.where(dl_ > 0.4, "filling",
np.where(lake_n > 0.85, "full",
np.where(dl_ < -0.25, "drawdown", "winter")))
PROG = { # 4 chords x 4 bars each (16-bar loops)
"winter": [0, 5, 3, 5], # i VI iv VI - rocking, patient
"filling": [0, 2, 3, 4], # i III iv v - bass climbs with the water
"full": [5, 2, 6, 0], # VI III VII i - broad; the VII->i 'Aeolian cadence'
"drawdown": [0, 6, 5, 4], # i VII VI v - bass walks down
}
bar_phase = W.phase_bars(phase_day, min_run=8)
chords = W.chord_track(bar_phase, PROG, bars_per_chord=4)
K = S.vsco_kit()
piano, bass, cello, cpizz, vln = K["piano"], K["bass"], K["cello"], K["cpizz"], K["vln"]
harp, glock, marimba, perc = K["harp"], K["glock"], K["marimba"], K["perc"]
n = g.n
B = {k: Bus(k, n, **kw) for k, kw in {
"piano": dict(gain=1.0, hp=70, send=0.4),
"drone": dict(gain=0.8, lp=800, send=0.15),
"pad": dict(gain=0.19, hp=220, send=0.6),
"fish": dict(gain=0.95, hp=120, send=0.4),
"heat": dict(gain=0.22, hp=400, send=0.7),
"precip": dict(gain=1.05, hp=250, send=0.55),
"frost": dict(gain=1.5, hp=500, send=0.5),
"perc": dict(gain=1.8, hp=60, send=0.18),
"birds": dict(gain=0.26, hp=700, send=0.4),
}.items()}
events = []
def mark(sec, text):
events.append((round(sec, 1), text))
def tones(b, lo, hi, prev=None, ext=(0, 2, 4, 8)):
return W.voice_lead(prev, W.chord(chords[b], modes[b], TONIC, ext), lo, hi)
# ---------------------------------------------------------- lake drone ----
F_PEDAL = 29 # F1 - the lake is a pedal point: the same low note under every chord
lake_bar, dlake_bar = per_bar(lake_n), per_bar(dl_)
for b in range(0, N_BARS, 4):
vel = 0.3 + 0.5 * lake_bar[b]
B["drone"].add(W.sustain(K["bass"], F_PEDAL + 12, vel, g.bar * 4, rng), g.s(b))
if dlake_bar[b] > 0.3: # filling: the cello adds the 5th
B["drone"].add(W.sustain(cello, F_PEDAL + 19, vel * 0.8, g.bar * 4, rng), g.s(b), gain=0.5, pan=0.2)
# -------------------------------------------------------- piano / outflow --
PATTERN = [0, 2, 3, 2, 1, 2, 3, 1]
DENSITY = {0: [0, 4], 1: [0, 2, 4, 6], 2: list(range(8)), 3: list(range(8))}
pen_bar = per_bar(pen_n)
# 'held' = the dam kept outflow within +-3% over the surrounding two weeks
rel_var = np.array([np.ptp(pen[max(0, i - 7):i + 8]) / pen[max(0, i - 7):i + 8].mean() for i in range(N_DAYS)])
held = per_bar(rel_var < 0.06)
dl_bar = per_bar(norm01(DL, 0, 100))
prev = None
fixed_vel = None
for b in range(N_BARS):
v = tones(b, 57, 79, prev)
prev = v
lvl = int(np.sqrt(pen_bar[b]) * 3.99)
strict = held[b] > 0.5
base = 0.26 + 0.3 * dl_bar[b] + 0.15 * pen_bar[b]
for step in DENSITY[lvl]:
if strict: # the dam's hand: identical, mechanical repetition
vel, jit = base + (0.08 if step in (0, 4) else 0), 0.0
else:
vel, jit = base + (0.1 if step in (0, 4) else 0) + rng.normal(0, 0.04), 9.0
B["piano"].add(piano.note(v[PATTERN[step]], float(np.clip(vel, 0.1, 0.9)), dur=g.beat * 1.4, release=1.0, rng=rng),
g.s(b, step * 0.5, jitter_ms=jit, rng=rng), pan=-0.1 + 0.06 * (PATTERN[step] - 1.5))
if b % 4 == 0:
root = W.voice_lead(None, [W.degree_to_midi(chords[b], modes[b], TONIC)], 41, 53)[0]
B["piano"].add(piano.note(root, 0.4, dur=g.bar * 2, release=1.5, rng=rng), g.s(b), pan=-0.2)
# mark the regulated flow's biggest step changes
steps = [i for i in range(7, N_DAYS - 7) if abs(pen[i + 3:i + 7].mean() - pen[i - 6:i - 2].mean()) > 2.5]
last = -99
for i in steps:
if i - last > 20:
mark(g.day_to_time(i), f"Penticton dam release changes {pen[i - 4]:.0f} -> {pen[i + 5]:.0f} m3/s ({DAYS[i]:%b %d})")
last = i
first_held = next((b for b in range(N_BARS) if held[b] > 0.5), None)
# ---------------------------------------------------------- string pad ----
gpp_bar = per_bar(gppn)
prev = None
for b in range(0, N_BARS, 2):
if gpp_bar[b] < 0.08:
prev = None
continue
v = tones(b, 65, 88, prev, ext=(2, 4, 8))
prev = v
for i, m in enumerate(v):
B["pad"].add(W.sustain(vln, m, 0.2 + 0.6 * gpp_bar[b], g.bar * 2, rng), g.s(b), gain=gpp_bar[b] ** 1.2,
pan=(-0.5, 0.1, 0.5)[i])
first_green = next((b for b in range(N_BARS) if gpp_bar[b] >= 0.08 and b > 40), None)
# -------------------------------------------------------------- sockeye ----
sock_bar = per_bar(sock, "sum")
cum = np.cumsum(sock_bar) / max(sock_bar.sum(), 1)
for b in range(N_BARS):
c = sock_bar[b]
if c < 5:
continue
k = int(min(12, round(1.7 * np.log1p(c / 40))))
if k == 0:
continue
# climb through the chord's own tones (never a passing note that could rub),
# starting higher as the run progresses: the season's fish moving upstream
arp = sorted({m + 12 * o for m in W.chord(chords[b], modes[b], TONIC) for o in range(-1, 3)})
arp = [m for m in arp if 60 <= m <= 96]
start = int(cum[b] * (len(arp) - 6))
slots = np.sort(rng.choice(16, size=min(k, 16), replace=False))
for j, sl in enumerate(slots):
m = arp[min(start + j, len(arp) - 1)]
B["fish"].add(marimba.note(m, float(np.clip(0.3 + 0.05 * k + rng.normal(0, 0.05), 0.2, 0.9)), rng=rng),
g.s(b, sl * 0.25, jitter_ms=4, rng=rng), pan=float(np.clip(-0.6 + 0.1 * j, -0.7, 0.7)))
pk = int(np.argmax(sock))
fs = next((i for i in range(N_DAYS) if sock[i] > 0 and DAYS[i].year == y1), None)
if fs is not None:
mark(g.day_to_time(fs), f"First sockeye of the {y1} run pass Wells Dam ({DAYS[fs]:%b %d}) - marimba")
mark(g.day_to_time(pk), f"Sockeye peak at Wells: {sock[pk]:,.0f} fish in one day ({DAYS[pk]:%b %d})")
# thermal stress: Wells water temperature (a real measurement) above 18 C
wt = np.nan_to_num(wells_T, nan=0.0)
stress = np.clip((smooth(wt, 5) - 18.0) / 2.5, 0, 1)
st_bar = per_bar(stress)
for b in range(0, N_BARS, 2):
if st_bar[b] <= 0.02:
continue
m = W.voice_lead(None, [W.degree_to_midi(8, modes[b], TONIC)], 84, 96)[0] # the 9th, held high
B["heat"].add(W.sustain(vln, m, 0.3 + 0.5 * st_bar[b], g.bar * 2, rng), g.s(b), gain=st_bar[b])
hs = next((i for i in range(N_DAYS) if stress[i] > 0.05), None)
if hs is not None:
mark(g.day_to_time(hs), f"Wells water passes 18 C ({DAYS[hs]:%b %d}) - held high violin: thermal stress")
wmax = int(np.nanargmax(wells_T))
mark(g.day_to_time(wmax), f"Warmest water at Wells, {wells_T[wmax]:.1f} C ({DAYS[wmax]:%b %d})")
# spawning near Oliver (simulated timing, scaled by the run that produced it)
# WY2024 keeps the counts it was rendered with on 09-24 (the cache now sums 2024 to 491,040, one more)
runs_size = {2023: 136956, 2024: 491039} if Y.label == "WY2024" else W.dart_run_totals(Y)
spawn = sum(runs_size[y] / runs_size[y1] * np.exp(-0.5 * (np.array([(d - W.dt.date(y, 10, 15)).days for d in DAYS]) / 9) ** 2)
for y in runs_size)
sp_bar = per_bar(spawn)
sp_fire = W.fire_bars(sp_bar, 2.2)
for b in range(N_BARS):
if b in sp_fire:
v = sorted(tones(b, 41, 57, ext=(0, 2, 4)), reverse=True)
for i, m in enumerate(v):
B["fish"].add(cpizz.note(m, 0.5, rng=rng), g.s(b, 1 + i * 0.75, jitter_ms=5, rng=rng), gain=0.8, pan=-0.35)
mark(0.0, f"{DAYS[0]:%b %Y}: the {y0} sockeye spawning near Oliver (timing simulated) - low cello pizzicato")
# ------------------------------------------------------- soft percussion --
active = np.clip((smooth(oli_n, 7) - 0.3) / 0.25, 0, 1)
drive = np.sqrt(np.clip(oli_n, 0, 1)) * active
drive_bar = per_bar(drive)
SW = 0.1
for b in range(N_BARS):
lv = drive_bar[b]
if lv < 0.12:
continue
for step in range(8):
beat = step * 0.5 + (SW * 0.5 if step % 2 else 0)
at = lambda: g.s(b, beat, jitter_ms=6, rng=rng)
if step in (0, 6) and lv > 0.12:
B["perc"].add(perc.hit("LogDrumLo", 0.35 + 0.35 * lv, rng), at(), pan=-0.3)
if lv > 0.35 and step in (3, 4):
B["perc"].add(perc.hit("LogDrumHi", 0.25 + 0.3 * lv, rng), at(), pan=0.3)
if lv > 0.5 and step % 2:
B["perc"].add(perc.hit("Tamb1-Shake", 0.12 + 0.2 * lv, rng), at(), pan=0.45)
if lv > 0.7 and step in (2, 7):
B["perc"].add(perc.hit("Conga-Tap1", 0.35, rng), at(), pan=0.15)
po = int(np.argmax(oli))
mark(g.day_to_time(po), f"Okanagan River at Oliver peaks, {oli[po]:.0f} m3/s ({DAYS[po]:%b %d}) - fullest percussion")
# ------------------------------------------------------ precip + frost ----
for d in range(N_DAYS):
b = min(d // 3, N_BARS - 1)
t0 = g.day_to_time(d)
if rain[d] > 0.2:
pool = tones(b, 67, 91)
for _ in range(min(9, int(round(1.6 * np.sqrt(rain[d]))))):
B["precip"].add(harp.note(int(rng.choice(pool)), rng.uniform(0.25, 0.55), rng=rng),
(t0 + rng.uniform(0, g.day_to_time(1))) * SR, pan=rng.uniform(-0.7, 0.7))
if snow[d] > 0.2:
pool = tones(b, 81, 98)
for _ in range(min(7, int(round(1.8 * np.sqrt(snow[d]))))):
B["precip"].add(glock.note(int(rng.choice(pool)), 0.22, dur=1.0, release=2.5, rng=rng),
(t0 + rng.uniform(0, g.day_to_time(1))) * SR, gain=0.6, pan=rng.uniform(-0.8, 0.8))
if T[d] < -8:
cold = np.clip((-T[d] - 8) / 12, 0.15, 1)
for _ in range(int(4 + 12 * cold)):
grain = ice_crinkle_grain(SR, dur=0.03, intensity=cold)[:, None].repeat(2, 1)
B["frost"].add(grain * 0.25, (t0 + rng.uniform(0, g.day_to_time(1))) * SR, pan=rng.uniform(-0.8, 0.8))
cold_d = int(np.argmin(T))
mark(g.day_to_time(cold_d), f"Coldest day, {T[cold_d]:.1f} C ({DAYS[cold_d]:%b %d}) - frost crackle")
# ------------------------------------------------------------------ birds --
in_window = lambda a, z: ((doy >= a) & (doy <= z)).astype(float)
BIRDS = [("western_meadowlark_a", (70, 197), 0.3, 0.9), ("western_meadowlark_b", (70, 197), 0.2, 0.9),
("sandhill_crane", (92, 116), 0.45, 0.8), ("sandhill_crane", (254, 285), 0.45, 0.8)]
for name, (a, z), p, gain in BIRDS:
phr = bird_phrases(name)
fire = W.fire_bars(per_bar(in_window(a, z)), p, phase=rng.uniform(0.3, 0.9))
for b in range(N_BARS):
if b in fire:
B["birds"].add(phr[rng.integers(len(phr))] * gain, g.s(b, rng.uniform(0, 3)), pan=rng.uniform(-0.7, 0.7))
mark(0.0, "Early Oct: sandhill cranes passing over (typical timing)")
mark(g.day_to_time(Y.doy_index(71)), "Western meadowlark returns (typical ~Mar 10)")
# ---------------------------------------------------------------- master --
humid = np.clip((smooth(hum, 7) - 20) / 50, 0, 1)
t_day = np.arange(n) / SR / g.day_to_time(1)
h_env = np.interp(t_day, np.arange(N_DAYS), humid)
irs = {"open": make_ir(4.0, bright=0.5, seed=3), "dry": make_ir(1.2, bright=0.8, seed=5)}
mix, stems = S.master(list(B.values()), irs, HEAVY / f"{NAME}.wav", stems_dir=HEAVY / f"{NAME}_stems",
ir_env={"open": h_env, "dry": 1 - h_env})
S.to_mp3(HEAVY / f"{NAME}.wav", OUT / f"{NAME}.mp3")
for b in range(0, N_BARS, 8):
if b == 0 or modes[b] != modes[b - 8]:
mark(g.t(b), f"Mode -> {W.MODES[modes[b]][0]} ({DAYS[b * 3]:%b %d})")
for b in range(1, N_BARS):
if bar_phase[b] != bar_phase[b - 1]:
mark(g.t(b), f"Lake phase -> {bar_phase[b]} ({DAYS[b * 3]:%b %d})")
if first_green is not None:
mark(g.t(first_green), f"Green-up: string pad enters ({DAYS[first_green * 3]:%b %d}; {GPP_SOURCE.split(' (')[0]})")
if first_held is not None:
mark(g.t(first_held), f"Outflow held flat by the dam ({DAYS[first_held * 3]:%b %d}) - piano turns strict")
lk = int(np.argmax(lake))
mark(g.day_to_time(lk), f"Lake peaks at gauge {lake[lk]:.2f} m ({DAYS[lk]:%b %d}); 2017 reached 3.01")
events.sort()
level = {k: float(20 * np.log10(np.sqrt((v ** 2).mean()) + 1e-12)) for k, v in stems.items()}
score = {
"name": NAME, "bpm": BPM, "tonic": "F", "bar_seconds": g.bar, "days_per_bar": 3,
"duration_s": g.dur, "gpp_source": GPP_SOURCE,
"days": [{"date": str(d), "t": round(g.day_to_time(i), 2), "lake_level": round(float(lake[i]), 3),
"lake_pct_normal": round(float(lake_pct[i]), 3), "penticton": round(float(pen[i]), 2),
"oliver": round(float(oli[i]), 2), "sockeye": float(sock[i]),
"wells_temp": None if np.isnan(wells_T[i]) else float(wells_T[i]),
"temp": round(float(T[i]), 1), "rain": float(rain[i]), "snow_we": float(snow[i]),
"humidity_min": round(float(hum[i]), 0), "daylength": round(float(DL[i]), 2),
"gpp_norm": round(float(gppn[i]), 3), "spawn_sim": round(float(spawn[i]), 3)}
for i, d in enumerate(DAYS)],
"bars": [{"bar": b, "t": round(g.t(b), 2), "mode": W.MODES[modes[b]][0], "degree": int(chords[b]),
"phase": bar_phase[b], "piano_density": int(np.sqrt(pen_bar[b]) * 3.99),
"held": bool(held[b] > 0.5), "perc_drive": round(float(drive_bar[b]), 3)} for b in range(N_BARS)],
"events": events, "stem_rms_db": level,
}
with open(OUT / f"{NAME}_score.json", "w") as fh:
json.dump(score, fh, indent=1)
print("\n".join(f"{W.fmt_time(t)} {e}" for t, e in events))
print("stem RMS dBFS:", {k: round(v, 1) for k, v in level.items()})
if __name__ == "__main__":
main()
Note
export_companion.py: data for the listening page
"""
Export a pilot's score + stem envelopes + audio into listen/ for the
listening-companion page (listen/index.html).
python export_companion.py # both pilots
Writes listen/<name>.json and listen/<name>.mp3 (160 kbps, encoded from the master WAV).
"""
import json, subprocess, sys
from pathlib import Path
import numpy as np
import soundfile as sf
sys.path.insert(0, str(Path(__file__).resolve().parent))
from dsp_core import ROOT, HEAVY_OUT
P = ROOT / "rendered" / "pilots"
HEAVY = HEAVY_OUT / "pilots" # master WAV + stems
OUT = ROOT / "listen"
ENV_HZ = 8
VOICES = { # stem -> (label, colour key)
"piano": ("Piano ostinato", "river"), "bass": ("Contrabass + cello", "deep"),
"drone": ("Lake drone", "deep"), "pad": ("Violins", "moss"), "ice": ("Ice", "ice"),
"frost": ("Frost", "ice"), "precip": ("Rain + snow", "rain"), "perc": ("Hand drums", "clay"),
"fish": ("Fish", "salmon"), "heat": ("Heat violin", "sun"), "birds": ("Birds", "sun"),
}
LANES = {
"granby_wy2024": [
# key, label, colour, transform, source, stem, what the sound does
("flow", "Granby River flow", "river", "log", "R", "piano", "Piano gets denser as flow rises; contrabass holds the root, cello adds a 5th at high flow."),
("flow_pct_normal", "Flow vs. 15-yr normal", "river2", "lin", "D", "piano", "Sets the harmony's brightness: below normal leans one mode darker, above normal one brighter."),
("burrell", "Burrell Creek flow", "clay", "log", "R", "perc", "Hand drums: tumba, conga, shaker, and a 3-3-2 log-drum figure at the peak."),
("ice", "Ice conditions (WSC flag)", "ice", "flag", "R", "ice", "Freeze: falling glockenspiel + bowed cymbal. Frozen: crackle and descending singing-ice chirps. Breakup: cymbal swell and rising harp."),
("temp", "Air temperature, Billings", "sun", "temp", "R", "ice", "Colder days make the ice crackle denser; also brightens the harmony across the season."),
("rain", "Rain", "rain", "sqrt", "R", "precip", "Each rainy day scatters harp notes from the current chord."),
("snow", "Snowfall", "snow", "sqrt", "R", "precip", "Snowfall scatters soft glockenspiel notes, high up."),
("snow_on_ground", "Snow on ground", "snow", "lin", "R", None, "Deep snow swaps the reverb to a dark, short room: snow absorbs sound."),
("gpp_norm", "Plant productivity (MODIS GPP)", "moss", "lin", "R", "pad", "Violins swell with photosynthesis."),
("trout_sim", "Redband trout spawning", "salmon", "lin", "S", "fish", "A short rising flute figure. Timing is modelled from temperature and season."),
("whitefish_sim", "Mountain whitefish spawning", "salmon2", "lin", "S", "fish", "A falling cello pizzicato figure. Timing is modelled."),
("daylength", "Day length", "sun2", "lin", "D", None, "With temperature, sets how bright the mode is."),
],
"okanagan_wy2024": [
("lake_level", "Okanagan Lake level, Kelowna", "deep", "lin", "R", "drone", "A low F held under everything, louder as the lake fills; cello adds a 5th while it rises."),
("lake_pct_normal", "Lake vs. 15-yr normal", "river2", "lin", "D", "drone", "Sets the harmony's brightness: below normal leans one mode darker."),
("penticton", "Outflow at Penticton dam", "river", "lin", "R", "piano", "Piano density follows the release. When the dam holds the flow flat, the piano plays strictly, with no variation."),
("oliver", "Okanagan River near Oliver", "clay", "log", "R", "perc", "Soft log drums, conga and shaker."),
("sockeye", "Sockeye passing Wells Dam", "salmon", "log", "R", "fish", "Marimba runs climb the chord; more fish means more notes, and the runs start higher as the season goes."),
("wells_temp", "Water temperature at Wells", "sun", "lin", "R", "heat", "Above 18 C a violin holds a high note: thermal stress for migrating sockeye."),
("spawn_sim", "Sockeye spawning near Oliver", "salmon2", "lin", "S", "fish", "Low, falling cello pizzicato around mid-October, scaled by run size. Timing is modelled."),
("temp", "Air temperature, Summerland", "sun2", "temp", "R", "frost", "Below -8 C: frost crackle."),
("rain", "Rain", "rain", "sqrt", "D", "precip", "Harp notes from the current chord. Rain vs. snow is split by temperature."),
("snow_we", "Snow", "snow", "sqrt", "D", "precip", "Soft, high glockenspiel."),
("humidity_min", "Humidity (daily min)", "ice", "lin", "R", None, "Humid days get a long hall reverb; dry days a short, bright room."),
("gpp_norm", "Plant productivity (MODIS GPP)", "moss", "lin", "R", "pad", "Violins swell with photosynthesis."),
],
}
# How the dashboard reads each record's value on a day: (unit, decimals). Special units:
# "pct" a percentile vs. the 2010-24 record, "flag" yes/no, "peak" a share of a
# typical peak day (MODIS 95th percentile, 2015-24), "sim" a modelled 0-1 window.
UNITS = {
"flow": ("m³/s", 1), "flow_pct_normal": ("pct", 0), "burrell": ("m³/s", 2), "ice": ("flag", 0),
"temp": ("°C", 1), "rain": ("mm", 1), "snow": ("cm", 1), "snow_on_ground": ("cm", 0),
"gpp_norm": ("peak", 2), "trout_sim": ("sim", 2), "whitefish_sim": ("sim", 2), "daylength": ("h", 1),
"lake_level": ("m", 3), "lake_pct_normal": ("pct", 0), "penticton": ("m³/s", 1), "oliver": ("m³/s", 1),
"sockeye": ("fish", 0), "wells_temp": ("°C", 1), "spawn_sim": ("sim", 2), "snow_we": ("mm", 1),
"humidity_min": ("%", 0),
}
TITLES = {
"granby_wy2024": dict(title="Granby River", sub="Grand Forks, BC · WSC 08NN002 + Burrell Creek 08NN023 · Billings climate",
key="D", year_label="WY 2024", blurb="A warm El Niño winter. The gauge flagged ice on only 28 days, in broken episodes, and the freshet came early and small: 147 m³/s, against 465 in the 2010–24 record."),
"okanagan_wy2024": dict(title="Okanagan", sub="Okanagan Lake 08NM083 · Penticton 08NM050 · Oliver 08NM085 · Summerland CS · Wells Dam",
key="F", year_label="WY 2024", blurb="A drought year on the lake, which rose only 0.7 m, that was also the biggest sockeye return in ten years of Wells Dam counts."),
}
def transform(v, how):
v = np.array([np.nan if x is None else float(x) for x in v])
if how == "flag":
return np.nan_to_num(v).clip(0, 1)
if how == "log":
v = np.log1p(np.nan_to_num(v) / max(np.nanpercentile(v, 5), 1e-3))
elif how == "sqrt":
v = np.sqrt(np.nan_to_num(v))
elif how == "temp":
lo, hi = np.nanmin(v), np.nanmax(v)
return np.round(np.nan_to_num((v - lo) / (hi - lo), nan=0), 3), float(-lo / (hi - lo))
hi, lo = np.nanpercentile(v, 99), np.nanmin(v)
return np.round(np.clip(np.nan_to_num((v - lo) / (hi - lo + 1e-12), nan=0), 0, 1), 3)
# Records whose daily values this site does not republish: the page plays them and draws their
# normalized shape (as the 09-25 release did), but its dashboard shows no numbers for them.
WITHHELD = {
"sockeye": "Columbia River DART counts are not republished here; you hear them, and the ring draws their shape.",
"wells_temp": "Columbia River DART temperatures are not republished here; you hear them, and the ring draws their shape.",
}
def main():
OUT.mkdir(exist_ok=True)
for name, lanes in LANES.items():
with open(P / f"{name}_score.json") as fh:
sc = json.load(fh)
days = sc["days"]
out_lanes = []
for key, label, col, how, src, stem, does in lanes:
r = transform([d[key] for d in days], how)
if isinstance(r, tuple): # "temp" also returns where 0 C sits; the page doesn't draw it
r = r[0]
unit, dp = UNITS[key]
raw = [None if d[key] is None else (float(d[key]) if unit != "flag" else int(bool(d[key])))
for d in days]
lane = dict(key=key, label=label, color=col, source=src, stem=stem, does=does,
v=[float(x) for x in r], raw=raw, unit=unit, dp=dp)
if key in WITHHELD: # heard, and drawn as a shape, but no values published
lane.update(raw=None, withheld=WITHHELD[key])
out_lanes.append(lane)
env = {}
for f in sorted((HEAVY / f"{name}_stems").glob("*.flac")):
x, sr = sf.read(f)
m = x.mean(axis=1) * 4
blk = sr // ENV_HZ
r = np.sqrt(np.array([(m[i:i + blk] ** 2).mean() for i in range(0, len(m) - blk, blk)]) + 1e-12)
db = 20 * np.log10(r)
env[f.stem] = [round(float(x), 2) for x in np.clip((db + 62) / 40, 0, 1)]
voices = [dict(stem=k, label=VOICES[k][0], color=VOICES[k][1]) for k in env]
with open(OUT / f"{name}.json", "w") as fh:
json.dump(dict(**TITLES[name], name=name, duration=sc["duration_s"], bar_seconds=sc["bar_seconds"],
n_days=len(days), dates=[d["date"] for d in days], t=[d["t"] for d in days], lanes=out_lanes,
bars=sc["bars"], events=sc["events"], env_hz=ENV_HZ, env=env, voices=voices,
gpp_source=sc["gpp_source"]),
fh, separators=(",", ":"))
subprocess.run(["ffmpeg", "-loglevel", "error", "-y", "-i", str(HEAVY / f"{name}.wav"), "-b:a", "160k",
str(OUT / f"{name}.mp3")], check=True)
print(name, (OUT / f"{name}.json").stat().st_size // 1024, "KB json,",
(OUT / f"{name}.mp3").stat().st_size // 1024 // 1024, "MB mp3")
if __name__ == "__main__":
main()
Note
dsp_core.py: paths, day length, ice grains
"""Core helpers shared by the water-year pieces: where files go, the sample
rate, astronomical day length, and the ice-crackle grain. Design notes live on
the methods page (methods.qmd).
Paths default to folders inside the repo:
rendered/ WAVs and stems (git-ignored) override: SONIFICATION_HEAVY
rendered/ MP3s and score JSON override: SONIFICATION_LIGHT
audio/ the site's copies of some tracks override: SONIFICATION_SITE_AUDIO
samples/ downloaded instrument audio override: SONIFICATION_SAMPLES
Set the environment variables to keep large audio somewhere else; renders.py
sets the first three to a scratch folder to re-render without overwriting."""
import os
from pathlib import Path
import numpy as np
from scipy.signal import lfilter
SR = 44100
ROOT = Path(__file__).resolve().parent.parent
HEAVY_OUT = Path(os.environ.get("SONIFICATION_HEAVY", ROOT / "rendered"))
SAMPLES = Path(os.environ.get("SONIFICATION_SAMPLES", ROOT / "samples"))
LIGHT_OUT = Path(os.environ.get("SONIFICATION_LIGHT", ROOT / "rendered"))
SITE_AUDIO = Path(os.environ.get("SONIFICATION_SITE_AUDIO", ROOT / "audio"))
def ice_crinkle_grain(sr=SR, dur=0.035, intensity=0.5):
"""One grain of an ice-forming 'crinkle' texture - brittle, glassy,
higher and more clicky than the shaker grain, built from several tiny
sub-clicks rather than one smooth burst so it reads as crackling
rather than hiss."""
n = int(dur * sr)
sig = np.zeros(n)
n_clicks = 3 + int(4 * intensity)
for _ in range(n_clicks):
pos = np.random.randint(0, max(n - 20, 1))
click_n = min(np.random.randint(6, 20), n - pos)
noise = np.random.uniform(-1, 1, click_n)
a = np.exp(-2 * np.pi * 6500 / sr)
hp = lfilter([1, -1], [1, -a], noise)
env = np.exp(-np.arange(click_n) / (click_n * 0.25))
sig[pos:pos + click_n] += hp * env
return sig / (np.abs(sig).max() + 1e-9)
def day_length_hours(day_of_year, latitude_deg=49.03):
"""Real astronomical day length (hours) from solar declination - a
calendar-derived seasonal signal that needs no data fetch and, unlike
snowpack, peaks at the summer solstice rather than lagging into
early summer, so it moves out of phase with the flow-driven layers."""
decl = np.radians(23.44) * np.sin(2 * np.pi * (284 + day_of_year) / 365.0)
lat = np.radians(latitude_deg)
cos_h = -np.tan(lat) * np.tan(decl)
cos_h = np.clip(cos_h, -1, 1)
hour_angle = np.arccos(cos_h)
return (24.0 / np.pi) * hour_angle
Note
fetch_granby.py, fetch_okanagan.py and fetch_samples.py: data and sample downloads
"""
Fetch the Granby pilot's source data into data/ (cached JSON; skipped if the
file is already there). Uses the same ECCC OGC API helpers as fetch_okanagan.py.
Hydrometric (ECCC OGC API, hydrometric-daily-mean, 2010-2024):
08NN002 Granby River at Grand Forks - DISCHARGE + ice-condition flags
08NN023 Burrell Creek above Gloucester Creek - DISCHARGE (a flashy tributary)
Climate (ECCC OGC API, climate-daily, 2010-2024):
STN_ID 1100 Billings (Grand Forks area)
The Granby valley MODIS GPP pull is in fetch_okanagan.py ("gpp").
"""
import sys
from pathlib import Path
sys.path.insert(0, str(Path(__file__).resolve().parent))
from fetch_okanagan import hydro, climate
if __name__ == "__main__":
hydro("08NN002", "granby")
hydro("08NN023", "burrell")
climate(1100, "billings")"""
Fetch the Okanagan pilot's source data into data/ (cached JSON, same shape as
the Granby pulls).
Hydrometric (ECCC OGC API, hydrometric-daily-mean, 2010-2024):
08NM083 Okanagan Lake at Kelowna - lake LEVEL (m), the slow regulated signal
08NM050 Okanagan River at Penticton - lake outflow DISCHARGE (dam-regulated)
08NM085 Okanagan River near Oliver - lower-river DISCHARGE (sockeye reach)
Climate (ECCC OGC API, climate-daily, 2010-2024):
STN_ID 979 Summerland CS (lakeside, central valley)
Satellite (ORNL DAAC MODIS subset REST API, 2015-2024, 8-day):
MOD17A2HGF GPP at a Summerland-bench point and the Granby valley point -
"plant metabolism" retrieved rather than simulated.
Fish (Columbia Basin Research DART, 2015-2024):
Wells Dam daily adult sockeye counts + water temperature. Wells is on the
Columbia above the Wenatchee and below the Okanagan confluence, so its
sockeye are overwhelmingly Okanagan-bound; it is not an Okanagan gauge.
The raw DART file is not committed to the repo; run `fetch_okanagan.py dart`
to pull it. Citation: Columbia River DART, Columbia Basin Research,
University of Washington. Adult Passage Daily Counts.
https://www.cbr.washington.edu/dart/query/adult_daily
"""
import json, sys, time, urllib.request, urllib.parse
from pathlib import Path
DATA = Path(__file__).resolve().parent.parent / "data"
OGC = "https://api.weather.gc.ca/collections"
def get(url, tries=4):
for i in range(tries):
try:
# ORNL returns HTTP 500 to urllib's default User-Agent; curl's works
req = urllib.request.Request(url, headers={"User-Agent": "curl/8.5", "Accept": "application/json"})
with urllib.request.urlopen(req, timeout=120) as r:
return r.read()
except Exception as e:
print(" retry", i, e, file=sys.stderr)
time.sleep(3 * (i + 1))
raise RuntimeError(url)
def ogc_all(coll, params):
out, offset = [], 0
while True:
q = dict(params, f="json", limit=10000, offset=offset)
d = json.loads(get(f"{OGC}/{coll}/items?" + urllib.parse.urlencode(q)))
out += d["features"]
if len(d["features"]) < 10000:
return out
offset += 10000
def hydro(stn, name):
p = DATA / f"{name}_{stn}_2010_2024.json"
if p.exists():
return
f = ogc_all("hydrometric-daily-mean",
{"STATION_NUMBER": stn, "datetime": "2010-01-01/2024-12-31"})
with open(p, "w") as fh:
json.dump(f, fh)
print(stn, len(f))
def climate(stn_id, name):
p = DATA / f"{name}_{stn_id}_climate_2010_2024.json"
if p.exists():
return
f = ogc_all("climate-daily", {"STN_ID": stn_id, "datetime": "2010-01-01/2024-12-31"})
with open(p, "w") as fh:
json.dump(f, fh)
print(stn_id, len(f))
def modis_gpp(lat, lon, name):
p = DATA / f"modis_gpp_{name}_2015_2024.json"
if p.exists():
return
base = "https://modis.ornl.gov/rst/api/v1/MOD17A2HGF"
dates = json.loads(get(f"{base}/dates?latitude={lat}&longitude={lon}"))["dates"]
dates = [d for d in dates if "2015-01-01" <= d["calendar_date"] <= "2024-12-31"]
rows = []
for i in range(0, len(dates), 10):
chunk = dates[i:i + 10]
q = urllib.parse.urlencode(dict(latitude=lat, longitude=lon,
startDate=chunk[0]["modis_date"], endDate=chunk[-1]["modis_date"],
kmAboveBelow=1, kmLeftRight=1))
d = json.loads(get(f"{base}/subset?{q}"))
for s in d["subset"]:
if s["band"] == "Gpp_500m":
rows.append({"date": s["calendar_date"], "data": s["data"]})
with open(p, "w") as fh:
json.dump({"lat": lat, "lon": lon, "product": "MOD17A2HGF", "band": "Gpp_500m",
"units": "kg C m-2 per 8 days, scale 0.0001; fill >= 32761",
"rows": rows}, fh)
print("gpp", name, len(rows))
def dart_sockeye():
p = DATA / "wells_dam_sockeye_2015_2024.json"
if p.exists():
return
out = {}
for yr in range(2015, 2025):
q = urllib.parse.urlencode({"sc": 1, "outputFormat": "csv", "year": yr, "proj": "WEL",
"span": "no", "startdate": "1/1", "enddate": "12/31"})
txt = get("https://www.cbr.washington.edu/dart/cs/php/rpt/adult_daily.php?" + q).decode()
out[yr] = txt
with open(p, "w") as fh:
json.dump(out, fh)
print("dart", {k: len(v) for k, v in out.items()})
if __name__ == "__main__":
what = sys.argv[1:] or ["hydro", "climate", "gpp", "dart"]
if "hydro" in what:
hydro("08NM083", "okanagan_lake_kelowna")
hydro("08NM050", "okanagan_river_penticton")
hydro("08NM085", "okanagan_river_oliver")
if "climate" in what:
climate(979, "summerland_cs")
if "gpp" in what:
modis_gpp(49.60, -119.68, "summerland")
modis_gpp(49.10, -118.45, "granby_valley")
if "dart" in what:
dart_sockeye()"""
Build/refresh the sample library and its license manifest.
Audio goes to samples/ in the repo (git-ignored; about 500 MB), or to
$SONIFICATION_SAMPLES if set. The manifest is committed so provenance is
browsable:
samples_manifest.csv (one row per source SET, not per file)
Sources:
- VSCO 2 Community Edition (CC0, github.com/sgossner/VSCO-2-CE) - pulled by
the scratch download list in the 2026-09-24 session; the instrument
subset is recorded below.
- Wikimedia Commons bird recordings (mostly xeno-canto mirrors, CC BY-SA),
fetched by exact title. Attribution is kept per file.
"""
import csv, json, time, urllib.request, urllib.parse
from pathlib import Path
from dsp_core import SAMPLES as ROOT
MANIFEST = Path(__file__).resolve().parent.parent / "samples_manifest.csv"
UA = {"User-Agent": "sonify-viz/0.1 (https://github.com/cordilleran/sonify-viz; personal research)"}
VSCO_SETS = [
("Keys/Upright Piano", "Upright piano (Ivy Audio / Simon Dalzell), 3 dynamics, sampled every major 3rd"),
("Strings/Harp", "Harp, mf"),
("Strings/Cello Section/susvib", "Cello section sustain w/ vibrato"),
("Strings/Cello Section/pizzT", "Cello section pizzicato"),
("Strings/Violin Section/susVib", "Violin section sustain w/ vibrato"),
("Strings/Violin Section/Pizz", "Violin section pizzicato"),
("Strings/Violin Section/Trem", "Violin section tremolo"),
("Strings/Solo Contrabass/SusNV", "Solo contrabass sustain, no vibrato"),
("Strings/Solo Contrabass/Pizz", "Solo contrabass pizzicato"),
("Woodwinds/Flute/susNV", "Flute sustain, no vibrato"),
("Percussion/Marimba", "Marimba"),
("Percussion/Glock", "Glockenspiel"),
("Percussion/Timpani", "Timpani hits and rolls (added 2026-09-25 for the climate pair's storm layer)"),
("Percussion/BDrumNewhit_*", "Concert bass drum, 7 dynamics (added 2026-09-25)"),
("Percussion/gongHit_*", "Tam-tam / gong hits p..fff (added 2026-09-25)"),
("Percussion/cymbal-crash1_* + susCymb1-hit_*", "Crash and suspended-cymbal hits (added 2026-09-25)"),
("Keys/Organ/Quiet", "Pipe organ, quiet manual, for drones (added 2026-09-25)"),
("Percussion/*", "Hand percussion + cymbal/gong/triangle one-shots (conga, tumba, quinto, log drum, tambourine shake, bowed/crescendo cymbal, gong scrape, bass-drum rub, bell tree)"),
]
BIRDS = [
# (commons title, local name, species, why it's here)
("File:Catharus ustulatus - Swainson's Thrush XC250599.mp3", "swainsons_thrush_a.mp3",
"Swainson's Thrush", "Granby riparian; late-May arrival, spiralling evening song"),
("File:Catharus ustulatus - Swainson's Thrush XC250601.mp3", "swainsons_thrush_b.mp3",
"Swainson's Thrush", "second take for variety"),
("File:Yellowstone sound library - American Dipper & Canada Geese - 001.mp3", "american_dipper.mp3",
"American Dipper", "the river's winter singer; sings through ice season"),
("File:Sturnella neglecta - Western Meadowlark - XC104522.ogg", "western_meadowlark_a.ogg",
"Western Meadowlark", "Okanagan grassland benches; March arrival"),
("File:Sturnella neglecta - Western Meadowlark XC254454.mp3", "western_meadowlark_b.mp3",
"Western Meadowlark", "second take"),
("File:Grus canadensis - Sandhill Crane - XC111491.ogg", "sandhill_crane.ogg",
"Sandhill Crane", "Okanagan flyway passage, April and September"),
("File:Riparia riparia - Sand Martin XC487832.mp3", "bank_swallow.mp3",
"Bank Swallow (Sand Martin, European recording, same species)", "Grand Forks cutbank colonies, May-Aug"),
("File:Gavia immer - Common Loon XC139388.mp3", "common_loon.mp3",
"Common Loon", "Lake Superior spring: loons return as the ice leaves (icecover piece)"),
("File:Branta canadensis - Canada Goose XC127624.ogg", "canada_goose.ogg",
"Canada Goose", "early-returning geese on still-frozen lake ice (icecover piece)"),
]
def commons(title):
q = urllib.parse.urlencode({"action": "query", "titles": title, "prop": "imageinfo",
"iiprop": "url|extmetadata", "format": "json"})
req = urllib.request.Request("https://commons.wikimedia.org/w/api.php?" + q, headers=UA)
page = next(iter(json.load(urllib.request.urlopen(req, timeout=60))["query"]["pages"].values()))
ii = page["imageinfo"][0]
m = ii["extmetadata"]
strip = lambda s: __import__("re").sub("<[^>]+>", "", s or "").strip()
return ii["url"], strip(m.get("LicenseShortName", {}).get("value")), strip(m.get("Artist", {}).get("value")), ii["descriptionurl"]
def fetch_vsco():
"""Download the VSCO 2 CE subset listed in vsco_files.txt (skips files present)."""
lst = Path(__file__).resolve().parent / "vsco_files.txt"
base = "https://raw.githubusercontent.com/sgossner/VSCO-2-CE/master/"
for rel in lst.read_text().split("\n"):
if not rel:
continue
dest = ROOT / "vsco2ce" / rel
if dest.exists() and dest.stat().st_size:
continue
dest.parent.mkdir(parents=True, exist_ok=True)
req = urllib.request.Request(base + urllib.parse.quote(rel), headers=UA)
dest.write_bytes(urllib.request.urlopen(req, timeout=300).read())
print("vsco", rel)
def main():
fetch_vsco()
rows = [["set_or_file", "local_path", "license", "author", "source_url", "notes"]]
for sub, note in VSCO_SETS:
rows.append([sub, f"samples/vsco2ce/{sub}", "CC0-1.0", "Versilian Studios (VSCO 2 CE)",
"https://github.com/sgossner/VSCO-2-CE", note])
bird_dir = ROOT / "birds"
bird_dir.mkdir(parents=True, exist_ok=True)
for title, local, species, why in BIRDS:
url, lic, author, page = commons(title)
dest = bird_dir / local
if not dest.exists():
time.sleep(8) # Commons 429s back-to-back media downloads
req = urllib.request.Request(url, headers=UA)
dest.write_bytes(urllib.request.urlopen(req, timeout=120).read())
rows.append([f"{species}", f"samples/birds/{local}", lic, author, page, why])
print(local, lic, author)
with open(MANIFEST, "w", newline="") as f:
csv.writer(f).writerows(rows)
if __name__ == "__main__":
main()Sources and licences
- Hydrometric and climate data: Environment and Climate Change Canada, Water Survey of Canada. Contains information licensed under the Open Government Licence – Canada.
- Gross primary productivity: Running, S. and Zhao, M. (2021). MODIS/Terra Gross Primary Productivity Gap-Filled 8-Day L4 Global 500m SIN Grid V061 (MOD17A2HGF). NASA Land Processes DAAC. doi:10.5067/MODIS/MOD17A2HGF.061. Retrieved through the ORNL DAAC MODIS/VIIRS subset service.
- Wells Dam counts: Columbia River DART, Columbia Basin Research, University of Washington. Adult Passage Daily Counts. https://www.cbr.washington.edu/dart/query/adult_daily. Data courtesy of Douglas County PUD (Wells Dam). The raw counts are not redistributed here;
fetch_okanagan.pyretrieves them. - Instruments: VSCO 2 Community Edition, Versilian Studios (CC0). Upright piano sampled by Simon Dalzell (Ivy Audio).
- Birds (Wikimedia Commons): Swainson’s thrush, Niels Krabbe (CC BY-SA 4.0); American dipper, NPS & MSU Acoustic Atlas / Jennifer Jerrett (public domain); western meadowlark and sandhill crane, Jonathon Jongsma (CC BY-SA 3.0 / 4.0); bank swallow, Pascal Christe (CC BY-SA 4.0).
- Licences for this work: code under the Mozilla Public License 2.0; the Water Year audio CC BY-SA 4.0 (it contains CC BY-SA recordings); the Climate Pair audio CC BY-NC-SA 4.0 (see the README).