brt tutorial - pca

[eeholmes]

describing the spectral shape with pca

Pre-process spectra

• Optionally log-transform: log10(Rrs) (helps a lot with dynamic range)
• Remove any very noisy bands if needed (e.g., extremes of the spectrum)
• Standardize across samples per band (this is useful before PCA):

For each wavelength: subtract mean, divide by std

• PCA cares about covariance structure; this is standard practice
• Compute PCA on the 128 bands
• Use a big set of spectra across your training region(s)

Get:

• Loadings: PCᵢ(λ) → “basis spectral shapes”
• Scores: for each sample, how much of each PC it has
• Choose K PCs

e.g. enough to explain 95–99% of variance, or just fix K=5–10 for teaching

So now instead of 128 bands, each spectrum becomes: PC1_score, PC2_score, …, PCK_score

• Use PC scores as predictors in your BRT or
• Generate ratios among PCs (or between PCs and a reference band).

So your model sees:

“Here is the overall brightness, here is the blue–green tilt, here is the curvature, plus SST and location; please predict CHL / bbp700.”

You never gave it a blue/green ratio explicitly—just the PCs, which encode those combinations.

[eeholmes]

use bands

You basically want:

A systematic way to generate a moderate number of candidate ratios so the model can ‘discover’ which ones matter.

And you’d like it to feel like real hyperspectral practice (band selection, dimensionality reduction), not a toy hack.

Here’s a strategy that ticks all those boxes.

1. Don’t ratio all 128 bands — compress first

Common real-world move with hyperspectral:
compress spectra into a small set of “super-bands” that still capture the main shape.

You can do this a couple of ways; for a tutorial I’d use something simple and transparent:

Option A: Fixed wavelength windows (my pick for teaching)

Pick K wavelength windows across your spectrum, e.g. K=10:

blue: 400–430, 430–460

cyan/green: 460–490, 490–520, 520–550

yellow/red/NIR: 550–580, 580–610, 610–640, 640–670, 670–700

For each pixel/sample, compute the mean Rrs in each window:

pseudo-ish

windows = [
(400, 430),
(430, 460),
(460, 490),
(490, 520),
(520, 550),
(550, 580),
(580, 610),
(610, 640),
(640, 670),
(670, 700),
]

for k, (wmin, wmax) in enumerate(windows):
mask = (wavelengths >= wmin) & (wavelengths < wmax)
X[f"Rrs_win{k}"] = Rrs[:, mask].mean(axis=1)

Now you’ve turned 128 bands into, say, 10 super-bands that roughly correspond to spectral regions people care about, and that feels very “hyperspectral-ish”.

Option B: Representative subsampling

Simpler variant:

Just take every Nth wavelength (e.g. every 10th band) → 10–12 representative discrete bands.

Works fine too, just a bit less “spectral-region-y” than the windows.

2. Now let the model “discover” ratios among those super-bands

With K=10 super-bands, the pairwise ratios are:

10 choose 2 = 45 ratios → totally manageable.

Generate them programmatically:

import itertools
import numpy as np

bands = [f"Rrs_win{k}" for k in range(10)]

for i, j in itertools.combinations(bands, 2):
Xi, Xj = X[i].values, X[j].values
X[f"{i}over{j}"] = Xi / np.maximum(Xj, 1e-6) # avoid div by zero

Then your BRT gets:

The 10 super-bands

The 45 ratios between them

You didn’t decide “blue/green is important”; you just gave it many possible ratios and said: you tell me which are useful.

Teaching win: you can show feature importance and say:

“Look, out of 45 possible ratios, the model strongly prefers things like Rrs_win1 / Rrs_win4 – which is basically a blue/green index, discovered by the model.”

3. Connect this to what people actually do with hyperspectral

This is very much in line with real workflows:

Band grouping / binning to reduce spectral resolution (e.g., OCI → sensor-agnostic bands, or “bio-optical bands”).

Band selection / feature selection where people:

define a library of candidate indices (ratios, differences, slopes),

then use ML to rank/select them.

Index discovery: several papers literally auto-generate a big set of candidates (like NDVI variants) and let ML pick the best.

Your tutorial version is just a tidy, reproducible version of that.

If you want to make it more explicitly “hyperspectral-y”, you can:

Define windows that correspond loosely to:

“chlorophyll absorption region”

“CDOM sensitive region”

“backscattering-dominated region”

Then note when the model loves ratios between those.

4. A nice teaching story structure

You can structure the tutorial like this:

Model 1: raw (or log) super-bands only

Features: Rrs_win0..9

Train BRT, show performance.

Model 2: super-bands + all ratios

Features: Rrs_win0..9 + 45 ratios

Train BRT, show improved performance (ideally especially in some regimes).

Interpretation

Plot feature importance.

Show that the top features are a mix of:

Certain windows (e.g., mid-green, red)

Certain ratios (e.g., blue/green, green/red)

Message:

“We didn’t tell it which ratio to use.
We just provided a lot of possible combinations and let the model discover which ones encode useful information.”

That hits the conceptual target you want: ML + feature space exploration → emergent “interesting indices.”

[eeholmes]

re whether pca is capturing useful ratios

1. No, PCA doesn’t guarantee a blue/green ratio

PCA just says:

“Here are orthogonal directions in 128-D space that explain the most variance.”

What those directions mean depends on the data:

PC1 is usually “overall brightness” (all loadings same sign).

PC2 is often something slope-like (blue vs green/red).

But if your training set is weird (e.g. dominated by one water type), PCs might look different.

So: you’re right — you can’t assume PC2 must be a blue/green ratio.

2. How to check if a PC encodes a blue/green ratio

This is where you can build a nice teaching moment:

Compute PCA on your spectra (after log-transform + standardization).

You get:

loadings[pc, wavelength]

scores[sample, pc]

For interpretation only (not fed to the model), compute a few “classic” ratios on the same samples, e.g.:

B/G = Rrs(440) / Rrs(555)

B2/G2 = Rrs(443) / Rrs(560)

maybe a red/blue ratio

Correlate each PC score with each ratio:

corr(PC1_score, B/G)

corr(PC2_score, B/G), etc.

If you see something like:

corr(PC2, B/G) = 0.95

corr(PC1, B/G) = 0.1

corr(PC3, B/G) = 0.2

…then you can very cleanly say in the tutorial:

“We never told the PCA what ‘blue/green’ is.
It discovered a component (PC2) whose scores are almost perfectly correlated with the blue/green ratio. So PC2 is essentially a data-driven blue/green index.”

You didn’t hard-code a ratio into the model; you’re just using the ratio afterward as a diagnostic to interpret what PCA learned.

That’s a great “ML can rediscover known ocean-color indices” story.

3. Increase the chance that a PC really is “shape, not brightness”

If you want PCA to pick up shape (like blue/green balance) instead of just “everything brighter or darker,” you can pre-condition your spectra a bit:

Log-transform: log10(Rrs)
(turns multiplicative relationships into additive; helps a lot.)

Normalize per spectrum (optional but useful):

Divide each spectrum by its mean or integral across λ:

This removes overall brightness, leaving mostly shape.

Then standardize each band across samples and do PCA.

Now the biggest modes of variance are much more likely to be:

“Is this spectrum bluer or greener than average?”

“Does it have a stronger dip/bump around the chlorophyll region?”

Rather than just “brighter everywhere”.

That makes it more likely that PC2 or PC3 will look like a blue/green tilt, which you can then verify with the correlations trick above.

4. How this plays in a teaching notebook

You can structure it like:

Exploratory PCA on spectra

Preprocess spectra (log, maybe normalize).

Run PCA, keep first 3–5 PCs.

Plot PC loadings vs wavelength

“Look, PC1 is roughly flat; PC2 is more positive in the blue, more negative in the green…”

Compute a couple of “classic” color ratios for interpretation only.

Show correlation table or scatter of PC2_score vs B/G.

Explain:

“We didn’t tell the algorithm anything about blue/green. It chose a direction in spectral space that corresponds very closely to that contrast — because that’s an important axis of variation in the data.”

Then, in the supervised part:

Feed PC scores + SST + lat/lon into your BRT for CHL/bbp700.

Show that PC2, PC3, etc. are among the most important features.

So you get both:

Unsup step: PCA “discovers” latent spectral patterns (one of which ≈ blue/green).

Sup step: BRT shows which of those patterns actually matter for CHL/bbp700.

5. TL;DR

You’re right: PCA doesn’t promise a blue/green ratio by fiat.

But you can:

preprocess spectra so PCA focuses on shape,

run PCA,

and then post-hoc check which PCs correlate strongly with known ratios.

That gives you exactly the teaching story you want:

“We let the algorithm discover spectral patterns, then we checked and saw one of those patterns is effectively the classic blue/green ratio.”

If you’d like, I can sketch a small “analysis block” (pseudo-Python) that:
PCA → compute a couple of ratios → print a correlation matrix → plot PC2 vs blue/green, ready to drop into a tutorial notebook.