Script radianza riflettanza per FigSpec F60-C

Seguendo i due post precedenti (1 e 2) di seguito e' riportato lo script per convertire l'immagine da DN a radianza e riflettanza usando empirical line e il pannello di bianco di riferimento (85%) in dotazione con lo strumento 

Per lanciare lo script si puo' usare 

python3 make_reflectance_cube.py figspec.hdr figspec.spe \
    --black figspec.figspecblack \
    --white figspec.figspecwhite \
    --panel Reflectance_Calibration_Cloth.figspecref \
    --out figspec_reflectance

 visto che le strisciate possono essere anche molto lunghe si puo' selezionare un intervallo di righe

 

python3 make_reflectance_cube.py figspec.hdr figspec.spe \
    --black figspec.figspecblack \
    --white figspec.figspecwhite \
    --panel Reflectance_Calibration_Cloth.figspecref \
    --out figspec_reflectance_5k_10k \
    --from-line 5000 --to-line 10000 

E qui iniziano i problemi

Se si mettono a confronto la riflettanza (prima immagine) e la radianza (seconda immagine) si vede che la riflettanza e' decisamente piu' rumorosa  ed ha valori superiori ad 1 (dovrebbero essere sempre inferiori ad 1)

 

 se poi vado ad estrarre lo spettro sul pannello di riferimento lo spettro non e' piatto (il dato e' clippato ma e' un artefatto del programma, se si vede l'immagine successiva si deve che i DN sono intorno a 3000-3100 quando il limite superiore e' 4095)

Temo che questo derivi dal metodo con cui lo strumento viene calibrato...a terra viene effettuato una immagine di dark current con l'otturatore chiuso ed una immagine di bianco tenendo in mano il drone ed inquadrando il bianco di riferimento..in pratica il bianco non e' ricavato direttamente dalla scena e possono passare anche diversi minuti tra il bianco a terra e la ripresa di tutta la scena...il che vuol dire che le condizioni di illuminazione possono essere cambiate
 

#!/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_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()