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.

---