Hyperspectral deep dive — reading minerals from real EnMAP spectra¶

Take-home notebook · Workshop 2 · IC-TWEP 2026 · QuestFeed¶

An imaging spectrometer records a full spectrum — about 200 usable bands from 0.42 to 2.45 µm — for every 30 m pixel. In this notebook you will work directly with a real EnMAP scene:

  1. look at the image and at the spectra of individual pixels;
  2. remove the continuum and measure the depth of diagnostic absorptions;
  3. classify every pixel by its best match in the USGS Spectral Library with the spectral angle mapper;
  4. compare the result with NASA EMIT's independent mineral identification.

Data: an EnMAP L2A surface-reflectance subset (6 × 6 km; © DLR) and the USGS Spectral Library 7, both from the workshop data store; EMIT L2B mineral families.

Ready.
In [2]:
# @title 1 · Load an EnMAP subset
SUBSET = "Mahd adh Dhahab — Saudi Arabia"  # @param ["Mahd adh Dhahab — Saudi Arabia", "Sohar–Lasail — Oman", "Bou Azzer — Morocco"]
subsets = {s["label"]: s for s in fetch_json("hyperspectral/subsets.json")}
S = subsets[SUBSET]
cube_q, c_tr, c_crs = fetch_raster(f"hyperspectral/{S['id']}_enmap.tif")
bands = fetch_csv(f"hyperspectral/{S['id']}_enmap_bands.csv")           # wavelength_nm, fwhm_nm
wl = bands["wavelength_nm"].values / 1000.0                               # µm
cube = np.where(cube_q == -32768, np.nan, cube_q / 10000.0).astype("float32")   # reflectance
B, H, W = cube.shape
def band(nm):
    return int(np.argmin(abs(bands["wavelength_nm"].values - nm)))
rgb = np.dstack([cube[band(650)], cube[band(560)], cube[band(480)]])
rgb = np.clip(rgb / np.nanpercentile(rgb, 99), 0, 1)
print(f"{S['name']}: {B} bands, {W} × {H} pixels at 30 m, EnMAP scene of {S['date']}")
fig, ax = plt.subplots(figsize=(7, 7)); ax.imshow(np.nan_to_num(rgb)); ax.set_title("EnMAP true colour"); ax.axis("off"); plt.show()

Mahd adh Dhahab district, Arabian Shield: 195 bands, 200 × 200 pixels at 30 m, EnMAP scene of 2026-04-10
No description has been provided for this image

Pixel spectra¶

Pick pixels by row and column (the image is 200 × 200). The defaults are spread across the subset; try your own after looking at the image. Gaps are atmospheric water-vapour bands, removed in processing.

In [3]:
# @title 2 · Plot spectra of chosen pixels
PIXELS = [(40, 40), (100, 100), (160, 60), (60, 160)]   # (row, column) — edit freely
fig = go.Figure()
for r, c in PIXELS:
    fig.add_trace(go.Scatter(x=wl, y=cube[:, r, c], mode="lines", name=f"pixel ({r}, {c})"))
for x0, x1 in [(0.85, 0.95), (2.15, 2.22), (2.30, 2.36)]:
    fig.add_vrect(x0=x0, x1=x1, fillcolor="#FDE68A", opacity=0.35, line_width=0)
fig.update_layout(title="EnMAP surface reflectance", xaxis_title="wavelength (µm)", yaxis_title="reflectance", height=430,
                  **{k: v for k, v in LAYOUT.items() if k != "hovermode"})
fig.show()

Continuum removal and absorption depth¶

Overall brightness varies with grain size, moisture and illumination; the shape of an absorption is what identifies a mineral. Dividing a spectrum by a straight-line continuum drawn between the shoulders of an absorption isolates the feature. Its depth = 1 − (reflectance at the centre ÷ continuum).

In [4]:
# @title 3 · Band-depth maps at 2.20 µm (Al-OH) and 2.33 µm (Mg-OH / carbonate)
def depth(left, centre, right):
    l, c, r = (np.nanmean([cube[band(w)] for w in ws], axis=0) for ws in (left, centre, right))
    wl_l, wl_c, wl_r = np.mean(left), np.mean(centre), np.mean(right)
    cont = l + (r - l) * (wl_c - wl_l) / (wl_r - wl_l)
    return 1 - c / cont
al = depth([2110, 2130], [2190, 2200, 2210], [2250, 2270])
mg = depth([2250, 2270], [2320, 2330, 2340], [2390, 2410])
fig, ax = plt.subplots(1, 2, figsize=(14, 6.2))
for a, d, t, cm in [(ax[0], al, "Al-OH depth at 2.20 µm (clays, white mica)", "magma"),
                    (ax[1], mg, "Mg-OH / carbonate depth at 2.33 µm", "viridis")]:
    im = a.imshow(np.clip(d, 0, np.nanpercentile(d, 99)), cmap=cm); a.set_title(t); a.axis("off")
    fig.colorbar(im, ax=a, shrink=0.75)
plt.tight_layout(); plt.show()
No description has been provided for this image

Spectral angle mapper against the USGS library¶

The library spectra are first resampled to EnMAP's bands with a Gaussian response of each band's width. Real surfaces are mixtures with sloping, weathered spectra, so we compare absorption shapes, not raw spectra: each spectrum is divided by a straight-line continuum across 2.10–2.40 µm, and the spectral angle is computed between the resulting absorption-depth profiles. Each pixel gets the mineral with the smallest angle — if that angle is below a threshold and the pixel shows a real absorption (depth > 2 %); otherwise it stays unclassified.

In [5]:
# @title 4 · Spectral angle mapper (SAM) on absorption shapes, 2.10–2.40 µm
lib = fetch_csv("spectra/usgs_mineral_spectra.csv")
lib = lib[lib.sensor == "lab"]
MINERALS = ["Kaolinite", "Alunite", "Muscovite (sericite)", "Montmorillonite", "Calcite", "Dolomite", "Chlorite",
            "Epidote", "Serpentine (antigorite)", "Talc", "Gypsum"]
sel = (wl >= 2.10) & (wl <= 2.40) & np.isfinite(cube).all(axis=(1, 2))
w_sel, fw = wl[sel], bands["fwhm_nm"].values[sel] / 1000.0
def resample(m):
    d = lib[lib.mineral == m].sort_values("wavelength_um")
    x, y = d.wavelength_um.values, d.reflectance.values
    g = np.exp(-0.5 * ((x[None, :] - w_sel[:, None]) / (fw[:, None] / 2.355)) ** 2)
    return (g * y).sum(1) / g.sum(1)
def depth_profile(spec):                       # 1 − spectrum / straight-line continuum across the window
    left, right = spec[..., :2].mean(-1, keepdims=True), spec[..., -2:].mean(-1, keepdims=True)
    t = (w_sel - w_sel[:2].mean()) / (w_sel[-2:].mean() - w_sel[:2].mean())
    return 1 - spec / (left + (right - left) * t)
X = depth_profile(cube[sel].reshape(sel.sum(), -1).T)          # pixels × bands
R = np.stack([depth_profile(resample(m)[None, :])[0] for m in MINERALS])
Xn = X / np.linalg.norm(X, axis=1, keepdims=True)
Rn = R / np.linalg.norm(R, axis=1, keepdims=True)
angles = np.degrees(np.arccos(np.clip(Xn @ Rn.T, -1, 1)))      # pixels × minerals
best = angles.argmin(1)
best_angle = angles.min(1)
has_feature = X.max(1) > 0.02
THRESHOLD_DEG = 15.0                                           # edit: smaller = stricter
label = np.where((best_angle <= THRESHOLD_DEG) & has_feature, best + 1, 0).reshape(H, W)
names = ["unclassified"] + MINERALS
counts = pd.Series(label.ravel()).value_counts()
display(pd.DataFrame({"mineral": [names[i] for i in counts.index], "pixels": counts.values,
                      "share (%)": np.round(100 * counts.values / label.size, 1)}).set_index("mineral"))
cmap = ListedColormap(["#F3F4F6"] + list(plt.get_cmap("tab20").colors[:len(MINERALS)]))
fig, ax = plt.subplots(figsize=(9, 7))
ax.imshow(label, cmap=cmap, vmin=0, vmax=len(MINERALS), interpolation="nearest"); ax.axis("off")
ax.set_title(f"Best-matching library mineral (angle ≤ {THRESHOLD_DEG}°)")
present = [i for i in counts.index if i > 0][:10]
ax.legend([plt.Rectangle((0, 0), 1, 1, color=cmap(i)) for i in present], [names[i] for i in present],
          loc="center left", bbox_to_anchor=(1.02, 0.5), frameon=False)
plt.show()
pixels share (%)
mineral
unclassified 34323 85.8
Epidote 4041 10.1
Calcite 1635 4.1
Chlorite 1 0.0
No description has been provided for this image

Does an independent instrument agree?¶

NASA's EMIT identifies minerals with a different instrument, resolution and algorithm (the Tetracorder expert system). Agreement between the two is evidence; disagreement is a prompt to look closer.

In [6]:
# @title 5 · Compare with EMIT's mineral families
emit, e_tr, e_crs = fetch_raster(f"districts/{S['district']}/emit_minerals.tif")
emit_on_enmap = np.zeros((H, W), np.uint16)
reproject(emit[0], emit_on_enmap, src_transform=e_tr, src_crs=e_crs, dst_transform=c_tr, dst_crs=c_crs,
          resampling=Resampling.nearest)
legend = fetch_json(f"districts/{S['district']}/emit_legend.json")["swir"]
family = {"Kaolinite": 4, "Alunite": 7, "Muscovite (sericite)": 6, "Montmorillonite": 5, "Calcite": 1, "Dolomite": 1,
          "Chlorite": 2, "Epidote": 2, "Serpentine (antigorite)": 3, "Talc": 3, "Gypsum": 7}
sam_family = np.zeros((H, W), np.uint16)
for i, m in enumerate(MINERALS, start=1):
    sam_family[label == i] = family[m]
both = (sam_family > 0) & (emit_on_enmap > 0)
agree = (sam_family == emit_on_enmap) & both
print(f"Pixels classified by both: {both.sum():,}; same mineral family: {100 * agree.sum() / max(both.sum(), 1):.1f} %")
# The trade-off between coverage and confidence: re-run the decision at several thresholds
rows = []
for thr in (10, 15, 20, 25):
    lab_t = np.where((best_angle <= thr) & has_feature, best + 1, 0).reshape(H, W)
    fam_t = np.zeros((H, W), np.uint16)
    for i, m in enumerate(MINERALS, start=1):
        fam_t[lab_t == i] = family[m]
    b_t = (fam_t > 0) & (emit_on_enmap > 0)
    rows.append({"threshold (°)": thr, "pixels classified (%)": round(100 * (lab_t > 0).mean(), 1),
                 "agreement with EMIT (%)": round(100 * ((fam_t == emit_on_enmap) & b_t).sum() / max(b_t.sum(), 1), 1)})
display(pd.DataFrame(rows).set_index("threshold (°)"))
tab = pd.crosstab(pd.Series([legend[str(v)] for v in emit_on_enmap[both]], name="EMIT family"),
                  pd.Series([legend[str(v)] for v in sam_family[both]], name="EnMAP SAM family"))
display(tab)

Pixels classified by both: 5,675; same mineral family: 73.5 %
pixels classified (%) agreement with EMIT (%)
threshold (°)
10 1.1 97.1
15 14.2 73.5
20 35.7 57.9
25 56.6 47.7
EnMAP SAM family Carbonates (calcite, dolomite) Chlorite / epidote
EMIT family
Carbonates (calcite, dolomite) 1591 1374
Chlorite / epidote 26 2581
White mica / illite (sericite) 18 85

Take it further¶

  • The table in Step 5 shows the core trade-off: a strict threshold classifies few pixels with high agreement; a loose one classifies many with less. Choose it for your purpose — targeting needs confidence, screening needs coverage.
  • Add minerals from the library (e.g. jarosite, pyrophyllite) and see what changes.
  • Restrict the angle to a single feature (e.g. 2.12–2.28 µm) to separate kaolinite from white mica by the shape of the 2.2 µm doublet.

Data: EnMAP L2A — contains modified EnMAP data © DLR · USGS Spectral Library Version 7 (Kokaly et al. 2017, public domain) · EMIT L2B mineral identification (NASA JPL / LP DAAC). Method: spectral angle mapper (Kruse et al., 1993).

QuestFeed · IC-TWEP 2026 take-home materials · questfeed.com