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