Stima Smile in FigSpec FS60-CL

Per stimare lo smile effect in FigSpec FS60-CL e' stato usato il file raw dei digital numbers cercando di isolare lo spostamento  dell minimo dell'assorbimento dell'Ossigeno a 760 nm

 

 

prima di tutto si vede che il grafico non e' simmetrico indice di problemi a livello ottico. In generale pero' lo smile e' stimato inferiore ad 1 banda (indicativamente 0.2) quindi e' trascurabile

L'aspetto da sottilineare (non legata allo smile) e' che sembra ci sia un offset di calibrazione spettrale di circa 2 nm (l'assorbimento reale e' di letteratura a 760 nm mentre il FigSpec lo vede a 762 nm) 

python3 detect_smile.py figspec.hdr figspec.spe --out smile_diagnostic.png

 

#!/usr/bin/env python3

"""

Estimates spectral smile by tracking the position of a known, narrow

atmospheric absorption feature (default: O2-A band ~760-762nm) across the

spatial (sample/across-track) dimension of an ENVI BIL hyperspectral cube.

If the sensor has no smile, the absorption minimum should sit at the same

wavelength for every column. If it drifts smoothly from one edge to the

other, that drift IS the smile curve.

Usage:

python3 detect_smile.py figspec_reflectance_rid_savgol.hdr \

figspec_reflectance_rid_savgol.img \

--out smile_diagnostic.png

Optional: --feature-nm 760 --search-window-nm 25 (O2-A, default)

--feature-nm 940 --search-window-nm 40 (water vapor, cross-check)

"""

import os

import re

import argparse

import numpy as np

def parse_hdr(hdr_path):

txt = open(hdr_path, "r", errors="ignore").read()

def get_int(key):

m = re.search(rf"{key}\s*=\s*(\d+)", txt)

return int(m.group(1))

samples = get_int("samples")

lines = get_int("lines")

bands = get_int("bands")

dtype_code = get_int("data type")

m = re.search(r"interleave\s*=\s*(\w+)", txt, re.IGNORECASE)

interleave = m.group(1).lower()

m = re.search(r"byte order\s*=\s*(\d+)", txt)

byte_order = int(m.group(1)) if m else 0

m = re.search(r"wavelength\s*=\s*\{([^}]+)\}", txt, re.IGNORECASE)

wavelengths = [float(x) for x in m.group(1).split(",")] if m else None

dtype_map = {1: np.uint8, 2: np.int16, 3: np.int32, 4: np.float32,

5: np.float64, 12: np.uint16, 13: np.uint32,

14: np.int64, 15: np.uint64}

np_dtype = np.dtype(dtype_map[dtype_code])

np_dtype = np_dtype.newbyteorder(">" if byte_order == 1 else "<")

return dict(samples=samples, lines=lines, bands=bands, dtype=np_dtype,

interleave=interleave, wavelengths=np.array(wavelengths))

def main():

ap = argparse.ArgumentParser()

ap.add_argument("hdr_path")

ap.add_argument("data_path")

ap.add_argument("--out", default="smile_diagnostic.png")

ap.add_argument("--feature-nm", type=float, default=760.0,

help="Center wavelength of the known absorption feature (default: O2-A 760nm)")

ap.add_argument("--search-window-nm", type=float, default=25.0,

help="Half-width of the search window around feature-nm")

ap.add_argument("--line-stride", type=int, default=5,

help="Use every Nth line when averaging (speed/memory tradeoff)")

args = ap.parse_args()

meta = parse_hdr(args.hdr_path)

samples, lines, bands = meta["samples"], meta["lines"], meta["bands"]

dtype = meta["dtype"]

interleave = meta["interleave"]

wavelengths = meta["wavelengths"]

if interleave != "bil":

raise ValueError("This script currently assumes BIL interleave.")

cube = np.memmap(args.data_path, dtype=dtype, mode="r",

shape=(lines, bands, samples))

print(f"Averaging every {args.line_stride}th line "

f"({lines // args.line_stride} lines used) to build a per-column mean spectrum...")

idx_lines = np.arange(0, lines, args.line_stride)

accum = np.zeros((bands, samples), dtype=np.float64)

for i, li in enumerate(idx_lines):

accum += cube[li, :, :]

if i % 200 == 0:

print(f" {i}/{len(idx_lines)}", end="\r")

mean_spectrum = (accum / len(idx_lines)).astype(np.float32) # (bands, samples)

print(f"\nDone averaging. Shape: {mean_spectrum.shape}")

lo_nm = args.feature_nm - args.search_window_nm

hi_nm = args.feature_nm + args.search_window_nm

band_mask = (wavelengths >= lo_nm) & (wavelengths <= hi_nm)

band_idxs = np.where(band_mask)[0]

if len(band_idxs) < 5:

raise ValueError(f"Too few bands in [{lo_nm},{hi_nm}]nm window; "

f"check --feature-nm/--search-window-nm vs your wavelength range.")

print(f"Search window: bands {band_idxs[0]}-{band_idxs[-1]} "

f"({wavelengths[band_idxs[0]]:.1f}-{wavelengths[band_idxs[-1]]:.1f} nm)")

sub_wavelengths = np.full(samples, np.nan)

for col in range(samples):

window_vals = mean_spectrum[band_idxs, col].astype(np.float64)

if np.all(window_vals <= 1e-6):

continue # degenerate/empty column

local_min = np.argmin(window_vals)

# sub-pixel parabolic refinement (skip if at window edge)

if 0 < local_min < len(window_vals) - 1:

y0, y1, y2 = window_vals[local_min - 1], window_vals[local_min], window_vals[local_min + 1]

denom = (y0 - 2 * y1 + y2)

delta = 0.5 * (y0 - y2) / denom if abs(denom) > 1e-9 else 0.0

delta = np.clip(delta, -1, 1)

else:

delta = 0.0

band_global = band_idxs[local_min]

# interpolate wavelength at sub-pixel band position

if 0 <= band_global + int(np.sign(delta)) < bands:

w0 = wavelengths[band_global]

w_next = wavelengths[band_global + (1 if delta >= 0 else -1)]

sub_wl = w0 + delta * (w_next - w0)

else:

sub_wl = wavelengths[band_global]

sub_wavelengths[col] = sub_wl

valid = ~np.isnan(sub_wavelengths)

cols = np.arange(samples)[valid]

vals = sub_wavelengths[valid]

if len(vals) < 10:

print("Too few valid columns to estimate smile reliably.")

return

# quadratic fit (typical smile shape) for a smooth summary curve

coeffs = np.polyfit(cols, vals, 2)

fit_curve = np.polyval(coeffs, cols)

ptp_measured = vals.max() - vals.min()

ptp_fit = fit_curve.max() - fit_curve.min()

print(f"\nAbsorption feature position across {samples} columns:")

print(f" raw measured range: {vals.min():.2f} - {vals.max():.2f} nm "

f"(peak-to-valley = {ptp_measured:.2f} nm)")

print(f" smooth quadratic fit range: {fit_curve.min():.2f} - {fit_curve.max():.2f} nm "

f"(peak-to-valley = {ptp_fit:.2f} nm)")

print(f" quadratic fit coefficients (a*x^2+b*x+c): {coeffs}")

band_spacing_nm = np.mean(np.diff(wavelengths))

print(f" approx. band spacing: {band_spacing_nm:.2f} nm/band")

print(f" smile in band units: ~{ptp_fit / band_spacing_nm:.2f} bands peak-to-valley")

try:

import matplotlib

matplotlib.use("Agg")

import matplotlib.pyplot as plt

plt.figure(figsize=(11, 5))

plt.scatter(cols, vals, s=6, alpha=0.4, label="measured (sub-pixel)", color="steelblue")

plt.plot(cols, fit_curve, color="darkorange", linewidth=2,

label=f"quadratic fit (p2v={ptp_fit:.2f} nm)")

plt.axhline(args.feature_nm, color="gray", linestyle="--", alpha=0.5,

label=f"nominal {args.feature_nm:.0f} nm")

plt.xlabel("Sample (across-track column)")

plt.ylabel(f"Detected absorption minimum (nm)")

plt.title("Spectral smile estimate")

plt.legend()

plt.tight_layout()

plt.savefig(args.out, dpi=120)

print(f"\nSaved plot to {args.out}")

except ImportError:

print("\n(matplotlib not installed - numeric results above still usable)")

if __name__ == "__main__":

main()