Random Network Analisys in GRASS

Questa è un'esercitazione introduttiva sulla Network Analisys che verrà condotta su un reticolo di linee e punti creati in modo casuale; su detto reticolo verranno applicati i diversi tools: creazione del grafo delle isodistanze dai punti, percorso minimo, problema del commesso viaggiatore, ecc.

Si farà uso del monitor per la rappresentazione dei risultati: (d.mon); alcune volte le operazioni usate potrebbero risultare di difficile comprensione a chi non avesse dimestichezza con la shell; si consiglia, in tal caso, di visualizzare i risultati nel Map Display nel modo usuale. L'intento è evidente: fare qualche semplice studio sui layouts per riutilizzarli in altri lavori.

Set di dati per l'esercitazione

Verrà usata, come al solito, location simpleloc; essendo la creazione del network completamente casuale, qualsiasi location può andare bene.

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

La location simpleloc può essere scaricata da qui.

Creazione del network

Un network è formato da archi e nodi tra loro connessi a formare un reticolo.

Il procedimento per creare il reticolo in modo casuale è il seguente:

1. si estraggono a sorte un certo numero di punti (v.random);
2. da ogni punto si ricavano le coordinate (v.out.ascii) e si memorizzano in un file di testo;
3. si leggono le coordinate dal file e si crea un vettoriale di linee (v.in.ascii);
4. si pulisce il vettoriale di linee (v.clean)
5. si estraggono nuovamente a sorte un certo numero di punti che andranno a costituire i nodi del network;
6. si crea la connessione (v.net) tra nodi e reticolo di linee (archi).

È consigliabile scrivere tutti i comandi in un file di testo ed eseguirlo come script; sarà interessante ripetere l'esercitazione variando il numero di punti casuali ed alcuni altri parametri.

Si riportano nel seguito i comandi impartiti, lasciando i relativi commenti.

# # # set region to default values
g.region -d
# # # set number of random points
numpoint="120"
# # # create a random points
v.random out=rndpoints n=$numpoint --overwrite
# # # create a temporary file to store points coordiantes
tmpascii="$(g.tempfile pid=$$)"
# # # first row in standard mode: (L)ine number of points
echo "L $numpoint" > ${tmpascii}
# # # append coordinates to ascii file; formatting by awk
v.out.ascii in=rndpoints fs=" "| awk '{print $1 " " $2}' >> ${tmpascii}
# # # inport ascii file; bad => before cleaning
v.in.ascii -n in=${tmpascii} out=rndlines_bad format=standard --overwrite

Il numero di punti è settato a 120. È particolarmente utile osservare come si possa creare il reticolo di linee importando le coordinate dal file ASCII generato dai punti: quale applicazione pratica si pensi alla creazione di linee partendo da punti rilevati con strumenti GPS.

L'output di v.out.ascii viene incanalato (pipe, |) attraverso awk per formattare l'output delle coordinate.

Si passa ora alla pulizia del vettoriale:

# # # v.clean: order of tools is significant:
# # # try to change it (also you have to change order of
# # # thresholds) and see how vector shape changes
v.clean input=rndlines_bad out=rndlines_nocat type=line \
 tool=break,snap,rmline,rmsa,rmdupl \
 thresh=0,0.4,0,0,0 --overwrite
v.build.polylines in=rndlines_nocat out=polylines --overwrite
# # # add category to vector lines
v.category in=polylines out=rndlines opt=add --overwrite
# # # clean up vectors
g.remove vect=rndpoints,rndlines_bad,rndlines_nocat,polylines
rm $tmpascii

v.clean è tra i comandi più difficili di GRASS: provare ad invertire l'ordine dei tools per capire come agiscono...

Il passaggio attraverso la creazione di polylines (v.build.polylines) è necessario per rimuovere gli pseudonodi.

Si noti che il vettoriale di linee non ha categoria e pertanto occorre aggiungerla con v.category.

Il reticolo di linee è pronto; si passa ai nodi, randomizzando nuovamente:

g.region -d
numpoint="5"
v.random out=netpoints n=$numpoint --overwrite

Ottenuto il reticolato di linee (=archi) e i nodi, si crea il network connettendoli. Per prima cosa aggiungere una tabella a netpoints:

# # # add table to points
v.db.addtable map=netpoints col="cat integer, dist double precision"

Il campo dist viene popolato con il comando v.distance che calcola, per ciscun elemento puntuale, la distanza alla più vicina linea.

# # # upload distances from points to the nearest line
v.distance from=netpoints to=rndlines upload=dist col=dist

L'opzione connect di v.net connette i punti al reticolo di linee entro una certa soglia (threshold); per determinare tale soglia si sceglie il valore massimo che assume il campo dist:

# # # get the maximun distnce from points to lines
max_dist=$(db.select -c sql="SELECT max(dist) FROM netpoints" )

...e lo si aumenta di una certa quantità per avere la sicurezza che tutti i punti vengano connessi: viene usato awk per eseguire il semplice conto.

max_dist=$(echo "$max_dist" | awk '{print $1+0.01}')

L'ultimo passo è la creazione del network vero e proprio con la connessione dei punti algi archi:

# # # create network by connecting points and arcs at given threshold
v.net in=rndlines points=netpoints \
 out=network op=connect \
 alayer=1 thresh=$max_dist --overwrite

L'output di v.net è un vettoriale con due layer: sul layer 1 vi sono le linee (archi); sul layer 2 vi sono i punti (nodi).

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

./sl_draw.sh
d.vect network layer=1 display=shape color=grey
d.vect network layer=2 display=shape icon=basic/circle \
 size=8 color=red fcolor=blue
d.vect network layer=2 display=cat llayer=2 \
 lsize=10 lcolor=blue

Il random-network creato

Network Analysis

La network analisys si basa sul concetto di costo necessario per percorrere un determinato arco o per attraversare un determinato nodo. In questa sede si assumono i valori di default:

1. il costo di viaggio lungo ciascun arco è pari alla lunghezza dell'arco medesimo;
2. il costo di attraversamento di ciascun nodo è pari a 0

Isodistanze (isocosti)

Il comando v.net.iso peremtte di definire delle zone (costs) di uguale costo di viaggio a partire dai nodi inidcati (ccats):

v.net.iso in=network out=isonet alayer=1 nlayer=2 \
 ccats=1-100 costs=2,3,5 --overwrite
v.db.addtable isonet
./sl_draw.sh
for catcolor in "1 red 2" "2 green 2" \
 "3 blue 2" "4 grey 1"; do
 # # # parsing parameters in group with set
 # # # $1=category $2=color $3=width
 # # # see http://tldp.org/LDP/abs/html/loops1.html
 set -- $catcolor
 d.vect isonet width=$3 display=shape \
 color=$2 where="cat=$1" 
done
d.vect network layer=2 display=cat \
 icon=basic/circle size=6 llayer=2 \
 bgcolor=white bcolor=orange \
 fcolor=black lcolor=black xref=center yref=center

Il grafo dello isodistanze

Percorso minimo

Con v.net.path vengono calcolati i percorsi minimi dal nodo di categoria 1 a tutti gli altri nodi.

Con g.tempfile viene creato un file di testo temporaneo che memorizza, per ogni riga, i seguenti elementi:

1. la categoria univoca del percorso che si va a cercare;
2. la categoria del punto di partenza (nell'esempio si cercano le distanze dal punto con cat=1);
3. la categoria del punto di arrivo (nell'esempio tutti gli altri punti creati).

Il file temporaneo viene usato come imput per v.net.path.

 # # distance from point with cat=1 to all others points
points_net=`g.tempfile pid=$$`
# # cat_short_path from_node_cat to_node_cat
echo "1 1 2" > $points_net
echo "2 1 3" >> $points_net
echo "3 1 4" >> $points_net
echo "4 1 5" >> $points_net
v.net.path input=network output=shortpath \
 file=$points_net --overwrite
rm $points_net
# # # display 
./sl_draw.sh
d.erase
# # # split monitor 
d.frame -c frame=1 at=51,99,1,49
d.frame -c frame=2 at=51,99,51,99
d.frame -c frame=3 at=1,49,1,49
d.frame -c frame=4 at=1,49,51,99
# # # loop; catpath is the shortest category path
# # # and the frame name (1,2,3,4)
for catpath in 1 2 3 4; do
 d.frame frame=$catpath -s
 d.vect network display=shape color=grey
 d.vect map=shortpath display=shape,cat \
 width=2 color=red cats=$catpath \
 bgcolor=white bcolor=blue \
 lcolor=blue
 d.vect network layer=2 display=shape,cat \
 icon=basic/circle size=6 llayer=2 \
 bgcolor=white bcolor=orange \
 fcolor=black lcolor=black xref=center yref=center
 string=$(v.db.select shortpath col=cost -c where="cat=$catpath")
 d.text at=2,94 text="cost=$string" color=blue
done

Merita notare le istruzioni impartite per leggere, per ogni percorso, la distanza tra i due punti; v.net.path crea un vettoriale con attributi:

v.db.select shortpath col=cost -c where="cat=$catpath"

I percorsi minimi dal nodo 1 a tutti gli altri nodi.

Il problema del commesso viaggiatore e l'albero di Steiner

Connettere una serie di punti attraverso un percorso più corto possibile toccando (attraversando) ciascun punto una sola volta: questo è il problema del commesso viaggiatore (v.net.salesman).

L'albero di Steiner (v.net.steiner) realizza la connessione tra i nodi dati con il percorso a minor costo (più corto possibile).

Per ricordare: se voglio pianificare una gita o una vacanza scelgo l'algoritmo del travel salesman; se voglio posare un tubo per portare acqua a diversi irrigatori uso l'albero di Steiner.

La sintassi dei due comandi è abbastanza simile:

v.net.salesman in=network \
 output=salesmanpath alayer=1 \
 nlayer=2 ccats=1,2,3,4,5 --overwrite
# # #
v.net.steiner input=network \
 output=steiner alayer=1 \
 nlayer=2 tcats=1-5 --overwrite

Si passa alla visualizzazione dei risultati:

# # # dispay
./sl_draw.sh; d.erase
# # # split monitor
d.frame -c frame=1 at=51,99,1,49
d.frame -c frame=2 at=51,99,51,99
d.frame -c frame=3 at=1,49,1,49
d.frame -c frame=4 at=1,49,51,99
# # # network in frame 1
d.frame -s frame=1
d.vect network display=shape color=grey
d.vect network layer=2 display=shape \
 icon=basic/circle size=6 fcolor=red
d.vect network layer=2 display=cat \
 llayer=2 lcolor=red lsize=10 \
 xref=left yref=bottom
d.text at=2,2 text="Network" color=black 
# # # salesmanpath in frame 2
d.frame -s frame=2
d.vect salesmanpath display=shape \
 color=magenta width=2
d.text at=2,2 text="Salesman path" color=black 
d.vect network layer=2 display=cat \
 bgcolor=white bcolor=orange llayer=2 \
 fcolor=black lcolor=black xref=center yref=center
# # # steiner tree in frame 3
d.frame -s frame=3
d.vect steiner display=shape \
 color=violet width=2
d.text at=2,2 text="Steiner tree" color=black 
d.vect network layer=2 display=cat \
 bgcolor=white bcolor=orange llayer=2 \
 fcolor=black lcolor=black xref=center yref=center
# # # which is the shortest ?
# # # steiner or salesman ?
# # # text in frame 4
d.frame -s frame=4
top=90
for vector in "salesmanpath" "steiner"; do
 # # # add table to vector
 v.db.addtable $vector col="length double precision"
 # # # upload vector length
 v.to.db $vector opt=length col=length
 # # # summarize length 
 length=$(db.select -c sql="SELECT sum(length) FROM $vector")
 # # # print 
 d.text.freetype -s text="$vector = $length" \
 size=10 at=2,$top color=black font=Vera
 # # # decrease text position in frame
 top=$(($top-8))
done

Il salesmann path e lo Steiner tree.

L'uso dei frames rende possibile la visualizzazione delle differenze; la stessa area geografica viene ripetuta in tre frames, ogni volta rappresentando qualcosa di diverso ma mantenendo costanti alcuni riferimenti (nel nostro caso: i nodi del network): non si vuole far vedere dove esattamente l'albero di Steiner e il Salesman path passano sul network; si vuole semplicemente rtrasmettere il messaggio che i due vettoriali realizzano connessioni con proprietà diverse.

Allocazione delle risorse

v.net.alloc crea una ripartizione del network tra i diversi nodi forniti in input:

v.net.alloc input=network output=alloc \
 alayer=1 nlayer=2 ccats=1,2,3,4,5 --overwrite

Il vettoriale di output non ha tabella attributi;se ne aggiunge una calcolando la lunghezza complessiva dei diversi subnets (uno per ciascuna riga della tabella) e quella totale di tutto il network:

# # # add table 
v.db.addtable alloc layer=1 col="length double precision"
# # # upload lenght
v.to.db alloc layer=1 col=length opt=length
# # # summarize total netwrok lenght
tot_length=$(db.select -c \
 sql="SELECT sum(length) FROM alloc")

Si procede ora alla visualizzazione:

# # # colorize
v.colors alloc range=1,5 color=bcyr column=cat

./sl_draw.sh; d.erase;
# # # split monitor into 4 frames
d.frame -c at=1,99,1,65 frame=1
d.frame -c at=67,99,66,99 frame=2
d.frame -c at=34,66,66,99 frame=3
d.frame -c at=1,33,66,99 frame=4

# # # overview in frame 1
d.frame -s frame=1
d.vect -a alloc
d.vect network layer=2 display=cat \
 bgcolor=white bcolor=orange llayer=2 \
 fcolor=black lcolor=black xref=center yref=center
 d.text.freetype -s at=2,2 text="Total lenght: ${tot_length}" \
 size=8 color=black
# # # first number = frame
# # # second number = subnet
for subnet in "2 2" "3 3" "4 5"; do
 set -- $subnet
 d.frame -s frame=$1
 d.vect alloc color=grey
 d.vect -a alloc where="cat=$2" width=2
 subnet_length=$(db.select -c \
 sql="SELECT sum(length) FROM alloc WHERE cat=$2")
 d.text.freetype -s at=2,2 text="subnet n. ${2}: ${subnet_length}" \
 size=8 color=black 
done

Il grafo dell'allocazione delle risorse.

Conclusioni

Non si è affrontato la Network Analisys in modo approfondito: giusto una presentazione dei principali tools che GRASS  mette a disposzione, senza considerare i versi di percorrenza degli archi ed i relativi costi.

Il codice di questa esercitazione si presta a molteplici varianti: creazione di vettoriali random, effetto dei tools di v.clean nel conferire un particolare aspetto finale al vettoriale, allestimento dei layouts con i soli comandi d.mon e d.vect.

---