Lesson 4.1.3: Training, validation and testing data preparation

flowchart TD
    A[Raw Data] --> B[Data Cleaning]
    B --> C[Feature Engineering]
    C --> D[Split Dataset]
    D --> E[Training Set]
    D --> F[Validation Set]
    D --> G[Test Set]
    E --> H[Model Training]
    F --> I[Hyper‑parameter Tuning]
    I --> J[Model Evaluation]
    G --> J

Conceptual Foundations of Data Splitting

When building predictive models for hyperspectral imagery, the way we partition the available data into training, validation and test sets is not merely a procedural formality – it is a statistical design choice that determines the reliability of our performance estimates and the generalisability of the final model. Hyperspectral datasets are characteristically high‑dimensional, highly correlated across spectral bands, and often exhibit spatial autocorrelation. These properties necessitate a careful approach to data splitting that guards against over‑optimistic bias, ensures representative sampling of the underlying spectral variability, and preserves the spatial integrity of the scenes.

Why split at all?

In supervised learning we optimise model parameters on a subset of data (the training set) and then evaluate the model on unseen data to estimate its performance. If the same data used for optimisation is also used for performance assessment, the estimate becomes optimistically biased: the model is effectively tested on data it has already “seen.” By reserving a separate set for validation, we can tune hyper‑parameters and select models in a way that mimics real‑world deployment, where predictions are made on data that has not been part of the training process.

Moreover, the test set should be completely untouched until the final evaluation. It acts as a safeguard against accidental leakage of information during model selection and provides an unbiased estimate of how the model will perform in practice. In hyperspectral remote sensing, where the cost of misclassification can be high (e.g. in mineral exploration or precision agriculture), this unbiased estimate is essential for risk assessment.

Statistical considerations in splitting

In conventional tabular data, random stratified sampling is often sufficient: we preserve the class distribution in each split to avoid class imbalance issues. However, hyperspectral imagery introduces two additional layers of complexity:

1. Spatial autocorrelation – Pixels that are geographically close tend to be spectrally similar. Randomly sampling individual pixels can lead to training and validation sets that are not truly independent because neighbouring pixels are represented in both sets. This inflates performance estimates.
2. Spectral variability – Within a single land cover type, the spectral signature can vary substantially due to illumination, atmospheric effects, or sensor noise. A split that inadvertently groups similar spectra together in the training set and dissimilar spectra in the validation set can produce misleading results.

To mitigate these issues, several strategies are commonly adopted:

• Block sampling – Instead of pixel‑level sampling, we sample contiguous spatial blocks (e.g. 10 × 10 px tiles). Each block is assigned to one of the splits, ensuring that neighbouring pixels are all contained within the same set. This preserves the spatial independence assumption.
• Stratified block sampling – We combine block sampling with stratification by class. Within each class, we partition the available blocks into training, validation and test groups, maintaining the overall class proportions.
• Leave‑one‑scene‑out – For multi‑scene datasets, we hold out entire scenes (or scenes from a specific location) as the test set. This tests the model’s ability to generalise across different geographic areas.
• Cross‑validation with spatial folds – K‑fold cross‑validation can be adapted to spatial data by defining folds that respect spatial boundaries.

Mathematical rationale for block sampling

Consider a hyperspectral image as a matrix \(X \in \mathbb{R}^{n imes p}\), where \(n\) is the number of pixels and \(p\) the number of spectral bands. We can view the pixel locations as points in a two‑dimensional grid. The spatial autocorrelation function \(C(h)\) describes the expected correlation between two pixels separated by a lag \(h\). If \(C(h)\) decays slowly, neighbouring pixels are highly correlated. By selecting blocks of size \(b imes b\), we effectively sample from a set of independent units if the block size exceeds the spatial correlation length. The probability that a pixel in the training set is correlated with a pixel in the validation set is then reduced to \(O(b^{-2})\), assuming a simple spatial decay model.

In practice, we estimate the empirical autocorrelation by computing the variogram of the spectral bands or by inspecting the spatial distribution of class labels. Once a suitable block size is chosen, we partition the image into non‑overlapping blocks, assign each block to a split, and then flatten the blocks back into pixel indices for model training.

Illustrative example: block sampling in Python

The following snippet demonstrates how to perform stratified block sampling on a hyperspectral raster stored as a NumPy array. The array shape is (rows, cols, bands). We also assume a corresponding label array of shape (rows, cols).

import numpy as np
from sklearn.model_selection import StratifiedShuffleSplit

Hyperspectral data and labels
hs_data = np.load('hyperspectral.npy')      # shape: (rows, cols, bands)
labels  = np.load('labels.npy')             # shape: (rows, cols)

Define block size (e.g. 20x20 pixels)
block_size = 20
rows, cols, _ = hs_data.shape
n_blocks_row = rows // block_size
n_blocks_col = cols // block_size

Create block indices and associated labels
block_indices = []
block_labels  = []

for i in range(n_blocks_row):
    for j in range(n_blocks_col):
        # Extract block slice
        r_start, r_end = i  block_size, (i + 1)  block_size
        c_start, c_end = j  block_size, (j + 1)  block_size
        block = hs_data[r_start:r_end, c_start:c_end, :]
        block_lab = labels[r_start:r_end, c_start:c_end]

        # Flatten block to a single sample
        block_flat = block.reshape(-1, block.shape[2])  # (block_size^2, bands)
        block_indices.append(block_flat)

        # Assign block label as the majority class within the block
        majority_label = np.bincount(block_lab.flatten()).argmax()
        block_labels.append(majority_label)

block_indices = np.array(block_indices)   # shape: (n_blocks, block_size^2, bands)
block_labels  = np.array(block_labels)    # shape: (n_blocks,)

Stratified split at block level
split = StratifiedShuffleSplit(n_splits=1, test_size=0.15, random_state=42)
train_idx, test_idx = next(split.split(block_indices, block_labels))

Further split training blocks into training and validation
split_val = StratifiedShuffleSplit(n_splits=1, test_size=0.15, random_state=42)
train_idx, val_idx = next(split_val.split(block_indices[train_idx], block_labels[train_idx]))

Assemble final datasets
train_blocks = block_indices[train_idx]
val_blocks   = block_indices[val_idx]
test_blocks  = block_indices[test_idx]

Optionally, reshape blocks back to pixel-wise arrays for model input
train_data = np.concatenate(train_blocks, axis=0)  # shape: (n_train_pixels, bands)
val_data   = np.concatenate(val_blocks, axis=0)
test_data  = np.concatenate(test_blocks, axis=0)

train_labels = block_labels[train_idx]
val_labels   = block_labels[val_idx]
test_labels  = block_labels[test_idx]

Line‑by‑line explanation:

1. Import NumPy and scikit‑learn’s StratifiedShuffleSplit for balanced sampling.
2. Load the hyperspectral data and the associated class labels.
3. Define the desired block size. The block size should be larger than the spatial correlation length to reduce leakage.
4. Compute how many blocks will fit along rows and columns.
5. Iterate over the grid, extracting each block and its corresponding label patch.
6. Flatten each block into a two‑dimensional array where each row is a pixel spectrum; this facilitates later stacking.
7. Determine the majority class within the block to assign a single label to the entire block.
8. Convert the list of block samples and labels into NumPy arrays.
9. Use StratifiedShuffleSplit to split blocks into training and test sets while preserving class proportions.
10. Repeat the split on the training blocks to carve out a validation set.
11. Concatenate the flattened blocks back into a single pixel‑wise dataset for each split.
12. Extract the corresponding block labels for each split.

Alternative approach: spatial cross‑validation in Python

When datasets contain multiple scenes or large spatial extents, a leave‑one‑scene‑out strategy can be more realistic. The following example shows how to implement spatial cross‑validation using scikit‑learn’s GroupKFold where each scene (or a pre‑defined spatial group) is treated as a single unit.

import numpy as np
from sklearn.model_selection import GroupKFold
from sklearn.metrics import accuracy_score

Load data and group identifiers (e.g. scene IDs)
hs_data = np.load('hyperspectral.npy')      # (rows, cols, bands)
labels  = np.load('labels.npy')             # (rows, cols)
groups  = np.load('scene_ids.npy')          # (rows, cols) – integer per scene

Flatten arrays for model input
X = hs_data.reshape(-1, hs_data.shape[2])   # (n_pixels, bands)
y = labels.flatten()                        # (n_pixels,)
g = groups.flatten()                        # (n_pixels,)

Define GroupKFold
gkf = GroupKFold(n_splits=5)

for fold, (train_idx, val_idx) in enumerate(gkf.split(X, y, groups=g)):
    X_train, X_val = X[train_idx], X[val_idx]
    y_train, y_val = y[train_idx], y[val_idx]

    # Train a simple classifier (e.g. RandomForest)
    from sklearn.ensemble import RandomForestClassifier
    clf = RandomForestClassifier(n_estimators=200, random_state=42)
    clf.fit(X_train, y_train)

    # Evaluate
    y_pred = clf.predict(X_val)
    acc = accuracy_score(y_val, y_pred)
    print(f'Fold {fold+1} accuracy: {acc:.4f}')

This approach ensures that all pixels from a given scene are kept together in either the training or validation set, thereby preventing spatial leakage across scenes.

Evaluation metrics and the role of the test set

After training a model on the training set and tuning hyper‑parameters on the validation set, we obtain a final model. The only honest performance estimate we can report is the one obtained on the test set. Common metrics in hyperspectral classification include overall accuracy, Cohen’s kappa, per‑class recall and precision, and the area under the receiver operating characteristic curve (AUC). Each of these metrics can be computed on the test set in the same way they were calculated on the validation set, providing a transparent comparison.

Because the test set was never used during model development, its performance reflects the model’s true generalisation ability. If the test accuracy is markedly lower than the validation accuracy, it may indicate that the validation set was too easy (e.g. due to residual spatial leakage) or that the model is over‑fitting to the training data. In such cases, revisiting the data splitting strategy (e.g. increasing block size or adopting leave‑one‑scene‑out validation) is warranted.

Practical tips for hyperspectral datasets

• Document block sizes and split ratios. Future reviewers should be able to reproduce the exact partitions.
• Check for class imbalance. If a class is under‑represented, consider oversampling its blocks or using focal loss during training.
• Preserve spectral fidelity. When normalising spectra, do so consistently across all splits to avoid introducing bias.
• Use random seeds. Fixed seeds ensure that the same split can be recreated, which is essential for reproducibility.
• Visualise the splits. Overlay training, validation and test blocks on the original image to inspect spatial coverage.

Comparison of splitting strategies

Conclusion

Data partitioning in hyperspectral remote sensing is a nuanced exercise that balances statistical rigour with practical constraints. By adopting block‑based, stratified sampling and carefully preserving spatial independence, we safeguard against over‑optimistic performance estimates and build models that truly generalise to new scenes. The code examples above provide concrete starting points; however, the specific parameters (block size, split ratios, stratification criteria) should be tuned to the characteristics of each dataset and the objectives of the analysis.

Practical Implementation and Code

In this section we translate the theoretical steps of training‑validation‑testing data preparation into executable Python code. The workflow emphasises reproducibility, stratified sampling to respect class distributions, and data augmentation to enrich the training set while preserving the integrity of the validation and test partitions. Throughout, we utilise the scikit‑learn library for splitting, encoding and scaling, and the spectral package for specialised spectral‑image transformations. The code snippets are deliberately compact yet fully commented so that each line’s purpose is crystal clear.

1. Loading and Inspecting the Dataset

Before any manipulation we load the data into a pandas DataFrame and perform a quick sanity check. It is vital to confirm that the target variable contains discrete labels and that no missing values are present.

import pandas as pd
from pathlib import Path

Path to the CSV file containing the spectral measurements
data_path = Path("data") / "spectral_dataset.csv"

Load the dataset
df = pd.read_csv(data_path)

Basic sanity checks
print(f"Dataset shape: {df.shape}")
print("Target class distribution:
", df['label'].value_counts())
print("Missing values per column:
", df.isnull().sum())

Explanation: We use Path for robust file path handling, pd.read_csv to read the tabular data, and value_counts to visualise class imbalance. The isnull().sum() call guarantees that downstream procedures will not fail due to missing entries.

2. Stratified Train‑Validation‑Test Split

To ensure that each split mirrors the overall class distribution, we perform a two‑stage stratified split: first separating the test set, then splitting the remaining data into training and validation sets. The train_test_split function from scikit‑learn accepts a stratify parameter that enforces this property.

from sklearn.model_selection import train_test_split

Separate features and target
X = df.drop(columns=['label'])
y = df['label']

1. Reserve 20% of the data for the final test set
X_temp, X_test, y_temp, y_test = train_test_split(
    X, y,
    test_size=0.20,
    random_state=42,
    stratify=y
)

2. Split the remaining 80% into training (70%) and validation (30%)
X_train, X_val, y_train, y_val = train_test_split(
    X_temp, y_temp,
    test_size=0.30,
    random_state=42,
    stratify=y_temp
)

print(f"Training set: {X_train.shape[0]} samples")
print(f"Validation set: {X_val.shape[0]} samples")
print(f"Test set: {X_test.shape[0]} samples")

Explanation: The first call allocates 20 % of the samples to the test set, guaranteeing that the test set remains untouched until model evaluation. The second call further splits the remaining data; the test_size=0.30 indicates that 30 % of the 80 % (i.e. 24 % of the full dataset) becomes the validation set. Both splits use a fixed random_state so that results are reproducible.

3. Feature Scaling with Standardisation

Neural networks and many other algorithms perform better when features are centred and scaled. StandardScaler computes the mean and standard deviation from the training data only, then applies the same transformation to validation and test sets to avoid data leakage.

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()

Fit on training data
X_train_scaled = scaler.fit_transform(X_train)

Apply the same transformation to validation and test data
X_val_scaled = scaler.transform(X_val)
X_test_scaled = scaler.transform(X_test)

Explanation: By fitting on the training set only, we preserve the statistical properties of unseen data. The transformed arrays are now ready for ingestion by downstream models.

4. One‑Hot Encoding of the Target Variable

For multi‑class classification, one‑hot encoding converts discrete class labels into binary vectors, which are compatible with categorical cross‑entropy loss functions.

from sklearn.preprocessing import OneHotEncoder
import numpy as np

encoder = OneHotEncoder(sparse_output=False)

Fit on training labels
y_train_enc = encoder.fit_transform(y_train.values.reshape(-1, 1))

Transform validation and test labels
y_val_enc = encoder.transform(y_val.values.reshape(-1, 1))
y_test_enc = encoder.transform(y_test.values.reshape(-1, 1))

Explanation: The reshape(-1, 1) call ensures a column vector shape, as required by the encoder. The sparse_output=False argument returns dense NumPy arrays, which are convenient for most frameworks such as TensorFlow or PyTorch.

5. Data Augmentation for Spectral Images

While tabular data can be augmented via simple transformations, spectral images benefit from domain‑specific operations such as spectral shifting, intensity scaling and random noise addition. The spectral library offers a convenient Augmenter class that can be composed of multiple Transform objects.

from spectral import Augmenter, transforms as T

Define a composition of augmentation transforms
augmenter = Augmenter(
    T.RandomSpectralShift(offset_range=(-0.05, 0.05)),
    T.RandomIntensityScaling(scale_range=(0.9, 1.1)),
    T.RandomGaussianNoise(noise_std=0.01)
)

Apply augmentations to the training set only
X_train_augmented = augmenter.apply(X_train_scaled)

Explanation: The RandomSpectralShift modifies the wavelength axis slightly, mimicking calibration drift. RandomIntensityScaling simulates varying illumination. Finally, RandomGaussianNoise adds sensor noise. Importantly, we apply these augmentations only to the training set; validation and test sets remain untouched to provide a fair assessment of model generalisation.

6. Packaging the Data into PyTorch Datasets

For deep learning workflows, it is common to wrap the data into torch.utils.data.Dataset objects. Below we construct a lightweight custom dataset that serves the scaled, augmented features and one‑hot labels.

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

class SpectralDataset(Dataset):
    def __init__(self, features, labels):
        self.features = torch.tensor(features, dtype=torch.float32)
        self.labels = torch.tensor(labels, dtype=torch.float32)

    def __len__(self):
        return self.features.shape[0]

    def __getitem__(self, idx):
        return self.features[idx], self.labels[idx]

Create datasets
train_dataset = SpectralDataset(X_train_augmented, y_train_enc)
val_dataset   = SpectralDataset(X_val_scaled,       y_val_enc)
test_dataset  = SpectralDataset(X_test_scaled,      y_test_enc)

DataLoaders for batching
batch_size = 64
train_loader = DataLoader(train_dataset, batch_size=batch_size, shuffle=True)
val_loader   = DataLoader(val_dataset,   batch_size=batch_size, shuffle=False)
test_loader  = DataLoader(test_dataset,  batch_size=batch_size, shuffle=False)

Explanation: Each dataset returns a tuple of feature tensor and label tensor. The training loader shuffles the data at every epoch, while validation and test loaders do not. This structure aligns with standard PyTorch training loops.

7. Summary Table of Data Splits

Replace the template placeholders with the actual numbers after execution.

8. Common Pitfalls and Best Practices

1. Data leakage: Never fit scalers, encoders or augmenters on the full dataset. Always fit exclusively on the training split.
2. Stratification consistency: Ensure that the target column used for stratification contains only the final, clean labels. Any preprocessing that modifies class frequencies must be applied after splitting.
3. Random state control: Use a fixed random_state for every stochastic operation (splitting, augmentations, shuffling) to guarantee reproducibility across runs.
4. Augmentation extent: Over‑aggressive augmentations can degrade performance. Validate the impact on a small validation subset before scaling up.
5. Memory considerations: Augmented datasets can be large. When memory is a constraint, consider on‑the‑fly augmentation within the __getitem__ method rather than precomputing the entire augmented set.

9. Extending to Distributed Training

When scaling to distributed environments, the data preparation pipeline should be serialisable and stateless. One approach is to serialise the fitted StandardScaler and OneHotEncoder objects using joblib.dump, then load them on each worker node. The SpectralDataset can be wrapped in a DistributedSampler to ensure each process receives a distinct portion of the data.

import joblib
from torch.utils.data.distributed import DistributedSampler

Serialise scalers and encoders
joblib.dump(scaler, "artifacts/scaler.joblib")
joblib.dump(encoder, "artifacts/encoder.joblib")

On each worker
scaler = joblib.load("artifacts/scaler.joblib")
encoder = joblib.load("artifacts/encoder.joblib")

Distributed sampler for training
train_sampler = DistributedSampler(train_dataset, num_replicas=args.world_size, rank=args.rank)
train_loader = DataLoader(train_dataset, batch_size=batch_size, sampler=train_sampler)

Explanation: The DistributedSampler guarantees that each process receives a unique, non‑overlapping subset of the training data, which is essential for efficient parallel training.

10. Final Remarks

The code presented above constitutes a robust, reproducible foundation for any spectral‑image classification project. By carefully separating concerns—splitting, scaling, encoding, augmenting, and batching—researchers can minimise pitfalls such as data leakage and class imbalance. The modular design also facilitates experimentation: swapping the augmentation pipeline, testing alternative scalers, or integrating new feature extraction methods can be achieved with minimal friction. Armed with this implementation, you are now ready to proceed to model definition, training loops and, ultimately, evaluation on the held‑out test set.

11. Additional Example: PCA Dimensionality Reduction

In some cases, reducing feature dimensionality before training can improve model efficiency and mitigate overfitting. The following snippet demonstrates how to apply Principal Component Analysis (PCA) to the scaled features, preserving 95 % of the variance.

from sklearn.decomposition import PCA

Fit PCA on training data
pca = PCA(n_components=0.95, random_state=42)
X_train_pca = pca.fit_transform(X_train_scaled)

Apply the same transformation to validation and test sets
X_val_pca = pca.transform(X_val_scaled)
X_test_pca = pca.transform(X_test_scaled)

print(f"PCA reduced dimensionality to {X_train_pca.shape[1]} components")

Explanation: By fitting PCA solely on the training data, we avoid leakage. The resulting lower‑dimensional representations can be fed into downstream models, often with a reduced computational burden.

12. Comparison Table: Feature Scaling Approaches

Worked Example: Urban Land‑Cover Classification

In this section we trace the entire pipeline that transforms raw hyperspectral imagery of an urban area into a tidy, model‑ready dataset. The example uses a publicly available hyperspectral cube collected over a metropolitan region, and demonstrates every step that a data scientist must perform: quality checks, georeferencing, spectral calibration, scene mosaic, cloud masking, ROI extraction, normalisation, and finally the split into training, validation and testing partitions. The goal is to give you a concrete, repeatable workflow that you can adapt to your own data.

1. Data Acquisition and Initial Inspection

We begin with the hyperspectral cube, a multi‑band raster where each band represents a narrow wavelength interval. The cube is stored in a NetCDF file named urban_cube.nc. The first thing to do is to inspect the file structure and metadata to confirm that the spatial resolution, spectral range and temporal information are as expected.

import xarray as xr

Load the NetCDF cube
cube = xr.open_dataset('urban_cube.nc')

Print global attributes and dimensions
print(cube.attrs)
print(cube.dims)

Line‑by‑line explanation:

• import xarray as xr – imports the Xarray library, which specialises in labelled multi‑dimensional arrays.
• cube = xr.open_dataset('urban_cube.nc') – opens the NetCDF file and creates an xarray.Dataset object.
• print(cube.attrs) – outputs global attributes such as sensor name, acquisition date and processing level.
• print(cube.dims) – shows the dimensions of the cube, typically (y, x, band) for spatial and spectral axes.

Checking the band coordinate confirms that the spectral axis spans 400 nm to 2500 nm with 200 bands, which is suitable for detailed surface material discrimination.

2. Radiometric Calibration and Spectral Normalisation

Raw sensor counts must be converted to reflectance. We use the sensor’s calibration constants stored in the NetCDF file. After calibration, we apply a simple continuum removal to reduce the influence of illumination variability.

# Assume calibration constants are stored in "calibration_factor"
cal_factor = cube['calibration_factor'].values

Convert digital numbers to reflectance
reflectance = cube['DN'] / cal_factor

Continuum removal: normalise each pixel spectrum
norm_reflectance = reflectance / reflectance.max(dim='band')

Explanation:

• cal_factor = cube['calibration_factor'].values – extracts the calibration factor array.
• reflectance = cube['DN'] / cal_factor – element‑wise division converts raw DN to reflectance.
• norm_reflectance = reflectance / reflectance.max(dim='band') – divides each pixel spectrum by its maximum value across bands, yielding a unitless, illumination‑independent representation.

3. Georeferencing and Masking

Hyperspectral cubes often lack a proper georeferencing grid. We align the cube to a reference orthophoto by applying the affine transform provided in the metadata. We also mask out the ocean and large water bodies that are irrelevant for urban land‑cover.

import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

Load the reference orthophoto
with rasterio.open('urban_ortho.tif') as ortho:
    transform, width, height = calculate_default_transform(
        ortho.crs, cube.attrs['crs'], ortho.width, ortho.height, *ortho.bounds)
    kwargs = ortho.meta.copy()
    kwargs.update({
        'crs': ortho.crs,
        'transform': transform,
        'width': width,
        'height': height
    })

Reproject the reflectance cube
reproj = xr.DataArray(
    data=np.empty((height, width, cube.sizes['band'])),
    dims=('y', 'x', 'band'),
    coords={
        'y': np.arange(height),
        'x': np.arange(width),
        'band': cube['band']
    }
)

for i in range(cube.sizes['band']):
    reproject(
        source=norm_reflectance.isel(band=i).values,
        destination=reproj.isel(band=i).data,
        src_transform=cube.attrs['transform'],
        src_crs=cube.attrs['crs'],
        dst_transform=transform,
        dst_crs=ortho.crs,
        resampling=Resampling.nearest
    )

After reprojection we generate a water mask from the orthophoto’s NDVI band: pixels with NDVI < 0.1 are considered water.

ndvi = (ortho.read(4) - ortho.read(3)) / (ortho.read(4) + ortho.read(3))
water_mask = ndvi < 0.1

We then apply this mask to the cube, setting water pixel spectra to NaN so they are excluded from later steps.

4. Region of Interest (ROI) Extraction

For supervised learning we need ground truth. Field crews have provided vector polygons for key urban classes: impervious surface, vegetation, water, building and road. We rasterise these polygons onto the cube’s grid and use the resulting labels to sample spectra.

import geopandas as gpd
from rasterio.features import rasterize

Load vector shapefile
shp = gpd.read_file('urban_labels.shp')

Create an empty label raster
label_raster = np.zeros((height, width), dtype=np.uint8)

Map class names to integer codes
class_map = {
    'impervious': 1,
    'vegetation': 2,
    'water': 3,
    'building': 4,
    'road': 5
}

Rasterise each class
for cls_name, cls_code in class_map.items():
    cls_polys = shp[shp['class'] == cls_name]
    shapes = ((geom, cls_code) for geom in cls_polys.geometry)
    rasterized = rasterize(
        shapes,
        out_shape=(height, width),
        transform=transform,
        fill=0,
        dtype=np.uint8
    )
    label_raster = np.maximum(label_raster, rasterized)

Apply water mask
label_raster[water_mask] = 0

Pixels labelled as 0 are ignored during sampling. The remaining labelled pixels form our candidate set for training.

5. Spectral Feature Engineering

While raw reflectance can be fed directly into a model, certain engineered features improve separability. We compute the following per‑pixel metrics:

• Mean reflectance in the visible band (400‑700 nm)
• Mean reflectance in the near‑infrared band (700‑1100 nm)
• Band ratio 800/700 nm (indicative of vegetation)
• Standard deviation across all bands (texture proxy)

# Define band indices for visible and NIR
vis_band_idx = slice(0, 50)   # 400‑700 nm
nir_band_idx = slice(50, 80)  # 700‑1100 nm

features = xr.Dataset({
    'vis_mean': norm_reflectance.isel(band=vis_band_idx).mean(dim='band'),
    'nir_mean': norm_reflectance.isel(band=nir_band_idx).mean(dim='band'),
    'ratio_800_700': norm_reflectance.isel(band=800) / norm_reflectance.isel(band=700),
    'spectral_sd': norm_reflectance.std(dim='band')
})

These engineered features are added as additional bands to the cube, resulting in a 4‑band feature cube ready for extraction.

6. Dataset Construction

We now pair each labelled pixel’s feature vector with its class label. To avoid bias from class imbalance we enforce a stratified sampling strategy: we randomly select an equal number of samples from each class, up to a maximum of 2000 per class.

import pandas as pd
import numpy as np

Flatten label raster and feature cube
labels_flat = label_raster.flatten()
features_flat = features.to_array().values.reshape(len(features), -1)

Create a DataFrame
df = pd.DataFrame(features_flat.T, columns=features.data_vars.keys())
df['label'] = labels_flat

Remove unlabeled pixels
df = df[df['label'] != 0]

Stratified sampling
sampled_df = df.groupby('label', group_keys=False).apply(
    lambda x: x.sample(n=min(2000, len(x)), random_state=42)
)

Shuffle
sampled_df = sampled_df.sample(frac=1, random_state=42).reset_index(drop=True)

7. Train‑Validate‑Test Split

We adopt an 70/15/15 split, ensuring that the spatial distribution of samples is preserved. To achieve this, we perform a spatial stratification: we first divide the study area into 5 equal‑area tiles, then split within each tile.

from sklearn.model_selection import train_test_split

Assign each pixel to a spatial tile
tile_size = 100  # pixels
tiles = (np.arange(height) // tile_size)[:, None] + (np.arange(width) // tile_size)

Flatten tile indices
tile_flat = tiles.flatten()

Merge tile info into DataFrame
sampled_df['tile'] = tile_flat[sampled_df.index]

Function to split within tiles
def split_by_tile(df, test_size=0.15, val_size=0.15):
    train_idx, temp_idx = [], []
    for tile in df['tile'].unique():
        tile_df = df[df['tile'] == tile]
        t, v = train_test_split(tile_df.index, test_size=test_size/(1-test_size), random_state=42)
        train_idx.extend(t)
        temp_idx.extend(v)
    val_idx, test_idx = train_test_split(temp_idx, test_size=val_size/(val_size+0.15), random_state=42)
    return df.loc[train_idx], df.loc[val_idx], df.loc[test_idx]

train_df, val_df, test_df = split_by_tile(sampled_df)

The three resulting DataFrames contain the feature vectors and labels for training, validation and testing respectively. They are now ready to be fed into a machine learning pipeline.

8. Persisting the Prepared Dataset

For reproducibility we store the splits as compressed Parquet files. Parquet preserves the columnar structure and is efficient for large tabular data.

train_df.to_parquet('urban_train.parquet', compression='snappy')
val_df.to_parquet('urban_val.parquet', compression='snappy')
test_df.to_parquet('urban_test.parquet', compression='snappy')

At this stage the dataset is fully prepared: every pixel in the training set is a 4‑dimensional feature vector with an associated class label; the validation and test sets are spatially balanced and disjoint. The next lesson will cover how to ingest this data into a deep learning model and evaluate its performance.

# Alternative approach: use a dedicated library for labelled dataset creation
Here we use the scikit‑learn library to create a custom splitter that preserves
spatial contiguity without manual tile assignment.

from sklearn.model_selection import GroupShuffleSplit

Define groups based on spatial coordinates (e.g., row index)
groups = sampled_df.index // tile_size  # simple grouping by row blocks

Initial train/validation split
gss1 = GroupShuffleSplit(n_splits=1, train_size=0.7, random_state=42)
train_idx, temp_idx = next(gss1.split(sampled_df, groups=groups))

Validation/test split from the remaining data
gss2 = GroupShuffleSplit(n_splits=1, train_size=0.5, random_state=42)
val_idx, test_idx = next(gss2.split(sampled_df.iloc[temp_idx], groups=groups[temp_idx]))

train_df = sampled_df.iloc[train_idx]
val_df = sampled_df.iloc[temp_idx][val_idx]
test_df = sampled_df.iloc[temp_idx][test_idx]

Common Pitfalls and Mitigation Strategies

Preparing data for training, validation and testing is often perceived as a routine step, yet it is a critical gatekeeper that determines the fidelity of any machine‑learning model. In this section we dissect the most prevalent errors that arise during data preparation, analyse why they occur, and provide concrete, best‑practice solutions that can be implemented immediately. The emphasis is on actionable advice that a data scientist or analyst can apply in a typical Python‑based workflow.

1. Inadequate Train‑Test Splitting

One of the most frequent mistakes is to split the dataset randomly without regard for temporal, spatial or categorical structure. This can lead to data leakage, where information from the test set inadvertently influences the training process. The result is an over‑optimistic estimate of model performance.

Mitigation:

• When dealing with time‑series data, always split chronologically: the earliest data for training, the most recent for validation and testing.
• For grouped data (e.g., patients, customers), perform a group‑wise split so that all records belonging to a particular group appear exclusively in one of the sets.
• Use stratified sampling for classification problems to preserve class proportions across splits.

Example: Splitting a dataset of customer transactions into training, validation and test sets while ensuring that each customer’s records stay together.

import pandas as pd
from sklearn.model_selection import GroupShuffleSplit

Load data
df = pd.read_csv('transactions.csv')

Treat each customer as a group
groups = df['customer_id']

Initialise splitters
train_val_split = GroupShuffleSplit(n_splits=1, test_size=0.3, random_state=42)
train_idx, temp_idx = next(train_val_split.split(df, groups=groups))

Second split: validation and test
val_test_split = GroupShuffleSplit(n_splits=1, test_size=0.5, random_state=42)
val_idx, test_idx = next(val_test_split.split(df.iloc[temp_idx], groups=groups.iloc[temp_idx]))

Create DataFrames
train_df = df.iloc[train_idx]
val_df   = df.iloc[val_idx]
test_df  = df.iloc[test_idx]

print(f'Train size: {len(train_df)}')
print(f'Validation size: {len(val_df)}')
print(f'Test size: {len(test_df)}')

Line by line explanation:

• import pandas as pd – Import the pandas library for data manipulation.
• from sklearn.model_selection import GroupShuffleSplit – Import the splitter that respects group boundaries.
• df = pd.read_csv('transactions.csv') – Read the data into a DataFrame.
• groups = df['customer_id'] – Define the grouping variable.
• train_val_split = GroupShuffleSplit(n_splits=1, test_size=0.3, random_state=42) – Create a splitter that will allocate 70 % of the groups to training and 30 % to a temporary set.
• train_idx, temp_idx = next(train_val_split.split(df, groups=groups)) – Execute the split; train_idx contains indices for training, temp_idx for the remaining 30 %.
• val_test_split = GroupShuffleSplit(n_splits=1, test_size=0.5, random_state=42) – Prepare a second split to divide the temporary set equally into validation and test.
• val_idx, test_idx = next(val_test_split.split(df.iloc[temp_idx], groups=groups.iloc[temp_idx])) – Perform the split on the subset of rows.
• train_df = df.iloc[train_idx] – Extract training rows.
• val_df = df.iloc[val_idx] – Extract validation rows.
• test_df = df.iloc[test_idx] – Extract test rows.
• print(...) – Output the sizes of each split for sanity checking.

2. Mis‑labelled Target Variables

In many real‑world datasets, especially those sourced from legacy systems or scraped from the web, the target column may contain inconsistencies: mixed data types, typos, or ambiguous categories. If left unchecked, the model will learn from noise rather than signal.

Mitigation:

• Standardise all categorical values: strip whitespace, normalise case, and map synonyms to a single label.
• Convert numeric targets to the correct type (e.g., float or integer) and handle missing values explicitly.
• Use a validation set to spot outliers that could indicate label corruption.

3. Feature Leakage

Feature leakage occurs when a feature contains information that would not be available at prediction time. For instance, including a future sales figure in a model predicting current sales will spuriously inflate performance metrics.

Mitigation:

• Analyse the provenance of each feature: ask whether the data would realistically be known at the time of prediction.
• Apply domain knowledge to exclude obviously leaky columns such as target‑derived aggregates.
• When in doubt, perform a sensitivity analysis: remove a suspect feature and observe the change in performance.

4. Ignoring Class Imbalance

Highly imbalanced datasets are common in fraud detection, medical diagnosis, and churn prediction. A naïve model may achieve high accuracy by always predicting the majority class, yet it would be useless in practice.

Mitigation:

• Visualise the class distribution with a bar plot or a table.
• Employ resampling strategies: oversample the minority class (e.g., SMOTE) or undersample the majority class.
• Use appropriate evaluation metrics such as the area‑under‑the‑receiver‑operating‑characteristic curve (AUC‑ROC), precision‑recall curves, or F1‑score.
• Incorporate class weights into the loss function to penalise misclassifications of the minority class more heavily.

5. Inconsistent Feature Engineering Across Splits

Feature engineering steps such as scaling, encoding or dimensionality reduction must be fitted only on the training data and then applied to validation and test sets. Fitting on the entire dataset before partitioning leads to information leakage.

Mitigation:

• Use scikit‑learn pipelines to chain preprocessing and modelling steps, ensuring that each step is correctly fitted.
• When manual feature creation is required, always fit the transformer (e.g., StandardScaler) on the training split and then transform the validation and test splits using the same fitted transformer.
• Document every transformation step in a reproducible script.

6. Over‑Engineering Features

Adding numerous handcrafted features can increase dimensionality and introduce multicollinearity, which may degrade model performance and interpretability, particularly for linear models.

Mitigation:

• Start with a modest feature set and incrementally add features, evaluating performance after each addition.
• Apply feature selection techniques such as recursive feature elimination, L1‑regularisation, or tree‑based importance scores.
• Use domain expertise to prioritise features that are both predictive and interpretable.

7. Poor Handling of Missing Data

Missing values are ubiquitous. Simply dropping rows or columns with missing values can severely reduce sample size and bias the dataset. Conversely, indiscriminate imputation may introduce artefacts.

Mitigation:

• Quantify the extent of missingness per feature.
• For numerical features, consider imputing with the median or a model‑based approach (e.g., KNN imputation).
• For categorical features, impute with the mode or a dedicated “missing” category.
• If the missingness is informative (i.e., not random), create a binary indicator feature that records whether the value was missing.

8. Neglecting Reproducibility

Data preparation is often performed interactively in notebooks, making it difficult to reproduce results later. This hampers collaboration and auditability.

Mitigation:

• Encapsulate the entire preprocessing pipeline in a script or a Jupyter notebook with clear version control.
• Store intermediate artefacts (e.g., transformed datasets, scaler objects) in a structured manner.
• Use random seeds consistently across all randomised operations to ensure deterministic behaviour.

9. Skipping Data Quality Checks

Low‑quality data—duplicates, outliers, inconsistent units—can silently degrade model performance. A common oversight is to assume that cleaning steps are unnecessary if the dataset appears large.

Mitigation:

• Perform a comprehensive data profiling step: summarise statistics, detect duplicates, and visualise distributions.
• Apply outlier detection methods (e.g., isolation forest) and decide on an appropriate treatment strategy (capping, removal, or transformation).
• Standardise units across features (e.g., metres vs. feet) before merging datasets.

10. Failing to Document Pre‑processing Choices

Without thorough documentation, future users cannot understand the rationale behind data transformations, leading to confusion when models are deployed or when new data arrives.

Mitigation:

• Maintain a data dictionary that records the source, type, transformation applied, and any assumptions made for each feature.
• Embed comments in code blocks explaining the purpose of each preprocessing step.
• Use versioning for both data and code to track changes over time.

Conclusion

Data preparation is an art that blends statistical rigour with domain insight. By recognising these ten common pitfalls—and applying the mitigation strategies outlined above—practitioners can minimise the risk of model degradation, ensure that performance metrics are trustworthy, and lay a solid foundation for downstream analytics. The key is to treat data preparation as a disciplined, reproducible process rather than a one‑off task. With careful attention to splitting, leakage, class balance, feature engineering, missing values, and documentation, the resulting models will be robust, interpretable, and ready for real‑world deployment.

Additional Example: Stratified K‑Fold Cross‑Validation for Imbalanced Data

When dealing with class imbalance, a common alternative to a single train‑validation split is to use stratified K‑fold cross‑validation. This ensures that each fold maintains the same class distribution as the overall dataset, providing a more reliable estimate of model performance.

import pandas as pd
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import f1_score
from sklearn.ensemble import RandomForestClassifier

Load dataset
df = pd.read_csv('fraud.csv')
X = df.drop('label', axis=1)
y = df['label']

Initialise stratified K‑fold
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

fold = 1
for train_idx, val_idx in skf.split(X, y):
    X_train, X_val = X.iloc[train_idx], X.iloc[val_idx]
    y_train, y_val = y.iloc[train_idx], y.iloc[val_idx]

    # Train model
    clf = RandomForestClassifier(n_estimators=200, random_state=42)
    clf.fit(X_train, y_train)

    # Evaluate
    preds = clf.predict(X_val)
    score = f1_score(y_val, preds)
    print(f'Fold {fold} F1‑score: {score:.4f}')
    fold += 1

Comparing Splitting Approaches

Further Reading

• J. R. Smith (2018) — Hyperspectral Image Analysis: Algorithms and Applications. Provides comprehensive algorithms for preprocessing hyperspectral data before machine learning.
• U.S. Geological Survey (USGS) (2020) — USGS Spectral Library. Official documentation for spectral data acquisition and calibration, essential for training data quality.
• G. I. Zhang & A. G. Smith (2021) — Deep Learning for Remote Sensing: A Practical Approach. Chapter on data splitting strategies and cross‑validation for hyperspectral datasets.
• International Organisation for Standardisation (ISO) (2019) — ISO 19115: Geographic Information — Metadata. Standardised metadata requirements critical for dataset reproducibility in AI workflows.
• J. M. Jones (2022) — Machine Learning in Earth Observation. Discusses best practices for preparing training, validation and test splits in high‑dimensional spectral data.
• European Space Agency (ESA) (2023) — Sentinel‑2 Data Processing Guidance. Provides guidelines on atmospheric correction and scaling, key for preparing clean training samples.