TLDR: volare sempre orientati con l'azimuth del Sole
l'azimuth del Sole cambia con il giorno e l'ora.All'alba ed al tramonto e' molto vicino a Est ed Ovest ma a mezzogiorno e' orientato verso Sud per la latitudini settentrionali
Frugando sulla chiavetta USB del software del sensore ho trovato un file denominato FS6261009镜头系数 che tradotto suona come "coefficiente dell'obbiettivo"....aprendolo si vede che si tratta di un file json in cui sono inseriti due array (uno monodimensionale di 480 valori ed uno bidimensionale 72x480)
il grafico assomiglia ad una correzione di vignettatura ma e' solo per una banda (sembra la 121)
Ho provato ad usare questa informazione per correggere la vignettatura ma nonostante la correzione l'effetto al bordo non e' stato eliminato ..l'ipotesi successiva e' che non sia qualcosa di relativo al sensore ma alla geometria di ripresa
Dai dati GPS si ha
Punti usati: 166
Lunghezza traccia: 134.1 m
Azimuth linea di volo: 251.86 deg (reciproco: 71.86 deg)
Direzione across-track (perpendicolare al volo): 341.86 deg / 161.86 deg
Al giorno e l'ora di ripresa il Sole aveva una elevazione di 66 gradi ed un azimuth, quindi la traiettoria di volo era praticamente parallela al direzione del Sole
Il sensore della camera ha una FOV di 24.75 gradi (IFOV 0.9 mrad × 480 canali spaziali)
Il fenomeno ai bordi quindi puo' essere non solo la vignettatura ma anche il forward scattering ed il backscattering che generano il fenomeno di hotspot BRDF
In pratica la direzione di ripresa e' tale che in una direzione l'ombra dell'oggetto e' coperta dall'oggetto stesso (back scattering, massima riflettanza hotspot propriamente detto, guardando dalla parte opposta dal Sole) mentre nella direzione opposta si osserva l'ombra piuttosto che l'oggetto stesso (forward scattering, guardando verso il Sole)
Il fenomeno affligge tutti i sensori iperspettrali sia da satellite che airborne ed e' particolarmente importante quanto piu' grande e' l'angolo di field of view. Un satellite che ha tipicamente questo problema e' Modis con una FOV di 55 gradi con una correzione basata sul modello a kernel RossThick-LiSparse. Lo stesso fenomeno si puo' verificare in analisi di immagini multitemporali anche di satelliti a bassa FOV (tipo Landsat 15 gradi) pur essendo eliosincroni
Altrimenti per una correzione piu' semplice viene utilizzato il modello RPV — Rahman-Pinty-Verstraete
#!/usr/bin/env python3
"""
brdf_hotspot_validation.py
Conferma quantitativa (non solo qualitativa) dell'ipotesi che l'asimmetria
sinistra/destra osservata nel profilo across-track del true-color sia
spiegabile con l'effetto hotspot/anti-hotspot BRDF, usando il modello
semi-empirico RPV (Rahman-Pinty-Verstraete) — lo stesso tipo di modello
usato in letteratura (es. MISR, POLDER) per descrivere la riflettanza
bidirezionale di superfici naturali includendo l'hotspot.
Idea generale
-------------
1. Dalla traccia GPS (o da un azimuth di volo fornito manualmente) si
ricava la direzione across-track (perpendicolare al volo).
2. Nota la posizione del sole (azimuth, elevazione) e il FOV del sensore,
si calcola per ogni colonna dello swath la geometria di vista (zenith e
azimuth relativo al sole) e si predice con RPV il fattore BRDF atteso
(normalizzato, cioe' solo la FORMA della variazione, non il valore
assoluto di riflettanza).
3. Dal true-color si estrae il profilo across-track empirico (luma media
per colonna), e se disponibile un file di calibrazione vignettatura
(vendor, vedi analyze_crosstrack_illumination.py --vignette-json) lo si
usa per "scorporare" la componente ottica (vignettatura) da quella di
scena (BRDF), isolando il residuo asimmetrico da confrontare col modello.
4. Si calcola la correlazione tra la forma prevista da RPV e il residuo
osservato: un'alta correlazione conferma quantitativamente che
l'asimmetria e' spiegabile con l'hotspot/anti-hotspot, non con un
artefatto strumentale.
Il modello RPV (riassunto)
---------------------------
R(theta_i, theta_r, phi) = rho0 * M(theta_i,theta_r) * F(g,Theta) * H(G)
M(theta_i,theta_r) = (cos(theta_i) cos(theta_r) (cos(theta_i)+cos(theta_r)))^(k-1)
-> termine di forma "bowl/dome"; k=1 lo rende neutro (default)
F(g,Theta) = (1-Theta^2) / (1 + 2*Theta*cos(g) + Theta^2)^1.5
-> funzione di fase di Henyey-Greenstein; g = angolo di scattering;
Theta<0 favorisce il backscattering (tipico canopy vegetale)
H(G) = 1 + (1-h) / (1 + G/h)
-> funzione hotspot; G = sqrt(tan^2(theta_i)+tan^2(theta_r)
-2*tan(theta_i)*tan(theta_r)*cos(phi)) e' la "distanza" in spazio
tangente tra le direzioni sole e sensore (G=0 esattamente
all'hotspot); h controlla quanto e' "stretto" il picco (tipico
0.05-0.3 per canopy vegetali)
rho0 viene posto a 1 perche' qui interessa solo la FORMA relativa, non il
valore assoluto di riflettanza.
Uso
---
python3 brdf_hotspot_validation.py true_color.png \\
--gps-csv figspec.gps \\
--sun-azimuth 154 --sun-elevation 66 \\
--vignette-json test.json \\
--out out_brdf
# oppure con azimuth di volo manuale invece del GPS:
python3 brdf_hotspot_validation.py true_color.png \\
--flight-azimuth 251.86 --sun-azimuth 154 --sun-elevation 66 --out out_brdf
"""
import argparse
import json
import sys
from pathlib import Path
import numpy as np
from PIL import Image
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
# ---------------------------------------------------------------------
# Geometria di volo (azimuth) dai punti GPS
# ---------------------------------------------------------------------
def flight_azimuth_from_gps(csv_path, lat_col="Latitude", lon_col="Longitude", delimiter=","):
import pandas as pd
df = pd.read_csv(csv_path, delimiter=delimiter)
lat = df[lat_col].values.astype(float)
lon = df[lon_col].values.astype(float)
lat0 = np.radians(lat.mean())
R = 6371000.0
x = np.radians(lon - lon.mean()) * R * np.cos(lat0)
y = np.radians(lat - lat.mean()) * R
pts = np.column_stack([x, y])
pts_centered = pts - pts.mean(axis=0)
cov = np.cov(pts_centered.T)
eigvals, eigvecs = np.linalg.eigh(cov)
main_dir = eigvecs[:, np.argmax(eigvals)]
az = np.degrees(np.arctan2(main_dir[0], main_dir[1])) % 360
return az
# ---------------------------------------------------------------------
# Modello RPV
# ---------------------------------------------------------------------
def rpv_reflectance(theta_i_deg, theta_r_deg, phi_deg, k=1.0, theta_hg=-0.15, h=0.1, rho0=1.0):
"""
Calcola il fattore di riflettanza bidirezionale RPV (non normalizzato
in valore assoluto, rho0 e' un fattore di scala arbitrario).
theta_i_deg : zenith solare (scalare, gradi)
theta_r_deg : zenith di vista (array, gradi)
phi_deg : azimuth relativo vista-sole (array, gradi) = az_vista - az_sole
k : parametro di forma bowl/dome (1.0 = neutro)
theta_hg : parametro di asimmetria della fase di Henyey-Greenstein (-1..1)
h : parametro di nitidezza dell'hotspot (piu' piccolo = picco piu' stretto)
rho0 : fattore di scala (irrilevante se si normalizza il risultato)
"""
ti = np.radians(theta_i_deg)
tr = np.radians(theta_r_deg)
phi = np.radians(phi_deg)
cos_ti, cos_tr = np.cos(ti), np.cos(tr)
sin_ti, sin_tr = np.sin(ti), np.sin(tr)
# angolo di scattering (fase)
cos_g = cos_ti * cos_tr + sin_ti * sin_tr * np.cos(phi)
cos_g = np.clip(cos_g, -1.0, 1.0)
# funzione di fase Henyey-Greenstein
F = (1 - theta_hg**2) / (1 + 2 * theta_hg * cos_g + theta_hg**2) ** 1.5
# funzione hotspot (spazio tangente)
tan_ti, tan_tr = np.tan(ti), np.tan(tr)
G = np.sqrt(np.clip(tan_ti**2 + tan_tr**2 - 2 * tan_ti * tan_tr * np.cos(phi), 0, None))
H = 1 + (1 - h) / (1 + G / h)
# termine di forma bowl/dome
M = (cos_ti * cos_tr * (cos_ti + cos_tr)) ** (k - 1)
R = rho0 * M * F * H
return R, np.degrees(np.arccos(cos_g)), G
# ---------------------------------------------------------------------
# Profilo across-track empirico dal true-color
# ---------------------------------------------------------------------
def empirical_column_profile(image_path, nodata_thresh=3.0):
img = Image.open(image_path).convert("RGB")
arr = np.asarray(img).astype(np.float64)
luma = arr.mean(axis=2)
mask = luma > nodata_thresh
cols = arr.shape[1]
profile = np.full(cols, np.nan)
for c in range(cols):
if mask[:, c].any():
profile[c] = luma[mask[:, c], c].mean()
return profile
def load_vignette_gain(vignette_json, band_key, cols):
with open(vignette_json, "r", encoding="utf-8") as f:
data = json.load(f)
coeffs = {e["Key"]: np.array(e["Value"], dtype=float) for e in data.get("Coefficients", [])}
if not coeffs:
raise ValueError(f"Nessun coefficiente trovato in {vignette_json}")
key = band_key if band_key is not None and band_key in coeffs else list(coeffs.keys())[0]
gain = coeffs[key]
if len(gain) != cols:
x_orig = np.linspace(0, 1, len(gain))
x_new = np.linspace(0, 1, cols)
gain = np.interp(x_new, x_orig, gain)
return gain, key
def fit_symmetric_vignetting(profile):
"""
Fallback se non e' disponibile un file di calibrazione vendor: stima
la componente di vignettatura come un fit quadratico del profilo (bella
approssimazione per una caduta di illuminazione a "campana" simmetrica),
e ritorna il "gain" equivalente (1/vignetting_shape_normalizzata) cosi'
l'interfaccia e' la stessa di load_vignette_gain.
"""
cols = len(profile)
x = np.arange(cols)
valid = np.isfinite(profile)
coeffs = np.polyfit(x[valid], profile[valid], deg=2)
fitted = np.polyval(coeffs, x)
fitted_norm = fitted / np.nanmean(fitted)
gain = 1.0 / fitted_norm
return gain, None
# ---------------------------------------------------------------------
# Analisi principale
# ---------------------------------------------------------------------
def analyze(image_path, out_dir, sun_azimuth, sun_elevation, fov_deg,
flight_azimuth=None, gps_csv=None, vignette_json=None, band_key=None,
rpv_k=1.0, rpv_theta=-0.15, rpv_h=0.1, nodata_thresh=3.0):
out_dir = Path(out_dir)
out_dir.mkdir(parents=True, exist_ok=True)
# --- 1. azimuth di volo ---
if flight_azimuth is None:
if gps_csv is None:
raise ValueError("Serve --flight-azimuth oppure --gps-csv per determinare la geometria di volo")
flight_azimuth = flight_azimuth_from_gps(gps_csv)
print(f"Azimuth di volo stimato dai punti GPS: {flight_azimuth:.2f} deg")
else:
print(f"Azimuth di volo (fornito manualmente): {flight_azimuth:.2f} deg")
across_1 = (flight_azimuth + 90) % 360
across_2 = (flight_azimuth - 90) % 360
print(f"Direzioni across-track: {across_1:.2f} deg / {across_2:.2f} deg")
theta_i = 90 - sun_elevation
print(f"Sole: azimuth={sun_azimuth:.1f} deg, elevazione={sun_elevation:.1f} deg (zenith={theta_i:.1f} deg)")
def fold180(a):
a = a % 360
return min(a, 360 - a)
d1 = fold180(across_1 - sun_azimuth)
d2 = fold180(across_2 - sun_azimuth)
print(f"Scarto tra direzioni across-track e azimuth solare: {d1:.1f} deg / {d2:.1f} deg")
# --- 2. profilo empirico dal true-color ---
profile = empirical_column_profile(image_path, nodata_thresh)
cols = len(profile)
x = np.arange(cols)
center = (cols - 1) / 2.0
# --- 3. geometria di vista per colonna ---
half_fov = fov_deg / 2.0
view_zenith = np.abs((x - center) / center) * half_fov # 0 al centro, half_fov ai bordi
# due possibili assegnazioni lato<->azimuth (non sappiamo a priori quale
# lato del sensore corrisponde a quale bearing fisico, dipende dal
# montaggio/scansione del sensore: testiamo entrambe e scegliamo quella
# che correla meglio con l'osservato, riportando comunque i risultati
# di entrambe per trasparenza)
offset = x - center
view_az_A = np.where(offset < 0, across_1, across_2)
view_az_B = np.where(offset < 0, across_2, across_1)
results = {}
for label, view_az in [("A", view_az_A), ("B", view_az_B)]:
phi = view_az - sun_azimuth
R, phase_angle, G = rpv_reflectance(theta_i, view_zenith, phi, k=rpv_k, theta_hg=rpv_theta, h=rpv_h)
results[label] = {
"view_az": view_az,
"R": R,
"phase_angle": phase_angle,
"G": G,
}
# --- 4. vignettatura (per isolare il residuo BRDF dall'empirico) ---
if vignette_json is not None:
gain, used_key = load_vignette_gain(vignette_json, band_key, cols)
vign_source = f"file vendor ({vignette_json}, Key={used_key})"
else:
gain, used_key = fit_symmetric_vignetting(profile)
vign_source = "fit quadratico simmetrico (fallback, nessun file vendor fornito)"
vignetting_shape = 1.0 / gain
vignetting_shape_norm = vignetting_shape / np.nanmean(vignetting_shape)
print(f"Componente di vignettatura stimata da: {vign_source}")
empirical_norm = profile / np.nanmean(profile)
residual_empirical = empirical_norm / vignetting_shape_norm
residual_empirical_norm = residual_empirical / np.nanmean(residual_empirical)
# --- 5. confronto modello vs osservato, per entrambe le assegnazioni ---
def compare(model_R):
model_norm = model_R / np.nanmean(model_R)
valid = np.isfinite(residual_empirical_norm) & np.isfinite(model_norm)
r = np.corrcoef(residual_empirical_norm[valid], model_norm[valid])[0, 1]
rmse = np.sqrt(np.nanmean((residual_empirical_norm[valid] - model_norm[valid]) ** 2))
# combinato: vignettatura * BRDF, confrontato con il profilo grezzo (non scorporato)
combined = model_norm * vignetting_shape_norm
combined_norm = combined / np.nanmean(combined)
r_combined = np.corrcoef(empirical_norm[valid], combined_norm[valid])[0, 1]
return model_norm, r, rmse, combined_norm, r_combined
lines = []
lines.append(f"Immagine: {image_path}")
lines.append(f"Colonne (samples): {cols}")
lines.append(f"FOV assunto: {fov_deg:.2f} deg (default: da datasheet FS-60, IFOV 0.9 mrad x 480 samples)")
lines.append(f"Azimuth di volo: {flight_azimuth:.2f} deg")
lines.append(f"Direzioni across-track: {across_1:.2f} / {across_2:.2f} deg")
lines.append(f"Sole: azimuth={sun_azimuth:.1f} deg, elevazione={sun_elevation:.1f} deg (zenith={theta_i:.1f} deg)")
lines.append(f"Scarto tra direzioni across-track e azimuth solare: {d1:.1f} deg / {d2:.1f} deg")
lines.append(f"Parametri RPV: k={rpv_k}, Theta={rpv_theta}, h={rpv_h}")
lines.append(f"Vignettatura stimata da: {vign_source}")
lines.append("")
best_label, best_r = None, -2
for label in ("A", "B"):
view_az = results[label]["view_az"]
model_norm, r, rmse, combined_norm, r_combined = compare(results[label]["R"])
results[label]["model_norm"] = model_norm
results[label]["r_residual"] = r
results[label]["rmse_residual"] = rmse
results[label]["combined_norm"] = combined_norm
results[label]["r_combined"] = r_combined
side_desc = f"colonne basse -> {view_az[0]:.1f} deg, colonne alte -> {view_az[-1]:.1f} deg"
lines.append(f"--- Assegnazione {label} ({side_desc}) ---")
lines.append(f" Correlazione modello BRDF vs residuo empirico (vignettatura scorporata): r={r:.3f} (RMSE={rmse:.4f})")
lines.append(f" Correlazione modello combinato (vignettatura x BRDF) vs profilo grezzo: r={r_combined:.3f}")
lines.append(f" Range angolo di fase coperto dal FOV: {results[label]['phase_angle'].min():.1f} - {results[label]['phase_angle'].max():.1f} deg")
lines.append(f" Range fattore hotspot G (0=hotspot esatto): {results[label]['G'].min():.3f} - {results[label]['G'].max():.3f}")
lines.append("")
if r > best_r:
best_r = r
best_label = label
lines.append(f"Assegnazione con correlazione migliore: {best_label} (r={best_r:.3f})")
hotspot_zenith_in_fov = theta_i <= half_fov
lines.append(f"Zenith solare ({theta_i:.1f} deg) {'RIENTRA' if hotspot_zenith_in_fov else 'NON rientra'} "
f"nel range di zenith di vista coperto dal FOV (0-{half_fov:.1f} deg): "
f"{'ci si aspetta un picco hotspot pieno vicino al bordo' if hotspot_zenith_in_fov else 'ci si aspetta solo un aumento progressivo verso il bordo, senza raggiungere il picco pieno'}")
report_text = "\n".join(lines)
print("\n" + report_text)
(out_dir / "report.txt").write_text(report_text, encoding="utf-8")
# --- Figura 1: confronto residuo empirico vs modello BRDF (entrambe le assegnazioni) ---
fig, ax = plt.subplots(figsize=(10, 5.5))
ax.plot(x, residual_empirical_norm, color="black", linewidth=1.8, label="residuo empirico (vignettatura scorporata)")
ax.plot(x, results["A"]["model_norm"], color="crimson", linestyle="--",
label=f"modello RPV - assegnazione A (r={results['A']['r_residual']:.3f})")
ax.plot(x, results["B"]["model_norm"], color="steelblue", linestyle="--",
label=f"modello RPV - assegnazione B (r={results['B']['r_residual']:.3f})")
ax.set_xlabel("Colonna (direzione across-track)")
ax.set_ylabel("Valore normalizzato (media=1)")
ax.set_title("Residuo empirico (dopo rimozione vignettatura) vs modello BRDF (RPV)")
ax.legend(fontsize=8)
ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig(out_dir / "residual_vs_brdf_model.png", dpi=150)
plt.close(fig)
# --- Figura 2: profilo grezzo vs modello combinato (vignettatura x BRDF) ---
fig, ax = plt.subplots(figsize=(10, 5.5))
ax.plot(x, empirical_norm, color="black", linewidth=1.8, label="profilo empirico grezzo (normalizzato)")
ax.plot(x, vignetting_shape_norm, color="gray", linestyle=":", label="sola componente vignettatura")
ax.plot(x, results[best_label]["combined_norm"], color="darkorange", linewidth=2,
label=f"modello combinato vignettatura x BRDF (assegn. {best_label}, r={results[best_label]['r_combined']:.3f})")
ax.set_xlabel("Colonna (direzione across-track)")
ax.set_ylabel("Valore normalizzato (media=1)")
ax.set_title("Profilo grezzo vs modello combinato (vignettatura ottica x effetto BRDF hotspot)")
ax.legend(fontsize=8)
ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig(out_dir / "raw_vs_combined_model.png", dpi=150)
plt.close(fig)
print(f"\nOutput salvati in: {out_dir.resolve()}")
print(" - report.txt")
print(" - residual_vs_brdf_model.png")
print(" - raw_vs_combined_model.png")
def main():
parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument("image", help="Immagine true-color (PNG/JPEG/TIFF)")
parser.add_argument("--out", default="out_brdf", help="Cartella di output")
parser.add_argument("--gps-csv", default=None, help="File CSV GPS per stimare l'azimuth di volo automaticamente")
parser.add_argument("--flight-azimuth", type=float, default=None, help="Azimuth di volo manuale (gradi), alternativa a --gps-csv")
parser.add_argument("--sun-azimuth", type=float, required=True, help="Azimuth solare (gradi da nord)")
parser.add_argument("--sun-elevation", type=float, required=True, help="Elevazione solare sopra l'orizzonte (gradi)")
parser.add_argument("--fov", type=float, default=24.75,
help="FOV across-track totale del sensore (gradi). Default 24.75 deg, "
"stimato da datasheet FS-60 (IFOV 0.9 mrad x 480 samples)")
parser.add_argument("--vignette-json", default=None, help="File JSON di calibrazione vignettatura vendor (opzionale)")
parser.add_argument("--band-key", type=int, default=None, help="Key della banda da usare dal file vignette-json")
parser.add_argument("--rpv-k", type=float, default=1.0, help="Parametro RPV k (forma bowl/dome, default 1.0=neutro)")
parser.add_argument("--rpv-theta", type=float, default=-0.15, help="Parametro RPV Theta (asimmetria fase, default -0.15)")
parser.add_argument("--rpv-h", type=float, default=0.1, help="Parametro RPV h (nitidezza hotspot, default 0.1)")
parser.add_argument("--nodata-thresh", type=float, default=3.0, help="Soglia luma per nodata (default 3.0)")
args = parser.parse_args()
analyze(args.image, args.out, args.sun_azimuth, args.sun_elevation, args.fov,
flight_azimuth=args.flight_azimuth, gps_csv=args.gps_csv,
vignette_json=args.vignette_json, band_key=args.band_key,
rpv_k=args.rpv_k, rpv_theta=args.rpv_theta, rpv_h=args.rpv_h,
nodata_thresh=args.nodata_thresh)
if __name__ == "__main__":
main()
Nessun commento:
Posta un commento