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