Machine Learning on Hyperspectral Scenes

# !pip install tanagerspec

# Standard library imports
from datetime import datetime               # For timestamping output folders
from pathlib import Path                    # For handling file and directory paths

# TanagerSpec package imports
from tanagerspec import download_scene      # Utility to download a Tanager scene
from tanagerspec import TanagerSpec         # Main class for Tanager hyperspectral data

# third party imports
import numpy as np
from scipy.stats import pearsonr
from sklearn.preprocessing import StandardScaler

# 1. Setup data and output directories
DATA_DIR = Path("data")
DATA_DIR.mkdir(parents=True, exist_ok=True)

timestamp = datetime.now().strftime("%Y%m%d_%H")   # unified format across the series
OUTPUT_DIR = Path(f"outputs_{timestamp}")
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

FIGURE_DIR = OUTPUT_DIR / "figures"
FIGURE_DIR.mkdir(parents=True, exist_ok=True)

# 2. Define scene parameters (clear agricultural scene, same as the Visualization lesson)
URL = "https://storage.googleapis.com/open-cogs/planet-stac/tanager1-release2-core-imagery/ortho_sr_hdf5/20250501_143138_87_4001_ortho_sr_hdf5.h5"
FILE_NAME = "ortho_sr_scene.h5"          # shared clear scene, reused across lessons
target_path = DATA_DIR / FILE_NAME

# 3. Smart download: only fetch if the file isn't already present
if target_path.exists():
    print(f"Scene already exists locally at: {target_path}. Skipping download.")
    TANAGER_SCENE_PATH = target_path
else:
    print("Local scene not found. Downloading from STAC...")
    TANAGER_SCENE_PATH = download_scene(URL, target_path)

# 4. Verify final path
if TANAGER_SCENE_PATH is not None and Path(TANAGER_SCENE_PATH).exists():
    print(f"Ready for analysis with scene: {TANAGER_SCENE_PATH}")
else:
    print("Error: Scene file is missing or download failed.")

# 5. Initialize TanagerSpec and prepare a clean cube for analysis
tanager_ortho_sr = TanagerSpec.from_file(TANAGER_SCENE_PATH)
tanager_ortho_sr.preprocess(masking=True, clipping=True)   # mask invalid/cloud pixels, clip to [0,1]

# 6. Denoise the scene (Optional)
tanager_ortho_sr.denoise(
    n_components=4,
)

print("TanagerSpec object initialized and preprocessed!")

pca = tanager_ortho_sr.analysis.dim_reduction(
    method="MNF",                              # Replace it with "PCA" or "ICA"
    n_components=2,                            # Define the number of components
    present_rgb="false_color_nir",             # Define the RGB preset
    plot=True,                                 # Define if the plot is shown
    scatter_color_by="index",                  # Define how the scatter is colored
    scatter_index="NDVI",                      # Define the index used for coloring
    scatter_cmap="RdYlGn",                     # Define the colormap used for coloring
    save_png=FIGURE_DIR / "dim_reduction_plot.png"
)

# --- Layer 1: a spectral index, straight from a TanagerSpec API ------------------
evi_array = tanager_ortho_sr.analysis.calculate_index(
    index_name="EVI",
    plot=False,
)

# --- Layer 2: the first reduced component from dim_reduction() --------------------
component1 = pca[:, :, 0]

# --- Standardize both arrays without flattening -----------------------------------
scaler1 = StandardScaler()
component1_scaled = scaler1.fit_transform(component1)

scaler2 = StandardScaler()
evi_scaled = scaler2.fit_transform(evi_array)

# --- Correlation (use only valid overlapping pixels) ------------------------------
# Mask out invalid (nan) values so correlation isn't affected
valid_mask = np.isfinite(component1_scaled) & np.isfinite(evi_scaled)
r, _ = pearsonr(component1_scaled[valid_mask], evi_scaled[valid_mask])
print(f"Pearson r between reduced component 1 and EVI: {r:.3f}")

# --- Visualization ----------------------------------------------------------------
comparison = tanager_ortho_sr.analysis.compare_layers(
    component1_scaled,
    evi_scaled,
    index1_name="Component 1 (scaled)",
    index2_name="EVI (scaled)",
    save_png=FIGURE_DIR / "comparison_PCA1_EVI.png",
)

clustered_array = tanager_ortho_sr.analysis.clustering(
    dr_method="PCA",                  # Dimension reduction: "PCA", "ICA", "MNF", or None for all bands
    n_components=3,                   # Number of reduced dimensions (if using dr_method)
    clustering_method="KMEANS",       # Clustering algorithm: "KMEANS" or "GMM"
    n_clusters=5,                     # Number of clusters to identify
    random_state=42,                  # For reproducible labels
    save_png=FIGURE_DIR / "clustering_plot.png",        # Save the cluster plot as PNG
    save_dpi=300,                                     # PNG resolution
    save_geotiff=OUTPUT_DIR / "ortho_clustering.tif",  # Export clusters as GeoTIFF
    # present_rgb="false_color_nir",                   # Optionally set RGB background
)

# run hunt_pixels to get the coordinates of the training pixels
# tanager_ortho_sr.plot.hunt_pixels()

library_stats, df_training = tanager_ortho_sr.build_spectral_library(
    # class name -> (col, row), or a list of (col, row) tuples for several training pixels.
    # NOTE: this is (col, row) = (x, y) -- the TRANSPOSE of hunt_pixels' (row, col).
    # Verify/adjust each point with tanager_ortho_sr.plot.hunt_pixels() and swap the
    # two numbers (hunt_pixels reports row, col) before pasting them here.
    targets={
        "veg1":     [(369, 364)],
        "veg2":     (188, 287),
        "building": (294, 509),
        "water":    (518, 451),
    },
    window_size=5,                # 5x5 window (25 pixels) averaged around each point
    export_csv=OUTPUT_DIR / "ml_spectral_training_set.csv",
)

# Prepare the mean spectra for each class, required for the SAM algorithm
pure_means = {k: v["mean"] for k, v in library_stats.items()}

# Set the SAM threshold: the maximum spectral angle (in radians) to accept a match; larger angles are left unclassified
sam_threshold = 0.1

# Classify every pixel in the scene using the Spectral Angle Mapper (SAM) method.
# This compares the spectrum at each pixel to the mean spectrum of every class, assigning the closest match.
# Outputs a GeoTIFF where each pixel is labelled with its predicted class.
sam_prediction = tanager_ortho_sr.analysis.classify_scene(
    method="SAM",                             # Use the SAM algorithm for classification
    library_means=pure_means,                 # Pass the class mean spectra
    sam_threshold=sam_threshold,              # Specify the SAM threshold
    present_rgb="false_color_nir",            # Use NIR false color as background for visualization
    save_geotiff=OUTPUT_DIR / "ortho_sam_classification.tif",  # Save the predicted label map as GeoTIFF
)

# Classify every pixel in the scene using the Random Forest (RF) algorithm.
# This method uses the training DataFrame produced earlier to train a collection of decision trees.
# Each pixel is labeled based on the class predictions of the trained RF model.
rf_prediction = tanager_ortho_sr.analysis.classify_scene(
    method="RF",                                    # Use the Random Forest algorithm
    df_training=df_training,                        # Training data (features + labels)
    rf_estimators=100,                              # Number of trees in the forest
    present_rgb="false_color_nir",                  # Use NIR false color as the background for visualization
    confidence_threshold=0.6,                       # Only classify pixels with at least 60% class probability
    save_geotiff=OUTPUT_DIR / "ortho_rf_classification.tif",  # Save the classification result as a GeoTIFF
)

# Classify every pixel in the scene using a Neural Network (NN) classifier.
# - This method uses a Multi-layer Perceptron (MLP) trained on the earlier training dataframe.
# - Each pixel’s spectrum is fed to the trained NN; pixels are labelled with their predicted class.
# - Only confidently classified pixels (≥ 50% probability) receive a class label, otherwise "Unclassified".
# - Two hidden layers are used (10, 10, 5 neurons), with up to 1000 training iterations.
# - The result is visualized over a NIR false color background,
#   and saved as both a PNG and a GeoTIFF.
nn_predictions = tanager_ortho_sr.analysis.classify_scene(
    method="NN",                                        # Use Neural Network classification
    df_training=df_training,                            # Training dataset (features + labels)
    nn_hidden_layer_sizes=(10, 10, 5),                  # Three hidden layers: 10, 10, and 5 neurons
    nn_max_iter=1000,                                   # Maximum training iterations
    present_rgb="false_color_nir",                      # NIR false color as visualization background
    confidence_threshold=0.5,                           # Minimum classification confidence (50%)
    nn_early_stopping=False,                            # Do not use early stopping in NN training
    plot=True,                                          # Display classification results
    save_png=FIGURE_DIR / "nn_classification.png",      # Save PNG visualization
    save_geotiff=OUTPUT_DIR / "nn_predictions.tif",     # Save predicted labels as GeoTIFF
)