Analisi MNF su spettri di riflettanza di plastica

Devo cerca di lavorare su spettri di riflettanza di plastica e la prima domanda e': quale sono le bande significative?

Sono partito dal dataset delle plastiche (36 in tutto) contenuti in USGS Spectral Library

Per estrarre le feature ho usato MNF Maximum Noise Fraction (al posto di PCA)  

Nel grafico sottostante sono plottate le prime 3 componenti MNF 

La terza componente (quella di solito in cui sono presenti i dati spettrali, la prima componente e' influenzata dalla ) e' stata utilizzata per addestrare un rete random tree per trovare quali siano le bande effettivamente significative

 Per finire un confronto tra gli spettri non elaborati e la selezione random forest 

 

 Le bande piu' significative estratte 

 ind           nm     peso

150          600     0.018671
1289        1739   0.018409
149          599     0.014763
143          593     0.013682
1998        2448   0.012601 

 

import numpy as np

import glob

import os

from scipy.signal import savgol_filter

from scipy.linalg import eigh

folder = "./plastica"

files = sorted(glob.glob(os.path.join(folder, "*.txt")))

arrays = []

for f in files:

    data = np.loadtxt(f, skiprows=1)

    if data.shape[0] > 2100:

        arrays.append(data[100:2100])

       

X = np.array(arrays)

X = np.where((X >= 0) & (X <= 1), X, 0)

print(f"Shape originale: {X.shape}")

# 1. Smoothing e calcolo rumore

X_smooth = savgol_filter(X, window_length=11, polyorder=2, axis=1)

noise_residue = X - X_smooth

# 2. Calcolo matrici di covarianza

noise_cov = np.cov(noise_residue, rowvar=False)

data_cov = np.cov(X, rowvar=False)

# Regolarizzazione (Ridge) alla diagonale della matrice di rumore.

# Questo garantisce che la matrice sia definita positiva.

eps = 1e-6 * np.trace(noise_cov) / noise_cov.shape[0]

noise_cov_reg = noise_cov + np.eye(noise_cov.shape[0]) * eps

# 3. Esecuzione MNF

eigenvalues, eigenvectors = eigh(data_cov, noise_cov_reg)

idx = np.argsort(eigenvalues)[::-1]

eigenvalues = eigenvalues[idx]

eigenvectors = eigenvectors[:, idx]

# 5. Proiezione

X_mnf = np.dot(X, eigenvectors)

print(f"Prime 5 componenti MNF (SNR più alto): {eigenvalues[:5]}")

import matplotlib.pyplot as plt

# L'indice 0 corrisponde alla prima colonna dopo l'ordinamento decrescente

first_mnf_loadings = eigenvectors[:, 0]

wavelengths = np.arange(350, 2500)

wl_tagliate = wavelengths[100:2100]

idx_max = np.argmax(np.abs(first_mnf_loadings))

print(f"La banda più importante è la numero {idx_max}, che corrisponde a {wl_tagliate[idx_max]} nm")

wls = np.arange(350, 2501)[100:2100]

fig, axs = plt.subplots(3, 1, figsize=(12, 15), sharex=True)

colors = ['#1f77b4', '#ff7f0e', '#2ca02c']

for i in range(3):

    # Estraiamo i loadings per la componente i-esima

    loadings = eigenvectors[:, i]

   

    axs[i].plot(wls, loadings, color=colors[i], lw=1.5)

    axs[i].set_title(f"Loadings MNF Componente {i+1} (Autovalore: {eigenvalues[i]:.2e})")

    axs[i].set_ylabel("Peso")

    axs[i].grid(True, alpha=0.3)

   

    # Evidenziamo l'area diagnostica dei carbonati (2300-2350 nm)

    axs[i].axvspan(2300, 2350, color='red', alpha=0.1, label='Zona Carbonati')

axs[2].set_xlabel("Lunghezza d'onda (nm)")

plt.tight_layout()

plt.show()

from sklearn.ensemble import RandomForestRegressor

import pandas as pd

y = X_mnf[:, 2]

rf = RandomForestRegressor(n_estimators=100, random_state=42)

rf.fit(X, y)

importances = rf.feature_importances_

wls = np.arange(350, 2501)[100:2100]

feature_results = pd.DataFrame({'Wavelength': wls, 'Importance': importances})

# Ordiniamo per importanza decrescente

top_features = feature_results.sort_values(by='Importance', ascending=False).head(10)

print("Le 10 lunghezze d'onda più influenti (Feature Importance):")

print(top_features)

# 5. Visualizzazione

plt.figure(figsize=(10, 5))

plt.plot(wls, importances, color='purple', label='Feature Importance')

plt.fill_between(wls, importances, color='purple', alpha=0.3)

plt.title("Rilevanza delle Bande Spettrali (Random Forest)")

plt.xlabel("Lunghezza d'onda (nm)")

plt.ylabel("Importanza Relativa")

plt.grid(True, alpha=0.3)

plt.show()

plt.figure(figsize=(12, 7))

plt.plot(wls, X.T, color='steelblue', alpha=0.3)

plt.title(f"Spettri di Riflettanza del Dataset ({X.shape[0]} campioni)")

plt.xlabel("Lunghezza d'onda (nm)")

plt.ylabel("Riflettanza (0-1)")

plt.grid(True, linestyle='--', alpha=0.5)

# Per evitare legende infinite, ne mettiamo una generica

plt.plot([], [], color='steelblue', label='Spettri campioni')

plt.legend()

plt.show()

# --- GRAFICO 1: SPETTRI DI RIFLETTANZA ---

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 12), sharex=True)

ax1.plot(wls, X.T, color='steelblue', alpha=0.2)

ax1.plot(wls, np.mean(X, axis=0), color='red', lw=2, label='Spettro Medio')

ax1.set_title("Spettri di Riflettanza (Dati Originali Puliti)")

ax1.set_ylabel("Riflettanza")

ax1.grid(True, alpha=0.3)

ax1.legend()

# --- GRAFICO 2: RILEVANZA DELLE BANDE (Feature Importance) ---

ax2.plot(wls, importances, color='darkgreen', lw=1.5)

ax2.fill_between(wls, importances, color='darkgreen', alpha=0.2)

ax2.set_title("Rilevanza delle Bande Spettrali (Random Forest Importance)")

ax2.set_ylabel("Importanza Relativa")

ax2.set_xlabel("Lunghezza d'onda (nm)")

ax2.grid(True, alpha=0.3)

plt.tight_layout()

plt.show()