Spettri di Prisma

Sono ritornato a lavorare sui dati Prisma, stavolta in maniera piu' convinta

Il primo problema e' trovare il software giusto. Se si aprono le immagini con Esa Snap non ci sono particolari problemi tranne quando si arriva ad aprire Spectrum View e si scopre che le bande dell'immagini sono divise in due VNIR e SWIR e quindi non e' possibile avere uno spettro completo da 400 a 2500 nm. Si devono aprire le due immagini ed a seconda di dove e' il cursore Spectrum View si spostera' tra i dati VNIR e SWIR 

 

Per avere uno spettro completo l'unica soluzione e' stata passare da uno script Python (compilato via Claude AI)

 

"""

PRISMA Hyperspectral Spectrum Plotter (L2C)

=============================================

Reads a PRISMA L2C HE5 file and plots the full VNIR+SWIR spectrum

for a user-selected pixel.

Confirmed HDF5 layout:

Cubes : HDFEOS/SWATHS/PRS_L2C_HCO/Data Fields/{VNIR,SWIR}_Cube

shape = (rows, bands, cols)

Wavelengths : root attrs List_Cw_Vnir / List_Cw_Swir (nm, length=n_bands)

Band flags : root attrs List_Cw_Vnir_Flags / List_Cw_Swir_Flags (1=valid)

Scale : root attrs L2ScaleVnirMin/Max, L2ScaleSwirMin/Max

DN → reflectance = DN / 65535 * (ScaleMax - ScaleMin) + ScaleMin

Usage

-----

python prisma_spectrum_plot.py <file.he5> [row] [col] [-o out.png]

Dependencies

------------

pip install h5py numpy matplotlib

"""

import argparse

import numpy as np

import h5py

import matplotlib.pyplot as plt

import matplotlib.ticker as ticker

from pathlib import Path

CUBE_VNIR = "HDFEOS/SWATHS/PRS_L2C_HCO/Data Fields/VNIR_Cube"

CUBE_SWIR = "HDFEOS/SWATHS/PRS_L2C_HCO/Data Fields/SWIR_Cube"

def load_prisma(filepath: str, row: int | None = None, col: int | None = None):

"""

Extract a single-pixel VNIR+SWIR spectrum from a PRISMA L2C HE5 file.

Parameters

----------

filepath : path to the .he5 file

row, col : pixel coordinates (0-based). Default = image centre.

Returns

-------

wavelengths : np.ndarray – wavelengths in nm (valid bands only, sorted)

spectrum : np.ndarray – reflectance [0–1]

meta : dict

"""

filepath = Path(filepath)

if not filepath.exists():

raise FileNotFoundError(filepath)

with h5py.File(filepath, "r") as f:

# ── wavelengths and valid-band flags from root attributes ─────────────

vnir_wl = np.array(f.attrs["List_Cw_Vnir"], dtype=np.float32)

swir_wl = np.array(f.attrs["List_Cw_Swir"], dtype=np.float32)

vnir_flags = np.array(f.attrs["List_Cw_Vnir_Flags"], dtype=np.int8)

swir_flags = np.array(f.attrs["List_Cw_Swir_Flags"], dtype=np.int8)

vnir_valid = vnir_flags == 1 # boolean mask for valid VNIR bands

swir_valid = swir_flags == 1 # boolean mask for valid SWIR bands

# ── scale factors: uint16 DN → reflectance ────────────────────────────

vnir_scale_min = float(f.attrs["L2ScaleVnirMin"])

vnir_scale_max = float(f.attrs["L2ScaleVnirMax"])

swir_scale_min = float(f.attrs["L2ScaleSwirMin"])

swir_scale_max = float(f.attrs["L2ScaleSwirMax"])

# ── cube shapes: (rows, bands, cols) ──────────────────────────────────

vnir_ds = f[CUBE_VNIR]

swir_ds = f[CUBE_SWIR]

n_rows, n_vnir_bands, n_cols = vnir_ds.shape

_, n_swir_bands, _ = swir_ds.shape

if row is None: row = n_rows // 2

if col is None: col = n_cols // 2

if not (0 <= row < n_rows and 0 <= col < n_cols):

raise ValueError(

f"Pixel ({row}, {col}) outside image bounds ({n_rows}x{n_cols})."

)

# ── read single pixel: cube[row, :, col] ─────────────────────────────

vnir_dn = vnir_ds[row, :, col].astype(np.float32)

swir_dn = swir_ds[row, :, col].astype(np.float32)

# ── DN → reflectance ──────────────────────────────────────────────────────

vnir_ref = vnir_dn / 65535.0 * (vnir_scale_max - vnir_scale_min) + vnir_scale_min

swir_ref = swir_dn / 65535.0 * (swir_scale_max - swir_scale_min) + swir_scale_min

# ── apply valid-band masks ────────────────────────────────────────────────

# Flags array length matches n_bands; trim to cube size just in case

vnir_mask = vnir_valid[:n_vnir_bands]

swir_mask = swir_valid[:n_swir_bands]

vnir_wl = vnir_wl[:n_vnir_bands][vnir_mask]

vnir_ref = vnir_ref[vnir_mask]

swir_wl = swir_wl[:n_swir_bands][swir_mask]

swir_ref = swir_ref[swir_mask]

# ── sort by wavelength (PRISMA stores bands longest-first) ───────────────

def sort_wl(wl, sp):

idx = np.argsort(wl)

return wl[idx], sp[idx]

vnir_wl, vnir_ref = sort_wl(vnir_wl, vnir_ref)

swir_wl, swir_ref = sort_wl(swir_wl, swir_ref)

# ── remove VNIR/SWIR overlap: keep VNIR below SWIR start ─────────────────

vnir_keep = vnir_wl < swir_wl[0]

wavelengths = np.concatenate([vnir_wl[vnir_keep], swir_wl])

spectrum = np.concatenate([vnir_ref[vnir_keep], swir_ref])

meta = dict(

row=row, col=col,

n_rows=n_rows, n_cols=n_cols,

n_vnir_valid=int(vnir_mask.sum()),

n_swir_valid=int(swir_mask.sum()),

filename=filepath.name,

)

return wavelengths, spectrum, meta

# ── plot ───────────────────────────────────────────────────────────────────────

COLOUR_REGIONS = [

(400, 700, "#e8f4e8", "VIS"),

(700, 1000, "#f5e8f5", "NIR"),

(1000, 1800, "#e8eef5", "SWIR-1"),

(1800, 2500, "#f5f0e8", "SWIR-2"),

]

ABSORPTION_BANDS = [

(1340, 1460, "H\u2082O"),

(1800, 1960, "H\u2082O"),

]

def plot_spectrum(wavelengths, spectrum, meta, output_path=None):

fig, ax = plt.subplots(figsize=(13, 5))

wl_min, wl_max = wavelengths[0], wavelengths[-1]

sp_max = np.nanmax(spectrum)

# Spectral region backgrounds

for lo, hi, colour, label in COLOUR_REGIONS:

lo_c, hi_c = max(lo, wl_min), min(hi, wl_max)

if lo_c < hi_c:

ax.axvspan(lo_c, hi_c, color=colour, alpha=0.45, zorder=0)

ax.text((lo_c + hi_c) / 2, sp_max * 1.02,

label, ha="center", va="bottom", fontsize=8,

color="#888888", zorder=2)

# Atmospheric absorption windows

for lo, hi, label in ABSORPTION_BANDS:

lo_c, hi_c = max(lo, wl_min), min(hi, wl_max)

if lo_c < hi_c:

ax.axvspan(lo_c, hi_c, color="#bbbbbb", alpha=0.55, zorder=1)

ax.text((lo_c + hi_c) / 2, sp_max * 0.96,

label, ha="center", va="top", fontsize=8,

color="#444444", zorder=3)

ax.plot(wavelengths, spectrum, color="#1a5276", linewidth=1.3, zorder=4)

ax.fill_between(wavelengths, spectrum, alpha=0.12, color="#1a5276", zorder=3)

ax.set_xlim(wl_min, wl_max)

ax.set_ylim(bottom=0)

ax.set_xlabel("Wavelength (nm)", fontsize=11)

ax.set_ylabel("Reflectance", fontsize=11)

ax.set_title(

f"PRISMA L2C Spectrum - pixel ({meta['row']}, {meta['col']}) "

f"[{meta['n_rows']}x{meta['n_cols']}]\n{meta['filename']}",

fontsize=10, pad=8

)

ax.xaxis.set_minor_locator(ticker.AutoMinorLocator(5))

ax.yaxis.set_minor_locator(ticker.AutoMinorLocator(4))

ax.tick_params(which="both", direction="in")

ax.grid(True, which="major", linestyle="--", alpha=0.35)

ax.annotate(

f"VNIR: {meta['n_vnir_valid']} bands | SWIR: {meta['n_swir_valid']} bands",

xy=(0.99, 0.02), xycoords="axes fraction",

ha="right", va="bottom", fontsize=8, color="#666666"

)

plt.tight_layout()

if output_path:

fig.savefig(output_path, dpi=150, bbox_inches="tight")

print(f"[OK] Saved: {output_path}")

else:

plt.show()

plt.close(fig)

# ── CLI ────────────────────────────────────────────────────────────────────────

def main():

parser = argparse.ArgumentParser(

description="Plot a VNIR+SWIR pixel spectrum from a PRISMA L2C HE5 file."

)

parser.add_argument("he5_file")

parser.add_argument("row", nargs="?", type=int, default=None,

help="Pixel row (default: centre)")

parser.add_argument("col", nargs="?", type=int, default=None,

help="Pixel col (default: centre)")

parser.add_argument("-o", "--output", default=None,

help="Save figure to file (e.g. spectrum.png)")

args = parser.parse_args()

print(f"[...] Reading {args.he5_file}")

wavelengths, spectrum, meta = load_prisma(args.he5_file, args.row, args.col)

print(

f"[OK] {meta['n_vnir_valid']} VNIR + {meta['n_swir_valid']} SWIR valid bands\n"

f" Image : {meta['n_rows']} rows x {meta['n_cols']} cols\n"

f" Pixel : row={meta['row']}, col={meta['col']}\n"

f" Range : {wavelengths[0]:.1f} - {wavelengths[-1]:.1f} nm "

f"({len(wavelengths)} merged bands)"

)

plot_spectrum(wavelengths, spectrum, meta, output_path=args.output)

if __name__ == "__main__":

main()

 

 

 

 

 

 

 

 

 

 

 

QGis e PythonPath

Ultimamente ho incasinato la Debian box e mi compare questo errore aprendo QGis...dopo funziona ma non disponibili i plugin 

la soluzione e' settare in modo esplicito la PythoPath

PYTHONPATH=/usr/lib/python3/dist-packages:$PYTHONPATH 

qgis

HyperCoast Python Library

Hypercoast e' una libreria Python per leggere dati da formati iperspettrali e visualizzare spettri

L'uso migliore e' quello all'interno di Jupyter Lab

 

import hypercoast

filepath = "ang20210604t105418_rfl_v2z1_img"

ds = hypercoast.read_aviris(filepath)

print(ds)

m = hypercoast.Map()

print(m)

m

m.add_aviris(ds, wavelengths=[1000, 700, 400], vmin=0, vmax=0.2)

m.add("spectral")

display(m)

 

Aviris-NG Grosseto

 

 

 

 

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