Conversione da punti a linee: come?

(english version: http://gisoftw.blogspot.com/2010/11/conversions-from-points-to-lines-which.html)

Uno dei primi concetti che vengono illustrati nei corsi GIS per descrivere i dati vettoriali è che esistono geometrie puntuali, lineari e poligonali. E che è possibile passare dall'una all'altra geometria. Tentare di metterlo in pratica può riservare sorprese, come si impara con l'esperienza.

Qual'è il caso specifico? Ho un livello poligonale, che rappresenta la geologia in un settore della Valnerina (Umbria, ad est di Foligno). Voglio isolare le faglie in un nuovo livello lineare, partendo dai contatti geologici che sono già presenti tra formazione e formazione geologica (trascuro le faglie che si sviluppano esclusivamente all'interno di una stessa formazione).

Convertendo il livello da poligonale a lineare mi presenta tutta una serie di problemi: i limiti possono avere significato sia di faglia sia di limite stratigrafico. I loro identificativi non sono coerenti. I versi delle linee non sono costanti lungo una stessa faglia. L'eliminazione dei limiti lineari che non interessano è macchinoso e soggetto ad errori.

In pratica, è più comodo retrocedere i dati a puntuali, cancellare tutti i punti inutili, attribuire un identificativo comune a tutti i punti che rappresentano una singola faglia, e poi trasformare tutti in lineare in base all'identificativo.

Ma non è così banale l'ultimo passaggio (quelli precedenti sono solo noiosi). Provo con alcuni strumenti di vari software, ed i risultati lasciano a desiderare. Linee ridondanti, che collegano punti che dovrebbero essere sconnessi, e simili amenità. Può dipendere dall'ordine, non sequenziale e con versi variabili, dei punti, ed anche dal fatto che purtroppo i punti si sono anche duplicati e triplicati nella conversione poligonale → lineare → puntuale.

Buona occasione per programmare con Python, tentando di creare uno strumento che risparmi ulteriori grattacapi effettuando correttamente la trasformazione punti → linee. A parte la cancellazione dei duplicati, serve una logica di conversione punto → linea, magari anche semplificata considerando che le tracce di faglie hanno uno sviluppo grosso modo rettilineo, senza particolari zig-zag interni. Può essere quella di trovare le due estremità della faglia, poi di partire da un estremo e man mano collegare ad ogni punto quello più vicino che non è ancora stato utilizzato.

Può essere sufficiente una strategia così semplice per dati di faglia?

E come trovare in maniera veloce, dal punto di vista dei tempi di esecuzione, i punti più distanti o viceversa quelli fra loro più vicini?

La seconda questione può essere risolta calcolando una matrice delle distanze, che conserva per ogni punto la distanza dai punti dell'intero set. Così si rintracciano subito i due punti più distanti fra loro, e nell'assunzione che rappresentino gli estremi della faglia, se ne sceglie uno e da questo si sceglie via via il punto più vicino, marcando come già usato il precedente (nella matrice delle distanze si settano i suoi valori in riga e in colonna ad infinito), fino ad esaurire i valori di distanza non infiniti nella matrice. Naturalmente, dopo avere scelto la distanza massima, i valori di distanza nulla di ogni punto con se stesso vengono resi inoffensivi trasformandoli in infinito. Il modulo Numpy di Python è adatto per questo genere di elaborazioni.

La risposta alla prima questione, almeno nel mio caso, è stata positiva. Una strategia così semplice permette di ottenere una geometria lineare corretta che rispetta gli andamenti delle faglie. Ovviamente su dati più complessi occorrerebbero strategie più sofisticate, per esempio per trattare i casi in cui non si selezionano correttamente gli estremi.

Lo script Python e i dati reali usati sono disponibili in questo file zippato. I moduli richiesti sono numpy e ogr: io ho effettuato le elaborazioni in Python(x,y), che incorpora già questi moduli. In FWTools numpy non è presente, ma solo il suo predecessore Numeric, e quindi questo modulo non funzionarebbe a meno di convertire la sintassi numpy in numeric. Altre distribuzioni osgeo con Python, Numpy e ogr dovrebbero invece andare bene.

Il programma si lancia dalla shell con la linea di comando:

python pt2line_05.py par.txt

dove pt2line_05.py è il nome dello script e par.txt è quello del file testuale che contiene i parametri (nome dello shapefile puntuale di input, dello shapefile lineare di output, del campo che contiene l'identificativo numerico di ogni linea), per esempio: 

flt_pts.shp

flt_lines.shp

flt_id