Calcolo azimuth ottimale di volo

 Script per il calcolo della direzione di volo per minimizzare l'asimmetria derivante da BRDF

 

pip install astral tzdata  

python3 optimal_flight_azimuth.py --lat 43.7696 --lon 11.2558 \
    --date 2026-07-22 --time 11:46 --tz Europe/Rome 

 

Data/ora locale: 2026-07-22 11:46:00+02:00
Posizione: lat=43.7696, lon=11.2558
Sole: azimuth=132.17 deg, elevazione=59.20 deg (zenith=30.80 deg)
FOV assunto: 24.75 deg  |  parametri RPV: k=1.0, Theta=-0.15, h=0.1

--- Risultato ---
Azimuth di volo OTTIMALE (minima asimmetria BRDF): 132.0 deg (reciproco: 312.0 deg)  ->  escursione prevista: 1.8%
Azimuth di volo PEGGIORE (massima asimmetria BRDF): 42.0 deg (reciproco: 222.0 deg)  ->  escursione prevista: 20.8%
 

il parametro teta indica l'asimmetria della funzione di fase ovvero quanto e' preferenziale il backscattering rispetto al forward scattering da parte della superficie (teta = 0 isotropia, valori minori di zero privilegia backascattering)

 

 il parametro H indica quanto e' largo il picco di hotspot o meglio quanto gradualmente il picco si presenta sul bordo

 

 

 

///////////////////////////////////////////////////////////////////////////////////////////// 

#!/usr/bin/env python3
"""
optimal_flight_azimuth.py

Calcola l'azimuth di volo che MINIMIZZA l'asimmetria across-track dovuta
all'effetto BRDF hotspot/anti-hotspot, data una data, un'ora e una
posizione geografica (lat/lon).

Principio
---------
L'asimmetria hotspot/anti-hotspot e' massima quando l'asse across-track
(perpendicolare alla linea di volo) giace nel piano principale solare,
cioe' quando la linea di volo e' perpendicolare all'azimuth del sole.
E' MINIMA quando la linea di volo e' invece allineata con l'azimuth del
sole (o il suo reciproco) - "vola con il sole davanti o dietro, non di
lato". Questo script:

1. Calcola la posizione del sole (azimuth, elevazione) per la data/ora/
   posizione fornite, con la libreria 'astral'.
2. Usa il modello BRDF semi-empirico RPV (lo stesso di
   brdf_hotspot_validation.py) per stimare quantitativamente l'ampiezza
   dell'asimmetria across-track attesa per OGNI possibile azimuth di volo
   (0-179 gradi, dato che una linea di volo e' bidirezionale e quindi il
   problema ha periodo 180 gradi), tenendo conto del FOV del sensore.
3. Riporta l'azimuth ottimale (minima asimmetria), quello peggiore
   (massima asimmetria) e un grafico dell'ampiezza attesa in funzione
   dell'azimuth di volo scelto, cosi' puoi vedere anche quanto e'
   'largo' il minimo (cioe' quanto puoi discostarti dall'ottimo restando
   comunque in una zona a basso impatto).

Uso
---
    python3 optimal_flight_azimuth.py --lat 43.7696 --lon 11.2558 \\
        --date 2026-07-22 --time 11:46 --tz Europe/Rome

    # con FOV/parametri RPV personalizzati:
    python3 optimal_flight_azimuth.py --lat 43.7243059 --lon 10.2791424 \\
        --date 2026-07-10 --time 11:47 --tz Europe/Rome \\
        --fov 24.75 --rpv-theta -0.15 --rpv-h 0.1
"""

import argparse
from datetime import datetime
from zoneinfo import ZoneInfo

import numpy as np
from astral import LocationInfo
from astral.sun import elevation as sun_elevation_fn, azimuth as sun_azimuth_fn
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt


# ---------------------------------------------------------------------
# Modello RPV (identico a brdf_hotspot_validation.py)
# ---------------------------------------------------------------------

def rpv_reflectance(theta_i_deg, theta_r_deg, phi_deg, k=1.0, theta_hg=-0.15, h=0.1, rho0=1.0):
    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)

    cos_g = cos_ti * cos_tr + sin_ti * sin_tr * np.cos(phi)
    cos_g = np.clip(cos_g, -1.0, 1.0)

    F = (1 - theta_hg**2) / (1 + 2 * theta_hg * cos_g + theta_hg**2) ** 1.5

    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)

    M = (cos_ti * cos_tr * (cos_ti + cos_tr)) ** (k - 1)

    return rho0 * M * F * H


def brdf_amplitude_for_flight_azimuth(flight_az, sun_az, sun_zenith, fov_deg,
                                       k=1.0, theta_hg=-0.15, h=0.1, n_samples=480):
    """
    Per un dato azimuth di volo, calcola l'ampiezza (escursione normalizzata
    max-min) del fattore BRDF atteso sullo swath, dato il FOV del sensore.
    """
    across_1 = (flight_az + 90) % 360
    across_2 = (flight_az - 90) % 360

    half_fov = fov_deg / 2.0
    col = np.arange(n_samples)
    center = (n_samples - 1) / 2.0
    view_zenith = np.abs((col - center) / center) * half_fov
    offset = col - center
    view_az = np.where(offset < 0, across_1, across_2)

    phi = view_az - sun_az
    R = rpv_reflectance(sun_zenith, view_zenith, phi, k=k, theta_hg=theta_hg, h=h)
    R_norm = R / np.mean(R)
    return R_norm.max() - R_norm.min(), R_norm


# ---------------------------------------------------------------------
# Posizione solare
# ---------------------------------------------------------------------

def get_sun_position(lat, lon, date_str, time_str, tz_str):
    tz = ZoneInfo(tz_str)
    hh, mm = [int(v) for v in time_str.split(":")]
    y, m, d = [int(v) for v in date_str.split("-")]
    dt = datetime(y, m, d, hh, mm, tzinfo=tz)

    loc = LocationInfo(latitude=lat, longitude=lon)
    el = sun_elevation_fn(loc.observer, dt)
    az = sun_azimuth_fn(loc.observer, dt)
    return dt, az, el


# ---------------------------------------------------------------------
# Analisi principale
# ---------------------------------------------------------------------

def analyze(lat, lon, date_str, time_str, tz_str, fov_deg, rpv_k, rpv_theta, rpv_h, out_path):
    dt, sun_az, sun_el = get_sun_position(lat, lon, date_str, time_str, tz_str)
    sun_zenith = 90 - sun_el

    print(f"Data/ora locale: {dt}")
    print(f"Posizione: lat={lat}, lon={lon}")
    print(f"Sole: azimuth={sun_az:.2f} deg, elevazione={sun_el:.2f} deg (zenith={sun_zenith:.2f} deg)")

    if sun_el <= 0:
        print("\n[attenzione] Il sole e' sotto l'orizzonte a quest'ora: nessun volo possibile/sensato.")
        return

    print(f"FOV assunto: {fov_deg:.2f} deg  |  parametri RPV: k={rpv_k}, Theta={rpv_theta}, h={rpv_h}")

    # scansione di tutti gli azimuth di volo possibili (0-179, periodo 180 gradi)
    azimuths = np.arange(0, 180, 0.5)
    amplitudes = np.array([
        brdf_amplitude_for_flight_azimuth(az, sun_az, sun_zenith, fov_deg, rpv_k, rpv_theta, rpv_h)[0]
        for az in azimuths
    ])

    best_idx = np.argmin(amplitudes)
    worst_idx = np.argmax(amplitudes)
    best_az = azimuths[best_idx]
    worst_az = azimuths[worst_idx]

    print("\n--- Risultato ---")
    print(f"Azimuth di volo OTTIMALE (minima asimmetria BRDF): {best_az:.1f} deg "
          f"(reciproco: {(best_az+180)%360:.1f} deg)  ->  escursione prevista: {100*amplitudes[best_idx]:.1f}%")
    print(f"Azimuth di volo PEGGIORE (massima asimmetria BRDF): {worst_az:.1f} deg "
          f"(reciproco: {(worst_az+180)%360:.1f} deg)  ->  escursione prevista: {100*amplitudes[worst_idx]:.1f}%")
    print(f"\n(per confronto: l'azimuth ottimale teorico e' semplicemente l'azimuth del sole stesso, "
          f"{sun_az:.1f} deg mod 180 = {sun_az % 180:.1f} deg - 'vola con il sole davanti o dietro')")

    # ampiezza per un eventuale azimuth di volo specifico gia' pianificato
    # (utile come confronto, mostrata solo nel grafico)

    fig, ax = plt.subplots(figsize=(9, 5.5))
    ax.plot(azimuths, 100 * amplitudes, color="darkred", linewidth=1.8)
    ax.axvline(best_az, color="green", linestyle="--", label=f"ottimale: {best_az:.1f} deg ({100*amplitudes[best_idx]:.1f}%)")
    ax.axvline(worst_az, color="crimson", linestyle="--", label=f"peggiore: {worst_az:.1f} deg ({100*amplitudes[worst_idx]:.1f}%)")
    ax.set_xlabel("Azimuth di volo (gradi, periodo 180°)")
    ax.set_ylabel("Escursione BRDF attesa sullo swath (%)")
    ax.set_title(f"Ampiezza asimmetria BRDF attesa in funzione dell'azimuth di volo\n"
                 f"{dt.strftime('%Y-%m-%d %H:%M %Z')} @ lat={lat}, lon={lon} "
                 f"(sole: az={sun_az:.1f}°, el={sun_el:.1f}°)")
    ax.legend(fontsize=9)
    ax.grid(alpha=0.3)
    fig.tight_layout()
    fig.savefig(out_path, dpi=150)
    plt.close(fig)
    print(f"\nGrafico salvato in: {out_path}")


def main():
    parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument("--lat", type=float, required=True, help="Latitudine (gradi decimali)")
    parser.add_argument("--lon", type=float, required=True, help="Longitudine (gradi decimali)")
    parser.add_argument("--date", required=True, help="Data locale, formato YYYY-MM-DD")
    parser.add_argument("--time", required=True, help="Ora locale, formato HH:MM")
    parser.add_argument("--tz", default="Europe/Rome", help="Timezone IANA (default: Europe/Rome)")
    parser.add_argument("--fov", type=float, default=24.75, help="FOV across-track del sensore, gradi (default: 24.75)")
    parser.add_argument("--rpv-k", type=float, default=1.0, help="Parametro RPV k (default 1.0)")
    parser.add_argument("--rpv-theta", type=float, default=-0.15, help="Parametro RPV Theta (default -0.15)")
    parser.add_argument("--rpv-h", type=float, default=0.1, help="Parametro RPV h (default 0.1)")
    parser.add_argument("--out", default="optimal_flight_azimuth.png", help="Percorso del grafico di output")
    args = parser.parse_args()

    analyze(args.lat, args.lon, args.date, args.time, args.tz,
             args.fov, args.rpv_k, args.rpv_theta, args.rpv_h, args.out)


if __name__ == "__main__":
    main() 

BRDF FigSpec FS60-C

Visto che nel volo ad Empoli era emersa l'importanza della BRDF nei dati acquisiti ho voluto provare con un volo precedente

Si tratta di un volo effettuato il 10 luglio 2026 ore 11:47 presso San Rossore (43.7243059,10.2791424) 

L'azimuth medio e' di 345.6 gradi 

Il Sole aveva azimuth 129.1° ed  elevazione 60.6°  

Lo scarto tra la track di volo ed il piano principale solare e' di 53.5 gradi (molto meglio rispetto ad Empoli l'angolo era di circa 8 gradi)

Usando il modello RPF si ha che tra i due lati dell'immagine si ha uno scarto di circa 11% contro il 23% di Empoli 

 Il problema in questo caso e' che si tratta di alternanza di tracce ascendenti e discendenti con stesso azimuth ...quindi si invertono ad ogni passate le condizioni di BRDF 

 

 

 

Ortorettifica FIgSpec FS-60C

Per georiferire il dato di FigSpec FS60-C si puo' usare il file figspec.gps

(da notare che il file gps e' un csv che include anche dei campi sull'assetto di volo del drone yaw,pitch and roll che sono valorizzati sempre a zero) 

Da un primo volo con una risoluzione spaziale a 480 px, 1 nm e quota di volo 120 m AGL si ha un GSD di 5 cm al Nadir 

Considerando che il volo ha ripreso 135.5  m x 68.8 m quindi circa 9322 metri quadri si ha una produzione di 88 Mb a metro quadro

python3 orthorectify_pushbroom.py true_color.png figspec.gps --flip-across-track --fov 12.375

 

#!/usr/bin/env python3

"""

orthorectify_pushbroom.py

Ortorettifica un'immagine true-color pushbroom usando la traccia GPS/IMU

(CSV) e il FOV noto del sensore, tramite georeferenziazione diretta (Direct

Georeferencing, DG) e resampling su una griglia regolare UTM.

Come funziona (in breve)

-------------------------

1. Interpola la posizione (lat, lon, altitudine) del velivolo per OGNI riga

dell'immagine, a partire dai fix GPS sparsi nel file CSV (che ne contiene

uno ogni N righe, non uno per riga).

2. Stima la prua/heading del volo dalla traccia GPS stessa (bearing tra punti

consecutivi, o un singolo heading globale se il volo e' una passata

dritta), dato che i campi Pitch/Roll/Yaw nel file GPS risultano spesso

non popolati (costanti a 0) nei log di alcuni sistemi.

3. Per ogni pixel (riga, colonna), calcola l'angolo di vista rispetto al

nadir (assumendo FOV a distribuzione angolare lineare/equiangolare sulle

colonne) e proietta il raggio di vista a terra assumendo TERRENO PIATTO

all'altitudine data (nessun DEM), con velivolo puntato a nadir (pitch=

roll=0, cioe' nessuna correzione di assetto oltre l'heading).

4. Converte le coordinate di terra cosi' ottenute (lat/lon per ogni pixel,

la "griglia di geolocazione") in coordinate UTM metriche, e ricampiona

(resample) i pixel sorgente su una griglia regolare in UTM per produrre

un'ortofoto vera e propria, salvata come GeoTIFF georeferenziato.

Assunzioni e limiti (importanti)

----------------------------------

- NESSUN modello del terreno (DEM): il terreno e' assunto piatto

all'altitudine di volo fornita. Su terreno con rilievo significativo

introduce "relief displacement" (spostamento radiale dal nadir

proporzionale al dislivello). Per un rilievo di pianura/campo agricolo

come questo, l'errore e' generalmente trascurabile.

- Nessuna correzione di boresight/misalignment tra IMU e camera (si assume

il sensore perfettamente allineato con l'asse di volo).

- Pitch e Roll sono assunti nulli (nadir stabilizzato) se non disponibili

nel file GPS o se costanti a zero (tipico quando quel canale non e'

popolato dal logger). Se il tuo file GPS ha valori di pitch/roll validi

(non tutti zero), lo script li usa automaticamente.

- Il FOV e' assunto equiangolare e lineare sulle colonne (buona

approssimazione per sensori pushbroom a stretto campo di vista).

- La direzione "sinistra/destra" del sensore rispetto alla direzione di

volo (cioe' se la colonna 0 e' a destra o sinistra del velivolo) NON e'

nota a priori: usa --flip-across-track se l'output risulta specchiato

rispetto alla realta' (confronta con una mappa/basemap nota).

Uso

---

python3 orthorectify_pushbroom.py true_color.png figspec.gps \\

--fov 24.75 --out ortho_true_color.tif

# se l'immagine esce specchiata lateralmente:

python3 orthorectify_pushbroom.py true_color.png figspec.gps \\

--fov 24.75 --flip-across-track --out ortho_true_color.tif

# con dimensione pixel di output forzata (metri):

python3 orthorectify_pushbroom.py true_color.png figspec.gps \\

--fov 24.75 --pixel-size 0.05 --out ortho_true_color.tif

"""

import argparse

import sys

from pathlib import Path

import numpy as np

import pandas as pd

from PIL import Image

import pyproj

import rasterio

from rasterio.transform import Affine

from scipy.ndimage import distance_transform_edt

import matplotlib

matplotlib.use("Agg")

import matplotlib.pyplot as plt

# ---------------------------------------------------------------------

# GPS: interpolazione posizione per riga + heading

# ---------------------------------------------------------------------

def load_and_interpolate_gps(csv_path, n_rows, line_col="Lines",

lat_col="Latitude", lon_col="Longitude", alt_col="Altitude"):

"""

Interpola lat/lon/altitudine per ogni riga dell'immagine (0..n_rows-1),

a partire dai fix GPS registrati (uno ogni N righe). Extrapolazione

costante (clamp) fuori dal range coperto dai fix, con avviso.

"""

df = pd.read_csv(csv_path).sort_values(line_col)

lines = df[line_col].values.astype(float)

lat = df[lat_col].values.astype(float)

lon = df[lon_col].values.astype(float)

alt = df[alt_col].values.astype(float)

row_idx = np.arange(n_rows, dtype=float)

if row_idx.min() < lines.min() or row_idx.max() > lines.max():

print(f"[attenzione] Le righe dell'immagine (0-{n_rows-1}) escono dal range coperto "

f"dai fix GPS ({lines.min():.0f}-{lines.max():.0f}); le righe fuori range "

f"useranno il valore del fix piu' vicino (estrapolazione costante).")

lat_i = np.interp(row_idx, lines, lat)

lon_i = np.interp(row_idx, lines, lon)

alt_i = np.interp(row_idx, lines, alt)

return lat_i, lon_i, alt_i, df

def heading_from_track(lat, lon, mode="global"):

"""

Stima l'heading (azimuth, gradi da nord) per ogni riga.

mode='global': un unico heading (PCA su tutta la traccia) per tutte le righe

- robusto, adatto a voli rettilinei (il caso tipico).

mode='local' : bearing punto-punto (differenza centrata), utile se il

volo non e' perfettamente dritto.

"""

n = len(lat)

lat0 = np.radians(np.mean(lat))

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

if mode == "global":

pts_c = pts - pts.mean(axis=0)

cov = np.cov(pts_c.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

# verifica il verso (potrebbe essere az o az+180): usa il verso del

# moto reale (dal primo all'ultimo punto) per orientarlo correttamente

overall = np.degrees(np.arctan2(x[-1] - x[0], y[-1] - y[0])) % 360

if min(abs(az - overall), 360 - abs(az - overall)) > 90:

az = (az + 180) % 360

return np.full(n, az)

# locale: bearing centrato, con gestione dei bordi

heading = np.zeros(n)

for i in range(n):

i0 = max(0, i - 1)

i1 = min(n - 1, i + 1)

dx = x[i1] - x[i0]

dy = y[i1] - y[i0]

heading[i] = np.degrees(np.arctan2(dx, dy)) % 360

return heading

# ---------------------------------------------------------------------

# Geometria pushbroom: da (riga, colonna) a (lat, lon) di terra

# ---------------------------------------------------------------------

def build_geolocation_grid(lat, lon, alt_agl, heading_deg, n_cols, fov_deg,

pitch_deg=None, roll_deg=None, flip_across_track=False):

"""

Calcola lat/lon di terra per ogni pixel (n_rows x n_cols), assumendo

terreno piatto e velivolo nadir-stabilizzato (salvo pitch/roll forniti).

Ritorna (lat_grid, lon_grid) shape (n_rows, n_cols).

"""

n_rows = len(lat)

half_fov = fov_deg / 2.0

col_idx = np.arange(n_cols)

center = (n_cols - 1) / 2.0

theta_c = (col_idx - center) / center * half_fov # gradi, - a sx, + a dx (o viceversa)

if flip_across_track:

theta_c = -theta_c

theta_c_rad = np.radians(theta_c)

if pitch_deg is None:

pitch_deg = np.zeros(n_rows)

if roll_deg is None:

roll_deg = np.zeros(n_rows)

R_earth = 6371000.0

lat_grid = np.empty((n_rows, n_cols))

lon_grid = np.empty((n_rows, n_cols))

psi = np.radians(heading_deg) # heading (yaw), rad

right_az = psi + np.pi / 2 # direzione "a destra" del volo (azimuth)

for r in range(n_rows):

h_agl = alt_agl[r]

if h_agl <= 0:

h_agl = np.nan # evita risultati assurdi se l'altitudine e' invalida

# NB: pitch/roll qui non applicati con rotazione completa (si assume

# nadir); un'estensione futura puo' aggiungere la rotazione 3D

# completa (roll attorno all'asse di avanzamento, pitch attorno

# all'asse trasversale) se il file GPS fornisce valori attendibili.

ground_offset = h_agl * np.tan(theta_c_rad) # metri, lungo l'asse across-track

north_off = ground_offset * np.cos(right_az[r])

east_off = ground_offset * np.sin(right_az[r])

dlat = (north_off / R_earth) * (180.0 / np.pi)

dlon = (east_off / (R_earth * np.cos(np.radians(lat[r])))) * (180.0 / np.pi)

lat_grid[r, :] = lat[r] + dlat

lon_grid[r, :] = lon[r] + dlon

return lat_grid, lon_grid

# ---------------------------------------------------------------------

# Proiezione UTM e resampling su griglia regolare

# ---------------------------------------------------------------------

def utm_epsg_from_lonlat(lon, lat):

zone = int((lon + 180) / 6) + 1

return (32600 + zone) if lat >= 0 else (32700 + zone)

def fill_small_gaps(image_uint8, filled_mask, max_fill_dist_px):

"""

Riempie i piccoli vuoti isolati (1-2 pixel) con il valore del pixel

valido piu' vicino, entro una distanza massima (in pixel di output).

Questo NON riempie le grandi zone nodata ai bordi/angoli del bounding

box (quelle sono corrette: lo swath e' un rettangolo ruotato rispetto

agli assi E/N della griglia UTM, quindi il bounding box allineato agli

assi contiene inevitabilmente due triangoli vuoti agli angoli - non e'

un difetto, e' la geometria della rotazione). Serve solo a colmare

micro-vuoti dovuti a piccola sotto-densita' locale di campionamento.

"""

empty = ~filled_mask

if not empty.any():

return image_uint8, filled_mask

dist, (idx_r, idx_c) = distance_transform_edt(empty, return_indices=True)

fillable = empty & (dist <= max_fill_dist_px)

out = image_uint8.copy()

new_mask = filled_mask.copy()

for ch in range(image_uint8.shape[2]):

channel = image_uint8[:, :, ch]

filled_channel = channel[idx_r, idx_c]

ch_out = out[:, :, ch]

ch_out[fillable] = filled_channel[fillable]

out[:, :, ch] = ch_out

new_mask[fillable] = True

return out, new_mask

def orthorectify(image_path, gps_csv, out_path, fov_deg, nodata_thresh=3.0,

pixel_size=None, flip_across_track=False, heading_mode="global",

agl_height_override=None, line_col="Lines", fill_gap_px=3):

img = Image.open(image_path).convert("RGB")

arr = np.asarray(img).astype(np.float64)

n_rows, n_cols = arr.shape[0], arr.shape[1]

print(f"Immagine caricata: {n_rows} righe x {n_cols} colonne")

luma = arr.mean(axis=2)

valid_mask = luma > nodata_thresh

lat, lon, alt, gps_df = load_and_interpolate_gps(gps_csv, n_rows, line_col=line_col)

if agl_height_override is not None:

alt = np.full(n_rows, agl_height_override)

print(f"Altitudine AGL forzata a {agl_height_override} m per tutte le righe")

else:

print(f"Altitudine AGL dal file GPS: min={alt.min():.1f} m, max={alt.max():.1f} m, media={alt.mean():.1f} m")

pitch = gps_df["Pitch"].values.astype(float) if "Pitch" in gps_df.columns else None

roll = gps_df["Roll"].values.astype(float) if "Roll" in gps_df.columns else None

if pitch is not None and np.allclose(pitch, 0) and roll is not None and np.allclose(roll, 0):

print("[info] Pitch/Roll nel file GPS sono costanti a zero (probabilmente non popolati): "

"assumo velivolo nadir-stabilizzato (nessuna correzione di assetto oltre l'heading).")

pitch_interp = roll_interp = None

elif pitch is not None and roll is not None:

row_idx = np.arange(n_rows, dtype=float)

lines = gps_df[line_col].values.astype(float)

pitch_interp = np.interp(row_idx, lines, pitch)

roll_interp = np.interp(row_idx, lines, roll)

print("[info] Uso i valori di Pitch/Roll forniti nel file GPS (non applicata rotazione 3D completa "

"in questa versione semplificata: solo heading applicato).")

else:

pitch_interp = roll_interp = None

heading = heading_from_track(lat, lon, mode=heading_mode)

print(f"Heading stimato: min={heading.min():.1f} deg, max={heading.max():.1f} deg "

f"({'costante (modo global)' if heading_mode == 'global' else 'variabile (modo local)'})")

print(f"FOV: {fov_deg:.2f} deg | flip_across_track={flip_across_track}")

lat_grid, lon_grid = build_geolocation_grid(

lat, lon, alt, heading, n_cols, fov_deg,

pitch_deg=pitch_interp, roll_deg=roll_interp, flip_across_track=flip_across_track

)

# --- proiezione in UTM ---

epsg = utm_epsg_from_lonlat(np.nanmean(lon_grid), np.nanmean(lat_grid))

print(f"CRS di output scelto: EPSG:{epsg} (UTM automatico dal baricentro dell'area)")

transformer = pyproj.Transformer.from_crs("EPSG:4326", f"EPSG:{epsg}", always_xy=True)

easting, northing = transformer.transform(lon_grid, lat_grid)

valid = valid_mask & np.isfinite(easting) & np.isfinite(northing)

if not valid.any():

raise ValueError("Nessun pixel valido da ortorettificare (controlla nodata_thresh o i dati GPS)")

e_min, e_max = easting[valid].min(), easting[valid].max()

n_min, n_max = northing[valid].min(), northing[valid].max()

print(f"Estensione ortofoto: E [{e_min:.1f}, {e_max:.1f}] m, N [{n_min:.1f}, {n_max:.1f}] m "

f"({e_max-e_min:.1f} x {n_max-n_min:.1f} m)")

if pixel_size is None:

# GSD approssimato al nadir: altitudine media * (FOV in rad / n_cols)

gsd = float(np.nanmean(alt)) * np.radians(fov_deg) / n_cols

pixel_size = max(gsd, 0.01)

print(f"Dimensione pixel di output (auto, da GSD al nadir): {pixel_size:.4f} m")

else:

print(f"Dimensione pixel di output (specificata): {pixel_size:.4f} m")

out_cols = int(np.ceil((e_max - e_min) / pixel_size)) + 1

out_rows = int(np.ceil((n_max - n_min) / pixel_size)) + 1

print(f"Dimensioni griglia di output: {out_rows} righe x {out_cols} colonne")

# --- binning (splat) dei pixel sorgente sulla griglia di output ---

col_out = ((easting - e_min) / pixel_size).astype(np.int64)

row_out = ((n_max - northing) / pixel_size).astype(np.int64) # riga 0 = Nord (northing massimo)

ok = valid & (col_out >= 0) & (col_out < out_cols) & (row_out >= 0) & (row_out < out_rows)

sum_rgb = np.zeros((out_rows, out_cols, 3), dtype=np.float64)

count = np.zeros((out_rows, out_cols), dtype=np.int64)

flat_row = row_out[ok]

flat_col = col_out[ok]

flat_idx = flat_row * out_cols + flat_col

for ch in range(3):

vals = arr[:, :, ch][ok]

sums = np.bincount(flat_idx, weights=vals, minlength=out_rows * out_cols)

sum_rgb[:, :, ch] = sums.reshape(out_rows, out_cols)

counts_flat = np.bincount(flat_idx, minlength=out_rows * out_cols)

count = counts_flat.reshape(out_rows, out_cols)

with np.errstate(invalid="ignore", divide="ignore"):

ortho = sum_rgb / count[:, :, np.newaxis]

nodata_out = count == 0

ortho[nodata_out] = 0

ortho_uint8 = np.clip(ortho, 0, 255).astype(np.uint8)

coverage_pct = 100 * (~nodata_out).sum() / nodata_out.size

print(f"Copertura griglia di output (prima del gap-fill): {coverage_pct:.1f}% delle celle riempite. "

f"NB: una copertura ben sotto il 100% e' normale e attesa: lo swath e' un rettangolo ruotato "

f"rispetto agli assi E/N della griglia UTM, quindi il bounding box allineato agli assi include "

f"inevitabilmente due triangoli vuoti agli angoli (non e' un artefatto di ricampionamento).")

if fill_gap_px > 0:

valid_out_mask = ~nodata_out

ortho_uint8, valid_out_mask = fill_small_gaps(ortho_uint8, valid_out_mask, fill_gap_px)

nodata_out = ~valid_out_mask

coverage_pct_after = 100 * valid_out_mask.sum() / valid_out_mask.size

print(f"Copertura dopo gap-fill dei soli micro-vuoti (raggio max {fill_gap_px} px): {coverage_pct_after:.1f}%")

# --- scrittura GeoTIFF ---

transform = Affine.translation(e_min, n_max) * Affine.scale(pixel_size, -pixel_size)

out_path = Path(out_path)

out_path.parent.mkdir(parents=True, exist_ok=True)

with rasterio.open(

out_path, "w",

driver="GTiff",

height=out_rows, width=out_cols, count=3,

dtype=np.uint8, crs=f"EPSG:{epsg}", transform=transform,

nodata=0,

compress="deflate",

) as dst:

for ch in range(3):

dst.write(ortho_uint8[:, :, ch], ch + 1)

print(f"\nOrtofoto GeoTIFF salvata in: {out_path}")

# --- anteprima PNG rapida ---

preview_path = out_path.with_suffix(".preview.png")

fig, ax = plt.subplots(figsize=(8, 8 * out_rows / max(out_cols, 1)))

ax.imshow(ortho_uint8, extent=[e_min, e_min + out_cols * pixel_size,

n_max - out_rows * pixel_size, n_max])

ax.set_xlabel(f"Easting (m, EPSG:{epsg})")

ax.set_ylabel("Northing (m)")

ax.set_title("Anteprima ortofoto")

fig.tight_layout()

fig.savefig(preview_path, dpi=150)

plt.close(fig)

print(f"Anteprima salvata in: {preview_path}")

return out_path

def main():

parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)

parser.add_argument("image", help="Immagine true-color pushbroom (PNG/JPEG/TIFF)")

parser.add_argument("gps_csv", help="File CSV GPS/IMU (con colonne Lines, Latitude, Longitude, Altitude, ...)")

parser.add_argument("--out", default="ortho.tif", help="Percorso del GeoTIFF di output")

parser.add_argument("--fov", type=float, required=True, help="FOV across-track totale del sensore (gradi)")

parser.add_argument("--nodata-thresh", type=float, default=3.0, help="Soglia luma per nodata (default 3.0)")

parser.add_argument("--pixel-size", type=float, default=None, help="Dimensione pixel output in metri (default: auto da GSD al nadir)")

parser.add_argument("--flip-across-track", action="store_true",

help="Inverte il lato sinistra/destra (usa se l'output risulta specchiato)")

parser.add_argument("--heading-mode", choices=["global", "local"], default="global",

help="'global': heading unico stimato via PCA su tutta la traccia (default, adatto a voli dritti). "

"'local': bearing punto-punto (per voli non rettilinei)")

parser.add_argument("--agl-height", type=float, default=None,

help="Forza l'altitudine AGL (m) per tutte le righe invece di usare la colonna Altitude del GPS")

parser.add_argument("--line-col", default="Lines", help="Nome della colonna indice-riga nel CSV GPS (default: Lines)")

parser.add_argument("--fill-gap-px", type=int, default=3,

help="Raggio massimo (in pixel di output) per il gap-fill nearest-neighbor "

"dei piccoli vuoti da aliasing griglia ruotata (default: 3, 0=disattiva)")

args = parser.parse_args()

orthorectify(args.image, args.gps_csv, args.out, args.fov,

nodata_thresh=args.nodata_thresh, pixel_size=args.pixel_size,

flip_across_track=args.flip_across_track, heading_mode=args.heading_mode,

agl_height_override=args.agl_height, line_col=args.line_col,

fill_gap_px=args.fill_gap_px)

if __name__ == "__main__":

main()