Advanced Examples

This page provides advanced usage examples for the HSI Mars package.

Batch Processing Multiple Images

Process multiple hyperspectral images in a loop:

from pathlib import Path
from hsimars import HSIMars

# Directory containing CRISM data
data_dir = Path("data/crism_images")
output_dir = Path("results/batch_analysis")
output_dir.mkdir(parents=True, exist_ok=True)

# Find all .hdr files
hdr_files = list(data_dir.glob("*.hdr"))

print(f"Processing {len(hdr_files)} images...")

for hdr_path in hdr_files:
    print(f"\nProcessing: {hdr_path.name}")

    # Load image
    hsi = HSIMars(hdr_path=hdr_path)
    img_data = hsi.get_img()

    # Save spectrum from center pixel
    center = [img_data.height // 2, img_data.width // 2]
    output_file = output_dir / f"{hdr_path.stem}_spectrum.png"

    hsi.plot_spectra(
        px=center,
        convex_hull=True,
        bands=True,
        output=output_file
    )

    print(f"  - Saved spectrum to {output_file}")

Region of Interest Analysis

Analyze a specific region of interest within an image:

import numpy as np
from hsimars import HSIMars

# Load image
hsi = HSIMars(hdr_path="data/sample.hdr")
img_data = hsi.get_img()

# Define region of interest (ROI)
roi_top_left = [100, 150]
roi_bottom_right = [200, 250]

# Extract all pixels in ROI
roi_pixels = []
for r in range(roi_top_left[0], roi_bottom_right[0] + 1):
    for c in range(roi_top_left[1], roi_bottom_right[1] + 1):
        roi_pixels.append([r, c])

# Plot average spectrum of ROI with standard deviation
hsi.plot_spectra(
    px=roi_pixels,
    convex_hull=True,
    bands=True,
    output="results/roi_spectrum.png"
)

# Calculate statistics
roi_spectra = img_data.hsi[
    roi_top_left[0]:roi_bottom_right[0]+1,
    roi_top_left[1]:roi_bottom_right[1]+1,
    :
]

print(f"ROI dimensions: {roi_spectra.shape[:2]}")
print(f"ROI pixels: {roi_spectra.shape[0] * roi_spectra.shape[1]}")
print(f"Mean intensity: {roi_spectra.mean():.2f}")
print(f"Std intensity: {roi_spectra.std():.2f}")

Class-Based Spectral Analysis

Analyze spectral signatures for different annotated classes:

import numpy as np
from hsimars import HSIMars

# Load data with annotations
hsi = HSIMars(
    hdr_path="data/sample.hdr",
    annotations="data/sample_labels.mat"
)

img_data, ann_data = hsi.data()

# Get unique class labels (excluding background = 0)
classes = np.unique(ann_data.labels)
classes = classes[classes > 0]

print(f"Found {len(classes)} classes")

# Analyze each class
for class_id in classes:
    # Find all pixels belonging to this class
    class_mask = ann_data.labels == class_id
    class_coords = np.argwhere(class_mask)

    # Convert to list of [row, col] pairs
    class_pixels = class_coords.tolist()

    print(f"\nClass {class_id}: {len(class_pixels)} pixels")

    # Plot average spectrum for this class
    if len(class_pixels) > 0:
        hsi.plot_spectra(
            px=class_pixels,
            convex_hull=True,
            bands=True,
            output=f"results/class_{class_id}_spectrum.png"
        )

Spectral Band Comparison

Compare specific spectral bands across the image:

import numpy as np
import matplotlib.pyplot as plt
from hsimars import HSIMars

# Load image
hsi = HSIMars(hdr_path="data/sample.hdr")
img_data = hsi.get_img()

# Select wavelengths of interest
wavelengths_of_interest = [700, 1000, 1500, 2000]  # nm

# Find closest band indices
band_indices = [
    np.argmin(np.abs(img_data.wavelength - wl))
    for wl in wavelengths_of_interest
]

# Create comparison plot
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes = axes.ravel()

for idx, (band_idx, wl) in enumerate(zip(band_indices, wavelengths_of_interest)):
    actual_wl = img_data.wavelength[band_idx]
    band_data = img_data.hsi[:, :, band_idx]

    im = axes[idx].imshow(band_data, cmap='gray')
    axes[idx].set_title(f'Band: {actual_wl:.1f} nm (target: {wl} nm)')
    axes[idx].axis('off')
    plt.colorbar(im, ax=axes[idx], fraction=0.046)

plt.tight_layout()
plt.savefig('results/band_comparison.png', dpi=150, bbox_inches='tight')
print("Saved band comparison to results/band_comparison.png")

Statistical Summary Generation

Generate comprehensive statistics for a dataset:

import numpy as np
from hsimars import HSIMars

def generate_statistics_report(hsi_path, ann_path=None, output_file="report.txt"):
    """Generate a comprehensive statistics report."""

    # Load data
    if ann_path:
        hsi = HSIMars(hdr_path=hsi_path, annotations=ann_path)
        img_data, ann_data = hsi.data()
    else:
        hsi = HSIMars(hdr_path=hsi_path)
        img_data = hsi.get_img()
        ann_data = None

    # Prepare report
    lines = []
    lines.append("="*60)
    lines.append("Hyperspectral Image Statistics Report")
    lines.append("="*60)
    lines.append(f"\nFile: {hsi_path}")
    lines.append(f"\nImage Properties:")
    lines.append(f"  - Dimensions: {img_data.height} x {img_data.width}")
    lines.append(f"  - Spectral bands: {img_data.channels}")
    lines.append(f"  - Data type: {img_data.dtype}")
    lines.append(f"  - Total pixels: {img_data.height * img_data.width:,}")

    lines.append(f"\nSpectral Information:")
    lines.append(f"  - Wavelength range: {img_data.wavelength.min():.2f} - {img_data.wavelength.max():.2f} nm")
    lines.append(f"  - Wavelength spacing: {np.diff(img_data.wavelength).mean():.2f} nm (mean)")

    lines.append(f"\nIntensity Statistics:")
    lines.append(f"  - Global mean: {img_data.hsi.mean():.4f}")
    lines.append(f"  - Global std: {img_data.hsi.std():.4f}")
    lines.append(f"  - Global min: {img_data.hsi.min():.4f}")
    lines.append(f"  - Global max: {img_data.hsi.max():.4f}")

    # Per-band statistics
    band_means = img_data.hsi.mean(axis=(0, 1))
    band_stds = img_data.hsi.std(axis=(0, 1))

    lines.append(f"\nPer-Band Statistics:")
    lines.append(f"  - Mean intensity range: {band_means.min():.4f} - {band_means.max():.4f}")
    lines.append(f"  - Highest mean at: {img_data.wavelength[band_means.argmax()]:.2f} nm")
    lines.append(f"  - Lowest mean at: {img_data.wavelength[band_means.argmin()]:.2f} nm")

    # Annotation statistics
    if ann_data is not None:
        unique_labels, counts = np.unique(ann_data.labels, return_counts=True)

        lines.append(f"\nAnnotation Statistics:")
        lines.append(f"  - Total classes: {len(unique_labels)}")
        lines.append(f"  - Class distribution:")

        for label, count in zip(unique_labels, counts):
            percentage = (count / ann_data.labels.size) * 100
            label_name = "Background" if label == 0 else f"Class {label}"
            lines.append(f"    * {label_name}: {count:,} pixels ({percentage:.2f}%)")

    lines.append("\n" + "="*60)

    # Write to file
    report_text = "\n".join(lines)
    with open(output_file, 'w') as f:
        f.write(report_text)

    print(report_text)
    print(f"\nReport saved to: {output_file}")

    return report_text

# Usage
generate_statistics_report(
    hsi_path="data/sample.hdr",
    ann_path="data/sample_labels.mat",
    output_file="results/statistics_report.txt"
)

Working with Specific Wavelength Ranges

Extract and analyze specific portions of the spectrum:

import numpy as np
from hsimars import HSIMars

# Load image
hsi = HSIMars(hdr_path="data/sample.hdr")
img_data = hsi.get_img()

# Define wavelength range of interest (e.g., SWIR: 1400-3000 nm)
wl_min, wl_max = 1400, 3000

# Find bands within range
mask = (img_data.wavelength >= wl_min) & (img_data.wavelength <= wl_max)
selected_bands = np.where(mask)[0]
selected_wavelengths = img_data.wavelength[mask]

print(f"Selected {len(selected_bands)} bands in range {wl_min}-{wl_max} nm")
print(f"Wavelength range: {selected_wavelengths.min():.1f} - {selected_wavelengths.max():.1f} nm")

# Extract subset
swir_data = img_data.hsi[:, :, selected_bands]

# Analyze average spectrum in this range
center = [img_data.height // 2, img_data.width // 2]
center_spectrum = swir_data[center[0], center[1], :]

# Plot
import matplotlib.pyplot as plt

plt.figure(figsize=(10, 6))
plt.plot(selected_wavelengths, center_spectrum, 'k-')
plt.xlabel('Wavelength (nm)')
plt.ylabel('Intensity (a.u.)')
plt.title(f'SWIR Spectrum ({wl_min}-{wl_max} nm) at Center Pixel')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('results/swir_spectrum.png', dpi=150, bbox_inches='tight')
print("Saved SWIR spectrum plot")

Custom Visualisation

Create custom visualisations using matplotlib:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from hsimars import HSIMars

# Load data
hsi = HSIMars(
    hdr_path="data/sample.hdr",
    annotations="data/sample_labels.mat"
)
img_data, ann_data = hsi.data()

# Create custom figure
fig = plt.figure(figsize=(15, 5))

# Panel 1: False-colour image
ax1 = plt.subplot(131)
# Use specific bands for RGB
false_colour = img_data.hsi[:, :, [50, 100, 150]]
false_colour = (false_colour - false_colour.min()) / (false_colour.max() - false_colour.min())
ax1.imshow(false_colour)
ax1.set_title('False Colour Image')
ax1.axis('off')

# Panel 2: Annotations
ax2 = plt.subplot(132)
# Create custom colourmap
n_classes = len(np.unique(ann_data.labels))
colours = plt.cm.tab10(np.linspace(0, 1, n_classes))
cmap = ListedColormap(colours)

im2 = ax2.imshow(ann_data.labels, cmap=cmap, interpolation='nearest')
ax2.set_title('Ground Truth Labels')
ax2.axis('off')
plt.colorbar(im2, ax=ax2, fraction=0.046)

# Panel 3: Single band
ax3 = plt.subplot(133)
band_idx = img_data.channels // 2  # Middle band
im3 = ax3.imshow(img_data.hsi[:, :, band_idx], cmap='gray')
ax3.set_title(f'Band {band_idx} ({img_data.wavelength[band_idx]:.1f} nm)')
ax3.axis('off')
plt.colorbar(im3, ax=ax3, fraction=0.046)

plt.tight_layout()
plt.savefig('results/custom_visualization.png', dpi=150, bbox_inches='tight')
print("Saved custom visualisation")

Performance Tips

For Large Datasets

# Process in chunks to manage memory
import numpy as np
from hsimars import HSIMars

hsi = HSIMars(hdr_path="data/large_image.hdr")
img_data = hsi.get_img()

# Process image in tiles
tile_size = 100
results = []

for i in range(0, img_data.height, tile_size):
    for j in range(0, img_data.width, tile_size):
        tile = img_data.hsi[
            i:min(i+tile_size, img_data.height),
            j:min(j+tile_size, img_data.width),
            :
        ]
        # Process tile
        result = tile.mean()  # Example operation
        results.append(result)

print(f"Processed {len(results)} tiles")

Parallel Processing

from pathlib import Path
from multiprocessing import Pool
from hsimars import HSIMars

def process_image(hdr_path):
    """Process a single image."""
    hsi = HSIMars(hdr_path=hdr_path)
    img_data = hsi.get_img()

    # Extract some metric
    mean_spectrum = img_data.hsi.mean(axis=(0, 1))

    return {
        'filename': hdr_path.name,
        'mean_intensity': mean_spectrum.mean(),
        'shape': img_data.shape
    }

# Find all images
data_dir = Path("data/crism_images")
hdr_files = list(data_dir.glob("*.hdr"))

# Process in parallel
with Pool(processes=4) as pool:
    results = pool.map(process_image, hdr_files)

# Print summary
for result in results:
    print(f"{result['filename']}: mean={result['mean_intensity']:.4f}")

Further Resources

• See the HSIMars for complete API documentation

• Check the GitHub repository for example notebooks

• Join discussions on GitHub Issues