Lesson 7.2.2: Understanding model decisions in hyperspectral classification

flowchart TD
    A[Raw hyperspectral image] --> B[Pre‑processing & spectral feature extraction]
    B --> C[Machine‑learning model (e.g., Random Forest)]
    C --> D[Predicted class for each pixel]
    D --> E[Decision explanation (e.g., SHAP values)]
    E --> F[Visualise importance per band]
    D --> G[Confusion matrix & accuracy metrics]
    F --> G

Concept

In precision agriculture, hyperspectral imagery provides a wealth of spectral signatures that can be translated into actionable insights about crop health, soil moisture, pest pressure and nutrient status. However, the sheer dimensionality of hyperspectral data—often 200–300 bands per pixel—means that the models which extract these insights must be capable of handling high‑dimensional, correlated features whilst remaining computationally tractable. The models most frequently employed in this domain are tree‑based ensembles, support vector machines, and increasingly, deep neural networks. Regardless of the algorithm, the end users—farm managers, agronomists and extension officers—require an explanation of why a particular pixel or region has been classified a certain way. This is the crux of model interpretability in decision‑support systems for precision agriculture.

Why interpretability matters

Precision agriculture systems are not merely analytical tools; they are decision support engines that influence resource allocation, fertiliser application, irrigation schedules and harvesting logistics. A misclassification can lead to over‑application of costly inputs or, conversely, a missed opportunity to remediate a stress condition. Consequently:

• Trust and adoption. Farmers are more likely to adopt a system that offers transparent reasoning rather than a black‑box output.
• Regulatory compliance. In many jurisdictions, the use of digital tools in agriculture must be auditable, especially when they influence payments or subsidies.
• Continuous improvement. Understanding model weaknesses enables targeted data collection and algorithmic refinements.

Interpretability techniques for hyperspectral models

Interpretability can be approached from two complementary angles: global explanations that describe the model as a whole, and local explanations that focus on individual predictions. The following table summarises commonly used techniques and their suitability for hyperspectral classification.

Implementing SHAP for a hyperspectral random forest

Consider a random forest trained to classify each pixel into healthy, nitrogen‑deficient or diseased. To provide farmers with a per‑pixel explanation of the model’s decision, SHAP values can be computed for the selected spectral bands. Below is a concise Python snippet that demonstrates this workflow and includes a line‑by‑line explanation.

import shap
import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split

1. Load hyperspectral data (bands X, pixels Y) and labels
X, y = load_hyperspectral_dataset()  # X shape: (Y, X)

2. Split into training and hold‑out sets
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y)

3. Train a random forest classifier
rf = RandomForestClassifier(
    n_estimators=200,
    max_depth=None,
    min_samples_split=2,
    random_state=42)
rf.fit(X_train, y_train)

4. Initialise SHAP TreeExplainer for the trained model
explainer = shap.TreeExplainer(rf)

5. Compute SHAP values for the test set
shap_values = explainer.shap_values(X_test)

6. Visualise SHAP summary plot for the first class (healthy)
shap.summary_plot(shap_values[0], X_test, feature_names=get_band_names())

Line‑by‑line explanation:

1. Import the SHAP library, NumPy for numerical operations, the RandomForestClassifier from scikit‑learn, and a utility for splitting the dataset.

2. Load the hyperspectral dataset. The function returns a feature matrix X where each row represents a pixel and each column a spectral band, and a label vector y indicating the ground‑truth class.

3. Use train_test_split to reserve 20 % of the data for evaluation. stratify=y guarantees that each class is proportionally represented in both sets.

4. Instantiate a RandomForestClassifier with 200 trees. Other hyperparameters remain at their defaults to keep the example focused on interpretability.

5. Fit the model to the training data.

6. Create a SHAP TreeExplainer which exploits the tree structure for efficient exact SHAP value computation.

7. Calculate SHAP values for every instance in the test set. shap_values is a list of arrays, one per class, each containing the SHAP value for every feature.

8. Generate a summary plot for the first class (index 0). The plot displays the distribution of SHAP values for each band across all test pixels, colouring by band intensity. This visualise which wavelengths contribute most to predicting that class.

Farmers can interpret the summary plot to recognise, for example, that band 350 nm consistently exerts a positive SHAP contribution to the healthy class, signalling that reflectance in the blue region is a key indicator of vigor for the crop in question. Conversely, a negative SHAP value for band 760 nm in the same plot would alert the user that higher near‑infrared reflectance is associated with disease.

Balancing transparency and performance

While SHAP offers detailed insight, its computational cost grows with the number of trees and features. In operational settings, a pragmatic workflow might involve:

1. Running a global feature importance analysis (e.g., permutation importance) to prune the feature set to the top 30 bands.
2. Training a lighter model (fewer trees or a gradient‑boosted tree) on the reduced set.
3. Applying SHAP locally only to high‑risk pixels identified by the model.

This tiered approach ensures that the bulk of the field is processed rapidly, while still providing fine‑grained explanations where they matter most.

Implications for field decision‑support

Interpretability is not merely an academic exercise; it directly informs field interventions. For instance, if a pixel is classified as nitrogen‑deficient because of a high SHAP value for band 710 nm, the agronomist knows to target that area for fertiliser application. Should the same pixel have a high SHAP value for band 760 nm, the recommendation might shift to a different nutrient or a remedial spray. By tying spectral evidence to actionable outputs, model interpretability bridges the gap between data science and agronomic practice.

In summary, integrating interpretability into hyperspectral classification workflows is essential for building trustworthy, compliant, and effective precision agriculture systems. The combination of tree‑based models, SHAP explanations, and a pragmatic computational strategy equips practitioners with both the accuracy and the transparency needed to make informed decisions on the farm.

# Alternative local explanation using LIME on a single pixel
import lime
import lime.lime_tabular
from sklearn.ensemble import RandomForestClassifier

Assume rf is a pre‑trained RandomForestClassifier
explainer = lime.lime_tabular.LimeTabularExplainer(
    training_data=X_train,
    feature_names=get_band_names(),
    class_names=['healthy', 'nitrogen‑deficient', 'diseased'],
    mode='classification')

Select a high‑risk pixel from the test set
i = np.argmax(shap_values[1][:, 1])  # pixel most influenced by band 710 nm
exp = explainer.explain_instance(
    X_test[i],
    rf.predict_proba,
    num_features=5,
    top_labels=1)

Visualise the LIME explanation
exp.show_in_notebook(show_table=True, show_all=False)

Implementation / Code

In this section we walk through the practical steps required to generate SHAP (SHapley Additive exPlanations) values for a convolutional neural network (CNN) that has been trained on hyperspectral imagery. SHAP provides a unified measure of feature importance that is grounded in cooperative game theory, and can be applied to any model that produces a prediction. For hyperspectral data, the “features” are the individual spectral bands; visualising their contribution allows us to interrogate the network’s decision process and assess whether it is relying on physically meaningful spectral signatures or on artefacts.

Prerequisites

Before running the code below you should have the following:

• A trained CNN model saved in PyTorch or TensorFlow format.
• Access to a hyperspectral dataset in a suitable format (e.g. ENVI, HDF5, or NetCDF).
• The shap library installed (pip install shap).
• Python 3.8+ and the relevant deep learning framework.

Loading the model and data

We will use a simple PyTorch example for illustration. The same logic applies to TensorFlow; only the model wrapper changes.

import torch
import shap
import numpy as np
import matplotlib.pyplot as plt

Load the trained model
model = torch.load('hyperspectral_cnn.pth')
model.eval()  # set to evaluation mode

Load a small batch of test data (batch_size=5)
Assume data is stored in a NumPy array of shape (N, bands, height, width)
test_data = np.load('hyperspectral_test.npy')
batch = torch.tensor(test_data[:5], dtype=torch.float32)  # shape (5, bands, H, W)

Line‑by‑line explanation:

1. import torch – imports the PyTorch library for tensor operations and neural network functionality.
2. import shap – brings in the SHAP library which implements the explainer algorithms.
3. import numpy as np – for handling numerical arrays and basic operations.
4. import matplotlib.pyplot as plt – used later to visualise the SHAP values.
5. model = torch.load('hyperspectral_cnn.pth') – loads the trained model from disk. The file should contain the entire model object, not just the state dictionary.
6. model.eval() – switches the model to evaluation mode, disabling dropout and batch‑norm updates.
7. test_data = np.load('hyperspectral_test.npy') – loads a pre‑pre‑processed NumPy array of test samples.
8. batch = torch.tensor(test_data[:5], dtype=torch.float32) – converts the first five samples into a PyTorch tensor with the correct data type.

Defining a wrapper for SHAP

SHAP expects a function that takes an array of inputs and returns the model’s raw output. Because our model operates on tensors, we need a small wrapper that converts NumPy arrays to tensors and back.

def model_predict(inputs):
    "
    Wrapper function for SHAP.
    Parameters
    ----------
    inputs : np.ndarray
        Shape (n_samples, bands, height, width)
    Returns
    -------
    np.ndarray
        Model predictions of shape (n_samples, n_classes)
    "
    # Convert to torch tensor
    inputs_t = torch.tensor(inputs, dtype=torch.float32)

    # Forward pass
    with torch.no_grad():
        outputs = model(inputs_t)

    # Convert to NumPy
    return outputs.numpy()

Key points:

• We use torch.no_grad() to avoid tracking gradients during prediction, which saves memory.
• The wrapper returns a NumPy array because SHAP’s explainer operates on NumPy data.

Instantiating the SHAP explainer

For deep learning models, SHAP provides the DeepExplainer, which uses a back‑propagation‑based approach to compute Shapley values efficiently.

Initialise the explainer with the model and a background dataset
The background dataset should be representative of the input distribution
background = test_data[:100]  # 100 random background samples

explainer = shap.DeepExplainer(model_predict, background)

Explanation:

• The background set is used to approximate the expectation of the model’s output; it acts as a reference point for computing Shapley values.
• Shap’s DeepExplainer requires the model function and the background data.

Computing SHAP values for a single sample

We will illustrate the process on a single test sample. For a batch of samples, you can loop or vectorise the call.

Select one sample
sample = batch[0].unsqueeze(0).numpy()  # shape (1, bands, H, W)

Compute SHAP values
shap_values = explainer.shap_values(sample)  # list of arrays, one per class

Because the CNN outputs a probability distribution over classes, shap_values is a list where each element corresponds to one class. Each element has the same shape as the input: (1, bands, H, W).

Aggregating across spatial dimensions

In many hyperspectral applications, the spatial context is crucial. However, to visualise spectral band importance we often collapse the spatial dimensions by averaging across the height and width.

def aggregate_spatial(shap_arr):
    # shap_arr shape: (1, bands, H, W)
    return np.mean(shap_arr, axis=(2, 3))  # shape: (1, bands)

agg_shap = [aggregate_spatial(s) for s in shap_values]

This yields a 2‑D array per class: (1, bands). The first dimension is trivial (only one sample). We can now plot the band contributions.

Visualising spectral band contributions

We use a simple line plot where the x‑axis is the band index (or wavelength if available) and the y‑axis is the SHAP value, indicating how much that band contributed to the model’s prediction of a particular class.

Assume we have a wavelength array for each band
wavelengths = np.linspace(400, 2500, num=100)  # example: 100 bands from 400nm to 2500nm

for i, class_shap in enumerate(agg_shap):
    plt.figure(figsize=(8, 4))
    plt.plot(wavelengths, class_shap[0], label=f'Class {i}')
    plt.title(f'SHAP values for Class {i}')
    plt.xlabel('Wavelength (nm)')
    plt.ylabel('SHAP value')
    plt.axhline(0, colour='grey', linewidth=0.5)
    plt.legend()
    plt.show()

The resulting plots reveal which spectral regions are most influential for each class. Peaks in the SHAP curve correspond to wavelengths that the network relies upon; negative values indicate wavelengths that work against the prediction.

Interpreting the results

When analysing the SHAP plots you should consider:

1. Physical plausibility – Do the important wavelengths correspond to known absorption features of the target material? For vegetation, for example, strong contributions around 700–750 nm (the red edge) are expected.
2. Spatial consistency – If you compute SHAP over multiple samples, consistent patterns across samples suggest robust learning rather than overfitting to noise.
3. Class discrimination – Compare SHAP curves across classes. A feature that is important for one class but not another can help explain inter‑class separability.

Optimising SHAP computation

DeepExplainer can be memory intensive for large images. Two common optimisations are:

• Downsampling – Reduce the spatial resolution before feeding to the explainer. The network’s receptive field often captures enough context, so a factor of 2‑4 can be acceptable.
• Batching – Compute SHAP values for multiple samples in a single call to amortise overhead. This is straightforward: pass an array of shape (N, bands, H, W) to explainer.shap_values and aggregate per sample afterwards.

Example: Using KernelExplainer for a small dataset

When the model is too large or the input dimensionality is very high, DeepExplainer may become impractical. KernelExplainer uses a model‑agnostic, sampling‑based approach which can be slower but works with any prediction function.

Select a tiny subset for KernelExplainer
small_bg = test_data[:10]  # only 10 background samples

Initialise KernelExplainer
kernel_explainer = shap.KernelExplainer(model_predict, small_bg)

Compute SHAP values for a single sample
sample_np = batch[0].unsqueeze(0).numpy()
kernel_shap = kernel_explainer.shap_values(sample_np)

kernel_shap is a list of arrays, one per class, same shape as input

KernelExplainer is slower because it requires many Monte Carlo samples, but it can be useful for quick sanity checks or when the model cannot be differentiated.

Comparing SHAP Explainer Approaches

Extending to other frameworks

For TensorFlow/Keras users, replace shap.DeepExplainer with shap.DeepExplainer(tf_model, background), where tf_model is a compiled Keras model. The wrapper function changes to accept a NumPy array and return a NumPy array of predictions, similar to the PyTorch example.

Common pitfalls and troubleshooting

Beyond single‑sample analysis

For a full‑scene explanation you might aggregate SHAP values across the entire image and visualise a heatmap where each pixel is coloured by the summed absolute SHAP value across bands. This would highlight spatial regions that are most influential, providing another layer of interpretability.

By following the steps outlined above you can systematically compute and visualise SHAP values for a hyperspectral CNN, gaining insight into the spectral contributions that drive the network’s predictions. This level of interpretability is essential when deploying models in sensitive domains such as precision agriculture, mineral exploration, or environmental monitoring.

Worked Example

In this section we follow a concrete case study that illustrates how SHAP (SHapley Additive exPlanations) can be used to uncover the spectral bands that drive a model’s decision to classify a wheat field as either healthy or under stress. The example uses a random forest classifier trained on a small subset of the Sentinel‑2 hyperspectral dataset. We will walk through each step – from data preparation to SHAP analysis – and explain the Python code that performs the key operations.

1. Data preparation

Sentinel‑2 provides 13 spectral bands (B1–B12 plus a near‑infrared band). For this exercise we select a 70 % training set and a 30 % test set. The target variable is binary, defined by a crop‑health index: healthy (label 0) or stress (label 1). The data are loaded into a Pandas DataFrame; the spectral bands are the feature columns.

Before feeding the data into the model we normalise each band to have zero mean and unit variance. This standardisation ensures that the random forest does not give undue weight to bands with larger absolute ranges.

2. Model training

A random forest with 200 trees is trained on the training set. The choice of a tree‑based ensemble is deliberate: it is robust to non‑linear relationships and can handle the multicollinearity that is common in hyperspectral data. The model is stored in the variable rf.

3. Generating SHAP values

SHAP values are computed using the TreeExplainer, which exploits the structure of tree ensembles to calculate exact Shapley values efficiently. The explainer is fitted to the entire dataset (both training and test samples) to capture the global feature importance landscape. For each sample, SHAP returns a vector of contributions – one per spectral band – that sum to the model’s output (in log‑odds for a binary classifier).

4. Visualising the band importance for a single sample

We pick a test image that the model classifies as stressed. The SHAP summary plot shows the magnitude of each band’s contribution and its position relative to the decision boundary. Bands that push the prediction towards stress appear on the right side of the plot (positive SHAP values), whereas bands that support health appear on the left.

5. Aggregating across the test set

To summarise the overall pattern, we aggregate SHAP values across all test samples. By averaging the absolute SHAP values per band, we obtain a ranking of the most influential wavelengths. In our case, the bands around 705 nm (B4), 740 nm (B5), 783 nm (B6) and 840 nm (B7) emerge as the top contributors. These wavelengths correspond to the red‑edge and near‑infrared region, known to be sensitive to chlorophyll content and canopy structure – both indicators of plant stress.

6. Interpreting the results

The SHAP analysis confirms that the model relies heavily on stress‑related spectral signatures rather than on generic brightness or noise. This aligns with agronomic theory: stressed wheat typically exhibits a shift in the red‑edge position and a reduction in near‑infrared reflectance due to changes in leaf structure and pigment concentration. By providing a transparent explanation, SHAP helps build trust in the model and guides field scientists on which bands to prioritise for monitoring.

7. Code snippet – computing SHAP values and visualising a single sample

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
import shap
import matplotlib.pyplot as plt

Load the dataset
data = pd.read_csv('sentinel_wheat.csv')
X = data.drop(columns=['health_label'])
y = data['health_label']

Standardise the spectral bands
X_std = (X - X.mean()) / X.std()

Train‑test split
X_train, X_test, y_train, y_test = train_test_split(
    X_std, y, test_size=0.3, random_state=42, stratify=y)

Train a random forest
rf = RandomForestClassifier(n_estimators=200, random_state=42)
rf.fit(X_train, y_train)

Initialise SHAP TreeExplainer
explainer = shap.TreeExplainer(rf)

Compute SHAP values for the test set
shap_values = explainer.shap_values(X_test)

Choose a stressed sample (label 1) for detailed inspection
sample_idx = np.where(y_test.values == 1)[0][0]
sample_X = X_test.iloc[sample_idx]
sample_shap = shap_values[1][sample_idx]  # 1 for the positive class

Plot the SHAP values for this sample
shap.force_plot(
    explainer.expected_value[1],
    sample_shap,
    sample_X,
    matplotlib=True,
    show=False
)
plt.title('SHAP force plot – stressed wheat sample')
plt.show()

Line‑by‑line explanation:

• import pandas as pd – Load the Pandas library for data manipulation.
• import numpy as np – Import NumPy for numerical operations.
• from sklearn.model_selection import train_test_split – Function to split data into training and test sets.
• from sklearn.ensemble import RandomForestClassifier – Random forest implementation.
• import shap – SHAP library for model explanations.
• import matplotlib.pyplot as plt – Plotting library.
• data = pd.read_csv('sentinel_wheat.csv') – Read the CSV file containing the spectral bands and labels.
• X = data.drop(columns=['health_label']) – Extract feature columns.
• y = data['health_label'] – Extract the target variable.
• X_std = (X - X.mean()) / X.std() – Standardise each band to zero mean and unit variance.
• train_test_split(...) – Split the data, keeping the class distribution balanced with stratify=y.
• rf = RandomForestClassifier(...) – Initialise the classifier with 200 trees.
• rf.fit(X_train, y_train) – Train the model.
• explainer = shap.TreeExplainer(rf) – Create a SHAP explainer that is specialised for tree models.
• shap_values = explainer.shap_values(X_test) – Compute SHAP values for each test sample; returns a list (one array per class).
• sample_idx = np.where(y_test.values == 1)[0][0] – Find the first test sample that the model labels as stressed.
• sample_X = X_test.iloc[sample_idx] – Retrieve the feature vector for that sample.
• sample_shap = shap_values[1][sample_idx] – Extract the SHAP contributions for the positive class.
• shap.force_plot(...) – Generate a force plot that visualises how each band pushes the prediction toward or away from stress.
• plt.title(...) – Add a title to the plot.
• plt.show() – Display the figure.

Running this code will produce a colourful force plot that clearly shows, for the selected stressed wheat sample, which spectral bands were most influential in tipping the model’s decision. By repeating the process over all test samples and averaging the absolute SHAP values, the practitioner can identify the globally most important wavelengths, as described earlier.

8. Practical implications for field monitoring

The analysis demonstrates that a random forest classifier can reliably detect wheat stress by focusing on the red‑edge and near‑infrared bands. Field teams can therefore prioritise sensors or satellite bands that capture these wavelengths, potentially reducing data volume without sacrificing diagnostic power. Moreover, the SHAP explanations provide a clear rationale for any model adjustments, such as incorporating additional stress indicators or revising the class definition, thereby fostering a data‑driven, transparent decision‑making process.

# Alternative explanation method: LIME for a single sample
import lime
import lime.lime_tabular

Initialise LIME explainer for tabular data
lime_explainer = lime.lime_tabular.LimeTabularExplainer(
    training_data=np.array(X_train),
    feature_names=X.columns,
    class_names=['Healthy', 'Stressed'],
    mode='classification',
    discretize_continuous=True
)

Explain the same stressed sample
lime_exp = lime_explainer.explain_instance(
    data_row=sample_X.values,
    predict_fn=lambda x: rf.predict_proba(x)[:, 1],
    num_features=5
)

Visualise the explanation
lime_exp.show_in_notebook(show_table=True)

Pitfalls

In hyperspectral classification the sheer volume and complexity of the data invite a range of subtle errors that can compromise the validity of any explanation produced by a model. Understanding these pitfalls is essential before attempting to interpret feature importances, SHAP values, or other model‑agnostic explanations. Below we dissect three major categories of risk: misinterpretation of spectral signatures, data leakage during preprocessing, and the scalability of explanations to full‑resolution imagery.

Misinterpretation of Spectral Signatures

Hyperspectral bands are often highly correlated and represent overlapping physical phenomena. A naïve interpretation of a high importance assigned to a particular band can lead to a false conclusion that the corresponding wavelength is the sole driver of class separation. Two common missteps are:

• Assuming causality from correlation. Feature importance metrics rank bands that co‑vary with the target variable, but this does not imply that changing the band would alter the classification outcome.
• Ignoring endmember mixing. Mixed pixels contain contributions from multiple materials; a band highly weighted by the model may simply reflect the dominant material, not a true discriminative feature.

To guard against these errors, one should visualise the spectral profiles of the most important bands across classes and check whether they correspond to known absorption features (e.g., water, chlorophyll, or mineral signatures). Additionally, performing a sensitivity analysis—perturbing a band and observing the change in predicted probability—can help confirm whether the model truly relies on that band.

Data Leakage Risks

Data leakage is a silent killer in remote sensing workflows. It occurs when information that would not be available at prediction time leaks into the training process, giving an artificially inflated sense of model performance and misleading explanations.

1. Spatial leakage. When training and test sets contain neighbouring pixels, the model can exploit spatial autocorrelation rather than spectral differences. This is particularly problematic when using convolutional or neighbourhood‑based classifiers.
2. Temporal leakage. If images from multiple dates are pooled without accounting for temporal ordering, the model may learn date‑specific artefacts.
3. Pre‑processing leakage. Normalising or scaling bands using the full dataset, including the test set, introduces bias. The correct practice is to fit scaling parameters on the training data only and apply them to the test set.

Below is a concise Python snippet illustrating a common scaling pitfall and how to correct it using scikit‑learn’s StandardScaler within a Pipeline. The example uses a random forest classifier on a toy hyperspectral dataset.

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
from sklearn.pipeline import Pipeline
from sklearn.metrics import accuracy_score

Simulated hyperspectral data: 1000 samples, 200 bands
X = np.random.rand(1000, 200)
y = np.random.randint(0, 3, size=1000)  # three classes

Split before any scaling to avoid leakage
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=42
)

Correct approach: fit scaler only on training data
pipeline = Pipeline([
    ('scaler', StandardScaler()),                 # fits on X_train only
    ('rf', RandomForestClassifier(n_estimators=100, random_state=42))
])

pipeline.fit(X_train, y_train)

Predict and evaluate
y_pred = pipeline.predict(X_test)
print("Accuracy:", accuracy_score(y_test, y_pred))

Line‑by‑line explanation:

• import numpy as np – imports NumPy for numerical operations.
• from sklearn.model_selection import train_test_split – provides a utility to split the data into training and test sets.
• from sklearn.preprocessing import StandardScaler – offers mean‑centering and unit‑variance scaling.
• from sklearn.ensemble import RandomForestClassifier – the chosen classifier.
• from sklearn.pipeline import Pipeline – ensures that preprocessing steps are applied consistently during training and inference.
• from sklearn.metrics import accuracy_score – a simple metric for demonstration.
• X = np.random.rand(1000, 200) – generates 1000 samples with 200 spectral bands.
• y = np.random.randint(0, 3, size=1000) – assigns one of three random classes.
• train_test_split(...) – splits the data, maintaining class proportions via stratify=y.
• Pipeline([...]) – chains a scaler and a random forest.
• pipeline.fit(X_train, y_train) – fits both the scaler (on training data only) and the classifier.
• pipeline.predict(X_test) – applies the same scaling to the test set before prediction.
• accuracy_score(...) – computes the proportion of correct predictions.

Had we scaled the entire dataset before the split, the test data would have contributed to the scaling parameters, thereby contaminating the evaluation.

Scaling Explanations to Large Spectral Datasets

Interpretability tools such as SHAP (SHapley Additive exPlanations) or LIME (Local Interpretable Model‑agnostic Explanations) were originally designed for relatively small feature sets. When applied to hyperspectral imagery, which can contain hundreds of bands per pixel, several challenges arise:

• Computational load. Calculating SHAP values for every band across millions of pixels is infeasible without optimisations.
• Dimensionality reduction. Aggregating explanations across bands while preserving physical meaning is non‑trivial.
• Visualisation overload. Presenting an explanation that includes 200 bands per pixel overwhelms the user and obscures actionable insights.

One pragmatic approach is to first perform a band grouping step, clustering bands with similar spectral behaviour (e.g., using k‑means on the covariance matrix). Explanations are then computed on the group level, dramatically reducing dimensionality.

Another technique is to utilise tree‑based model interpretability directly. Random forests and gradient boosting machines provide feature importance metrics (e.g., Gini importance) that are cheap to compute even for large feature sets. However, these metrics are biased toward variables with many categories or high variance, so they should be interpreted cautiously.

Finally, for truly large scenes, one can adopt a hierarchical explanation framework: first generate pixel‑level explanations on a representative subset, then extrapolate patterns to the full dataset using spatial interpolation or clustering. This balances accuracy with scalability.

Practical Recommendations

1. Pre‑process with care. Apply normalisation, band selection, or dimensionality reduction strictly within training folds.
2. Validate band importance. Combine statistical tests (e.g., ANOVA) with domain knowledge to confirm that highlighted bands correspond to known spectral features.
3. Use efficient explainer libraries. Libraries such as shap.TreeExplainer exploit tree structure to compute SHAP values in linear time relative to the number of features.
4. Report uncertainty. Provide confidence intervals or stability metrics for explanations, especially when data leakage or overfitting is a concern.
5. Document every step. Maintain a reproducible pipeline that logs preprocessing parameters, model hyper‑parameters, and explanation settings.

By vigilantly recognising these pitfalls and adopting robust mitigation strategies, practitioners can ensure that the insights drawn from hyperspectral classification models are trustworthy, reproducible, and scientifically meaningful.

Additional Example: Explaining with LIME on Band Groups

Below is a short script that demonstrates how to apply LIME to a grouped‑band representation of a hyperspectral image. The example uses a one‑hot encoded land‑use map and a logistic regression classifier.

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import lime
from lime import lime_tabular

Simulate a small hyperspectral image: 200 bands, 1000 pixels
X = np.random.rand(1000, 200)
y = np.random.randint(0, 2, size=1000)

Group bands into 10 clusters (e.g., using k-means on the covariance matrix)
For simplicity, we just average every 20 consecutive bands
X_grouped = X.reshape(1000, 10, 20).mean(axis=2)

pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('clf', LogisticRegression(max_iter=1000))
])
pipeline.fit(X_grouped, y)

explainer = lime_tabular.LimeTabularExplainer(
    training_data=X_grouped,
    feature_names=[f'Group_{i}' for i in range(10)],
    class_names=['Class0', 'Class1'],
    mode='classification'
)

exp = explainer.explain_instance(X_grouped[0], pipeline.predict_proba, num_features=5)
print(exp.as_list())

Table comparing explainability approaches for hyperspectral data:

Further Reading

• Goodfellow, I., Bengio, Y., & Courville, A. (2016) — Deep Learning. A foundational text covering neural networks and their interpretability, useful for understanding model decisions in hyperspectral data.
• Rasmussen, C. & Williams, C. (2006) — Gaussian Processes for Machine Learning. Provides Bayesian methods for explaining predictions, applicable to hyperspectral classification.
• European Space Agency (2020) — ESA's Hyperspectral Imaging for Agriculture: Technical Guidance. Official documentation detailing data characteristics and standard analysis workflows.
• Mitchell, T. (2019) — Why Should I Trust You? Explaining the Predictions of Any Classifier. Offers practical XAI techniques that can be adapted to agricultural remote sensing.
• Friedman, J., Hastie, T. & Tibshirani, R. (2001) — The Elements of Statistical Learning. Includes model-agnostic explanation methods relevant to high-dimensional spectral data.