martedì 6 gennaio 2009

Full planar topology in GRASS (Seconda parte)

Node Topology

Occorre ottenere un vettoriale di nodi; per fare ciò si può usare il tool v.to.points:

v.to.points -n input=arcs out=multiNode \
llayer=2 type=boundary --overwrite

I nodi vengono estratti dal layer 2 di arcs; volendo mantenere una relazione tra arco e nodo si usa l'opzione llayer che indica quale layer utilizzare (il $2$) per memorizzare la categoria dell'arco (=boundary) negli attributi dei nodi. v.to.points genera due layer: il primo con gli archi originali, il secondo con i nodi:

 v.db.select multiNode layer=2
cat|lcat|along
1|1|0
2|1|6.24264068711928
3|2|0
4|2|3
5|3|0
6|3|6.70820393249937
[...]
23|12|0
24|12|2.41421356237309
25|13|0
26|13|7.07106781186548

Il campo lcat memorizza la categoria dell'arco (layer 2 di arc) sul quale veine creato il nodo.

Sul layer 2 di multiNode si sono ottenuti i punti (=nodi) iniziali e finali di ciascun arco: va da sé che i nodi (finale di un arco ed iniziale dell'arco successivo) sono sovrapposti uno sull'altro.

Questo esercizio è una simulazione della Full Planar Topology: la tabella Node Topology necessita che i nodi siano unici ed univocamente individuati. Si procede allora a una pulizia dei nodi sovrapposti con v.clean tool=rmdupl:

v.clean in=multiNode out=nodeSimple \
type=point tool=rmdupl --overwrite

Al vettoriale risultato si attribuisce un valore di categoria univoco:

v.category nodeSimple out=node \
layer=3 opt=add --overwrite

v.clean tool=rmdupl ha rimosso i nodi sovrapposti da multiNode; il meccanismo di memorizzazione dei vettori implementato in GRASS  rende però possibile il recupero dei valori di categoria per ciascun layer grazie al mantenimento di una relazione uno a molti: interrogando dal layer 3 (parte uno) di nodeSimpel il layer 2 (parte molti), è possibile conoscere quanti e quali valori di categoria sono presenti in quest'ultimo layer. Non solo: con lo stesso meccanismo si possono popolare i campi del layer 3 con altri attributi del layer 2.

Per prima cosa occorre aggiungere alcuni campi al layer 3:

v.db.addtable node layer=3 \
columns="cat integer, countnodes integer, \
countarcs integer, nodeslist char(12), \
arcslist char(12), arcslistnames char(20)"

Il tool che consente di realizzare quanto sopra descritto è v.to.db opt=query : per ogni geometria esegue la funzione contenuta in qcolumn interrogando il layer qlayer e popolando il campo definito in column. I commenti lasciati fra le righe di comando servono da spiegazione:

# # # numero di nodi sovrapposti nel layer 2
v.to.db node option=query layer=3 \
qlayer=2 column=countnodes qcolumn="count(cat)"
# # # numero di archi collegati ai nodi
v.to.db node option=query layer=3 \
qlayer=2 column=countarcs qcolumn="count(lcat)"
# # # elenco delle categorie dei nodi sovrapposti sul layer 2
v.to.db node option=query layer=3 \
qlayer=2 column=nodeslist qcolumn="group_concat(cat)"
# # # elenco delle categorie delgi archi collegati
v.to.db node option=query layer=3 \
qlayer=2 column=arcslist qcolumn="group_concat(lcat)"
# # # elenco dei nomi degli archi collegati
v.to.db node option=query layer=3 \
qlayer=2 column=arcslistnames qcolumn="group_concat(arcLabel)"

Il risultato è il seguente:

v.db.select node layer=3
cat|countnodes|countarcs|nodeslist|arcslist|arcslistnames
1|2|2|1,2|1,1|a1,a1
2|3|3|6,9,17|3,5,9|a3,a5,a9
3|3|3|15,18,19|8,9,10|a8,a9,a10
4|3|3|4,5,20|2,3,10|a2,a3,a10
5|4|4|7,10,11,22|4,5,6,11|a4,a5,a6,a11
6|4|4|8,12,14,23|4,6,7,12|a4,a6,a7,a12
7|3|3|21,24,25|11,12,13|a11,a12,a13
8|4|4|3,13,16,26|2,7,8,13|a2,a7,a8,a13

Si può verificare la situazione delle categorie nei diversi layers:

v.category node opt=report
Layer/table: 1/node_1
type count min max
point 25 1 13
line 0 0 0
boundary 0 0 0
centroid 0 0 0
area 0 0 0
all 25 1 13
Layer/table: 2/node_2
type count min max
point 26 1 26
line 0 0 0
boundary 0 0 0
centroid 0 0 0
area 0 0 0
all 26 1 26
Layer/table: 3/node_3
type count min max
point 8 1 8
line 0 0 0
boundary 0 0 0
centroid 0 0 0
area 0 0 0
all 8 1 8

La tabella Node Topology è semplicemente:

 v.db.select node layer=3 col=cat,arcslistnames
cat|arcslistnames
1|a1,a1
2|a3,a5,a9
3|a8,a9,a10
4|a2,a3,a10
5|a4,a5,a6,a11
6|a4,a6,a7,a12
7|a11,a12,a13
8|a2,a7,a8,a13

L'arco a1 si chiude su se stesso nel nodo di categoria 1. La tabella prodotta presenta in corrispondenza del nodo 1 la ripetizione del nome dell'arco: ciò costituisce una differenza fra quanto, sulla base del modello Full Planar Topology, dovrebbe essere.

Si crea ora una visualizzazione del vettoriale node :

./sl_draw.sh
d.font.freetype font=timesbd
d.vect arcs layer=2 display=shape,attr attrcol=arcLabel \
type=boundary llayer=2 \
xref=center yref=center color=grey \
lsize=10 bgcolor=white bcolor=black
d.vect node layer=3 display=shape,cat llayer=3 \
size=5 icon=basic/triangle lcolor=aqua \
xref=center yref=center \
lsize=10 bgcolor=white bcolor=aqua

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).

Relazioni tra i vettoriali node e arcs.

Polygon Topology

La terza ed ultima tabella Polygon Topology viene ottenuta combinando gli attributi sui due layers del vettoriale arcs.

Nel campo arcList sarà memorizzato l'elenco degli archi che delimitano ciascuna area:

 v.db.addcol arcs col="arcList char(15)"

Si provi dapprima questa istruzione SQL:

echo "select a.cat as cat,a.label as label, \
group_concat(b.arcLabel) as arcList \
from arcs a, arcs_2 b \
WHERE (a.cat==b.catPolyLeft) OR \
(a.cat==b.catPolyRight) \
GROUP BY a.cat" | db.select

Ancora una volta si fa uso della funzione di aggregazione comando di aggregazione group_concat supportato nelle recenti versioni di SQLite; inoltre le diverse tabelle sono indicate da semplici lettere.

Si crea una tabella temporanea per memorizzare il risultato dell'istruzione SELECT e subito dopo la si interroga per aggiornare il valore del campo arcList:

echo "CREATE TEMPORARY TABLE foo AS \
select arcs.cat as cat, group_concat(arcs_2.arcLabel) as List \
from arcs, arcs_2 \
WHERE (arcs.cat==arcs_2.catPolyLeft) OR \
(arcs.cat==arcs_2.catPolyRight) \
GROUP BY arcs.cat;
UPDATE arcs SET arcList=(select List FROM \
foo WHERE arcs.cat=foo.cat)" | db.execute

Notare che:

  1. si è usata la notazione consueta: nomeTabella.nomeCampo anzichè le lettere per indicare le tabelle;
  2. si sono concatenate le due istruzioni SQL con il ';';
  3. viene creata una tabella temporanea; la sintassi completa del comando CREATE TABLE è riportato manuale di SQLite;
  4. la tabella temporanea verrà distrutta appena la connessione al database cesserà.

Conclusioni

Arrivati al fondo di questa tortuosa esercitazione merita puntualizzare alcune cose:

  1. l'esercitazione condotta è stata una simulazione: GRASS  implementa il modello topologico in modo ben più efficiente!
  2. la procedura si presta a numerose varianti: si provi ad esempio ad usare l'opzione -v (anziché -n) del tool v.to.points; la via seguita è quella che si è ritenuto più didattica...
  3. mentre le tabelle Node Topology e Arc Topology sono state ottenute interrogando le geometrie, la terza tabella, Polygon Topology è stata ottenuta ricombinando altre tabelle; si provi ad ottenere Polygon Topology direttamente dalle geometrie...
  4. per i debuttanti di GRASS: non scoraggiarsi; alcune cose del tipo
     v.category nodeSimple out=node layer=3 opt=add --overwrite
    non si capiscono al volo (a meno di persone particolarmente portate). Questo è normale: rifare alcune volte un esercizio non è una perdita di tempo.

Tutti i comandi usati in questa esercitazione sono contenuti in uno script scaricabile da da questo link.


Nessun commento:

Posta un commento