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()
