Truecolor FigSpec FS60-CL

Un piccolo script per avere un truecolor dal sensore iperspettrale FigSpex F60-CL

L'unica elaborazione e' uno stretch delle bande dei colori RGB 

 

Red band:   idx 123 -> 649.58 nm
Green band: idx 75 -> 549.73 nm
Blue band:  idx 28 -> 450.00 nm

#!/usr/bin/env python3

"""

Extract a true-color (RGB) composite from a large ENVI BIL hyperspectral cube

without loading the whole file into memory.

Usage:

python3 make_true_color.py figspec.hdr figspec.img true_color.png

If your data file has a different extension (.dat, .raw, .bin, or no

extension at all), just pass its actual path as the second argument.

"""

import sys

import os

import re

import argparse

import numpy as np

from PIL import Image

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

# ENVI data type codes -> numpy dtype

dtype_map = {

1: np.uint8,

2: np.int16,

3: np.int32,

4: np.float32,

5: np.float64,

6: np.complex64,

9: np.complex128,

12: np.uint16,

13: np.uint32,

14: np.int64,

15: np.uint64,

}

np_dtype = np.dtype(dtype_map[dtype_code])

if byte_order == 1:

np_dtype = np_dtype.newbyteorder(">")

else:

np_dtype = np_dtype.newbyteorder("<")

return {

"samples": samples,

"lines": lines,

"bands": bands,

"dtype": np_dtype,

"interleave": interleave,

"wavelengths": wavelengths,

}

def closest_band(wavelengths, target):

return min(range(len(wavelengths)), key=lambda i: abs(wavelengths[i] - target))

def stretch(band, low_pct=2, high_pct=98):

"""Percentile contrast stretch to 0-255, robust to outliers/no-data."""

band = band.astype(np.float32)

valid = band[np.isfinite(band)]

lo, hi = np.percentile(valid, [low_pct, high_pct])

if hi <= lo:

hi = lo + 1

out = np.clip((band - lo) / (hi - lo), 0, 1)

return (out * 255).astype(np.uint8)

def main():

ap = argparse.ArgumentParser()

ap.add_argument("hdr_path")

ap.add_argument("data_path")

ap.add_argument("out_path", nargs="?", default="true_color.png")

ap.add_argument("--offset", type=int, default=None,

help="Byte offset into data file where the cube starts "

"(overrides auto-detection / header offset)")

ap.add_argument("--lines", type=int, default=None,

help="Override number of lines (in case header value is wrong)")

args = ap.parse_args()

hdr_path = args.hdr_path

data_path = args.data_path

out_path = args.out_path

meta = parse_hdr(hdr_path)

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

dtype = meta["dtype"]

interleave = meta["interleave"]

wavelengths = meta["wavelengths"]

if args.lines is not None:

lines = args.lines

itemsize = dtype.itemsize

per_line_bytes = bands * samples * itemsize

file_size = os.path.getsize(data_path)

if args.offset is not None:

offset = args.offset

else:

expected = per_line_bytes * lines

if file_size == expected:

offset = 0

elif file_size > expected:

# extra trailing bytes; ignore them, assume no leading offset

offset = 0

else:

# file smaller than expected: leftover bytes at the START are

# likely a header the .hdr file didn't account for.

remainder = file_size % per_line_bytes

offset = remainder

lines = file_size // per_line_bytes

print(f"NOTE: file size ({file_size}) doesn't match header dims. "

f"Auto-detected a {offset}-byte leading offset and adjusted "

f"lines to {lines}. Override with --offset/--lines if wrong.")

print(f"Cube: {lines} lines x {bands} bands x {samples} samples, "

f"dtype={dtype}, interleave={interleave}, offset={offset} bytes")

if wavelengths:

r_idx = closest_band(wavelengths, 650)

g_idx = closest_band(wavelengths, 550)

b_idx = closest_band(wavelengths, 450)

print(f"Red band: idx {r_idx} -> {wavelengths[r_idx]:.2f} nm")

print(f"Green band: idx {g_idx} -> {wavelengths[g_idx]:.2f} nm")

print(f"Blue band: idx {b_idx} -> {wavelengths[b_idx]:.2f} nm")

else:

raise ValueError("No wavelength info in header; specify band indices manually.")

# Memory-map the file so we only touch the pages we actually read.

if interleave == "bil":

# file layout: [line][band][sample]

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

shape=(lines, bands, samples))

red = np.array(cube[:, r_idx, :])

green = np.array(cube[:, g_idx, :])

blue = np.array(cube[:, b_idx, :])

elif interleave == "bsq":

# file layout: [band][line][sample]

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

shape=(bands, lines, samples))

red = np.array(cube[r_idx, :, :])

green = np.array(cube[g_idx, :, :])

blue = np.array(cube[b_idx, :, :])

elif interleave == "bip":

# file layout: [line][sample][band]

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

shape=(lines, samples, bands))

red = np.array(cube[:, :, r_idx])

green = np.array(cube[:, :, g_idx])

blue = np.array(cube[:, :, b_idx])

else:

raise ValueError(f"Unsupported interleave: {interleave}")

print("Bands extracted, stretching contrast...")

r8 = stretch(red)

g8 = stretch(green)

b8 = stretch(blue)

rgb = np.dstack([r8, g8, b8])

img = Image.fromarray(rgb, mode="RGB")

img.save(out_path)

print(f"Saved true-color image to {out_path} ({img.size[0]}x{img.size[1]})")

if __name__ == "__main__":

main()

 

 

 

 

 

Empirical line on image reference FigSpec

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

 

Script radianza riflettanza per FigSpec F60-CL

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