Monday 29 November 2010

Adding a GUI to a Python script

We can quite easily add a GUI to a Python script by using the PyQt module (http://wiki.python.org/moin/PyQt ).

As an example, I use a command-line script described in a previous post (in Italian): La pendenza topografica lungo direzioni variabili (http://gisoftw.blogspot.com/2010/06/la-pendenza-topografica-lungo-direzioni.html). This script calculates the directional slope of a DEM along orientations defined in an accompanying grid. These orientations typically vary in space and can represent for instance wind directions or glacial flows.

A more recent version of this script has also a GUI, obtained by changing just a few tens of code lines and without the need of Qt Designer (http://qt.nokia.com). In this last version, the GUI manages the definition of the input files and calculates the output slope grids according to two methods: a central differences method and Horn (1981) method. One can look at http://www.malg.eu/opensoftware/dem/directionalslope.php for additional details, download and a case study.

A valuable book for learning GUI scripting in Python is Rapid GUI Programming with Python and Qt, by Mark Summerfield (http://www.qtrac.eu/pyqtbook.html ). As a development environment, Python(x,y) (http://code.google.com/p/pythonxy) has all the required dependencies, so that one can start coding GUI without additional installations.

Saturday 27 November 2010

'Copiapó', i.e. QuantumGIS 1.6.0

The developers of QuantumGIS have just announced the release of QuantumGIS 1.6.0, named 'Copiapó. The previous one was released less than 4 months ago.
For more info see:
 http://blog.qgis.org/node/146

Friday 12 November 2010

Conversions from points to lines: which methods?

(english version of post Conversione da punti a linee: come?”, published 2010-11-09)


One of the first concepts presented in GIS courses is that vector data have point, line and polygon geometries. And we can switch between them without significant problems. However, such changes can be problematic in practice.


The actual case? I have a polygon layer of the geology in an area of the Nera Valley (Umbria, east of Foligno). I want to extract faults into a new line layer, using the contacts that already exist between the geological formation and geological formation (I ignore the faults developed entirely inside a formation). By converting polygon to line geometries, I face two main drawbacks. Line editing to delete unwanted stratigraphic limits is error prone. The latter is that line orientations are often inconsistent along a same fault.


In practice, it is more convenient to convert the lines into points, clear all unnecessary points, assign a common identifier to all the points representing a single fault, and then convert the geometry from point to line, based on fault identifiers.





While point editing is just boring, the last step can be tricky. I try with some software tools, and the result is not the expected one. There is an abundance of redundant lines, connecting points that should not be connected. It may depend on the inconsequential order of the original lines and also the presence of duplicated points generated by the polygon-to-line conversion.





Good opportunity for programming with Python. I try to create a tool ensuring proper processing ofpoint ? line conversion. Apart from removing duplicated points, we need a successful logic of conversion point ? line. It could be quite straightforward, considering that fault traces have an almost straight alignment. Therefore, a tentative algorithm for a single fault could consist in finding the two extremes of the fault trace, then, starting from an extreme, connect to each point the closest, unused one.


Can a straightforward strategy as the one depicted convert effectively points to elementary linear geometries? How to find quickly, from the execution time point of view, the most distant points, or vice versa the closest ones?


The second question can be approached via a distance matrix, which stores the distances between each point couple. So the largest distance can be easily picked, or alternatively the shortest distance between a point and a set of other points. Assuming that the largest distance is due to the two fault ends, we can find them. From one of the ends, we can then select the shortest distance with another point and connect the two points. Then, we mark the starting point as used (i.e., set all its distances to infinity) and re-apply the same logic until we use all the points. For this processing, the numpy Python module is suitable.


The answer to the first question is positive, at least in this case. This simple strategy allows an accurate conversion to line geometry. Obviously for more complex and general cases, a more sophisticated strategy would be needed.






The Python script and the used data are available in this zipped file. Input and output data format is shapefile. Required Python modules are numpy and ogr: I used Python (x, y), which already incorporates these modules. The FWTools software has the older Numeric module instead of the recent numpy. Therefore, the script should not work in FWTools, or you should modify from numpy to Numeric terminology. Other Python distributions with numpy ogr and could work without problems. The program run from the shell with the following command line:


python pt2line_05.py par.txt


where pt2line_05.py is the script name and par.txt is the text file that contains the parameters (name of the input point shapefile, name of the output line shapefile and name of the field that contains the numeric identifier of each line). An example is: 
flt_pts.shp
flt_lines.shp
flt_id

Tuesday 9 November 2010

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











Monday 8 November 2010

Due eventi "open gis" a breve in Italia

Sono:
GFOSS DAY 2010 - Foligno (Pg) - 18 e 19 Novembre 2010: http://www.gfoss.it/drupal/gfossday2010
PgDay 2010 - Roma - 10 Dicembre 2010: http://2010.pgday.it/.

Le registrazioni al secondo evento, dedicato a PostgreSQL, si sono aperte oggi, 8 Novembre.

Wednesday 3 November 2010

Antartide e GIS

E' un continente quasi sempre lontano dai riflettori della cronaca. Viene citato talvolta quando qualche grande iceberg si distacca dalle sue piattaforme glaciali, evento magari presentato come "diretta conseguenza" del riscaldamento globale. La scienza invece ha gli occhi puntati su questo continente da più di un secolo. La ricerca  ha bisogno di un sistema di riferimento cartografico solido, oltre che dati per studiarlo, per esempio immagini satellitari e dati cartografici.
Quali sono i sistemi cartografici usati e quali i dati di base disponibili?
Tra i ricercatori il sistema più usato, soprattutto se si studia l'intero continente, è il Polare Stereografico. Il parallelo standard è il -71°, ovvero 71°S. Il datum, come al solito, è il WGS84. Ovviamente questo stesso sistema può essere usato specularmente per le zone polari settentrionali.
Questo sistema è quello usato per la più comprensiva raccolta di dati cartografici a scala dell'intera Antartide,  l'Antartic Digital Database: http://www.add.scar.org:8080/add/. Questo database cartografico è open, salvo registrazione preliminare. E' stato supervisionato dal British Antarctic Survey (BAS) e raccoglie, suddivisi in una serie di tiles (vista la vasta superficie del continente), una serie di livelli che vanno dalle curve di livello, alle colonie di animali (pinguini, petrelle e simili), alle linee di flusso glaciali. E' di qualità buona, anche se localmente possono essere presenti alcune incongruenze o i dati possono non essere aggiornati, per esempio dove le lingue glaciali si sono modificate recentemente. Le fonti dei dati sono infatti i prodotti cartografici di svariate nazioni e su lunghi periodi di tempo, il che giustifica le variazioni di qualità ed aggiornamento nonostante il tentativo di uniformazione praticato sotto la guida del BAS.
Accoppiabili a questi dati vi sono i mosaici satellitari a scala continentale: due classici sono il MOA (MODIS Mosaic of Antartica): http://nsidc.org/data/moa/index.html,  ed il RAMP Radarsat: http://nsidc.org/data/nsidc-0103.html. Entrambi sono monobanda a scala di grigi, a 125 m di risoluzione nominale e già georiferiti in Polare Stereografica. Mentre il MOA deriva da sensori nell'ottico, il RAMP si basa su un sensore radar (quindi con le situazioni di spostamento apparente delle zone rilevate tipiche di questa metodologia).

Nelle due immagini sottostanti è visualizzata la lingua glaciale Drygalski, lunga circa 100 km, prima del suo impatto con l'iceberg B15A (avvenuto ad Aprile 2005, http://www.esa.int/esaCP/SEMEGLW797E_index_0.html). La prima immagine è RAMP-Radarsat, la seconda MOA.




Più recente è il Landsat Image Mosaic of Antarctica (LIMA): http://lima.usgs.gov/, che deriva da immagini Landsat con una risoluzione nominale di 30 m e viene redistribuito ricampionato a 240 m di risoluzione in formato jp2.
Oltre a questo è possibile anche scaricare le immagini Landsat. Ma con la recente liberalizzazione dei dati Landsat, forse è più semplice e ricco scegliere le immagini da scaricare dal sito glovis: http://glovis.usgs.gov/.