Combinazione di linee con punti

Questa esercitazione riguarda la combinazione di primitive vettoriali di tipo linea e punti. Le operazioni effettuate sono:

Calcolo della distanza di punti e linee e determinazione del tipo di linea più vicina a punti

Si vuole determinare la distanza dei siti archeologici contenuti nella mappa archsites dalla rete stradale (mappa roads):

Si imposta la regione su quella di default con

> g.region -d

e si visualizzano la rete stradale in nero ed i siti archeologici in rosso

> d.erase

> d.vect map=roads@PERMANENT

> d.vect map=archsites@PERMANENT color=red size=6 icon=basic/pushpin


Siti archeologici di Spearfish estrade
Si calcolano le distanze dalle strade dei punti di archsites e le si visualizza sul terminale, il parametro upload vale dist per indicare che si vuole calcolare la distanza, il parametro column indica il nome (dist in questo caso) per la colonna che conterrà la distanza, la flag -p indica che l'output deve essere mandato sul terminale:

> v.distance from=archsites@PERMANENT to=roads@PERMANENT from_type=point to_type=point,line,area from_layer=1 to_layer=1 dmax=-1 upload=dist column=dist -p

Si visualizza:

from_cat|dist
1|58.862398
2|521.737468
3|52.485691
4|62.672967
5|3.133039
6|25.666562
7|0.392152
8|53.277011
9|80.370538
10|19.351228
11|77.802276
12|884.668624
13|222.098356
14|256.744119
15|50.195620
16|105.516975
17|305.790858
18|146.890604
19|578.305319
20|183.446934
21|596.708712
22|935.814755
23|846.492376
24|2.499839
25|39.683481
Statistics:
0 categories with more than 1 feature in 'from'
0 categories - no nearest feature found

dove la prima colonna indica la categoria del punto della mappa archsites e la seconda la distanza dalla strada più vicina nella mappa roads.
In modo simile si visualizza il tipo di strada più vicino a ciascun punto, il parametro upload vale to_attr per indicare che si vuole determinare l'attributo delle linea più vicina, il parametro to_column vale label per indicare la colonna della tabella roads da cui si vuole estrarre l'attributo, il parametro column indica il nome (road_type in questo caso) per la colonna che conterrà il tipo di strada, la flag -p indica che l'output deve essere mandato sul terminale:

> v.distance from=archsites@PERMANENT to=roads@PERMANENT from_type=point to_type=point,line,area from_layer=1 to_layer=1 dmax=-1 upload=to_attr column=road_type to_column=label -p

Viene visualizzato:

from_cat|road_type
1|light-duty road, improved surface
2|unimproved road
3|unimproved road
4|secondary highway, hard surface
5|interstate
6|primary highway, hard surface
7|light-duty road, improved surface
8|interstate
9|secondary highway, hard surface
10|unimproved road
11|secondary highway, hard surface
12|unimproved road
13|unimproved road
14|unimproved road
15|light-duty road, improved surface
16|primary highway, hard surface
17|secondary highway, hard surface
18|unimproved road
19|light-duty road, improved surface
20|light-duty road, improved surface
21|secondary highway, hard surface
22|unimproved road
23|light-duty road, improved surface
24|secondary highway, hard surface
25|unimproved road
Statistics:
0 categories with more than 1 feature in 'from'
0 categories - no nearest feature found

Se si vuole visualizzare i segmenti di minima distanza tra i siti archeologici contenuti nella mappa archsites dalla rete stradale (mappa roads) si aggiunge il parametro output in cui si indica il nome della mappa vettoriale in cui salvare i segmenti di minima distanza

> v.distance from=archsites@PERMANENT to=roads@PERMANENT output=archistes_to_roads from_type=point to_type=point,line,area from_layer=1 to_layer=1 dmax=-1 upload=dist column=dist --overwrite

e si visualizzano i segmenti creati in verde, senza cancellare prima il monitor:

> d.vect archistes_to_roads color=green


Distanze tra i siti archeologici di Spearfish e le strade

Calcolo della distanza di punti e linee e determinazione del tipo di linea più vicina a punti e aggiunta al database

Le distanze e il tipo di strada più vicino determinate nelle procedure precedenti possono essere aggiunti alla tabella dei punti, invece che visualizzati come nei casi precedenti. Si devono aggiungere alla tabella due colonne in cui inserire la distanza e il tipo di strada più vicino: poichè la mappa archsites è nella location PERMANENT essa è in sola lettura e la sua tabella non può essere modificata. Per questo motivo le procedure di seguito sono applicate ai soli siti archeolgici all'interno della Black Hills Natl. Forest, contenuti nella mappa archsites_natfor creata nella esercitazione Combinazione di aree con punti, la quale, essendo nel mapset corrente, è modificabile.

Si devono innnanzitutto aggiungere le colonne alla tabella archsites_natfor, se esistono già (ad es. perchè si è già fatta questa procedura) db.execute da errore, ma si può ignorare:

> echo "ALTER TABLE archsites_natfor ADD COLUMN road_type varchar(30)" | db.execute

> echo "ALTER TABLE archsites_natfor ADD COLUMN dist double" | db.execute

la colonna road_type è di tipo varchar e conterrà la stringa che identifica il tipo di strada, la colonna dist è di tipo double e conterrà le distanze.
Si aggiunge alla tabella degli archsites_natfor nella colonna road_type il tipo di strada più vicina, il parametro upload vale to_attr per indicare che si vuole determinare l'attributo delle linea più vicina, il parametro to_column vale label per indicare la colonna della tabella roads da cui si vuole estrarre l'attributo, il parametro column indica il nome (road_type in questo caso) per la colonna che conterrà il tipo di strada, a differenze della procedura precedente manca la flag -p che indica che l'output deve essere mandato sul terminale:

> v.distance from=archsites_natfor to=roads@PERMANENT from_type=point to_type=point,line,area from_layer=1 to_layer=1 dmax=-1 upload=to_attr column=road_type to_column=label

Si aggiunge poi alla tabella degli archsites_natfor nella colonna "dist" la distanza dalla strada più vicina, il parametro upload vale dist per indicare che si vuole calcolare la distanza, il parametro column indica il nome (dist in questo caso) per la colonna che conterrà la distanza:

> v.distance from=archsites_natfor to=roads@PERMANENT from_type=point to_type=point,line,area from_layer=1 to_layer=1 dmax=-1 upload=dist column=dist

Si visualizza quindi la nuova tabella della mappa archsites_natfor con:

> echo "select * from archsites_natfor"|db.select

che riporta

cat|str1|road_type|dist
1|Signature Rock|light-duty road, improved surf|58.862398
2|No Name|unimproved road|521.737468
3|Canyon Station|unimproved road|52.485691
4|Spearfish Creek|secondary highway, hard surfac|62.672967
10|Slaughterhouse Gulch|unimproved road|19.351228
13|No Name|unimproved road|222.098356
14|No Name|unimproved road|256.744119
16|Boulder Creek Cabin|primary highway, hard surface|105.516975
18|Cole Creek Mine|unimproved road|146.890604
19|No Name|light-duty road, improved surf|578.305319
20|No Name|light-duty road, improved surf|183.446934
22|Bob Miller|unimproved road|935.814755
23|No Name|light-duty road, improved surf|846.492376
25|No Name|unimproved road|39.683481

Ad ogni sito archeologico è quindi associato il tipo di strada più vicina e la distanza da essa.

Creazione di punti lungo una linea a intervalli regolari

Si crea una nuova mappa vettoriale che contiene un insieme di punti a distanza massima data lungo la ferrovia nella mappa railroads. La distanza indicata è quella massima perchè viene comunque inserito un ulteriore punto se il punto alla distanza indicata corrisponde ad ogni cambio di direzione: se si vuole inserire punti ad una distanza esatta si deve usare v.segment.

Si visualizza la ferrovia, cancellando prima il monitor, con

> d.erase

> d.vect map=railroads

e si creano i punti nella mappa railroads_points a distanza massima di 1000 m (indicata con l'attributo dmax) lungo la ferrovia nella mappa railroads.

> v.to.points input=railroads@PERMANENT type=point,line,boundary,centroid output=railroads_points llayer=1 dmax=1000 --overwrite

e si visualizza il risultato, senza cancellare il monitor:

> d.vect map=railroads_points color=red size=5 icon=basic/circle


Punti a distanza fissata sul tracciato della ferrovia
La nuova mappa ha la tabella della mappa di partenza collegata sul primo layer, in questo caso la mappa railroads non ha tabelle collegate e quindi nemmeno la mappa railroads_points ha tabelle collegate sul primo layer, ed una nuova tabella contente per ogni punto creato una colonna lcat che indica la categoria del segmento della mappa di partenza (railroads in questo caso) a cui il punto è sovrapposto ed una colonna along che contiene la distanza progressiva dal primo nodo. Questa tabella ha nome uguale alla mappa creata più il suffisso "_2".
Si visualizza quest'ultima tabella con:

> echo "select * from railroads_points_2"|db.select

che scrive sul terminale (per brevità solo le prime righe sono riportate)

cat|lcat|along
1|-1|0.000000
2|-1|991.740216
3|-1|1983.480431
4|-1|2975.220647
5|-1|3966.960862

E' possibile creare punti solo sui vertici (cambi di direzione) di una linea con l'opzione -v

> v.to.points input=railroads@PERMANENT type=point,line,boundary,centroid output=railroads_points llayer=1 dmax=1000 -v --overwrite

e si visualizza con:

> d.erase

> d.vect railroads

> d.vect map=railroads_points color=red size=5 icon=basic/circle


Cambio di direzione sul tracciato della ferrovia

Creazione di punti e linee lungo una linea a distanza date lungo la linea stessa

Questa procedura crea punti o segmenti lungo una linea di categoria data, nel caso di punti si deve indicare la distanza lungo la linea stessa dal primo vertice della linea, nel caso di segmenti si deve indicare la distanza di inizio e fine del segmento stesso. Si deve ovviamente sapere quale è il primo vertice della linea lungo cui si creano i punti e/o i segmenti, questo si può facilmente visualizzare indicando l'opzione dir nel parametro display di d.vect, e le categorie delle linee, visualizzabili con l'opzione cat nel parametro display di d.vect:

> d.erase

> d.vect railroads

Si crea prima una mappa railroads_segment_point contenente un punto con categoria 1 lungo la linea con categoria 1 a distanza 1000 dall'inizio della linea stessa

> echo "P 1 1 1000"|v.segment input=railroads@PERMANENT output=railroads_segment_point llayer=1 --overwrite

Si crea poi una mappa railroads_segment_line contenente un segmento con categoria 2 lungo la linea con categoria 1 con inizio a distanza 2000 dall'inizio della linea e fine a distanza 3000

> echo "L 2 1 2000 3000"|v.segment input=railroads@PERMANENT output=railroads_segment_line llayer=1 --overwrite

e si visualizza il risultato, senza cancellare il monitor; il punto in rosso

> d.vect map=railroads_segment_point color=red size=5 icon=basic/circle

e il segmento in verde

> d.vect map=railroads_segment_line color=green

Per avere come output una sola mappa che contiene sia punti che segmenti è sufficiente unire le due mappe appena create con v.patch:

> v.patch input=railroads_segment_line,railroads_segment_point output=railroads_segment --overwrite

e si visualizza, ora tutti sia punti che linee sono rossi (si possono visualizzare con colori diversi usando il parametro type del modulo d.vect)

> d.erase

> d.vect railroads

> d.vect map=railroads_segment color=red size=5 icon=basic/circle

Divisione di una linea in segmenti

E' possibile dividere una linea in segmenti di lunghezza massima data con il modulo v.split. In questa procedura si dividono le linee nella mappa railroads in segmenti di 100 metri, creando la mappa railroads_split_lenght:

> v.split input=railroads@PERMANENT output=railroads_split_lenght length=100 --overwrite

e si visualizza

> d.erase

> d.vect railroads

> d.vect map=railroads_split_lenght color=red

le linee sono ovviamente sovrappposte, ma si può verificare la differenza di lunghezza nelle due mappe con

> d.what.vect

che restituisce per le due mappe, ad es.

La lunghezza del segmento è di circa 3000 metri nella mappa railroads e di circa 100 metri nella mappa railroads_split_lenght.
Si può verificare anche che è cambiato il numero di segmenti per ogni categoria, nella mappa originale railroads

> v.to.db map=railroads type=point,line,boundary,centroid layer=1 qlayer=1 option=count units=meters column=lenght -p

restituisce

cat|count
-1|7
1|17
2|8

e nella mappa railroads_split_lenght

> v.to.db map=railroads_split_lenght type=point,line,boundary,centroid layer=1 qlayer=1 option=count units=meters column=lenght -p

cat|count
-1|1914
1|826
2|139

il numero di segmenti è quindi molto aumentato.

E' possible inoltre dividere una linea in segmenti con un numero massimo di vertici, si indica cioè il numero di vertici consecutivi della linea di partenza che ogni segmento deve contenere (2=solo segmenti rettilinei), ovviamente quindi diminuendo il numero massimo di vertici il numero di segmenti aumenta e viceversa. Si vuole dividere la mappa railroads in segmenti con al più 10 vertici, creando la mappa railroads_split_num:

> v.split input=railroads@PERMANENT output=railroads_split_num vertices=10 --overwrite

e si visualizza

> d.erase

> d.vect railroads

> d.vect map=railroads_split_num color=red

le linee sono ovviamente sovrappposte, ma come in precedenza si può verificare la differenza di lunghezza nelle due mappe con

> d.what.vect

che restituisce per le due mappe, ad es.

Ovviamente, solo le linee con più di 10 vertici nella mappa di partenza sono divise. Come prima si può verificare anche che è cambiato il numero di segmenti per ogni categoria, nella mappa originale railroads

> v.to.db map=railroads type=point,line,boundary,centroid layer=1 qlayer=1 option=count units=meters column=lenght -p

restituisce come prima

cat|count
-1|7
1|17
2|8

e nella mappa railroads_split_num

> v.to.db map=railroads_split_num type=point,line,boundary,centroid layer=1 qlayer=1 option=count units=meters column=lenght -p

cat|count
-1|7
1|72
2|17

In alternativa si possono assegnare categorie diverse ad ogni segmento sul secondo layer e visualizzare con colori diversi i segmenti (che ora hanno categoria diversa). Si aggiunge un secondo layer con categorie da 1 al numero di segmenti, con incremento 1:

> v.category input=railroads_split_num output=railroads_split_num_cat type=point,line,boundary,centroid,area option=add cat=1 layer=2 step=1 --overwrite

e si visualizzano i segmenti 3, 4 e 5 rispettivamente in rosso, verde e blu

> d.erase

> d.vect map=railroads_split_num_cat layer=2 llayer=2 cats=3 color=red

> d.vect map=railroads_split_num_cat layer=2 llayer=2 cats=4 color=green

> d.vect map=railroads_split_num_cat layer=2 llayer=2 cats=5 color=blue

Lo script che realizza queste operazioni è geoprocessing_spearfish_linee_punti.sh.