Per ovviare al problema della calibrazione visto che nel precedente post (calibrazione da manuale del sensore) ho estratto la firma in radianza di 4 punti del pannello di bianco di riferimento direttamente dall'immagine (ho avuto l'accortezza di inserirlo nella scena) usando Snap ed estrando 4 spettri in csv
qui i risultati sono decimente migliori...e' corretto che il pannello di bianco di riferimento sia a 0.85 perche' e' calibrato all'80%
 |
| Calibrazione di fabbrica del riferimento di bianco |
Confrontando con gli assorbimenti atmosferici
O₂ A-band~758-770 (minimo a ~760-762nm) Ossigeno Forte
H₂O~810-840 Vapore acqueo Debole-moderata
H₂O (complesso 940nm)~890-990 (minimo tipico ~935-945, ma variabile)
Confrontando con gli spettri misurati
python3 make_reflectance_cube.py /media/disco/Arpat/FigSpec_260522_171651/figspec.hdr /media/disco/Arpat/FigSpec_260522_171651/figspec.spe --black /media/disco/Arpat/FigSpec_260522_171651/figspec.figspecblack --white-csv bianco.csv --panel ./Reflectance_Calibration_Cloth.figspecref --out figspec_reflectance --from-line 33000 --to-line 40000
#!/usr/bin/env python3
"""
Computes a full calibrated reflectance cube (all bands) from an ENVI BIL
hyperspectral cube using dark-current and white-reference correction files
plus the actual spectral reflectance curve of the calibration panel, and
saves the result as a standard ENVI file (float32, BIL).
Usage:
python3 make_reflectance_cube.py figspec.hdr figspec.spe \
--black figspec.figspecblack \
--white figspec.figspecwhite \
--panel Reflectance_Calibration_Cloth.figspecref \
--out reflectance_cube
Produces:
reflectance_cube.hdr
reflectance_cube.img (float32, BIL, same samples/lines/bands as input)
Reflectance is computed per band, per spatial column as:
reflectance = (raw - dark) / (white - dark + eps) * panel_reflectance(wavelength)
Processes the cube in chunks of lines to keep memory usage low regardless
of file size.
"""
import os
import re
import argparse
import numpy as np
from scipy.signal import savgol_filter
def auto_offset_and_lines(path, bands, samples, itemsize, forced_lines=None,
forced_offset=None):
per_line_bytes = bands * samples * itemsize
file_size = os.path.getsize(path)
if forced_offset is not None and forced_lines is not None:
return forced_offset, forced_lines
expected_lines = forced_lines if forced_lines is not None else file_size // per_line_bytes
expected = per_line_bytes * expected_lines
if file_size == expected:
return 0, expected_lines
remainder = file_size % per_line_bytes
lines = file_size // per_line_bytes
print(f" [{os.path.basename(path)}] size={file_size} -> "
f"auto offset={remainder} bytes, lines={lines}")
return remainder, lines
def parse_correction_file(path, n_bands):
"""Parse a .figspecblack/.figspecwhite correction file -> (bands, spatial_bins)."""
data = open(path, "rb").read()
text = data.decode("latin1")
matches = list(re.finditer(r"[0-9a-fA-F]{200,}", text))
if not matches:
raise ValueError(f"No hex-encoded correction block found in {path}")
m = matches[0]
hexstr = data[m.start():m.end()].decode("ascii")
raw = bytes.fromhex(hexstr)
arr = np.frombuffer(raw, dtype="<f4")
if arr.size % n_bands != 0:
raise ValueError(f"{path}: {arr.size} floats doesn't divide evenly by {n_bands} bands")
spatial_bins = arr.size // n_bands
return arr.reshape(n_bands, spatial_bins)
def parse_white_csv(path, n_bands, wavelengths, wl_tolerance=0.1):
"""Parse a CSV of raw radiance spectra extracted directly from the image at the
calibration panel's location (e.g. several 'Pin N' columns from ENVI Spectrum View).
Averages all numeric columns together and returns a (bands, 1) map (spatially
uniform - no across-track/vignetting correction, since it's just point samples).
"""
with open(path, "r", errors="ignore") as f:
header = f.readline().strip().split("\t")
if len(header) == 1: # try comma if tab didn't split
header = header[0].split(",")
sep = ","
else:
sep = "\t"
rows = []
for line in f:
line = line.strip()
if not line:
continue
parts = line.split(sep)
rows.append([float(p) for p in parts])
rows = np.array(rows) # (bands, 1+n_pins) : col0 = wavelength
csv_wl = rows[:, 0]
csv_vals = rows[:, 1:] # (bands, n_pins)
if len(csv_wl) != n_bands:
raise ValueError(f"{path}: has {len(csv_wl)} rows but cube has {n_bands} bands")
max_diff = np.max(np.abs(csv_wl - wavelengths))
if max_diff > wl_tolerance:
print(f" WARNING: wavelength mismatch between {path} and cube header "
f"(max diff {max_diff:.3f} nm) - check band alignment!")
avg_spectrum = csv_vals.mean(axis=1) # (bands,)
n_pins = csv_vals.shape[1]
print(f" Parsed {path}: {n_pins} replicate spectra, averaged -> single (bands,) profile")
return avg_spectrum.reshape(n_bands, 1).astype(np.float32) # (bands, spatial_bins=1)
def parse_panel_reflectance(path):
"""Parse the calibration cloth CSV -> sorted (wavelengths_nm, reflectance_fraction) arrays."""
data = open(path, "rb").read().decode("latin1")
rows = []
for line in data.splitlines()[1:]:
line = line.strip().strip(",")
if not line:
continue
parts = [p.strip() for p in line.split(",") if p.strip() != ""]
if len(parts) < 2:
continue
try:
wl = float(parts[0])
refl_pct = float(parts[1])
rows.append((wl, refl_pct))
except ValueError:
continue
rows.sort()
wl = np.array([r[0] for r in rows])
refl = np.array([r[1] for r in rows]) / 100.0
return wl, refl
def upsample_spatial_2d(profile_2d, target_len):
"""Upsample a (bands, spatial_bins) map to (bands, target_len) via nearest-neighbor repeat."""
bins = profile_2d.shape[1]
factor = target_len / bins
if float(factor).is_integer():
return np.repeat(profile_2d, int(factor), axis=1)
x_old = np.linspace(0, 1, bins)
x_new = np.linspace(0, 1, target_len)
out = np.empty((profile_2d.shape[0], target_len), dtype=np.float32)
for b in range(profile_2d.shape[0]):
out[b] = np.interp(x_new, x_old, profile_2d[b])
return out
def write_envi_hdr(out_hdr_path, samples, lines, bands, wavelengths, wl_units,
interleave="bil", description="Calibrated reflectance cube"):
with open(out_hdr_path, "w") as f:
f.write("ENVI\n")
f.write(f"description = {{{description}}}\n")
f.write(f"samples = {samples}\n")
f.write(f"lines = {lines}\n")
f.write(f"bands = {bands}\n")
f.write("header offset = 0\n")
f.write("file type = ENVI Standard\n")
f.write("data type = 4\n") # float32
f.write(f"interleave = {interleave}\n")
f.write("byte order = 0\n")
f.write(f"wavelength units = {wl_units}\n")
wl_str = ", ".join(f"{w:.4f}" for w in wavelengths)
f.write(f"wavelength = {{{wl_str}}}\n")
f.write("reflectance scale factor = 1.0\n")
if __name__ == "__main__":
main()
Nessun commento:
Posta un commento