domenica 28 dicembre 2008

Trovare i contorni duplicati con GRASS

Obbiettivi

Con questa esercitazione si vuole dimostrare come GRASS  possa individuare quali e quanti sono, in una copertura poligonale non topologica, i contorni (boundaries) duplicati.

Si farà ricorso a moduli di alto livello, senza dover accedere ai vettoriali a basso livello, come nel caso del'impiego di v.build  (alto e basso livello vengono qui usati per riferirsi al livello di interfacciabilità verso l'utente).

Set di dati per l'esercitazione

La copertura poligonale scelta per l'esercitazione è polyNoTopology, contenuta nella location simpleloc.

Questa location è stata creata appositamente per fini didattici; per avere maggiori notizie vedere questo  precedente blog.

La location simpleloc può essere scaricata da qui.

polyNoTopology è un vettoriale non topologico (modello CAD o modello spaghetti) con le seguenti caratteristiche:

  • i contorni dei poligoni adiacenti sono duplicati;
  • i contorni, se si intersecano, non generano nodi.

Prime verifiche

Avviare GRASS  ed entrare nella location simpleloc; si dovrebbero trovare i seguenti vettoriali:

GRASS 6.4.svn (simpleloc):~ > g.list vect
----------------------------------------------
vector files available in mapset <user1>:
node_grid polyBreakCat polyNoTopology punti xylabels
polyBoundCat polyClean polyTopology vgrid
polyBreak polyFromArc poly_arc vrastergrid

Si può verificare che vi sono lati duplicati controllando le coordinate dei vertici dei contorni semplicemente con v.out.ascii:

v.out.ascii in=polyNoTopology \
format=standard | sed -n '/^B/,/^C/p' | sed -n '/^C/!p'

Leggendo le coordinate in output si può verificare che i lati in comune sono digitati due volte. Il pipe ('|') verso sed non è strettamente necessario e serve solo per filtrare i dati non utili in questo contesto.

Cercare i contorni duplicati

I contorni di polyNoTopology, oltre ad essere duplicati in corrispondenza delle adiacenze tra poligoni, non sono spezzati (break).

La procedura per trovare i duplicati può essere riassunta nei seguenti passi:

  1. spezzare in boundaries (= contorni dei poligoni), inserendo un nodo nei punti di incrocio;
  2. numerare in modo univoco i boundaries;
  3. rimuovere i boundaries duplicati;
  4. attribuire nuovamente una numerazione univoca ai boundaries senza duplicati
  5. interrogare, dallo strato dei boundaries senza duplicati, le geometrie dello strato con i duplicati: per ciascun elemento geometrico del primo ottenere quanti e quali elementi del secondo esistono: se si ottiene un valore pari ad 1 il boundaries in quel tratto non è duplicato; se si ottiene un valore maggiore di 1 il boundaries è duplicato.

Per spezzare i boundaries si usa il modulo  v.clean con il tool break:

v.clean in= vectInput out=polyBreak tool=break type=boundary --overwrite

Si può facilmente verificare che i boundaries non hanno categoria; nell'output di v.category opt=report, infatti, vengono riportate le categorie per i soli centroidi:

v.category in=polyBreak opt=report
Layer/table: 1/polyBreak
type count min max
point 0 0 0
line 0 0 0
boundary 0 0 0
centroid 6 1 6
area 0 0 0
all 6 1 6

Si aggiunge allora la categoria ai boundaries spezzati (break) creando il layer 2 ; considerato che i centroidi hanno categoria da 1 a 6 , per non confondere i valori, si decide di numerare i boundaries da 100 in poi:

v.category in=polyBreak type=boundary out=polyBreakCat \
layer=2 opt=add cat=100 --overwrite

La situazione dei layers e delle categorie può essere evidenziata ancora una volta con v.category opt=report:

v.category in=polyBreakCat opt=report
Layer: 2
type count min max
point 0 0 0
line 0 0 0
boundary 18 100 117
centroid 0 0 0
area 0 0 0
all 18 100 117
Layer/table: 1/polyBreakCat
type count min max
point 0 0 0
line 0 0 0
boundary 0 0 0
centroid 6 1 6
area 0 0 0
all 6 1 6

La situazione è quindi la seguente: sul layer 1 di polyBreakCat vi sono i centroidi con categorie da 1 a 6 ; sul layer 2 vi sono i buondaries numerati, in modo univoco, da 100 a 117 .

Occorre ora rimuovere i duplicati da polyBreakCat:

v.clean in=polyBreakCat out=polyClean \
type=boundary tool=rmdupl --overwrite

GRASS  ha rimosso i pezzi di boundaries duplicati; il vettoriale ottenuto (polyClean) conserva sul layer 1 le categorie dei centroidi e sul layer 2 le categorie di tutti i boundaries, compresi quelli duplicati; ed è proprio quest'ultima caratterisitca il segreto di questo eserzio.

Se infatti ora si rinumerano i boundaries, verrano rinumerati solo ed esclusivamente i boundaries senza i duplicati ma esisterà una relazione uno a molti tra il layer 2 e il layer 3 .

v.category in=polyClean out=polyBoundCat \
type=boundary layer=3 opt=add cat=100 --overwrite
v.category in=polyBoundCat opt=report
Layer/table: 2/polyBoundCat_2
type count min max
point 0 0 0
line 0 0 0
boundary 18 100 117
centroid 0 0 0
area 0 0 0
all 18 100 117
Layer/table: 3/polyBoundCat_3
type count min max
point 0 0 0
line 0 0 0
boundary 13 100 112
centroid 0 0 0
area 0 0 0
all 13 100 112
Layer/table: 1/polyBoundCat
type count min max
point 0 0 0
line 0 0 0
boundary 0 0 0
centroid 6 1 6
area 0 0 0
all 6 1 6

Riassumendo:

  • sul layer 1 : i centroidi (categorie da 1 a 6 );
  • sul layer 2 : i boundaries spezzati (break) ma con i duplicati (categorie da 100 a 117 );
  • sul layer 3 : i boundaries spezzati (break) e senza dupicati (categorie da 100 a 112 ).

Si aggiungono ora le tabelle al layer 2 e 3 di polyBoundCat:

v.db.addtable polyBoundCat layer=2

v.db.addtable polyBoundCat layer=3 \
col="cat integer, count integer,listarcs varchar(12)"

GRASS  offre un potente modulo per interrogare le geometrie e caricare i dati dell'interrogazione nella tabella degli attributi: v.to.db.

Nel modello vettoriale di GRASS, per ogni primitiva geometrica possono venir memorizzati i numeri dei layers e, per ciascun numero di layer, il numero della categoria (difficilmente si capirà il meccanismo di questo esercizio se non si ha chiaro questo concetto; vedere il manuale di GRASS).

Viene ora il momento di interrogare le geometrie del layer 2 , popolando i due campi count e listarcs:

v.to.db map=polyBoundCat option=query type=boundary \
layer=3 qlayer=2 column=count qcolumn="count(*)"
v.to.db map=polyBoundCat option=query type=boundary \
layer=3 qlayer=2 column=listarcs qcolumn="group_concat(cat)"

Nella seconda operazione si usa il comando di aggregazione group_concat supportato nelle recenti versioni di SQLite: questo comando concatena i valori trovati, creando così l'elenco delle categorie dei boundaries.

Si passa ora alla visualizzazione della tabella:

v.db.select polyBoundCat layer=3
cat|count|listarcs
100|1|100
101|1|105
102|1|106
103|1|107
104|1|111
105|1|112
106|1|115
107|1|117
108|2|101,103
109|2|102,104
110|2|108,109
111|2|110,113
112|2|114,116

GRASS  ha correttamente indviduato quali e quanti boundaries erano duplicati:

  • nella colonna count, per ogni record (=per ogni geometria) del layer 3 sono stati contati quanti elementi sono presenti nel layer 2 ;
  • nella colonna listarcs, per ogni record (=per ogni geometria) del layer 3 sono stati elencati quali elementi sono presenti nel layer 2 ; in questo caso sono state elencate le categorie corrispondenti ai boundaries contenuti nel layer 2 .

Una visualizzazione grafica renderà più chiare le idee; si può usare la linea di comando:

# # # disegno
./sl_draw.sh
d.vect polyBoundCat layer=3 display=shape,cat color=grey llayer=3 \
xref=center yref=center lcolor=aqua \
lsize=10 bgcolor=white bcolor=black
# # # output del file png
d.out.file out=polyBoundCat_l3 format=png --overwrite
./sl_draw.sh
d.vect polyBoundCat layer=2 display=shape,cat color=grey llayer=2 \
xref=center yref=center lcolor=red \
lsize=10 bgcolor=white bcolor=black
# # # output del file png
d.out.file out=polyBoundCat_l2 format=png --overwrite

Lo script sl_draw.sh può essere scaricato da questo indirizzo; deve essere ovviamente scomapttato, reso eseguibile (chmod ugoa+x sl_draw.sh) e posto in una directory in modo che l'interprete dei comandi possa trovarlo (controllare con echo $PATH).

Tutto il codice utilizzato in questo esercizio può esere scaricato da qui

Il risultato delle ultime righe di comando dovrebbe essere la visualizzazione dei due layer di polyBoundCat in due separati monitors:

È pure possibile usare la comoda interfaccia grafica wxpython di GRASS:

Il caso dell'import di shapefiles

Tutta la procedura vista può essere utilmente applicata al caso dell'import di shapefiles in GRASS.

Lo shapefile è una struttura dati di tipo non topologico. Si provi allora ad esportare in shapefile la copertura poligonale (topologicamente corretta) della location simpleloc denominata polyTopology con v.out.ogr:

v.out.ogr -c input=polyTopology type=area dsn=~/tmp

L'opzione -c fa in modo che nell'export siano comprese solo le features geometriche con categorie (nel caso di polyTopology tutte le aree). Si può ora reimportare in GRASS  lo shapefile:

v.in ogr dsn=~/tmp layer=polyTopology out=polyFromShp snap=-1 --o

L'esercizio visto può essere ripetuto con il vettoriale ottenuto dall'import (polyFromShp) addivenendo ai medesimi risultati.

Il modulo v.clean offre dei tools appositi per la pulizia dei vettoriali importati in GRASS  da strutture dati (come gli shapefiles) non topologici: rmdupl e bpol; il primo rimuove i duplicati; il secondo spezza (=break) in modo topologico il contorno dei poligoni.

Conclusioni

Questa esercitazione permette di capire quali attenzioni devono essere poste quando si importa in un GIS  topologico come GRASS  una struttura dati non topologica.

Il modulo v.clean rappresenta il passaggio obbligato attraverso il quale pulire i vettoriali importati in GRASS: si è vista la possibilità di rimozione dei duplicati (rmdupl), l'inserimento di nodi nei punti di intersezione (break); si è pure accennato all'apposito tool bpol (per i polgioni).

Il modulo v.to.db consente di interrogare le geometrie tra layers diversi dello stesso vettoriale con l'opzione query.

Sono stati trattati solo due dei diversi tipi di errori di topologia. Per altre anomalie, quali pseudonodi, dangles, weird polygons, sliver polygons, GRASS  offre dei rimedi più o meno efficaci, che si esamineranno in dettaglio in prossime esercitazioni.



sabato 13 dicembre 2008

Piani Forestali e wget

Visto che per me è ormai diventata ormai una fissazione, mi appunto il comando per scaricare, tutti insieme e in un colpo solo, i Piani Forestali della Regione Piemonte (tutti i comuni e tutte le aree forestali) :



wget -r -nd -nH --level=3 --accept=zip \

http://www.regione.piemonte.it/montagna/foreste/pianifor/forestale.htm



Ricordarsi di registrarsi... prima.

ogr2ogr, codici epsg e Piani Forestali del Piemonte

In un precedente post ho usato ogr2ogr per convertire gli shapefiles dei Piani Forestali della Regione Piemonte da ED50 UTM 32N a WGS84 UTM 32N.
Quando si parla di conversione tra datum diversi divento diffidente: forse a causa della scarsa conoscenza della materia e delle ripetute fregature avute in passato...e nel presente!
In ogni caso voglio semplicemente provare ad usare ogr2ogr con diversi settaggi, tanto per chiarirmi le idee.
La sintassi (semplificata) di ogr2ogr per la conversione tra sistemi di riferimento diversi (di shapefiles) è la seguente:
ogr2ogr -s_srs (descrizione del sistema di riferimento dell'input)
-t_srs (descrizione del sistema di riferimento dell'output)
(nome shapefile di output) (nome shapefile di input)

Per quanto riguarda la descrizione del sistema di riferimento dell'imput (nel mio caso: ED50 UTM 32N) le possibilità che ho preso in considerazione sono:

  1. una variabile che immagazzina una lunga stringa wkt, come quella riportata qui:

    wkt_ed50='PROJCS["UTM Zone 32, Northern Hemisphere", GEOGCS["international",
    DATUM["European_Datum_1950",
    SPHEROID["International_1924",6378388,297],
    TOWGS84[-87.000,-98.000,-121.000]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],UNIT["meter",1]]'

    Per ottenere questo wkt ho semplicemente creato una Location di GRASSS in ED50 UTM32N e da shell ho digitato g.proj -w. Lo stesso risultato potrebbe essere ottenuto da qui (fare attenzione ai parametri di traslazione TOWGS84...);
  2. il codice epsg:
    '+init=epsg:23032'
  3. il codice epsg (23032) e i parametri di traslazione:
    '+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000'
  4. le stesse informazioni di cui sopra ma in formato proj:
    '+proj=utm +zone=32 +ellps=intl +units=m +towgs84=-87.000,-98.000,-121.000'
Per quanto riguarda l'output in WGS84 UTM 32N ho semplicemente settato al codice epsg 32632.
Ho quindi scritto questo semplice script per testare il tutto su uno shapefile dei Piani Forestali della Regione Piemonte:

#!/bin/sh
dir_work="$HOME/pft"
shpfile="004017"
suffix="_wkt"
ogr2ogr -s_srs "$wkt_ed50" \
    -t_srs '+init=epsg:32632' \
    "$dir_work/${shpfile}${suffix}.shp" \
    "$dir_work/${shpfile}.SHP" suffix="_epsg"

ogr2ogr -s_srs '+init=epsg:23032'\
    -t_srs '+init=epsg:32632' \
    "$dir_work/${shpfile}${suffix}.shp" \
    "$dir_work/${shpfile}.SHP" suffix="_epsgtowg"

ogr2ogr -s_srs '+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000' \
    -t_srs '+init=epsg:32632' \
    "$dir_work/${shpfile}${suffix}.shp" \
    "$dir_work/${shpfile}.SHP" suffix="_defepsg"

ogr2ogr -s_srs '+proj=utm +zone=32 +ellps=intl +units=m +towgs84=-87.000,-98.000,-121.000' \
    -t_srs '+init=epsg:32632' \
    "$dir_work/${shpfile}${suffix}.shp" \
    "$dir_work/${shpfile}.SHP" suffix="traspunto"

wine /usr/local/traspunto/traspunto.exe shp pia pia \
    "${dir_work}/${shpfile}" \
    "${dir_work}/${shpfile}${suffix}" ED50 WGS84 32 32

dove nella variabile $wkt_ed50 è contenuto il file wkt visto sopra.

Da notare che uso anche traspunto, disponibile qui, con wine.

Questa è una prova (empirica) nella quale apprezzo visivamente i risultati; non ho qui l'interesse di fare alcuna analisi statistica del tipo: esiste una differenza statistica (e se si, quale?)nel posizionamento spaziale di una serie di punti trasformati con traspunto rispetto a ogr2ogr?

Una volta eseguito lo script visualizzo il tutto in OpenJUMP (ma, ovviamente, vanno bene anche Quantum Gis, gvSIG, udig, ecc).

Il primo confronto interessante è quello tra traspunto (linee in rosso) e ogr2ogr con '+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000' (linee in verde).

Per apprezzare la differenza, nella figura ho settato una griglia di sfondo con una maglia di 10 metri.

Il secondo confronto riguarda i due settagggi di ogr2ogr: il primo con il solo codice epsg (in blu) (-s_srs '+init=epsg:23032'); il secondo (in verde) con il codice e i parametri di traslazione ('+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000'). Ho mantenuto la griglia a 10 m.




Com'era prevedibile non esistono differenze tra il risultato della conversione effettuata con il wkt, con '+init=epsg:23032 +towgs84=-87.000,-98.000,-121.000' e con
-s_srs '+proj=utm +zone=32 +ellps=intl +units=m +towgs84=-87.000,-98.000,-121.000'.

Conclusioni

Una prova rigorosa su modalità diverse di conversione può solo essere condotta con una serie di punti dei quali si conoscono a priori con esattezza le coordinate nei diversi sistemi di interesse.
Questa è una semplice prova per fissare alcune cose:
  • non fermarsi al codice epsg: cercare, se vi sono, i parametri di rototraslazione relativi alla zona di interesse;
  • nella definizione di una nuova Location porsi il problema di settare i parametri di rototraslazione (vedere g.proj);
  • l'uso di traspunto può essere un'alternativa;
  • l'uso del software Verto unitamente ai grigliati (in vendita presso l'Istituto Geografico Militare ) potrebbe costituire la migliore tra le alternative (?) ...

giovedì 11 dicembre 2008

Piani Forestali, ogr2ogr e conversione di coordinate

Da qualche tempo lavoro con i Piani Forestali della Regione Piemonte che possono essere scaricati da qui previa registrazione.
È possibile scaricare per singoli Comuni oppure per aree forestali.
Lavorare con un gran numero di shapefiles è abbastanza semplice, quando si capisce la possibilità di iterare all'interno di una directory contenente gli shapefiles medesimi.
Ecco i primi passaggi che faccio nell'ipotesi di aver scaricato un'area forestale.

L'iterazione ha sempre la stessa struttura base: for files in directory; do qualcosa; done

Ad esempio:

laptop@pcgis:~/pft$ for f in *.*; do echo $f; done

004002.ZIP

004006.ZIP

004025.ZIP

004079.ZIP

004093.ZIP

... etc ...



Appena scaricati, sul mio pc Linux-Ubuntu i files zippati non hanno permessi di scrittura-lettura


laptop@pcgis:~/pft$ for f in *.*; do echo $f; chmod ugoa+wr $f; done

004002.ZIP

004006.ZIP

004025.ZIP

004079.ZIP

004093.ZIP

... etc ...



A questo punto scompatto gli zip:

laptop@pcgis:~/pft$ for f in *.*; do echo $f;unzip $f; done004002.ZIP

Archive: 004002.ZIP

inflating: 004002.TXT

inflating: 004002.DBF

inflating: 004002.SHX

inflating: 004002.SHP

004006.ZIP

Archive: 004006.ZIP

inflating: 004006.TXT

inflating: 004006.DBF

inflating: 004006.SHX

inflating: 004006.SHP

004025.ZIP

Archive: 004025.ZIP

inflating: 004025.TXT

inflating: 004025.DBF

inflating: 004025.SHX

inflating: 004025.SHP

... etc ...




Gli shapefiles dei PFT della Regione Piemonte sono in ED50 UTM32 Nord che corrisponde ad un codice epsg pari a 23032 (parametri per la rototraslazione in WGS84: TOWGS84=-87.000,-98.000,-121.000); siccome la direttiva che ho ricevuto è quella di lavorare in WGS 84 UTM 32 Nord, devo fare una conversione tra datum

laptop@pcgis:~/pft$ suffix="_wgs84"



laptop@pcgis:~/pft$ for f in *.SHP; do name=${f%.*}; ogr2ogr \

-s_srs '+proj=utm +zone=32 +ellps=intl +units=m +towgs84=-87.000,-98.000,-121.000' \

-t_srs '+init=epsg:32632' "${name}${suffix}.shp" $f; done



Con poche righe di comando ho ottenuto nella cartella pft una lunga serie di shapefiles nominati con il vecchio nome dello shapefile in ED50 UTM 32N più un suffisso _wgs84. Tutti gli shapefiles generati da ogr2ogr hanno il file prj che descrive il sistema di riferimento (WGS84 UTM 32N)

In figura si può vedere il medesimo shapefile riferito al sistema ED50 UTM 32 N (in marrone) e il risultato della traslazione effettuata con ogr2ogr (in viola).



From Grass, gis e dintorni


Ovvio che il tutto può essere inglobato in uno script, per migliorarne la leggibilità e la riusabilità

Per saperne di più su ogr2ogr si può consultare il suo manuale qui.

I comandi impartiti richiedono una spiegazione che sarà oggetto di un prossimo post.

Smanettare in GRASS con simpleloc

Confesso di essere uno smanettone di GRASS. È un colpa?
Da quando ho scoperto la shell di GRASS non ho più potuto farne a meno. L'iniziazione è stata opera dell'ottimo libro Open Source Gis. A GRASS Gis Approach di M. Neteler e H. Mitasova (vedere la pagina web ); a darmi il colpo di grazia è stato simpleloc.
Quando si smanetta con GRASS la cosa più noiosa è la visualizzazione dei risultati. Scartata a priori l'ipotesi di usare la nuova e ottima interfaccia grafica wxPython e quella (vecchia ma funzionale) tck/tk, non resta che d.mon.

Si perde tuttavia del tempo ad aprire un monitor, creare un minimo di base per la visualizzazione, e disegnare ciò che si vuole.
Mi sono allora proposto di creare un semplice script che faccia tutto questo, senza dovermi preoccupare se ho ancora un monitor da aprire (nel caso siano già tutti aperti).
Condivido volentieri questo script: sl_draw.zip. Da notare: è frutto di un principiante, saprete senz'altro trovare codice ridondante e miglioramenti: fatemeli sapere. Grazie.
Scaricare, scompattare e rendere eseguibile lo script, nella medesima directory ove è già stato scaricato lo script per la creazione dei dati di simpleloc e il file di dati per la creazione dei vettoriali (entrambi i file possono essere scaricati da qui)

Per creare simpleloc vedere qui.
Con due semplici linee di comando è possibile creare una visualizzazione decente...

mercoledì 10 dicembre 2008

simpleloc

Premessa: datasets e principianti

Esistono diversi set di dati per fare esercizi con GRASSe per sperimentarne i moduli. I più usati sono senz'altro questi due:

A fianco di numerosi e innegabili vantaggi, questi datasets presentano alcuni difetti, per lo meno quando si parla di diattica e di vettoriali, che lasciano talvolta il principiante (come me) un po' disorientato:

troppo estesi:
coprono territori troppo ampi e spesso non si riesce a cogliere al volo il risultato di un modulo (necessità di zoom, di ridurre il campo di ricerca, ecc.); inoltre le coordinate, essendo reali, sono espresse con numeri grandi, che distolgono l'attenzione (facendo un esercizio, un principiante è portato a riporre una maggior attenzione nel digitare le coordinate rispetto al contenuto sostanziale dell'operazione medesima);
troppo dettagliati:
l'occhio si perde nei dettagli lasciandosi sfuggire l'essenziale (il risultato dell'operazione);
troppo complessi:
intendendo con complessità l'elevato numero di strati informativi (vettoriali e raster) che compongono i datasets; da ciò consegue un buon numero di nomi di mappe (che confonde il principiante);
troppo ...reali:
il principiante (come me) deve imparare ad astrarre la realtà ragionando direttamente con i termini di punti, linee(archi), poligoni(aree) e non di hospitals, roads, landocover, ecc.
troppo grossi:
vorrei avere un set di dati leggero, che si possa inviare facilmente con un'e-mail.

Per quanto sopra, mi sono proposto di creare una location veramente semplice: location simpleloc!

location simpleloc

location simpleloc  è una location XY, non georiferita. La regione è quadrata e le coordinate vanno da 0 a 12; per creare location simpleloc  lanciare GRASS; indipendentemente dall'interfaccia (tcl/tk, wxPython, testuale) e dalla lingua, i passi da fare sono simili:

  1. location simpleloc  si chiama esattamente così, questo è il suo nome;
  2. secgliere la voce che maggiormente assomigli a ``creare una nuova location'';
  3. alla richiesta di specificare il sistema di coordinate scegliere ``XY'' (o qualcosa di simile come, ad es.: use arbitrary non-earth coordinate system (XY) ;
  4. alla richiesta di settare la regione di default, rispondere ``si (yes)'' e inserire questi dati:
    • north=12, south=0;
    • east=12, west=0;
    • n-s resolution=1;
    • e-w resolution=1

Per popolare location simpleloc  sono sufficienti tre file:

  1. due file di dati (sl_poly.dat e sl_points.dat);
  2. uno script (sl_makesimploloc.sh).

scaricabili insieme da qui.

Scompattare il file; rendere lo script eseguibile; entrare in GRASS  nella location location simpleloc  (in PERMANENT) e lanciare lo script: location simpleloc  è pronta.

Se si vuole location simpleloc  già pronta all'uso, la si può scarcare da questo link.


mtk iBlue e Quantum Gis

L’mtk iBlue 747 è un GPS data logger; le sue caratteristiche sono riportate qui.
Ho sempre avuto qualche difficoltà nella gesione del gps.
Esistono molti programmmi per scaricare tracce e punti dal gps, ma devo confessare che gli unici che hanno funzionato al volo sono mtkbabel (grazie N. Rigacci) e GPSBabel.

Entrambi sono installabili in Ubuntu Intrepid Ibex con

sudo apt-get install mtkbabel gpsbabel

Nei link sopra riportati esiste una buona documentazione.

In Quantum GIS è possibile impostare l’import delle tracce e dei waypoints. Ho provato nel seguente modo:

1. attivare il plugin GPS:

2. aprire la finestra Strumenti GPS


3. Nella finestra Edit devices aggiungere una nuova periferica (che potremmo chiamare mtk-usb inserendo i seguenti comandi nelle caselle scaricawaypoints e scarica track (ricordarsi di premere il pulsante Aggiornaperiferica):

%babel -w -i mtk -o gpx -f /dev/ttyUSB0 %in %out

%babel -t -i mtk -o gpx -f /dev/ttyUSB0 %in %out



A questo punto collegando l’IBlue 747 (in posizione di LOG) ad una prota USB è possibile scaricare i waypoint e le tracce direttamente in Quantum GIS:

domenica 7 dicembre 2008

Compilare GDAL e GRASS in Ubuntu Intrepid Ibex

Resoconto di una compilazione portata a termine da un principiante.

Ho faticato molto per installare GRASS e GDAL , soprattutto perchè non avevo letto attentamente quanto segue:

  1. la pagina apposita del wiki di GRASS;

  2. la pagina dedicata al plugin GRASS-GDAL ;

  3. questa utilissima guida;

  4. il file README che spiega i requisiti per l'interfaccia wxPython di GRASS ;

  5. il file README che accompagna il sorgente del plugin GRASS-GDAL .

In pratica questi appunti sono un sintetico riassunto di quanto sopra nonchè una semplice relazione di cosa ho fatto. Nel dubbio fare affidamento sulla documentazione sopra citata... (questi sono appunti scritti da un novellino, possono servire come traccia (vaga) ma non come regola).

Requisiti

Lavorando in Ubuntu 8.10 ho dovuto per prima cosa installato una certa quantità di pacchetti che risultavano mancanti, tra i quali molti dev:

$ sudo apt-get update && sudo apt-get upgrade
$ sudo apt-get install bison build-essential byacc \
fftw3 fftw3-dev flex freeglut3-dev libavcodec-dev \
libfreetype6-dev libgcc1 libgdal1-dev libgl1-mesa-dev \
libjpeg-dev libjpeg62-dev libncurses5-dev libpng12-dev \
libpq-dev libreadline5 libreadline5-dev libsqlite3-dev \
libtiff4-dev libxmu-dev libxmu-dev tk8.5-dev zlib1g-dev

Inoltre ho installato i seguenti pacchetti:

$ sudo apt-get install libwxbase2.8-0 libwxbase2.8-dbg libwxbase2.8-dev  \
libwxgtk2.8-0 libwxgtk2.8-dbg libwxgtk2.8-dev \
python-wxgtk2.8 wx2.8-doc wx2.8-examples \
wx2.8-headers wx2.8-i18n wx-common

Occorre inoltre controllare che GEOS e PROJ siano installati (controllare con Synaptic l'esistenza di pacchetti dev).

Ulteriori pacchetti dev potrebbero essere richiesti (vedere dopo).

I sorgenti

Ho deciso di avere tutti i sorgenti in /usr/local/src/.

Per prima cosa ho creato la cartella e ho fatto in modo che anche l'utente normale (non root) vi potesse scrivere:

$ sudo mkdir /usr/local/src
$ sudo chown mionomeutente:miogruppo /usr/local/src/

Mi sono posizionato in tale cartella e ho scaricato i sorgenti di GDAL con svn (che può essere installato con Synaptic):

$ cd /usr/local/src
$ svn checkout https://svn.osgeo.org/gdal/trunk/gdal gdal

viene scaricata la versione in sviluppo (questo può avere vantaggi e qualche svantaggio...) all'interno di una cartella denominata gdal.

La stessa operazione per GRASS :

$ cd /usr/local/src
$ svn checkout https://svn.osgeo.org/grass/grass/branches/develbranch_6 grass6_devel

viene creata la cartella /usr/local/src/grass6_devel con tutti i sorgenti di GRASS .

Infine ho scaricato da qui il file gdal-grass-1.4.3.tar.gz in /usr/local/src ed ho scompattato tale file (tar xvfz gdal-grass-1.4.3.tar.gz).

Compilare

La compilazione deve avvenire secondo un ordine preciso:

  1. GDAL

  2. GRASS

  3. plugin GRASS-GDAL

Regola generale: si configura (configure), si compila (make) e si installa (make install).

GDAL

Per configurare GDAL :

$ cd /usr/local/src/gdal
$ ./configure \
--with-pg=/usr/bin/pg_config \
--prefix=/usr/local \
--with-python \
--with-ogdi \
--with-sqlite \
--with-mysql=/usr/bin/mysql_config \
--with-libtiff=internal \
--without-grass \
--with-libtiff=internal

Non essendo farina del mio sacco rimando alla documentazione richiamata all'inizio per comprendere il significato delle diverse opzioni, nonchè all'help che può essere visualizzato digitando:

$ ./configure --help

Quello che davvero è importante è che GDAL deve esssere compilato senza il supporto di GRASS (--without-grass).

Il risultato di configure è il seguente:

GDAL is now configured for i686-pc-linux-gnu

Installation directory: /usr/local
C compiler: gcc -g -O2
C++ compiler: g++ -g -O2

LIBTOOL support: yes

LIBZ support: external
GRASS support: no
CFITSIO support: no
PCRaster support: internal
NetCDF support: yes
LIBPNG support: external
LIBTIFF support: internal (BigTIFF=yes)
LIBGEOTIFF support: internal
LIBJPEG support: external
LIBGIF support: external
OGDI support: no
HDF4 support: yes
HDF5 support: yes
Kakadu support: no
JasPer support: yes (GeoJP2=yes)
ECW support: no
MrSID support: no
MSG support: no
GRIB support: yes
cURL support (wms/wcs/...):yes
PostgreSQL support: yes
MySQL support: yes
Ingres support: no
Xerces-C support: yes
NAS support: no
Expat support: yes
ODBC support: yes
PGeo support: yes
OCI support: no
GEORASTER support: no
SDE support: no
DODS support: no
SQLite support: yes
DWGdirect support no
INFORMIX DataBlade support:no
GEOS support: yes


Old-gen python no
SWIG Bindings: python

Statically link PROJ.4: no
enable OGR building: yes
enable pthread support: no
hide internal symbols: no

Controllare tutte le voci del resoconto di configure; ci possono essere difficoltà legate al fatto che non si hanno abilitate librerire per alcune tipologie di file. In questo caso si può:

  1. cercare le librerie in Synaptic (pacchetti deb);

  2. cercare indicazioni nel wiki e negli archivi della mailing list;

  3. cercare indicazioni nella mailing list di GRASS ;

  4. cercare in Google;

  5. postare un quesito in qualche mailing list.

Qualora non si proceda alla compilazione ma, dopo modifiche, si vuole riconfigurare (cioè impartire un'altra volta il comando ./configure...), occorre sempre cancellare i file di configurazione creati con il precedente ./configure.

Ricordare:

make clean:
cancella ogni file creato da make;
make distclean:
cancella ogni file creato da make e ogni files creato dal ./configure.

Quindi, prima di ri-eseguire ./configure occorre make distclean.

A questo punto inizia la compilazione:

$ make

Se non vengono segnalati errori dopo la compilazione si procede all'installazione

$ sudo make install 

Se non specificato diversamente durante la fase di configurazione:

i binari di GDAL (ogr2ogr, gdalwarp, ectc.): vengono installati in /usr/local/bin. Per poterli usare questa cartella deve essere contenuta nei percorsi di ricerca dei comandi (vedi qui per capire cos'è un percorso di ricerca di comandi e cos'è la variabile PATH); se necessario quindi occorre modificare la variabile di ambiente PATH; tra gli infiniti esempi di come farlo vedere qui;

le librerie di GDAL (es.: libgdal.so.1.13.0): vengono posizionate in /usr/local/lib; anche queste devono poter essere viste (condivise). Si vedrà dopo cosa fare al riguardo.

GRASS

Leggere i files README e INSTALL nella directory /usr/local/src/grass6_devel/.

Molte delle cose dette al riguardo della compilazione di GDAL valgono anche per GRASS .


$ cd /usr/local/src/grass6_devel/
$ CFLAGS="-g -Wall" LDFLAGS="-s" ./configure \
--with-tcltk-includes=/usr/include/tcl8.5 \
--with-postgres=yes \
--with-postgres-includes=/usr/include/postgresql \
--with-sqlite \
--with-cxx \
--with-blas \
--with-lapack \
--with-cairo \
--with-fftw \
--with-freetype=yes \
--with-freetype-includes=/usr/include/freetype2 \
--with-readline \
--with-opengl-includes=/usr/include/GL \
--with-mysql \
--with-mysql-includes=/usr/include/mysql \
--with-python=/usr/bin/python2.5-config \
--with-wxwidgets=/usr/bin/wx-config \
--enable-largefile=yes \
--with-proj-includes=/usr/include \
--with-proj-share=/usr/share/proj \
--with-proj-libs=/usr/lib

Controllare che sul sistema in uso i percorsi richiamati esistano realmente. Per il supporto blas/lapack ho fatto ricorso a Synaptic installando alcuni pacchetti.

Più o meno quanto sopra riportato è preso da qui, dove, tra le molte altre cose interessanti, si legge:

Think twice before using this script. Some users experienced
problems such as disabled XGL etc.
...

Al termine della procedura di configurazione ottengo un resoconto come questo:

GRASS is now configured for:  i686-pc-linux-gnu

Source directory: /usr/local/src/grass6_devel
Build directory: /usr/local/src/grass6_devel
Installation directory: ${prefix}/grass-6.4.svn
Startup script in directory: ${exec_prefix}/bin
C compiler: gcc -g -Wall
C++ compiler: c++ -g -O2
Building shared libraries: yes
64bit support: no
OpenGL platform: X11
MacOSX application: no

NVIZ: yes

BLAS support: yes
C++ support: yes
Cairo support: yes
DWG support: no
FFMPEG support: no
FFTW support: yes
FreeType support: yes
GDAL support: yes
GLw support: no
JPEG support: yes
LAPACK support: yes
Large File support (LFS): yes
Motif support: no
MySQL support: yes
NLS support: no
ODBC support: no
OGR support: yes
OpenGL support: yes
PNG support: yes
PostgreSQL support: yes
Python support: yes
Readline support: yes
SQLite support: yes
Tcl/Tk support: yes
wxWidgets support: yes
TIFF support: yes
X11 support: yes

Valgono le stesse considerazioni fatte per GDAL ; se non si vede il supporto per un qualcosa di proprio interesse:

  • provare ad installare il qualcosa con Synaptic (installare anche il corrispondente pacchetto dev);

  • si modifichino le opzioni dello script configure in accordo con l'help (./configure --help) in funzione delle proprie esigenze (ricordare: make distclean cancella ogni file creato da make e ogni files creato dal ./configure);

  • cercare nella mailing list di GRASS , nel wiki, ecc.

È ovvio che per avere GRASS con la nuova interfaccia occorre abilitare il supporto a Python e a wxPython (--with-python e --with-wxwidgets); al riguardo può essere letto questo documento.

Se la configurazione fornisce i risultati desiderati si passa alla compilazione:

$ make

che richiede abbastanza tempo. Al termine dovrebbe esserci un messaggio tipo questo:

GRASS GIS compilation log
-------------------------
Started compilation: sab dic 6 12:20:39 CET 2008
--
Errors in:
No errors detected.
--
Finished compilation: sab dic 6 12:30:58 CET 2008

Si passa alla compilazione di vdigit come spiegato qui:

  1. cerco le librerie _gdi_.so:

    $ locate _gdi_.so
  2. creo un link simbolico con le librerie trovate

    $ sudo ln -s /usr/lib/python2.4/site-packages/wx-2.8-gtk2-unicode/wx/_gdi_.so \
    /usr/local/lib/libgdi.so
  3. compilo vdigit:

    $ cd gui/wxpython/vdigit
    $ make

ritorno nella directory /usr/local/src/grass6_devel e

$ sudo make install 

Affinchè le librerire condivise possano essere realmente tali (cioè...condivise!) ho modficato il file etc/ld.so.conf in questo modo:

$ cat /etc/ld.so.conf
include /etc/ld.so.conf.d/*.conf
/usr/local/lib
/usr/local/grass-6.4.svn/lib

Tale file deve essere editato con i privilegi di superutente (es.: gksu gedit /etc/ld.so.conf).

Ricordarsi di dare il comando sudo ldconfig dopo aver aggiornato il file etc/ld.so.conf.

Plugin GRASS /GDAL

GRASS dipende da GDAL ; GDAL dipende da GRASS . Questo plugin risolve il problema dell'interdipendenza. Leggere questo:

$ cd /usr/local/src/gdal-grass-1.4.3/
$ cat README | less

Seguendo questa guida, ho proseguito così:

$ cd /usr/local/lib/ 
$ sudo ln -s /usr/local/grass-6.4.svn/lib/*.so .
$ cd /usr/local/src/gdal-grass-1.4.3/
./configure --with-gdal=/usr/local/bin/gdal-config \
--with-grass=/usr/local/grass-6.4.svn
$ make
$ sudo mkdir /usr/local/lib/gdalplugins
$ sudo cp *.so /usr/local/lib/gdalplugins/

Nella directory /usr/local/lib/gdalplugins/ vengono copiati due files: gdal_GRASS.so e ogr_GRASS.so.

Compilazioni successive

GDAL

sudo apt-get update && sudo apt-get upgrade
make distclean
svn up
U ogr/ogrsf_frmts/grass/ogrgrassdatasource.cpp
U frmts/grass/grass57dataset.cpp
U frmts/grass/GNUmakefile
Aggiornato alla revisione 15902.

GRASS

sudo apt-get update && sudo apt-get upgrade
cd /usr/local/src/grass6_devel/
make distclean
svn up
U scripts/r.colors.stddev/r.colors.stddev
U vector/v.category/main.c
U lib/driver/Polygon.c
U ps/ps.map/description.html
U ps/ps.map/ps_vlegend.c
U gui/wxpython/gui_modules/menudata.py
U gui/tcltk/gis.m/gmmenu.tcl
U gui/tcltk/gis.m/Makefile

U gui/tcltk/d.m/Makefile
U gui/tcltk/d.m/menu.tcl
Aggiornato alla revisione 34740.

Plugin GDAL-GRASS

Controllare se qui esiste una nuova versione. Ripetere le operazioni già esposte