Lesson 4.1.1: Structuring hyperspectral data for model training

flowchart TD
    A[Raw hyperspectral images] --> B[Pre‑processing]
    B --> C[Calibration & Noise reduction]
    B --> D[Band selection / Dimensionality reduction]
    C --> E[Normalisation]
    D --> E
    E --> F[Data labelling]
    F --> G[Train‑test split]
    G --> H[Model training]
    H --> I[Evaluation & Validation]

Additional Example: Spectral Normalisation and Augmentation

In many pipelines, normalising each band to a common scale and applying simple augmentations can improve model robustness. Below is a quick example that demonstrates per‑band min‑max scaling followed by random spectral shifts.

import numpy as np
import rasterio
import random
from pathlib import Path

Normalise each band to [0,1] using min-max per band
def minmax_normalise(arr):
    mins = arr.min(axis=(1,2), keepdims=True)
    maxs = arr.max(axis=(1,2), keepdims=True)
    return (arr - mins) / (maxs - mins + 1e-8)

Randomly shift spectral bands to simulate illumination changes
def spectral_shift(arr, shift_range=(-0.05, 0.05)):
    shift = np.random.uniform(*shift_range, size=(arr.shape[0], 1, 1))
    return arr + shift

def augment_scene(scene_path, out_path):
    with rasterio.open(scene_path) as src:
        data = src.read()  # shape (bands, H, W)
        data = minmax_normalise(data)
        data = spectral_shift(data)
        meta = src.meta.copy()
        meta.update(dtype=rasterio.float32)

        with rasterio.open(out_path, 'w', **meta) as dst:
            dst.write(data.astype(rasterio.float32))

if __name__ == "__main__":
    scene = Path("hyperspectral_project/raw/sensorA/scene_001/scene_001.tif")
    out = Path("hyperspectral_project/processed/normalised") / "scene_001_norm.tif"
    out.parent.mkdir(parents=True, exist_ok=True)
    augment_scene(scene, out)

This snippet can be plugged into the preprocess.py workflow to generate a second set of tiles that are normalised and augmented for training.

Comparing Organisation Approaches

Implementation/Code: Preprocessing Pipelines

When preparing hyperspectral data for machine‑learning models, the preprocessing stage is arguably the most critical. The raw measurements arrive as large, multi‑dimensional arrays that are rarely ready for ingestion by a neural network. Normalisation, resampling, and efficient batching must be orchestrated so that each pipeline stage respects the physical meaning of the data while keeping memory usage under control. Below we walk through a typical workflow in Python, emphasising the rationale behind each step and providing annotated code snippets that you can adapt to your own projects.

1. Normalising the Spectral Cubes

Hyperspectral sensors record radiance or reflectance values that span many orders of magnitude. Models trained on such raw values often converge slowly or become numerically unstable. A common remedy is to apply per‑band normalisation, where each spectral channel is centred and scaled independently. This preserves inter‑band relationships while reducing dynamic range.

We illustrate a simple yet powerful approach using scikit‑learn's StandardScaler. The example assumes the data is stored in a numpy array of shape (N, H, W, B), where N is the number of samples, H and W are spatial dimensions, and B is the spectral dimension.

import numpy as np
from sklearn.preprocessing import StandardScaler

Dummy hyperspectral dataset: 100 samples, 64x64 pixels, 200 spectral bands
X = np.random.rand(100, 64, 64, 200).astype(np.float32)

Reshape to (NHW, B) so each pixel becomes a sample for scaling
pixels = X.reshape(-1, X.shape[-1])

Initialise the scaler
scaler = StandardScaler()

Fit on the entire dataset (or on a validation set if you prefer)
scaler.fit(pixels)

Apply the transformation
X_norm = scaler.transform(pixels)

Reshape back to the original shape
X_norm = X_norm.reshape(X.shape)

Line‑by‑line explanation:

1. import numpy as np and from sklearn.preprocessing import StandardScaler bring in the required libraries.
2. The random array X mimics a real hyperspectral stack; in practice you would load your data from a .tiff or .mat file.
3. pixels = X.reshape(-1, X.shape[-1]) flattens the spatial dimensions. Each row now represents a single pixel’s spectral signature.
4. We instantiate StandardScaler(), which will compute the mean and standard deviation for each spectral band.
5. scaler.fit(pixels) calculates these statistics across all pixels in the dataset.
6. scaler.transform(pixels) centres and scales every band to unit variance.
7. Finally, reshape restores the four‑dimensional shape, ready for downstream processing.

Alternative strategies include min‑max scaling, logarithmic compression, or even spectral‑band‑wise clipping to a specified percentile range. The choice depends on sensor characteristics and the model’s sensitivity to outliers.

# Alternative: Min‑Max scaling to [0, 1] per band
from sklearn.preprocessing import MinMaxScaler

minmax_scaler = MinMaxScaler()
pixels_minmax = minmax_scaler.fit_transform(pixels)
X_minmax = pixels_minmax.reshape(X.shape)

2. Resampling to a Common Grid

Hyperspectral data acquired from different platforms often exhibit varying spectral resolutions. To train a single model, you need all inputs on a shared wavelength grid. Resampling can be performed via interpolation or spectral binning. Below we demonstrate cubic spline interpolation using scipy.interpolate.

from scipy.interpolate import interp1d

Original wavelengths (nm), e.g., 400–1000 nm sampled every 3 nm
orig_wl = np.arange(400, 1001, 3)

Desired common grid (e.g., 400–1000 nm every 5 nm)
common_wl = np.arange(400, 1001, 5)

Function to resample a single spectral cube
def resample_cube(cube, orig_wl, common_wl):
    # cube shape: (H, W, B)
    H, W, B = cube.shape
    # Create an array to hold resampled data
    resampled = np.empty((H, W, len(common_wl)), dtype=cube.dtype)
    # Resample each pixel individually
    for i in range(H):
        for j in range(W):
            f = interp1d(orig_wl, cube[i, j, :], kind='cubic')
            resampled[i, j, :] = f(common_wl)
    return resampled

Example usage on a single sample
sample_cube = X_norm[0]  # shape (64, 64, 200)
sample_resampled = resample_cube(sample_cube, orig_wl, common_wl)

This routine loops over every pixel, constructs a cubic spline interpolator, and evaluates it on the new grid. For large datasets, vectorised or GPU‑accelerated resampling (e.g., torch.nn.functional.interpolate with custom kernels) are preferable.

3. Constructing Data Loaders with PyTorch

Deep learning frameworks typically expect data to be streamed in batches. PyTorch’s DataLoader class, coupled with a custom Dataset, provides a flexible way to feed normalised, resampled cubes into a neural network. Below we define a straightforward dataset that encapsulates the three preprocessing steps.

import torch
from torch.utils.data import Dataset, DataLoader

class HyperspectralDataset(Dataset):
    def __init__(self, cubes, labels, scaler, orig_wl, common_wl):
        ""
        cubes: list or array of shape (N, H, W, B)
        labels: array of shape (N,)
        scaler: fitted StandardScaler instance
        orig_wl, common_wl: wavelength arrays for resampling
        ""
        self.cubes = cubes
        self.labels = labels
        self.scaler = scaler
        self.orig_wl = orig_wl
        self.common_wl = common_wl

    def __len__(self):
        return len(self.cubes)

    def __getitem__(self, idx):
        # Normalise
        cube = self.cubes[idx]
        pixels = cube.reshape(-1, cube.shape[-1])
        cube_norm = self.scaler.transform(pixels).reshape(cube.shape)

        # Resample
        cube_resampled = resample_cube(cube_norm, self.orig_wl, self.common_wl)

        # Convert to torch tensor
        tensor = torch.from_numpy(cube_resampled).permute(2, 0, 1)  # (B, H, W)
        label = torch.tensor(self.labels[idx], dtype=torch.long)
        return tensor, label

Instantiate dataset
dataset = HyperspectralDataset(X_norm, np.random.randint(0, 5, size=100),
                               scaler, orig_wl, common_wl)

Create DataLoader
loader = DataLoader(dataset, batch_size=8, shuffle=True, num_workers=4)

Key points:

• __getitem__ normalises each cube on the fly, ensuring that the scaler’s parameters remain fixed and that no data leakage occurs between training and validation splits.
• Resampling is performed per sample; if the target grid is the same for all samples, you could pre‑compute this step to speed up training.
• Reordering the dimensions to (B, H, W) matches the expectation of convolutional layers that treat the spectral axis as channels.
• Using num_workers=4 offloads data loading to separate processes, keeping the GPU busy.

4. Benchmarking Memory Footprint

Preprocessing can be memory‑intensive, especially when working with gigabyte‑scale datasets. It is prudent to monitor peak memory usage and optimise accordingly. The following table summarises typical memory consumption for a batch of 8 samples, each 64×64×200, under the pipeline above.

Strategies to reduce memory:

1. In‑place operations: normalise directly on the original array.
2. Half‑precision: cast tensors to torch.float16 if the model tolerates it.
3. On‑the‑fly resampling: perform interpolation only for the batch being processed, discarding the intermediate result after conversion to tensor.
4. Chunking the dataset: split the dataset into smaller files and load them sequentially.

5. Integrating with a Convolutional Neural Network

Once the data is in a DataLoader, you can feed it into a CNN. A minimal architecture for hyperspectral classification might look like this:

import torch.nn as nn
import torch.nn.functional as F

class SimpleHSNet(nn.Module):
    def __init__(self, in_channels, num_classes):
        super().__init__()
        self.conv1 = nn.Conv2d(in_channels, 32, kernel_size=3, padding=1)
        self.conv2 = nn.Conv2d(32, 64, kernel_size=3, padding=1)
        self.pool = nn.MaxPool2d(2)
        self.fc1 = nn.Linear(64  16  16, 128)
        self.fc2 = nn.Linear(128, num_classes)

    def forward(self, x):
        x = F.relu(self.conv1(x))
        x = self.pool(x)
        x = F.relu(self.conv2(x))
        x = self.pool(x)
        x = x.view(x.size(0), -1)
        x = F.relu(self.fc1(x))
        return self.fc2(x)

Here, in_channels equals the number of spectral bands after resampling (len(common_wl)). The pooling layers reduce spatial dimensions, allowing the fully connected layers to operate on a manageable feature vector.

6. Practical Tips and Common Pitfalls

• Sensor‑specific corrections: before normalisation, apply radiometric calibration, dark‑current subtraction, and atmospheric correction if available.
• Batch‑wise vs dataset‑wide scaling: fitting the scaler on the entire training set preserves global statistics but may leak information into validation folds if not carefully split.
• Interpolation artefacts: cubic interpolation can introduce ringing near sharp spectral features; consider windowed sinc or linear interpolation for smoother results.
• GPU memory vs CPU memory: keep the preprocessing on the CPU until tensors are ready; otherwise you risk saturating GPU memory with intermediate buffers.
• Reproducibility: set random seeds for all libraries (np.random.seed, torch.manual_seed, etc.) to guarantee deterministic behaviour.

By following the steps outlined above, you create a robust preprocessing pipeline that normalises, resamples, and batches hyperspectral data efficiently. This foundation enables the subsequent training of deep learning models that can leverage the rich spectral information inherent in your measurements.

9. Alternative Normalisation Strategy: Per‑Pixel Min‑Max Scaling

Some applications benefit from normalising each pixel independently so that every spectrum spans the full dynamic range. This can be useful when the absolute reflectance values are less informative than the relative spectral shape.

# Per‑pixel min‑max scaling to range [0, 1]
pixel_min = samples.min(axis=1, keepdims=True)
pixel_max = samples.max(axis=1, keepdims=True)
samples_minmax = (samples - pixel_min) / (pixel_max - pixel_min + 1e-8)

print("First pixel before scaling:", samples[0, :5])
print("First pixel after scaling:", samples_minmax[0, :5])

Data Preparation Checklist

Before embarking on model training, confirm that your dataset satisfies the following criteria. A quick audit against these items can save hours of debugging later.

• All scenes are labelled consistently, with a clear mapping from class names to numeric identifiers.

• Spectral metadata (wavelengths, sensor type, acquisition date) is stored in a machine‑readable format such as JSON or CSV and is accessible to all preprocessing scripts.

• The raw data has been verified for integrity (checksum, file size) and is free of corrupted tiles that could crash loaders.

• The chosen hold‑out split preserves spatial contiguity where necessary, avoiding leakage of neighbouring pixels between training and test sets.

Workflow Overview

The end‑to‑end pipeline can be summarised in three high‑level stages, each encapsulating a set of scripts and artefacts. Visualising the flow helps identify bottlenecks and ensures reproducibility across teams.

1. Data Ingestion – Load raw files, validate integrity, and store a canonical copy in a version‑controlled repository (e.g., Git LFS or a dedicated data store).

2. Pre‑processing – Apply calibration, normalisation, resampling, and augmentation in a deterministic, reproducible script. Output should be a set of normalised, resampled tiles ready for model consumption.

3. Training & Evaluation – Use the prepared tiles to train a model, monitor training dynamics, and generate evaluation artefacts (confusion matrices, ROC curves). Store checkpoints and final metrics in a structured experiment log.

Pitfalls: Common Issues and Fixes

When organising hyperspectral imagery for supervised model training, a handful of routine missteps can erode model performance more quickly than any hyper‑parameter tweak. Below we catalogue the most frequent errors—spectral leakage, class imbalance, and over‑fitting—followed by concrete strategies to mitigate each. The discussion is deliberately technical, aimed at practitioners who already understand the basics of spectral pre‑processing and machine learning pipelines.

1. Spectral Leakage: The Invisible Noise Source

Spectral leakage occurs when the discrete Fourier transform (DFT) of a signal is contaminated by energy from neighbouring wavelengths. In hyperspectral data this is often a consequence of inadequate windowing or of mixing pixels that have been coarsely resampled. Leakage can manifest as subtle smears in the reflectance curve, leading to mis‑classification in models that rely on fine spectral features.

Typical symptoms:

• Unexpected spikes in the spectral signature of a class that do not match ground truth.
• Poor separability in Principal Component Analysis (PCA) or Linear Discriminant Analysis (LDA) plots.
• High variance in cross‑validation scores across neighbouring bands.

Fixes:

1. Apply a Hamming or Hann window before computing the DFT if you use frequency‑domain augmentation.
2. Use spectral band selection to exclude bands that are known to be noisy (e.g., water absorption bands).
3. When resampling, adopt cubic convolution or Lanczos interpolation rather than nearest‑neighbour to preserve spectral continuity.

Illustration of a simple windowing step in Python:

import numpy as np
from scipy.signal import windows

Assume `spectrum` is a 1‑D numpy array of reflectance values
def apply_hamming_window(spectrum):
    ""
    Multiply the input spectrum by a Hamming window to reduce spectral leakage.
    ""
    # 1. Generate a Hamming window of the same length as the spectrum
    window = windows.hamming(len(spectrum))
    # 2. Element‑wise multiplication preserves the overall amplitude
    return spectrum * window

Example usage:
spectrum = np.array([...])  # Replace with real data
cleaned_spectrum = apply_hamming_window(spectrum)

Line‑by‑line explanation:

• import numpy as np imports the numerical library.
• from scipy.signal import windows brings in the window functions.
• Function apply_hamming_window takes a 1‑D array.
• Inside the function, windows.hamming(len(spectrum)) creates a Hamming window matching the spectrum length.
• The window is multiplied element‑wise with the spectrum, attenuating edge effects.
• The cleaned spectrum is returned for downstream processing.

2. Class Imbalance: When Rare Signals Get Overlooked

Hyperspectral datasets frequently exhibit skewed class distributions—e.g., a sparse canopy class versus a vast background. If left uncorrected, models will bias toward the majority class, yielding deceptively high overall accuracy but poor recall for minority classes.

Common indicators:

• Accuracy > 90 % but F1‑score for minority class < 0.3.
• Confusion matrix dominated by diagonal entries for the majority class.
• Training loss plateauing while validation loss diverges.

Mitigation tactics:

• Resampling: either oversample minority classes (e.g., SMOTE) or undersample the majority class.
• Class‑weighted loss functions: assign higher penalties to minority errors.
• Use stratified splits during cross‑validation to preserve class ratios.
• Implement focal loss for deep learning models to focus learning on hard examples.

Python example of a class‑weighted Random Forest:

from sklearn.ensemble import RandomForestClassifier
from sklearn.utils.class_weight import compute_class_weight
import numpy as np

X_train: feature matrix, y_train: class labels
class_weights = compute_class_weight(
    class_weight='balanced',
    classes=np.unique(y_train),
    y=y_train
)
weight_dict = dict(zip(np.unique(y_train), class_weights))

rf = RandomForestClassifier(
    n_estimators=200,
    class_weight=weight_dict,
    random_state=42
)
rf.fit(X_train, y_train)

Explanation:

1. Import the RandomForestClassifier and a util to compute class weights.
2. Calculate weights inversely proportional to class frequencies.
3. Map each class to its weight.
4. Instantiate the classifier with class_weight=weight_dict, ensuring minority classes exert greater influence.
5. Fit the model on the training data.

3. Over‑Fitting: The Model Learns Noise, Not Signal

Over‑fitting is the classic pitfall where a model captures idiosyncratic patterns in the training set that do not generalise. In hyperspectral work, high dimensionality combined with limited samples can exacerbate this.

Symptoms:

• Training accuracy close to 100 % while validation accuracy lags significantly.
• Large standard deviation in validation metrics across folds.
• Model predictions fluctuate wildly when presented with spectrally similar but unseen samples.

Prevention strategies:

1. Dimensionality reduction: PCA, t‑SNE, or autoencoders to summarise spectra.
2. Regularisation: L1/L2 penalties in linear models, dropout in neural nets.
3. Use early stopping based on validation loss.
4. Apply cross‑validation with a high number of folds to obtain a robust estimate of generalisation.
5. Limit model complexity: fewer layers/neurons or trees in ensemble methods.

Code snippet for PCA‑based feature reduction followed by a Gradient Boosting model:

from sklearn.decomposition import PCA
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

pipeline = Pipeline([
    ('scaler', StandardScaler()),      # Zero‑mean, unit‑variance
    ('pca', PCA(n_components=0.95)),   # Retain 95 % variance
    ('gb', GradientBoostingClassifier(
        n_estimators=300,
        learning_rate=0.05,
        max_depth=5,
        random_state=123
    ))
])

pipeline.fit(X_train, y_train)

Line‑by‑line:

1. Import PCA, GradientBoostingClassifier, Pipeline, and StandardScaler.
2. Instantiate a pipeline that first standardises data.
3. Apply PCA, automatically selecting enough components to explain 95 % of variance.
4. Fit a Gradient Boosting model with a moderate learning rate and depth.
5. Train the entire pipeline on the training set.

4. Neglecting Spatial Context: Treating Pixels as Isolated

Many novices treat each pixel as an independent sample, discarding the valuable spatial relationships inherent in imagery. This can lead to high intra‑class variance and over‑fitting.

Remedies:

• Incorporate spatial‑spectral fusion by concatenating local patches or using 3‑D convolutions.
• Use Markov Random Fields or conditional random fields to post‑process predictions.
• Apply morphological opening/closing to smooth prediction maps.

5. Inadequate Validation Design

Validation protocols that ignore the spatial autocorrelation inherent in remote sensing data can inflate performance estimates. Classic random splits may inadvertently place adjacent pixels in training and test sets, leading to overly optimistic results.

Best practice:

• Implement stratified block cross‑validation where each fold is a spatially contiguous block.
• Use a leave‑one‑scene-out strategy when multiple scenes are available.
• Report metrics per class and per scene to expose hidden biases.

6. Feature Scaling Missteps

In hyperspectral work, reflectance values are typically already normalised (0–1). However, when adding derived indices (NDVI, SAVI, etc.) or external covariates (DEM, soil maps), inconsistent scaling can distort model learning.

Solution:

1. Apply robust scaling (median and IQR) to outlier‑prone features.
2. When mixing data types, standardise all features to zero mean and unit variance.
3. Document scaling choices in the metadata for reproducibility.

7. Forgetting to Monitor Training Dynamics

Without visualising loss curves or per‑epoch metrics, subtle over‑fitting may go unnoticed.

Recommendation:

• Plot training and validation loss after each epoch.
• Track additional metrics such as confusion matrix evolution.
• Use TensorBoard or similar tools for deep learning workflows.

Summative Checklist