Il sensore FigSpec FS60-C mostra problemi di striping con bande non disturbate solo da 450 nm fino a 850 nm
Questo problema sembra derivare da due fattori 1) un basso rapporto segnale/rumore (stimato in 23 dB) 2) in problemi di elettronica legati probabilmente al tempo di integrazione
Per una correzione esclusivamente software si puo' implementare local moment matching con baseline a filtro mediano
python3 analyze_crosstrack_cube.py cubo.hdr --out out_cube --save-corrected
#!/usr/bin/env python3
"""
analyze_crosstrack_cube.py
Analizza un cubo iperspettrale pushbroom (formato ENVI: .hdr + .bil/.bsq/.bip,
oppure un array numpy .npy/.npz) per individuare e quantificare differenze di
illuminazione lungo la direzione across-track (le colonne del sensore),
banda per banda.
Perche' questo script serve
---------------------------
Nel true-color renderizzato avevi notato una fascia destra piu' scura della
sinistra. Su un solo composito RGB non si puo' distinguere se il problema
nasce:
(a) nel sensore/ottica (vignettatura radiometrica, gia' presente nella
riflettanza calibrata), oppure
(b) nel rendering RGB (stretch, white balance, ecc.), oppure
(c) da un vero effetto fisico di scena (BRDF/anisotropia, ombre, angolo
di illuminazione solare rispetto alla direzione di volo).
Lavorando sul cubo (idealmente gia' calibrato in riflettanza con dark/white/
pannello) e guardando il profilo across-track banda per banda, si puo' capire
se il gradiente e' costante su tutte le lunghezze d'onda (probabile
vignettatura/ottica del sensore, o dark/white non uniformi lungo lo swath)
oppure cambia con la banda (piu' probabile un effetto fisico/spettrale di
scena o un residuo di calibrazione spettrale).
Cosa fa lo script
-----------------
1. Carica il cubo (rows=along-track, cols=across-track/swath, bands=spettrale).
Auto-rileva l'interleave ENVI (BIL tipico per pushbroom) tramite 'spectral'.
2. Maschera i pixel nodata (es. tutti zero o NaN, o sotto una soglia).
3. Per ogni banda calcola il profilo medio/mediano per colonna (across-track).
4. Quantifica il gradiente sinistra/destra per banda con un fit lineare e
con il confronto diretto delle due meta'.
5. Produce:
- una heatmap (banda x colonna) del profilo across-track normalizzato,
per vedere a colpo d'occhio se il gradiente e' presente su tutte le bande
o solo in alcune (utile per capire la causa fisica),
- un grafico dello slope (pendenza) del gradiente in funzione della
lunghezza d'onda/banda,
- un grafico di alcune bande rappresentative (es. blu/verde/rosso/NIR se
riconoscibili, altrimenti prima/centro/ultima banda),
- un report testuale banda per banda,
- opzionalmente un cubo corretto (flat-field per colonna, per banda) in
formato .npy, utile per verificare se il problema sparisce.
Uso
---
# Cubo ENVI (serve il file .hdr, il dato binario associato viene trovato
# automaticamente in base al campo 'file' nell'header o stesso nome)
python3 analyze_crosstrack_cube.py cubo.hdr --out out_cube
# Cubo numpy gia' in memoria come array (righe, colonne, bande)
python3 analyze_crosstrack_cube.py cubo.npy --out out_cube --shape rows,cols,bands
# Salva anche il cubo corretto
python3 analyze_crosstrack_cube.py cubo.hdr --out out_cube --save-corrected
"""
import argparse
import sys
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
# --------------------------------------------------------------------------
# Caricamento cubo
# --------------------------------------------------------------------------
def _parse_envi_header_datafile(hdr_path):
"""
Legge il campo 'file =' / 'data file =' dall'header ENVI, se presente,
per capire come si chiama il file binario associato. Ritorna None se non
trovato o non parsabile.
"""
try:
text = Path(hdr_path).read_text(encoding="utf-8", errors="ignore")
except Exception:
return None
for line in text.splitlines():
low = line.strip().lower()
if low.startswith("file =") or low.startswith("data file ="):
value = line.split("=", 1)[1].strip()
value = value.strip("{}").strip()
if value:
return value
return None
def _find_envi_data_file(hdr_path):
"""
Cerca il file binario associato a un header ENVI quando 'spectral' non
riesce a trovarlo automaticamente (lo fa solo per estensioni .img, .dat,
.sli o nessuna estensione). Strategie, in ordine:
1. Campo 'file =' dentro l'header stesso.
2. Qualunque file nella stessa cartella con lo stesso nome base
(stesso nome del .hdr senza estensione) e estensione diversa da
.hdr, .png, .jpg, .txt, .csv, .py (per evitare falsi positivi).
Ritorna il Path del file trovato, o None.
"""
hdr_path = Path(hdr_path)
base = hdr_path.with_suffix("") # rimuove .hdr
declared = _parse_envi_header_datafile(hdr_path)
if declared:
candidate = (hdr_path.parent / declared)
if candidate.exists():
return candidate
candidate2 = hdr_path.parent / Path(declared).name
if candidate2.exists():
return candidate2
exclude_suffixes = {".hdr", ".png", ".jpg", ".jpeg", ".txt", ".csv", ".py", ".md"}
candidates = [p for p in hdr_path.parent.glob(base.name + "*")
if p.is_file() and p.suffix.lower() not in exclude_suffixes and p != hdr_path]
if candidates:
candidates.sort(key=lambda p: p.stat().st_size, reverse=True)
return candidates[0]
return None
def load_cube(path, data_file_override=None):
"""
Carica un cubo iperspettrale da:
- file ENVI header (.hdr) -> usa la libreria 'spectral'
- file binario ENVI (qualunque estensione, es. .spe, .bil, .bsq, .bip,
.raw, .dat...) che ha un .hdr associato nella stessa cartella
- file numpy (.npy / .npz)
Ritorna array numpy (rows, cols, bands) in float64, e un dict di metadati
(wavelengths se disponibili, nomi banda, ecc.)
"""
path = Path(path)
meta = {"wavelengths": None, "band_names": None, "source": str(path)}
if path.suffix.lower() == ".npy":
cube = np.load(str(path))
return cube.astype(np.float64), meta
if path.suffix.lower() == ".npz":
data = np.load(str(path))
for key in ("cube", "data", "arr_0"):
if key in data:
return data[key].astype(np.float64), meta
first_key = list(data.keys())[0]
return data[first_key].astype(np.float64), meta
hdr_path = None
data_path = Path(data_file_override) if data_file_override else None
if path.suffix.lower() == ".hdr":
hdr_path = path
else:
candidate_hdrs = [
path.with_suffix(".hdr"),
Path(str(path) + ".hdr"),
]
for c in candidate_hdrs:
if c.exists():
hdr_path = c
if data_path is None:
data_path = path
break
if hdr_path is None or not hdr_path.exists():
raise ValueError(
f"Formato non riconosciuto o header ENVI non trovato per {path}. "
f"Usa un file .hdr (ENVI), il file binario con .hdr associato, oppure .npy/.npz"
)
import spectral
if data_path is None:
try:
img = spectral.envi.open(str(hdr_path))
except spectral.io.envi.EnviDataFileNotFoundError:
found = _find_envi_data_file(hdr_path)
if found is None:
raise ValueError(
f"Impossibile trovare il file dati binario associato a {hdr_path}. "
f"Passa esplicitamente il file dati con --data-file, oppure rinomina/"
f"copia il binario accanto all'header con estensione .img/.dat/.bin."
)
print(f"[info] File dati non auto-rilevato da 'spectral'; uso: {found}")
img = spectral.envi.open(str(hdr_path), image=str(found))
else:
img = spectral.envi.open(str(hdr_path), image=str(data_path))
cube = img.load().astype(np.float64)
wl = getattr(img, "bands", None)
if wl is not None and getattr(wl, "centers", None) is not None:
meta["wavelengths"] = np.array(wl.centers, dtype=float)
# conserviamo l'header originale (dict) e l'interleave, utili per
# riesportare il cubo corretto in formato ENVI mantenendo wavelength,
# unita', sensor type, ecc.
meta["envi_header"] = dict(getattr(img, "metadata", {}) or {})
meta["interleave"] = getattr(img, "interleave", None)
return np.asarray(cube), meta
def ensure_rows_cols_bands(cube, shape_hint=None):
"""
Se il cubo non e' gia' (rows, cols, bands), e l'utente ha fornito
--shape rows,cols,bands, lo reshape. Altrimenti assume che sia gia'
nell'ordine corretto (che e' il default per 'spectral' e per la
convenzione usata nello script true-color della sessione precedente).
"""
if shape_hint is not None:
rows, cols, bands = shape_hint
cube = cube.reshape(rows, cols, bands)
return cube
# --------------------------------------------------------------------------
# Maschera nodata
# --------------------------------------------------------------------------
def nodata_mask(cube, thresh=1e-6):
"""
True dove il pixel e' valido. Un pixel e' nodata se e' NaN su tutte le
bande oppure se il valore assoluto medio su tutte le bande e' sotto
soglia (tipico bordo nero da rotazione/georeferenziazione).
"""
finite = np.isfinite(cube).all(axis=2)
mean_abs = np.nanmean(np.abs(cube), axis=2)
valid = finite & (mean_abs > thresh)
return valid
# --------------------------------------------------------------------------
# Profilo across-track per banda
# --------------------------------------------------------------------------
def crosstrack_profile_per_band(cube, mask):
"""
Ritorna:
mean_profile: (cols, bands) media per colonna e banda (NaN se nessun
pixel valido in quella colonna)
valid_count: (cols,) numero di pixel validi per colonna
"""
rows, cols, bands = cube.shape
mean_profile = np.full((cols, bands), np.nan)
valid_count = mask.sum(axis=0)
for c in range(cols):
col_mask = mask[:, c]
if col_mask.sum() == 0:
continue
mean_profile[c, :] = cube[col_mask, c, :].mean(axis=0)
return mean_profile, valid_count
def linear_fit(x, y):
good = np.isfinite(y)
if good.sum() < 2:
return np.nan, np.nan, np.nan
xs, ys = x[good], y[good]
A = np.vstack([xs, np.ones_like(xs)]).T
slope, intercept = np.linalg.lstsq(A, ys, rcond=None)[0]
pred = slope * xs + intercept
ss_res = np.sum((ys - pred) ** 2)
ss_tot = np.sum((ys - ys.mean()) ** 2)
r2 = 1 - ss_res / ss_tot if ss_tot > 0 else np.nan
return slope, intercept, r2
def export_envi_cube(out_path, cube, meta, interleave="bil"):
"""
Salva un cubo (rows, cols, bands) in formato ENVI (coppia .hdr + file dati),
riusando i metadati originali (wavelength, unita', sensor type, ecc.) se
disponibili in meta['envi_header'], cosi' il cubo esportato resta
apribile in ENVI/QGIS/altri tool con le stesse informazioni spettrali
dell'originale.
`out_path` puo' essere passato con o senza estensione: verra' sempre
scritto come coppia <out_path>.hdr + <out_path>.dat (interleave BIL di
default, tipico per dati pushbroom).
"""
import spectral
out_path = Path(out_path)
if out_path.suffix.lower() == ".hdr":
out_path = out_path.with_suffix("")
orig_header = dict(meta.get("envi_header") or {})
# Ripulisce campi che 'spectral' ricalcola da solo o che diventerebbero
# inconsistenti (dimensioni, nome file, offset) per evitare header corrotti
for key in ("samples", "lines", "bands", "header offset", "file type",
"byte order", "data type", "interleave"):
orig_header.pop(key, None)
metadata = orig_header # contiene ad es. wavelength, wavelength units, sensor type, map info...
spectral.envi.save_image(
str(out_path.with_suffix(".hdr")),
cube.astype(np.float32),
dtype=np.float32,
force=True,
interleave=interleave,
metadata=metadata,
)
return out_path.with_suffix(".hdr")
# --------------------------------------------------------------------------
# Analisi principale
# --------------------------------------------------------------------------
def analyze(cube_path, out_dir, shape_hint=None, nodata_thresh=1e-6,
split_fraction=0.5, save_corrected=False, data_file_override=None,
output_format="envi"):
cube, meta = load_cube(cube_path, data_file_override=data_file_override)
cube = ensure_rows_cols_bands(cube, shape_hint)
rows, cols, bands = cube.shape
print(f"Cubo caricato: {rows} righe x {cols} colonne x {bands} bande")
mask = nodata_mask(cube, nodata_thresh)
mean_profile, valid_count = crosstrack_profile_per_band(cube, mask) # (cols, bands)
x = np.arange(cols)
wavelengths = meta.get("wavelengths")
if wavelengths is not None and len(wavelengths) == bands:
band_axis = wavelengths
band_label = "Lunghezza d'onda (nm)"
else:
band_axis = np.arange(bands)
band_label = "Indice banda"
# --- Fit lineare per banda: slope, intercept, r2 ---
slopes = np.full(bands, np.nan)
r2s = np.full(bands, np.nan)
for b in range(bands):
s, i, r2 = linear_fit(x, mean_profile[:, b])
slopes[b] = s
r2s[b] = r2
# variazione percentuale bordo-bordo per banda, rispetto al valore medio globale della banda
global_mean_per_band = np.nanmean(mean_profile, axis=0)
variation_total = slopes * (cols - 1)
variation_pct = 100 * variation_total / np.where(global_mean_per_band != 0, global_mean_per_band, np.nan)
# --- Confronto meta' sinistra vs destra per banda ---
split_col = int(cols * split_fraction)
left_mask = mask.copy()
left_mask[:, split_col:] = False
right_mask = mask.copy()
right_mask[:, :split_col] = False
left_mean = np.array([cube[left_mask, b].mean() if left_mask.any() else np.nan for b in range(bands)])
right_mean = np.array([cube[right_mask, b].mean() if right_mask.any() else np.nan for b in range(bands)])
diff_pct = 100 * (right_mean - left_mean) / np.where(left_mean != 0, left_mean, np.nan)
# --- Report testuale ---
out_dir = Path(out_dir)
out_dir.mkdir(parents=True, exist_ok=True)
lines = []
lines.append(f"Cubo: {cube_path}")
lines.append(f"Dimensioni: {rows} righe x {cols} colonne x {bands} bande")
lines.append(f"Soglia nodata: {nodata_thresh}")
lines.append("")
lines.append("Riepilogo globale (mediato su tutte le bande):")
lines.append(f" Slope medio: {np.nanmean(slopes):.6f} per colonna")
lines.append(f" Variazione media bordo-bordo: {np.nanmean(variation_total):+.4f} ({np.nanmean(variation_pct):+.2f}% rispetto alla media di banda)")
lines.append(f" Differenza media destra-sinistra: {np.nanmean(diff_pct):+.2f}%")
lines.append("")
n_bands_show = min(bands, 30)
lines.append(f"Dettaglio per banda (prime {n_bands_show} mostrate nel report; tutte le bande nel file .csv):")
lines.append(f"{'banda':>6} {'asse':>10} {'slope':>12} {'R^2':>8} {'var_tot':>10} {'var_%':>8} {'diff_dx-sx_%':>14}")
for b in range(n_bands_show):
lines.append(f"{b:6d} {band_axis[b]:10.2f} {slopes[b]:12.6f} {r2s[b]:8.3f} {variation_total[b]:10.3f} {variation_pct[b]:8.2f} {diff_pct[b]:14.2f}")
if bands > n_bands_show:
lines.append(f"... (+{bands - n_bands_show} bande, vedi crosstrack_per_band_stats.csv)")
report_text = "\n".join(lines)
print(report_text)
(out_dir / "report.txt").write_text(report_text, encoding="utf-8")
# CSV completo per tutte le bande
csv_path = out_dir / "crosstrack_per_band_stats.csv"
with open(csv_path, "w", encoding="utf-8") as f:
f.write("band_index,band_axis,slope,r2,variation_total,variation_pct,diff_right_minus_left_pct,global_mean\n")
for b in range(bands):
f.write(f"{b},{band_axis[b]},{slopes[b]},{r2s[b]},{variation_total[b]},{variation_pct[b]},{diff_pct[b]},{global_mean_per_band[b]}\n")
print(f"CSV per-banda salvato in {csv_path}")
# ---------------------------------------------------------------
# Figura 1: heatmap banda x colonna, normalizzata per banda
# (ogni riga della heatmap divisa per la propria media, cosi' il
# confronto tra bande con intensita' molto diverse resta leggibile)
# ---------------------------------------------------------------
norm_profile = mean_profile / np.where(global_mean_per_band != 0, global_mean_per_band, np.nan)[np.newaxis, :]
fig, ax = plt.subplots(figsize=(10, 6))
im = ax.imshow(
norm_profile.T,
aspect="auto",
origin="lower",
extent=[0, cols, band_axis[0], band_axis[-1]],
cmap="RdBu_r",
vmin=0.7, vmax=1.3,
)
ax.set_xlabel("Colonna (direzione across-track)")
ax.set_ylabel(band_label)
ax.set_title("Profilo across-track normalizzato per banda\n(rosso=piu' luminoso della media di banda, blu=piu' scuro)")
fig.colorbar(im, ax=ax, label="valore / media di banda")
fig.tight_layout()
fig.savefig(out_dir / "heatmap_band_vs_column.png", dpi=150)
plt.close(fig)
# ---------------------------------------------------------------
# Figura 2: slope (gradiente) in funzione della banda
# ---------------------------------------------------------------
fig, axes = plt.subplots(2, 1, figsize=(9, 7), sharex=True)
axes[0].plot(band_axis, variation_pct, color="darkred")
axes[0].axhline(0, color="gray", linestyle=":")
axes[0].set_ylabel("Variazione bordo-bordo (%)")
axes[0].set_title("Gradiente across-track in funzione della banda")
axes[0].grid(alpha=0.3)
axes[1].plot(band_axis, diff_pct, color="darkblue")
axes[1].axhline(0, color="gray", linestyle=":")
axes[1].set_ylabel("Diff. destra-sinistra (%)")
axes[1].set_xlabel(band_label)
axes[1].grid(alpha=0.3)
fig.tight_layout()
fig.savefig(out_dir / "gradient_vs_wavelength.png", dpi=150)
plt.close(fig)
# ---------------------------------------------------------------
# Figura 3: profili across-track per alcune bande rappresentative
# ---------------------------------------------------------------
sample_idx = sorted(set([0, bands // 4, bands // 2, (3 * bands) // 4, bands - 1]))
fig, ax = plt.subplots(figsize=(9, 5))
cmap = plt.get_cmap("viridis")
for k, b in enumerate(sample_idx):
ax.plot(x, mean_profile[:, b], color=cmap(k / max(1, len(sample_idx) - 1)),
label=f"banda {b} ({band_axis[b]:.0f}{'nm' if band_label.startswith('Lunghezza') else ''})")
ax.axvline(split_col, color="gray", linestyle=":", label="split 50%")
ax.set_xlabel("Colonna (direzione across-track)")
ax.set_ylabel("Valore medio (DN o riflettanza)")
ax.set_title("Profilo across-track per bande rappresentative")
ax.legend(fontsize=8)
ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig(out_dir / "sample_bands_profile.png", dpi=150)
plt.close(fig)
# ---------------------------------------------------------------
# Correzione flat-field per colonna, per banda (opzionale)
# ---------------------------------------------------------------
if save_corrected:
print("Applico correzione flat-field per colonna (per banda)...")
correction = np.where(mean_profile > 1e-9, global_mean_per_band[np.newaxis, :] / mean_profile, 1.0)
correction = np.clip(correction, 0.3, 3.0) # evita amplificazioni estreme ai bordi rumorosi
corrected = cube * correction[np.newaxis, :, :]
corrected[~mask] = 0
if output_format == "envi":
corrected_path = export_envi_cube(out_dir / "cube_flatfield_corrected", corrected, meta,
interleave="bil")
print(f"Cubo corretto salvato in formato ENVI: {corrected_path} (+ file dati associato)")
else:
corrected_path = out_dir / "cube_flatfield_corrected.npy"
np.save(corrected_path, corrected)
print(f"Cubo corretto salvato in formato numpy: {corrected_path}")
# Verifica: ricalcola il profilo sul cubo corretto e stampa il nuovo slope medio
mean_profile_corr, _ = crosstrack_profile_per_band(corrected, mask)
slopes_corr = np.full(bands, np.nan)
for b in range(bands):
s, i, r2 = linear_fit(x, mean_profile_corr[:, b])
slopes_corr[b] = s
print(f"Slope medio PRIMA della correzione: {np.nanmean(slopes):.6f}")
print(f"Slope medio DOPO la correzione: {np.nanmean(slopes_corr):.6f}")
print(f"\nOutput salvati in: {out_dir.resolve()}")
print(" - report.txt")
print(" - crosstrack_per_band_stats.csv")
print(" - heatmap_band_vs_column.png")
print(" - gradient_vs_wavelength.png")
print(" - sample_bands_profile.png")
if save_corrected:
if output_format == "envi":
print(" - cube_flatfield_corrected.hdr (+ file dati binario associato)")
else:
print(" - cube_flatfield_corrected.npy")
def main():
parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument("input", help="Cubo iperspettrale: file .hdr (ENVI), .npy o .npz")
parser.add_argument("--out", default="crosstrack_cube_analysis", help="Cartella di output")
parser.add_argument("--shape", default=None,
help="Solo se il .npy non e' gia' (rows,cols,bands): 'rows,cols,bands'")
parser.add_argument("--nodata-thresh", type=float, default=1e-6,
help="Soglia sotto la quale un pixel e' nodata (default 1e-6)")
parser.add_argument("--split-fraction", type=float, default=0.5,
help="Frazione di colonne che definisce il confine sinistra/destra")
parser.add_argument("--save-corrected", action="store_true",
help="Salva anche il cubo corretto (flat-field per colonna, per banda) come .npy")
parser.add_argument("--data-file", default=None,
help="Percorso esplicito del file binario ENVI associato all'header, "
"da usare se l'auto-rilevamento fallisce (es. estensioni non standard come .spe)")
parser.add_argument("--output-format", choices=["envi", "npy"], default="envi",
help="Formato del cubo corretto salvato con --save-corrected (default: envi)")
args = parser.parse_args()
shape_hint = None
if args.shape:
shape_hint = tuple(int(v) for v in args.shape.split(","))
if len(shape_hint) != 3:
print("ERRORE: --shape deve essere 'rows,cols,bands'", file=sys.stderr)
sys.exit(1)
analyze(args.input, args.out, shape_hint, args.nodata_thresh,
args.split_fraction, args.save_corrected, data_file_override=args.data_file,
output_format=args.output_format)
if __name__ == "__main__":
main()
#!/usr/bin/env python3
"""
analyze_crosstrack_cube.py
Analizza un cubo iperspettrale pushbroom (formato ENVI: .hdr + .bil/.bsq/.bip,
oppure un array numpy .npy/.npz) per individuare e quantificare differenze di
illuminazione lungo la direzione across-track (le colonne del sensore),
banda per banda.
Perche' questo script serve
---------------------------
Nel true-color renderizzato avevi notato una fascia destra piu' scura della
sinistra. Su un solo composito RGB non si puo' distinguere se il problema
nasce:
(a) nel sensore/ottica (vignettatura radiometrica, gia' presente nella
riflettanza calibrata), oppure
(b) nel rendering RGB (stretch, white balance, ecc.), oppure
(c) da un vero effetto fisico di scena (BRDF/anisotropia, ombre, angolo
di illuminazione solare rispetto alla direzione di volo).
Lavorando sul cubo (idealmente gia' calibrato in riflettanza con dark/white/
pannello) e guardando il profilo across-track banda per banda, si puo' capire
se il gradiente e' costante su tutte le lunghezze d'onda (probabile
vignettatura/ottica del sensore, o dark/white non uniformi lungo lo swath)
oppure cambia con la banda (piu' probabile un effetto fisico/spettrale di
scena o un residuo di calibrazione spettrale).
Cosa fa lo script
-----------------
1. Carica il cubo (rows=along-track, cols=across-track/swath, bands=spettrale).
Auto-rileva l'interleave ENVI (BIL tipico per pushbroom) tramite 'spectral'.
2. Maschera i pixel nodata (es. tutti zero o NaN, o sotto una soglia).
3. Per ogni banda calcola il profilo medio/mediano per colonna (across-track).
4. Quantifica il gradiente sinistra/destra per banda con un fit lineare e
con il confronto diretto delle due meta'.
5. Produce:
- una heatmap (banda x colonna) del profilo across-track normalizzato,
per vedere a colpo d'occhio se il gradiente e' presente su tutte le bande
o solo in alcune (utile per capire la causa fisica),
- un grafico dello slope (pendenza) del gradiente in funzione della
lunghezza d'onda/banda,
- un grafico di alcune bande rappresentative (es. blu/verde/rosso/NIR se
riconoscibili, altrimenti prima/centro/ultima banda),
- un report testuale banda per banda,
- opzionalmente un cubo corretto (flat-field per colonna, per banda) in
formato .npy, utile per verificare se il problema sparisce.
Uso
---
# Cubo ENVI (serve il file .hdr, il dato binario associato viene trovato
# automaticamente in base al campo 'file' nell'header o stesso nome)
python3 analyze_crosstrack_cube.py cubo.hdr --out out_cube
# Cubo numpy gia' in memoria come array (righe, colonne, bande)
python3 analyze_crosstrack_cube.py cubo.npy --out out_cube --shape rows,cols,bands
# Salva anche il cubo corretto
python3 analyze_crosstrack_cube.py cubo.hdr --out out_cube --save-corrected
"""
import argparse
import sys
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
# --------------------------------------------------------------------------
# Caricamento cubo
# --------------------------------------------------------------------------
def _parse_envi_header_datafile(hdr_path):
"""
Legge il campo 'file =' / 'data file =' dall'header ENVI, se presente,
per capire come si chiama il file binario associato. Ritorna None se non
trovato o non parsabile.
"""
try:
text = Path(hdr_path).read_text(encoding="utf-8", errors="ignore")
except Exception:
return None
for line in text.splitlines():
low = line.strip().lower()
if low.startswith("file =") or low.startswith("data file ="):
value = line.split("=", 1)[1].strip()
value = value.strip("{}").strip()
if value:
return value
return None
def _find_envi_data_file(hdr_path):
"""
Cerca il file binario associato a un header ENVI quando 'spectral' non
riesce a trovarlo automaticamente (lo fa solo per estensioni .img, .dat,
.sli o nessuna estensione). Strategie, in ordine:
1. Campo 'file =' dentro l'header stesso.
2. Qualunque file nella stessa cartella con lo stesso nome base
(stesso nome del .hdr senza estensione) e estensione diversa da
.hdr, .png, .jpg, .txt, .csv, .py (per evitare falsi positivi).
Ritorna il Path del file trovato, o None.
"""
hdr_path = Path(hdr_path)
base = hdr_path.with_suffix("") # rimuove .hdr
declared = _parse_envi_header_datafile(hdr_path)
if declared:
candidate = (hdr_path.parent / declared)
if candidate.exists():
return candidate
candidate2 = hdr_path.parent / Path(declared).name
if candidate2.exists():
return candidate2
exclude_suffixes = {".hdr", ".png", ".jpg", ".jpeg", ".txt", ".csv", ".py", ".md"}
candidates = [p for p in hdr_path.parent.glob(base.name + "*")
if p.is_file() and p.suffix.lower() not in exclude_suffixes and p != hdr_path]
if candidates:
candidates.sort(key=lambda p: p.stat().st_size, reverse=True)
return candidates[0]
return None
def load_cube(path, data_file_override=None):
"""
Carica un cubo iperspettrale da:
- file ENVI header (.hdr) -> usa la libreria 'spectral'
- file binario ENVI (qualunque estensione, es. .spe, .bil, .bsq, .bip,
.raw, .dat...) che ha un .hdr associato nella stessa cartella
- file numpy (.npy / .npz)
Ritorna array numpy (rows, cols, bands) in float64, e un dict di metadati
(wavelengths se disponibili, nomi banda, ecc.)
"""
path = Path(path)
meta = {"wavelengths": None, "band_names": None, "source": str(path)}
if path.suffix.lower() == ".npy":
cube = np.load(str(path))
return cube.astype(np.float64), meta
if path.suffix.lower() == ".npz":
data = np.load(str(path))
for key in ("cube", "data", "arr_0"):
if key in data:
return data[key].astype(np.float64), meta
first_key = list(data.keys())[0]
return data[first_key].astype(np.float64), meta
hdr_path = None
data_path = Path(data_file_override) if data_file_override else None
if path.suffix.lower() == ".hdr":
hdr_path = path
else:
candidate_hdrs = [
path.with_suffix(".hdr"),
Path(str(path) + ".hdr"),
]
for c in candidate_hdrs:
if c.exists():
hdr_path = c
if data_path is None:
data_path = path
break
if hdr_path is None or not hdr_path.exists():
raise ValueError(
f"Formato non riconosciuto o header ENVI non trovato per {path}. "
f"Usa un file .hdr (ENVI), il file binario con .hdr associato, oppure .npy/.npz"
)
import spectral
if data_path is None:
try:
img = spectral.envi.open(str(hdr_path))
except spectral.io.envi.EnviDataFileNotFoundError:
found = _find_envi_data_file(hdr_path)
if found is None:
raise ValueError(
f"Impossibile trovare il file dati binario associato a {hdr_path}. "
f"Passa esplicitamente il file dati con --data-file, oppure rinomina/"
f"copia il binario accanto all'header con estensione .img/.dat/.bin."
)
print(f"[info] File dati non auto-rilevato da 'spectral'; uso: {found}")
img = spectral.envi.open(str(hdr_path), image=str(found))
else:
img = spectral.envi.open(str(hdr_path), image=str(data_path))
cube = img.load().astype(np.float64)
wl = getattr(img, "bands", None)
if wl is not None and getattr(wl, "centers", None) is not None:
meta["wavelengths"] = np.array(wl.centers, dtype=float)
# conserviamo l'header originale (dict) e l'interleave, utili per
# riesportare il cubo corretto in formato ENVI mantenendo wavelength,
# unita', sensor type, ecc.
meta["envi_header"] = dict(getattr(img, "metadata", {}) or {})
meta["interleave"] = getattr(img, "interleave", None)
return np.asarray(cube), meta
def ensure_rows_cols_bands(cube, shape_hint=None):
"""
Se il cubo non e' gia' (rows, cols, bands), e l'utente ha fornito
--shape rows,cols,bands, lo reshape. Altrimenti assume che sia gia'
nell'ordine corretto (che e' il default per 'spectral' e per la
convenzione usata nello script true-color della sessione precedente).
"""
if shape_hint is not None:
rows, cols, bands = shape_hint
cube = cube.reshape(rows, cols, bands)
return cube
# --------------------------------------------------------------------------
# Maschera nodata
# --------------------------------------------------------------------------
def nodata_mask(cube, thresh=1e-6):
"""
True dove il pixel e' valido. Un pixel e' nodata se e' NaN su tutte le
bande oppure se il valore assoluto medio su tutte le bande e' sotto
soglia (tipico bordo nero da rotazione/georeferenziazione).
"""
finite = np.isfinite(cube).all(axis=2)
mean_abs = np.nanmean(np.abs(cube), axis=2)
valid = finite & (mean_abs > thresh)
return valid
# --------------------------------------------------------------------------
# Profilo across-track per banda
# --------------------------------------------------------------------------
def crosstrack_profile_per_band(cube, mask):
"""
Ritorna:
mean_profile: (cols, bands) media per colonna e banda (NaN se nessun
pixel valido in quella colonna)
valid_count: (cols,) numero di pixel validi per colonna
"""
rows, cols, bands = cube.shape
mean_profile = np.full((cols, bands), np.nan)
valid_count = mask.sum(axis=0)
for c in range(cols):
col_mask = mask[:, c]
if col_mask.sum() == 0:
continue
mean_profile[c, :] = cube[col_mask, c, :].mean(axis=0)
return mean_profile, valid_count
def linear_fit(x, y):
good = np.isfinite(y)
if good.sum() < 2:
return np.nan, np.nan, np.nan
xs, ys = x[good], y[good]
A = np.vstack([xs, np.ones_like(xs)]).T
slope, intercept = np.linalg.lstsq(A, ys, rcond=None)[0]
pred = slope * xs + intercept
ss_res = np.sum((ys - pred) ** 2)
ss_tot = np.sum((ys - ys.mean()) ** 2)
r2 = 1 - ss_res / ss_tot if ss_tot > 0 else np.nan
return slope, intercept, r2
def export_envi_cube(out_path, cube, meta, interleave="bil"):
"""
Salva un cubo (rows, cols, bands) in formato ENVI (coppia .hdr + file dati),
riusando i metadati originali (wavelength, unita', sensor type, ecc.) se
disponibili in meta['envi_header'], cosi' il cubo esportato resta
apribile in ENVI/QGIS/altri tool con le stesse informazioni spettrali
dell'originale.
`out_path` puo' essere passato con o senza estensione: verra' sempre
scritto come coppia <out_path>.hdr + <out_path>.dat (interleave BIL di
default, tipico per dati pushbroom).
"""
import spectral
out_path = Path(out_path)
if out_path.suffix.lower() == ".hdr":
out_path = out_path.with_suffix("")
orig_header = dict(meta.get("envi_header") or {})
# Ripulisce campi che 'spectral' ricalcola da solo o che diventerebbero
# inconsistenti (dimensioni, nome file, offset) per evitare header corrotti
for key in ("samples", "lines", "bands", "header offset", "file type",
"byte order", "data type", "interleave"):
orig_header.pop(key, None)
metadata = orig_header # contiene ad es. wavelength, wavelength units, sensor type, map info...
spectral.envi.save_image(
str(out_path.with_suffix(".hdr")),
cube.astype(np.float32),
dtype=np.float32,
force=True,
interleave=interleave,
metadata=metadata,
)
return out_path.with_suffix(".hdr")
# --------------------------------------------------------------------------
# Analisi principale
# --------------------------------------------------------------------------
def analyze(cube_path, out_dir, shape_hint=None, nodata_thresh=1e-6,
split_fraction=0.5, save_corrected=False, data_file_override=None,
output_format="envi"):
cube, meta = load_cube(cube_path, data_file_override=data_file_override)
cube = ensure_rows_cols_bands(cube, shape_hint)
rows, cols, bands = cube.shape
print(f"Cubo caricato: {rows} righe x {cols} colonne x {bands} bande")
mask = nodata_mask(cube, nodata_thresh)
mean_profile, valid_count = crosstrack_profile_per_band(cube, mask) # (cols, bands)
x = np.arange(cols)
wavelengths = meta.get("wavelengths")
if wavelengths is not None and len(wavelengths) == bands:
band_axis = wavelengths
band_label = "Lunghezza d'onda (nm)"
else:
band_axis = np.arange(bands)
band_label = "Indice banda"
# --- Fit lineare per banda: slope, intercept, r2 ---
slopes = np.full(bands, np.nan)
r2s = np.full(bands, np.nan)
for b in range(bands):
s, i, r2 = linear_fit(x, mean_profile[:, b])
slopes[b] = s
r2s[b] = r2
# variazione percentuale bordo-bordo per banda, rispetto al valore medio globale della banda
global_mean_per_band = np.nanmean(mean_profile, axis=0)
variation_total = slopes * (cols - 1)
variation_pct = 100 * variation_total / np.where(global_mean_per_band != 0, global_mean_per_band, np.nan)
# --- Confronto meta' sinistra vs destra per banda ---
split_col = int(cols * split_fraction)
left_mask = mask.copy()
left_mask[:, split_col:] = False
right_mask = mask.copy()
right_mask[:, :split_col] = False
left_mean = np.array([cube[left_mask, b].mean() if left_mask.any() else np.nan for b in range(bands)])
right_mean = np.array([cube[right_mask, b].mean() if right_mask.any() else np.nan for b in range(bands)])
diff_pct = 100 * (right_mean - left_mean) / np.where(left_mean != 0, left_mean, np.nan)
# --- Report testuale ---
out_dir = Path(out_dir)
out_dir.mkdir(parents=True, exist_ok=True)
lines = []
lines.append(f"Cubo: {cube_path}")
lines.append(f"Dimensioni: {rows} righe x {cols} colonne x {bands} bande")
lines.append(f"Soglia nodata: {nodata_thresh}")
lines.append("")
lines.append("Riepilogo globale (mediato su tutte le bande):")
lines.append(f" Slope medio: {np.nanmean(slopes):.6f} per colonna")
lines.append(f" Variazione media bordo-bordo: {np.nanmean(variation_total):+.4f} ({np.nanmean(variation_pct):+.2f}% rispetto alla media di banda)")
lines.append(f" Differenza media destra-sinistra: {np.nanmean(diff_pct):+.2f}%")
lines.append("")
n_bands_show = min(bands, 30)
lines.append(f"Dettaglio per banda (prime {n_bands_show} mostrate nel report; tutte le bande nel file .csv):")
lines.append(f"{'banda':>6} {'asse':>10} {'slope':>12} {'R^2':>8} {'var_tot':>10} {'var_%':>8} {'diff_dx-sx_%':>14}")
for b in range(n_bands_show):
lines.append(f"{b:6d} {band_axis[b]:10.2f} {slopes[b]:12.6f} {r2s[b]:8.3f} {variation_total[b]:10.3f} {variation_pct[b]:8.2f} {diff_pct[b]:14.2f}")
if bands > n_bands_show:
lines.append(f"... (+{bands - n_bands_show} bande, vedi crosstrack_per_band_stats.csv)")
report_text = "\n".join(lines)
print(report_text)
(out_dir / "report.txt").write_text(report_text, encoding="utf-8")
# CSV completo per tutte le bande
csv_path = out_dir / "crosstrack_per_band_stats.csv"
with open(csv_path, "w", encoding="utf-8") as f:
f.write("band_index,band_axis,slope,r2,variation_total,variation_pct,diff_right_minus_left_pct,global_mean\n")
for b in range(bands):
f.write(f"{b},{band_axis[b]},{slopes[b]},{r2s[b]},{variation_total[b]},{variation_pct[b]},{diff_pct[b]},{global_mean_per_band[b]}\n")
print(f"CSV per-banda salvato in {csv_path}")
# ---------------------------------------------------------------
# Figura 1: heatmap banda x colonna, normalizzata per banda
# (ogni riga della heatmap divisa per la propria media, cosi' il
# confronto tra bande con intensita' molto diverse resta leggibile)
# ---------------------------------------------------------------
norm_profile = mean_profile / np.where(global_mean_per_band != 0, global_mean_per_band, np.nan)[np.newaxis, :]
fig, ax = plt.subplots(figsize=(10, 6))
im = ax.imshow(
norm_profile.T,
aspect="auto",
origin="lower",
extent=[0, cols, band_axis[0], band_axis[-1]],
cmap="RdBu_r",
vmin=0.7, vmax=1.3,
)
ax.set_xlabel("Colonna (direzione across-track)")
ax.set_ylabel(band_label)
ax.set_title("Profilo across-track normalizzato per banda\n(rosso=piu' luminoso della media di banda, blu=piu' scuro)")
fig.colorbar(im, ax=ax, label="valore / media di banda")
fig.tight_layout()
fig.savefig(out_dir / "heatmap_band_vs_column.png", dpi=150)
plt.close(fig)
# ---------------------------------------------------------------
# Figura 2: slope (gradiente) in funzione della banda
# ---------------------------------------------------------------
fig, axes = plt.subplots(2, 1, figsize=(9, 7), sharex=True)
axes[0].plot(band_axis, variation_pct, color="darkred")
axes[0].axhline(0, color="gray", linestyle=":")
axes[0].set_ylabel("Variazione bordo-bordo (%)")
axes[0].set_title("Gradiente across-track in funzione della banda")
axes[0].grid(alpha=0.3)
axes[1].plot(band_axis, diff_pct, color="darkblue")
axes[1].axhline(0, color="gray", linestyle=":")
axes[1].set_ylabel("Diff. destra-sinistra (%)")
axes[1].set_xlabel(band_label)
axes[1].grid(alpha=0.3)
fig.tight_layout()
fig.savefig(out_dir / "gradient_vs_wavelength.png", dpi=150)
plt.close(fig)
# ---------------------------------------------------------------
# Figura 3: profili across-track per alcune bande rappresentative
# ---------------------------------------------------------------
sample_idx = sorted(set([0, bands // 4, bands // 2, (3 * bands) // 4, bands - 1]))
fig, ax = plt.subplots(figsize=(9, 5))
cmap = plt.get_cmap("viridis")
for k, b in enumerate(sample_idx):
ax.plot(x, mean_profile[:, b], color=cmap(k / max(1, len(sample_idx) - 1)),
label=f"banda {b} ({band_axis[b]:.0f}{'nm' if band_label.startswith('Lunghezza') else ''})")
ax.axvline(split_col, color="gray", linestyle=":", label="split 50%")
ax.set_xlabel("Colonna (direzione across-track)")
ax.set_ylabel("Valore medio (DN o riflettanza)")
ax.set_title("Profilo across-track per bande rappresentative")
ax.legend(fontsize=8)
ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig(out_dir / "sample_bands_profile.png", dpi=150)
plt.close(fig)
# ---------------------------------------------------------------
# Correzione flat-field per colonna, per banda (opzionale)
# ---------------------------------------------------------------
if save_corrected:
print("Applico correzione flat-field per colonna (per banda)...")
correction = np.where(mean_profile > 1e-9, global_mean_per_band[np.newaxis, :] / mean_profile, 1.0)
correction = np.clip(correction, 0.3, 3.0) # evita amplificazioni estreme ai bordi rumorosi
corrected = cube * correction[np.newaxis, :, :]
corrected[~mask] = 0
if output_format == "envi":
corrected_path = export_envi_cube(out_dir / "cube_flatfield_corrected", corrected, meta,
interleave="bil")
print(f"Cubo corretto salvato in formato ENVI: {corrected_path} (+ file dati associato)")
else:
corrected_path = out_dir / "cube_flatfield_corrected.npy"
np.save(corrected_path, corrected)
print(f"Cubo corretto salvato in formato numpy: {corrected_path}")
# Verifica: ricalcola il profilo sul cubo corretto e stampa il nuovo slope medio
mean_profile_corr, _ = crosstrack_profile_per_band(corrected, mask)
slopes_corr = np.full(bands, np.nan)
for b in range(bands):
s, i, r2 = linear_fit(x, mean_profile_corr[:, b])
slopes_corr[b] = s
print(f"Slope medio PRIMA della correzione: {np.nanmean(slopes):.6f}")
print(f"Slope medio DOPO la correzione: {np.nanmean(slopes_corr):.6f}")
print(f"\nOutput salvati in: {out_dir.resolve()}")
print(" - report.txt")
print(" - crosstrack_per_band_stats.csv")
print(" - heatmap_band_vs_column.png")
print(" - gradient_vs_wavelength.png")
print(" - sample_bands_profile.png")
if save_corrected:
if output_format == "envi":
print(" - cube_flatfield_corrected.hdr (+ file dati binario associato)")
else:
print(" - cube_flatfield_corrected.npy")
def main():
parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument("input", help="Cubo iperspettrale: file .hdr (ENVI), .npy o .npz")
parser.add_argument("--out", default="crosstrack_cube_analysis", help="Cartella di output")
parser.add_argument("--shape", default=None,
help="Solo se il .npy non e' gia' (rows,cols,bands): 'rows,cols,bands'")
parser.add_argument("--nodata-thresh", type=float, default=1e-6,
help="Soglia sotto la quale un pixel e' nodata (default 1e-6)")
parser.add_argument("--split-fraction", type=float, default=0.5,
help="Frazione di colonne che definisce il confine sinistra/destra")
parser.add_argument("--save-corrected", action="store_true",
help="Salva anche il cubo corretto (flat-field per colonna, per banda) come .npy")
parser.add_argument("--data-file", default=None,
help="Percorso esplicito del file binario ENVI associato all'header, "
"da usare se l'auto-rilevamento fallisce (es. estensioni non standard come .spe)")
parser.add_argument("--output-format", choices=["envi", "npy"], default="envi",
help="Formato del cubo corretto salvato con --save-corrected (default: envi)")
args = parser.parse_args()
shape_hint = None
if args.shape:
shape_hint = tuple(int(v) for v in args.shape.split(","))
if len(shape_hint) != 3:
print("ERRORE: --shape deve essere 'rows,cols,bands'", file=sys.stderr)
sys.exit(1)
analyze(args.input, args.out, shape_hint, args.nodata_thresh,
args.split_fraction, args.save_corrected, data_file_override=args.data_file,
output_format=args.output_format)
if __name__ == "__main__":
main()
