Tesi di Laurea
La fotogrammetria digitale multi-temporale d'archivio per l'analisi delle variazioni planimetriche costiere

Costituzione di un nuovo database geografico.



Il presente paragrafo nasce come risposta ad una necessità, non prevista all’inizio di questo studio, che potremmo definire come il bisogno di unificare l’informazione geografica sotto un unico sistema di rappresentazione e gestione. Come si è già scritto, la copertura informativa dell’area deltizia è costituita da quattro campagne di rilievo aereofotogrammetrico, di recente produzione, quindi oggetto di nuova sistemazione informatica. Ma nell’ambito del presente studio, si è appalesata la possibilità/necessità di recuperare anche tutti i lavori antecedenti, e di inserirli all’interno di un unico database geografico di nuova costituzione. Onde evitare confusione tra le fonti a disposizione, vecchie e nuove che siano, si riportano nelle seguenti due tabelle, i livelli informativi geografici trattati e disponibili, con le loro caratteristiche peculiari.

Di primo acchito la questione potrebbe essere affrontata importando direttamente in ArcMap tutti i vari dati, ma purtroppo questo si scontra con un problema annoso, costituito dal diverso sistema di riferimento utilizzato, difatti l’area del Delta è nella zona di sovrapposizione dei due fusi, che possono così essere usati entrambi, ed è esattamente quello che accaduto; i nuovi dati si trovano sul fuso Ovest a differenza dei vecchi. La particolare condizione richiede una certa attenzione sulle deformazioni lineari presenti, considerando che l’area di analisi può essere inscritta in un rettangolo di (40×20)Km, . Dalla geodesia sappiamo che il Nord di rete è fissato dalla tangente al meridiano centrale del fuso, nel punto d’emanazione del sistema di riferimento cartografico, ebbene, avendo a che fare con due fusi diversi, le tangenti saranno certamente differenti nella loro direzione, perciò la trasformazione non sarà una semplice traslazione, ma avremmo pure una rotazione, accompagnata da deformazioni lineari dovute al diverso cilindro di proiezione utilizzato.  

Cosicché dopo aver provato metodi manuali che prevedevano l’estrazione diretta delle coordinate dai contenuti vettoriali, ma senza risultati soddisfacenti, date le difficoltà di ripristinare il contenuto vettoriale medesimo, si è deciso di usare il GIS della Autodesk: Autocad Map 3D. Questo software consente di assegnare un sistema di riferimento cartografico, o geografico, ai contenuti vettoriali memorizzanti in normali file .dwg, quindi è permessa una perfetta inter-operatività con gli altri software della Autodesk. La possibilità di georiferire i file vettoriali, ne consente pure il passaggio da un sistema di riferimento ad un altro, con semplici procedure (Fig. 5.6-5.7). 

Fig. 5.6 – I passaggi necessari per assegnare il sistema di riferimento al dwg contente la linea di riva 2.008, evidentemente deve essere coerente con quelle che sono le coordinate metriche del file. Per far questo si agisce sull’icona nella parte inferiore dell’area di lavoro [1], che consente di selezionare il sistema di coordinate desiderato [2]. Si ponga attenzione che qui il sistema Gauss-Boaga è espresso in termini internazionali, quindi il Fuso Est, corrisponde alla Zona 2. Fatto questo si procede al salvataggio del file di disegno, che ora sarà correttamente georiferito.

Fig. 5.7 – Il percorso logico continua con la creazione di un nuovo progetto, al quale andrà assegnato il sistema di riferimento di destinazione [3]-[4]; ora possiamo importare nuove carte nell’area di lavoro, cioè altri file di disegno esterni [5]. Il vero processo di trasformazione avviene con l’uso delle query [6], difatti dopo aver impostato tutte le opzioni come al passo [7], il programma procede alla fase di disegno, ove tutte le cartografie associate vengono ri-proiettate nel nuovo sistema di riferimento indicato all’inizio. Ora non ci rimane altro che salvare il lavoro svolto.

Le fasi operative appena esposte consentono di assegnare un nuovo sistema di riferimento ai dati di natura vettoriale, che potranno essere importati direttamente all’interno del nuovo database geografico in ArcMap, ed organizzati in opportuni layers, assieme a tutti gli altri livelli cartografici. Non è certamente possibile fare altrettanto con i contenuti raster (le immagini): il processo di ri-proiezione comportano delle deformazioni non accettabili, che si manifestano in spostamenti diversi dalla semplice roto-traslazione rigida, caso che potrebbe ancora essere affrontato con una certa facilità. Si tratta tecnicamente di una trasformazione geometrica tra due sistemi di riferimento piani, con la successiva riassegnazione del contenuto radiometrico alle nuove posizioni (pixel) definite dalla trasformazione.  


LA TRASFORMAZIONE GEOMETRICA PIANA.



Ogni trasformazione si definisce come quella serie di operazioni grafiche o numeriche, che fanno corrispondere biunivocamente ad un insieme piano di punti, un altro insieme piano avente lo stesso numero di elementi. Si noti che il numero di elementi trasformati è pari al numero delle entità puntuali presenti nell’insieme d’origine, difatti le procedure analitiche che ne consentano l’incremento sono già state affrontate nella sezione dedicata alla fotogrammetria, e sono le tecniche ricampionamento, che nulla hanno a che spartire con le trasformazioni geometriche propriamente definite. Anche se è pur vero che data la natura discreta dell’immagine digitale, spesso l’elemento trasformato, in termini di contenuto radiometrico, non può essere direttamente assegnato alla nuova posizione, per cui si utilizzano delle tecniche di riassegnazione, che sono ancora una volta i metodi di ricampionamento, ma sono invero pur sempre dei metodi analitici armonizzanti eseguiti successivamente alla trasformazione geometrica. Qualsiasi trasformazione piana trae origine, in termini cartesiani, da un polinomio di ordine n, che pone in relazione due sistemi di riferimento, [O;ξ,η] e [O';ξ',η' ]:

Da queste relazioni, del tutto generali, discendono tutte le possibili trasformazioni piane, che si possono classificare in due gruppi.

• Trasformazioni globali. 
• Trasformazioni locali. 

Le prime riferiscono la loro azione all’intero piano immagine, quindi i parametri trasformatori avranno validità globale, mentre le trasformazioni locali, suddividono il dominio in sottodomini d’azione, quindi l’immagine viene ad essere suddivisa in aree, ed i parametri saranno validi solo per la loro specifica area di riferimento. In entrambi i casi, la determinazione dei parametri della trasformazione avviene di norma conoscendo le coordinate di un certo numero di punti di controllo a terra (GCP), sulla base dei quali il processo determina il miglior adattamento tra i due piani immagine. Il numero minimo di punti di controllo richiesto è legato al tipo di trasformazione piana adottata, in quanto dipende dal numero di parametri della trasformazione stessa. Evidentemente è sempre opportuno lavorare con un numero di punti superiore allo stretto necessario, poiché solo in tal modo è consentita una valutazione della qualità della trasformazione, ad esempio mediante l’analisi dei residui di una compensazione ai minimi quadrati, come vedremmo. Le possibili trasformazioni di tipo globale si differenziano tra loro dal numero dei parametri coinvolti, e dunque per il tipo di deformazione che sono in grado di correggere. 

• Conforme o di Helmert.

Trattasi di una trasformazione a quattro parametri, due traslazioni nelle direzioni principali (ξ(0)',η(0)'), un fattore di scala isotropo m, ed una rotazione.

• affine.

Corrisponde ad una trasformazione polinomiale del primo grado, essa è composta da cinque/sei parametri: le due traslazioni, una variazione di scala anisotropa, una o due rotazioni (α per la rotazione delle direzioni orizzontali e β per quelle verticali, le due rotazioni possono essere anche identiche).

• Proiettiva od omografica.

Trasformazione ad otto parametri molto utilizzata in ambito fotogrammetrico, direttamente derivata dalle equazioni di collinearità nel caso di oggetti ripresi piani, per cui dai nove parametri iniziali del caso più generale, ci si riduce a otto.

• Polinomiale di ordine n.

Il numero di parametri è in funzione dell’ordine, ed è sostanzialmente il caso più generale, rappresentato dalla relazione più generale.

Tutte queste trasformazioni globali si possono classificare in lineari ed in non lineari, le prime si esplicano in una traslazione, una rotazione globale, o rotazione con scorrimento (affine a sei parametri), ed una variazione di scala (isotropa o anisotropa). Mentre per una trasformazione omografica il risultato è un’immagine prospettica, risultato della proiezione di un piano su di un altro. Tutte queste sono delle correzioni lineari, che trovano applicazione nell’affrontare gli effetti di scala, offest, rotazione e riflessione; una trasformazione globale non lineare invece corregge distorsioni non lineari, in questo caso i risultati sono profondamente legati al numero dei punti di controllo e alla loro distribuzione. Il numero minimo di GCP per una trasformazione polinomiale dipende dall’ordine n.

Dall’altro lato una trasformazione locale è una trasformazione i cui parametri vengono determinati per aree specifiche, le quali dovranno essere le più piccole possibili. Il principale vantaggio di questo metodo è la non dipendenza dal tipo di deformazione presente nel processo, la cui conoscenza non sempre è nota, ma sarebbe richiesta se stessimo per applicare una trasformazione globale, il cui adattamento non può prescindere dalla natura dei dati d’origine. Il numero di punti di controllo richiesto dipende dal modello matematico applicato, che usualmente è lo stesso delle trasformazioni globali, ma viene circoscritto ad aree limitate. La trasformazione locale può essere condotta con vari metodi, nel presente studio si è usato il metodo della scomposizione del dominio in elementi finti, più comunemente conosciuto con il nome di foglio di gomma (rubber-sheeting lineare o non lineare). In altri termini il dominio, il piano immagine, viene scomposto in una mesh triangolare o quadrangolare, spesso il metodo di meshing utilizza il criterio di Delauney. Per ogni singolo elemento areale viene applicata la stessa modalità di trasformazione geometrica, che potrebbe essere lineare, come una trasformazione affine, oppure non lineare, quindi una trasformazione polinomiale di grado n. I parametri calcolati per ogni sub-area avranno quindi validità limitata a tale sottodominio, e solo qui potranno essere utilizzati.  

Un difetto di questo tipo di trasformazione è il fatto che non viene garantita la continuità spaziale dell’immagine ricampionata, con conseguenti possibili problemi di taglio: punti perimetrali di aree adiacenti, caratterizzate da parametri di trasformazione differenti, potrebbero subire spostamenti differenti, in disaccordo tra di loro. Ovviamente, per ottenere risultati migliori in termini di precisione, oltre che qualitativi, occorre adottare sub-aree più piccole possibili, vale a dire un maggiore numero di punti di controllo GCP. La tecnica del foglio di gomma è nata per scopi cartografici, e largamente impiegata nei software tecnici, ed è spesso usata per ri-georiferire a datum attuali la cartografia storica, quindi si è deciso di mutuarla al nostro caso di studio.   


L'ASPETTO APPLICATIVO.



L’incapacità di Autocad 3D Map di gestire direttamente la trasformazione geometrica dei contenuti raster, ci obbliga ad introdurre un nuovo software di casa Autodesk: Autocad Raster Design, esso non è una vera e propria applicazione a se stante, ma si presenta come un componente aggiuntivo per le altre suite. Il programma porta con se capacità di vettorializzazione, ovverosia un modulo per procedere all’estrazione al tratto di elementi vettoriali, ma non è certo questo il suo scopo principale, difatti esso fornisce la capacità di procedere alla georeferenziazione delle immagini, che possono essere salvate con opportuni file di supporto, che ne consentano l’importazione georiferita in altri software. Ma il cuore del programma è il modulo dedicato alla procedura di correlazione grazie alla tecnica del rubber-sheeting, che consente di riproiettare contenuti raster su nuovi sistemi di riferimento, come si è visto al precedete paragrafo. Il procedimento prevede la definizione di un grigliato vettoriale (Fig. 5.8) che consenta l’individuazione univoca di un certo numero di punti d’appoggio GCP; con una ottantina di punti possiamo eseguire una trasformazione polinomiale locale del sesto ordine. La posizione di questi GCP è fissata da un normale file vettoriale, che potrà essere a sua volta proiettato nel nuovo sistema di rifermento con i metodi visti.

Fig. 5.8 – Sfruttando la barra di Raster Design si importa l’immagine [1] da georiferire, la quale deve essere coerentemente posizionata [2]-[3]. Si notino i due grigliati d’appoggio, quello riferito al fuso Ovest nella sua posizione naturale, mentre il grigliato afferente al fuso Est è in una posizione del tutto casuale. Operando in ambito non lineare, le posizioni reciproche dei due grigliati d’appoggio non sono ininfluenti ai termini del risultato finale, ma non ne inficiano la validità. In realtà la soluzione trovata è una possibile tra le infinte altre, spetterà all’operatore controllarne la qualità e la utilizzabilità

Fig. 5.9 – Rubber-Sheeting: il percorso operativo.

Una volta completato il processo di georeferenziazione, l’immagine può essere importata come nuovo livello cartografico in ArcMap: la sua posizione è fissata grazie alle informazioni contenute nel file ASCII (.tfw) di supporto, che rappresentano i sei parametri di una generica trasformazione affine (Fig. 5.10).

Fig. 5.10 – Il file .tfw di supporto al dato raster.

Una volta costituito il nuovo geodatabase, si può procedere ad una sua prima valutazione, anche qualitativa, facilitata dalla presenza in ArcMap di funzioni che consentono di rendere trasparente il colore di sfondo delle immagini, ma anche queste, possono essere visualizzate in semitrasparenza (Fig. 5.11).

Fig. 5.11 – Delta del Po, sovrapposizione tra l’ortofoto del 2008 e la cartografia storica (1924).

È doveroso distingue fra l’errore portato dalla trasformazione geometrica, e quello che è l’errore per i limiti della restituzione. L’errore da trattamento geometrico deve essere tale da garantire il mantenimento della qualità, intesa come precisione del prodotto. Ciò consiste, in termini ragionevoli, imporre che lo scarto quadratico medio massimo della trasformazione sia inferire a 0,5 GSD, considerata come la minor dimensione identificabile da parte del restitutore.


L’ANALISI DEI RESIDUI.



Come scritto, la valutazione della qualità di una trasformazione avviene normalmente attraverso l’analisi dei residui, in uscita dal calcolo di compensazione ai minimi quadrati. Nel nostro caso particolare, dopo aver fissato tutti i GCP, il programma ci fornisce una tabella dei risultati (passo 5 della Fig. 5.9), ove per ogni punto d’appoggio è riportato il suo errore di posizionamento, che si traduce analiticamente in una pitagorica.

Scomponendo l’errore lungo le due direzioni principali, e considerando la totalità dei punti, si perviene alle nota relazione che esprime lo scarto quadratico medio totale RMS (senza la correzione di Bessel).

Si faccia attenzione che il valor medio, in questo caso, è stato sostituito dal valor della posizione che teoricamente il punto trasformato dovrebbe occupare (ξ(i),η(i)), dato che questo è noto in modo aprioristico, essendo fornito dalla trasformazione geometrica del reticolo d’appoggio. Per una stima sufficientemente affidabile della trasformazione è stato disposto un numero abbastanza elevato di punti (81), valore che si è cercato con dei primi tentativi pilota. Ora, al di là dei singoli valori si vuole, in questo paragrafo, mostrare alcune tecniche di rappresentazione dei residui, in modo tale da consentirne l’immediata e corretta interpretazione. Ma prima di questo risulta utile riassumere nella tabella seguente i valori notevoli ricavati dall’elaborazione.

Grazie ad un programma di riconoscimento dei caratteri (ABBYY Fine Reader OCR) è stato estratto il listato dei valori dei residui, questo ci consente di operare alcune considerazioni sulla loro rappresentazione. Il tipo più semplice di visualizzazione è un comunissimo grafico a barre (Fig. 5.12), ove gli RMSp sono riportati secondo l’ordine della loro estrazione in Raster Design, in tutti i casi dall’angolo in alto a sinistra a quello in basso, a destra. Tale rappresentazione consente di individuare il valore massimo, e di confrontare le trasformazioni tra loro.

La rappresentazione appena ottenuta può essere utile anche in sede di controllo, dato che ci consente di individuare quei GCP caratterizzati da un elevato residuo. Ad ogni modo questa rappresentazione risente di un limite fondamentale dovuto alla sua monodimensionalità, cioè dall’incapacità di farci comprendere la distribuzione dei residui sullo spazio bidimensionale del piano cartografico. Ma se consideriamo i singoli RMSp come le ipotetiche altezze discrete di una superficie continua, potremmo ricostruirci la terza dimensione, dalla cui analisi ci potremmo pure ricavare le curve di livello di egual residuo, e sovrapporle all’immagine, dandoci la facoltà di comprendere quali zone sono maggiormente esposte all’errore di trasformazione geometrica. Il grande vantaggio di una simile rappresentazione è il fatto di estendere l’errore su tutto il piano dell’immagine; è ovviamente il risultato di una interpolazione, e come tale deve essere considerata, quindi per sua natura potremmo avere massimi e minimi “rilassati” dal processo di adattamento statistico. Altro aspetto alquanto antipatico è la possibile comparsa di valori negativi, del tutto privi di significato analitico, questo succede soprattutto nelle zone di bordo, ed in assenza di punti d’appoggio: tali valori semplicemente non dovranno essere considerati, essendo essi il frutto di una lacunosa informazione d’appoggio.

Le tecniche che si utilizzano non sono né più, né meno, quelle utilizzate per l’estrazione dei modelli di elevazione (DEM). Le tecniche vanno sotto il nome di gridding methods, e tutte hanno il fine di produrre un grigliato di elevazione regolare, partendo dalla conoscenza di un certo numero di nodi d’appoggio, anche distribuito in modo irregolare. Nei software di analisi statistica spesso sono implementati molti metodi di gridding, e certamente non tutti possono fornire le risposte cercate, dato che essi, spesso, sono pensati per agire in determinati ambiti. In termini più generali i gridding methods si possono classificare in due grandi categorie.

• Exact Interpolators. 
• Smoothing Interpolators. 

Nella prima categoria troviamo tutti quegli algoritmi che producono una superficie la quale contiene i nodi d’appoggio, spesso i risultati appaiono “spigolosi”, con una grande variabilità locale. Gli smoothing interpolators consentono di formare una superficie con una certa rigidità di curvatura, caratteristica che può essere diversamente fissata in ragione dei diversi fattori di smoothing adottati. Nel caso in esame si è utilizzato ancora un volta Surfer, programma di analisi che mette a disposizione diversi metodi di interpolazione spaziale: una prima scrematura dei vari metodi è avvenuta analizzando visivamente il risultato ottenuto, quindi si è trattato di operare una mera analisi qualitativa dell’output. I parametri di confronto oggettivo sono contenuti nel report che il programma fornisce alla fine dell’operazione di interpolazione, in tal caso si è posta particolare attenzione ai valori estremi dei nodi d’appoggio, raffrontati con quelli che sono i valori omologhi della superficie interpolante. Dopo vari tentativi si è convenuto che il metodo migliore per il nostro caso di studio è il Radial Basis Function del tipo Multiquadratic. Questo è un algoritmo exact interpolator, al quale però si possono aggiungere dei parametri per lo smoothing, e quindi assumere caratteristiche proprie degli algoritmi smoothing interpolators: 

Ove h è la distanza tra il punto ed il nodo, ed R^2 è il fattore di smoothing, il cui valore può essere impostato dall’utente, o assegnato direttamene dal software dopo aver analizzato il passo griglia, perciò B(h) rappresenta il peso da assegnare ai punti diversi dai nodi d’appoggio. I risultati sono riportati nelle pagine seguenti.

Fig. 5.13 – Impostazioni usate in Surfer 10 per l’implementazione dell’algoritmo interpolante Radial Basis Function, il fattore di smoothing R^2 è stato determinato dal programma.