Calculating the intersections between planes and DEM: a Python implementation
Dato un punto nello spazio 3D, con una orientazione associata (espressa p.e. come inclinazione ed immersione, secondo la convenzione geologica), qual'è l'intersezione della superficie col DEM, ovvero l'andamento della traccia della superficie sul DEM?
Questo problema è di interesse in geologia strutturale, in quanto spesso disponiamo di misure di superfici in stazioni strutturali, misurate come orientazione di un piano.
Possiamo essere interessati a conoscere quale sia la “proiezione” del piano sul DEM, nell'assunzione che la sua orientazione rimanga immutata in un certo intorno della misura.
Se noi non conosciamo l'andamento effettivo della superficie, questa “proiezione” può essere utile per ipotizzare la posizione della traccia nelle vicinanze della misura. Se invece, grazie ad osservazioni di campagna, ne conosciamo l'andamento, possiamo verificare quanto la nostra misura planare sia rappresentativa dell'andamento complessivo della superficie e se nel variarla riusciamo a ottenere un maggiore accordo tra misura e dato di campagna.
Oltre a queste finalità immediate, ce ne possono essere altre: questa determinazione può essere il primo passo per segmentare un volume roccioso, la cui superficie superiore è espressa dal DEM, in più sottoelementi, con superfici di delimitazione rappresentate da più piani o superfici più complesse (p.e. funzioni trigonometriche).
Inoltre, la determinazione delle intersezioni per superfici teoriche permette di disporre di dataset teorici di test per algoritmi che effettuano il procedimento inverso, cioé stimano in maniera (semi)automatica l'andamento di superfici naturali a partire da osservazioni puntuali su una superficie topografica.
Modello analitico
La misura planare rappresenta il caso geometrico più semplice, quindi più facilmente implementabile in un algoritmo e concettualmente più semplice. Questo aiuta la comprensione dell'andamento reale della superficie investigata.
Determinare l'intersezione tra DEM e superficie non è immediato, in quanto di tratta di un problema con dimensionalità superiore a 2, soglia oltre la quale i software GIS attuali sono poco attrezzati.
Il modello concettuale qui usato per calcolare le intersezioni tra DEM e piano nell'algoritmo implementato in Python si basa su semplici considerazioni di geometria analitica.
L'equazione di un piano passante per un punto (x0, y0, z0) e perpendicolare ad una retta con coseni direttori (l, m, n) è (Zwirner, 1983, p. 198):
l (x – x0) + m ( y – y0 ) + n ( z – z0 ) = 0
sviluppabile come:
z = - ( l / n ) (x – x0 ) – ( m / n ) ( y – y0 ) + z0
I coseni direttori di una retta con inclinazione δ e orientazione θ sono (Groshong, 1999, pp. 67-8; vedi anche Cox e Hart, p. 167, eq. 5.1-3):
l = sin θ cos δ
m = cos θ cos δ
n = - sin δ
Considerando le relazioni tra un piano con inclinazione δP e immersione θP, e la retta ad esso perpendicolare e che punta verso l'alto, si ha:
θ = θP
δ = δP - 90°
per cui i coseni direttori predetti, espressi rispetto alla giacitura del piano, hanno valori:
l = sin θP sin δP
m = cos θP sin δP
n = cos δP
Ritornando all'equazione del piano, abbiamo quindi:
z = a ( x – x0 ) + b ( y – y0 ) + z0
con:
a = - l / n = - sin θP tan δP
b = - m / n = - cos θP tan δP
I due coefficienti a, b rappresentano le derivate parziali di z rispetto a x e a y.
Essendo θP e δP costanti nel caso del piano, anche a e b lo sono.
Queste sono le basi analitiche per esprimere un piano in base alla sua giacitura geologica. Ma qual'è la metodologia per il calcolo dell'intersezione tra DEM e superficie planare?
Una procedura semplice può essere quella di calcolare, per ogni cella (e precisamente il suo punto centrale), l'intersezione tra due linee che rappresentano l'andamento locale del DEM e della superficie, sia in direzione parallela all'asse delle x, sia in quella dell'asse delle y.
Consideriamo il caso delle intersezioni lungo sezioni parallele all'asse delle x. Per quelle parallele all'asse delle y le considerazioni sono analoghe.
L'algoritmo determina le equazioni dei segmenti che esprimono l'andamento locale del DEM in base alle coordinate ( xn, zn ) e ( xn+1, zn+1 ).
Se i segmenti relativi al DEM e al piano lungo l'asse delle x hanno coefficienti angolari ed intercette ( m1, q1 ) e ( m2, q2 ), l'intersezione ha valore:
xinters = ( q2 – q1 ) / ( m1 – m2 )
mentre il valore y è quello della particolare riga considerata. Da ( x, y ) si determina poi z.
Un valore calcolato di intersezione rappresenta una intersezione reale se il suo valore delle x ricade nell'intervallo tra il centro della cella considerata e quello della cella successiva lungo la direzione dell'asse delle x.
Implementazione
Questa prima versione è stata implementata con Python e richiede i moduli numpy e gdal. I DEM corrispondono a matrici di elevazione, il che permette di elaborarli agevolmente in Python con numpy.
L'algoritmo creato non è ottimizzato ed effettua i calcoli necessari per la determinazione dei valori delle intersezioni descritti di seguito per ogni cella del DEM, anche quando questa non presenta nessuna intersezione con la superficie planare. Ma anche in questa versione preliminare, i tempi di esecuzione sono comunque ridotti grazie alla efficiente gestione delle operazioni su array da parte di numpy.
La versione attuale dell'algoritmo è a linee di codice, senza interfaccia grafica, ed ancora in fase di testing e miglioramento. Verrà rilasciata come software open source ad uno stadio successivo.
Esempio di applicazione su una faglia normale pliocenica della Valnerina (Umbria)
Presento un esempio di applicazione di questo modulo su dati geologici della zona della Valnerina, Umbria (mappati da me nella tesi di dottorato).
La sequenza stratigrafica affiorante è quella umbro-marchigiana, con sedimenti marini che vanno dal Dogger-Malm (Calcari Diasprini) al Miocene inferiore (Bisciaro). Nel Miocene-Pliocene questi sedimenti vengono piegati e fagliati, in corrispondenza della fase compressiva che ha originato la catena appenninica. Nel tardo Pliocene-Pleistocene inizia una fase distensiva, con formazione di faglie normali ad orientazione NW-SE. E' proprio su una di queste faglie che valutiamo la corrispondenza tra misurazioni di terreno e andamento cartografico, cioé traccia della superficie con la topografia.
Consideriamo una faglia normale per la quale si dispone di una misura strutturale, con immersione 227° ed inclinazione 58°. L'andamento mappato della faglia è evidenziato in figura.
Le strie circa E-W nello shaded relief sono dovute ad artefatti tipico del DEM Aster.
Il risultato con i valori della misura non concorda con la traccia mappata della faglia. L'immersione complessiva sembra essere di circa 10° inferiore a quella nel sito.
Testando una nuova intersezione, usando una immersione di 217° e mantenendo l'inclinazione a 58°, si ottiene una traccia che ricalca maggiormente l'andamento complessivo, anche se sembra avere una variabilità locale superiore a quella mappata. Come è noto, tanto più un piano è verticale, tanto più la sua traccia sulla superficie topografica sarà rettilinea. E' possibile quindi testare un piano con la stessa immersione ma inclinazione maggiore.
Con valori di immersione di 217° ed inclinazione 70° si ottiene una traccia della superficie che rispetta in maniera soddisfacente l'andamento mappato della faglia. Discordanze locali possono essere dovute a variazioni della superficie, o anche a interpolazioni erronee della traccia stessa in zone con assenza di dati. Nella porzione sud-orientale della traccia è comunque evidente che si ha una discordanza di circa 20° nell'immersione. Si può quindi ipotizzare che in questo settore la faglia varii la sua orientazione complessiva, oppure che si tratti di un segmento indipendente dalla faglia principale. Questo suggerisce le zone in cui possono essere utili ulteriori controlli di terreno.
Cox, A., Hart, R.B., 1990. La tettocnica delle palcche. Meccanismi e modalità. Zanichelli, Bologna. 383 pp.
Groshong, R. H. jr., 1999. 3-D Structural Geology. Springer Verlag, Berlin. 324 pp.
Zwirner, G. 1983. Istituzioni di matematiche. Parte seconda. Ed. CEDAM, Padova. 486 pp.
