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.



Nessun commento:

Posta un commento