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.


Full planar topology in GRASS (Prima parte)

Un breve richiamo: la Full Planar Topology

La Full Planar Topology è un modello topologico poligonale che soggiace ad alcune regole fondamentali:

Connettività:
gli archi sono tra loro connessi con un nodo in comune; ogni arco inizia e termina con un nodo (ogni arco ha due nodi); tutte le intersezioni tra archi sono nodi.
Contiguità:
poligoni adiacenti hanno un solo arco in comune; ogni arco divide due aree.
Contenimento:
un poligono chiuso da archi ha un'area misurabile; serie di archi definiscono perimetri e aree chiuse. Ogni area è circondata da archi e nodi. Ogni nodo è contenuto in aree.

Poste queste regole, il vettoriale polyTopology della location simpleloc rappresentato qui di seguito, dovrebbe avere le tabelle riportate sotto:

Le entità geometriche fondamentali: nodi, archi, poligoni. Il verso degli archi si riferisce al verso di digitazione.

    

Id (cat) nodo Archi
1 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
Node topology.

    

Id arco L_Poly R_Poly
a1 polyA polyC
a2 polyA polyB
a3 polyA polyD
a4 universe polyE
a5 polyA universe
a6 polyA universe
a7 polyF polyA
a8 universe polyB
a9 universe polyD
a10 polyB polyD
a11 universe polyE
a12 polyF polyE
a13 polyF universe
Arc Topology

    

Id (cat) poligono label Archi
1 polyA a1,a2,a3,a5,a6,a7
2 polyB a2,a8,a10
3 polyC a1
4 polyD a3,a9,a10
5 polyE a4,a11,a12
6 polyF a7,a12,a13
Polygon Topology

Obbiettivi

GRASS  è un GIS  topologico; le 3 tabelle fondamentali viste sopra devono quindi essere facilmente ottenibili.

Il comando

 v.build map=polyTopology opt=dump

consente la visualizzazione di tutti gli elementi della topologia del vettoriale polyTopology.

Lo scopo di questo esercizio è didattico: nella costruzione delle tre tabelle topologiche, si condurrà una simulazione senza ricorrere a v.build; questo consentirà di esplorare al meglio alcuni tools che interagiscono con le entità geometriche fondamentali (nodi,archi,poligoni).

Set di dati per l'esercitazione

La copertura poligonale scelta per l'esercitazione è polyTopology, 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.

Procedura

Si creerà per prima la tabella Arc Topology, quindi la Node Topology e da ultimo la Polygon Topology.

Le fasi logiche sono le seguenti:

  1. creazione di un layer di archi (=boundaries) numerato univocamente;
  2. agli attributi degli archi vengono aggiunti gli identificativi (=cat) dei poligoni posti a destra e a sinistra di ciascun elemento;
  3. creazione di un layer di nodi (=punti) numerati univocamente;
  4. agli attributi dei nodi vengono aggiunti gli identificativi degli archi ai quali sono uniti;
  5. dalla combinazione delle tabelle dei nodi e degli archi si ottiene la Polygon Topology.

Arc Topology

I boundaries per default non hanno una categoria. Si crea un nuovo layer aggiungendo un valore di categoria univoco:

v.category in=polyTopology type=boundary \
layer=2 out=arcs option=add cat=1 --overwrite

Agli attributi del layer appena creato si aggiunge una serie di campi: alcuni conterrano solo testo e serviranno per la visualizzazione. Inoltre alcuni campi testuali offorno l'occasione per fare esercizio anche con l'SQL.

v.db.addtable arcs layer=2 columns="cat integer, \
arcLabel char(3), \
catPolyLeft integer, \
polyLeftLabel char(5), \
catPolyRight int, \
polyRightLabel char(5), \
arcLength double precision"

Si inizia a popolare il campo arcLabel:

v.db.update map=arcs layer=2 column=arcLabel qcolumn="'a'||cat"

L'operatore SQL || concatena le stringhe; in questo caso concatena la stringa a con il contenuto del campo cat.

Con il tool v.to.db opt=sides si ottengono le categorie dei poligoni rispettivamente a sinistra e a destra di ciascun arco (boundary):

v.to.db arcs layer=2 option=sides col=catPolyLeft,catPolyRight

Si calcola inoltre la lunghezza di ciascun arco:

v.to.db arcs layer=2 option=length col=arcLength

Viene ora il momento di usare l'SQL: sulla base delle categorie delle aree ottenute con v.to.db opt=sides si popolano i due campi polyLeftLabel e polyRightLabel:

echo "UPDATE arcs_2 SET polyLeftLabel=IFNULL( \
(SELECT arcs.label FROM arcs \
WHERE arcs.cat=arcs_2.catPolyLeft),'universe')" | db.execute

echo "UPDATE arcs_2 SET polyRightLabel=IFNULL( \
(SELECT arcs.label FROM arcs \
WHERE arcs.cat=arcs_2.catPolyRight),'universe')" | db.execute

L'operazione effettuata merita qualche spiegazione; sul layer 1 del vettoriale arcs vi sono gli attributi delle aree:

v.db.select arcs layer=1
cat|label|GRASSRGB
1|polyA|141:211:199
3|polyC|190:186:218
2|polyB|255:255:179
4|polyD|251:128:114
5|polyE|128:177:211
6|polyF|253:180:98

Nei campi catPolyLeft e catPolyRight del layer2 sono state memorizzate le categorie delle aree del layer 1: esiste quindi una relazione biunivoca fra i due campi del layer 2 e il campo cat del layer1. Grazie a questa relazione si può stabilire un join fra le due tabelle imponendo la condizione WHERE arcs.cat=arcs_2.catPolyLeft e WHERE arcs.cat=arcs_2.catPolyRight.

Se tali condizioni non vengono soddisfatte significa che il valore della cat del layer 1 non ha corrispondenze nel layer 2: è il caso del poligono universe (valore -1), l'area che, nel modello topologico Full Planar Topology avvolge idealmente ogni altro oggetto.

Dopo queste operazioni gli attributi sul layer 2 sono:

v.db.select arcs layer=2
cat|arcLabel|catPolyLeft|polyLeftLabel|catPolyRight|polyRightLabel|arcLength
1|a1|3|polyC|1|polyA|6.242641
2|a2|2|polyB|1|polyA|3
3|a3|4|polyD|1|polyA|6.708204
4|a4|5|polyE|-1|universe|3.414214
5|a5|-1|universe|1|polyA|5.019765
6|a6|-1|universe|1|polyA|2.828427
7|a7|1|polyA|6|polyF|4.162278
8|a8|2|polyB|-1|universe|5.162278
9|a9|4|polyD|-1|universe|10.261297
10|a10|4|polyD|2|polyB|3
11|a11|5|polyE|-1|universe|5.828427
12|a12|5|polyE|6|polyF|2.414214
13|a13|-1|universe|6|polyF|7.071068

Si verifichi anche la situazione delle categorie sui due layers:

v.category arcs opt=report
Layer/table: 2/arcs_2
type count min max
point 0 0 0
line 0 0 0
boundary 13 1 13
centroid 0 0 0
area 0 0 0
all 13 1 13
Layer/table: 1/arcs
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

A questo punto la costruzione della tabella Arc Topology è conclusa:

v.db.select arcs layer=2 \
col=arcLabel,polyRightLabel,polyLeftLabel
arcLabel|polyRightLabel|polyLeftLabel
a1|polyA|polyC
a2|polyA|polyB
a3|polyA|polyD
a4|universe|polyE
a5|polyA|universe
a6|polyA|universe
a7|polyF|polyA
a8|universe|polyB
a9|universe|polyD
a10|polyB|polyD
a11|universe|polyE
a12|polyF|polyE
a13|polyF|universe

Si procede ora alla visualizzazione :

./sl_draw.sh
d.font.freetype font=timesbd
d.vect -a arcs layer=1 display=shape
d.vect arcs layer=1 display=attr attrcol=Label llayer=1 \
color=black lsize=10 xref=center yref=bottom
d.vect arcs layer=2 display=shape,attr,dir attrcol=arcLabel \
type=boundary llayer=2 \
xref=center yref=center color=grey \
lsize=10 bgcolor=white bcolor=black

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

Il vettoriale arcs: sul layer 1 mantiene le aree; sul layer 2 gli archi (=boundaries) con gli attributi calcolati.