mercoledì 30 dicembre 2009

Qgis e la sorgente dati non valida.

Raramente uso Qgis per l'elaborazione di raster memorizzati in mapset di GRASS; ho quindi, per lungo tempo, sopportato passivamente l'errore che ricevevo quando tentavo di caricare un raster di GRASS in Qgis: sorgente dati non valida.
L'ultima versione di Qgis ("Enceladus") è però molto ricca  e l'uso della shell di GRASS nel plugin per Qgis pare uno strumento efficace; ho fatto allora  qualche tentativo (positivo) per risolvere il problema dei raster.

Ecco gli appunti (non è... oro colato) di una soluzione poco ortodossa e temporanea:

 Linux Ubuntu Karmic; nel file sources.list ho, tra gli altri, questi repositories:

deb http://ppa.launchpad.net/ubuntugis/ubuntugis-unstable/ubuntu karmic main
deb-src http://ppa.launchpad.net/ubuntugis/ubuntugis-unstable/ubuntu karmic main

Per il gioco delle dipendenze, nel sistema sono installati questi files:

/usr/lib/gdal15plugins/gdal_GRASS.so
/usr/lib/gdal15plugins/ogr_GRASS.so

/usr/lib/gdal16plugins/ogr_GRASS.so
/usr/lib/gdal16plugins/gdal_GRASS.so

Sospetto che i due pacchetti gdal vadano in conflitto; rinomino il pacchetto meno aggiornato e creo un link:

pcgis@laptop:/usr/lib/gdal15plugins$ sudo mv gdal_GRASS.so copia_di_gdal_GRASS.so
pcgis@laptop:/usr/lib/gdal15plugins$ sudo ln -s /usr/lib/gdal16plugins/gdal_GRASS.so  gdal_GRASS.so

Qgis 1.4 ora carica senza problemi  i raster di  GRASS:

From Grass, gis e dintorni





lunedì 21 dicembre 2009

GPS Garmin in Windows: i due (+1) indispensabili

Dare gli strumenti indispensabili, facili da installare e liberi,  per  lavorare in Windows con apparecchi  GPS in... un'ottica GIS...: brevi istruzioni per l'uso rivolte a nuovi  utenti (entry level).

Map Source

Non essendo libero e facile da installare non dovrebbe essere menzionato in questo post; un accenno, però, occorre farlo. Considero Map Source utiile nel solo caso si debbano gestire mappe acquistate da  Garmin: pare che questo software  si  preoccupi in modo quasi maniacale di tutto ciò che riguarda autenticazione, sblocco e registrazione  delle mappe acquistate tralasciando di fornire all'utente un prodotto... utile e fruibile in tutte le occasioni di normale utilizzo (download  e upload di waypoitns, tracce, rotte).  Per quanto l'esperienza di chi scrive, l'utilità di Map Source inizia e finisce nel momento in cui, dopo una lunga e cervellotica  fase di autenticazione nella quale vengono richiesti 3-4 codici di attivazione, le mappe acquistate sono sbloccate e possono essere caricate nel GPS.
  

Quantum Gis

Qgis può essere scaricato da qui: consiglio di scegliere il file denominato QGIS-1.3.0-3-No-GrassSetup.exe (o un analogo file con un numero più recente di versione). Scaricarlo ed installarlo accettando la licenza ed ogni altra preimpostazione (non è necessario scaricare i set di dati di esempio proposti durante l'installazione).

Una guida introduttiva a  Qgis può essere trovata qui.

Il plugin GPS di Qgis sfrutta gpsbabel, un altro programma  che deve essere presente  (installato) nel sistema.

Gpsbabel

A giudizio di chi scrive, gpsbabel è uno tra i migliori software per lavorare con i  GPS.
Scaricare da qui il file denominato gpsbabel-1.3.6.zip (o un analogo file con un numero più recente di versione). Per l'installazione:
  1. creare una cartella sull'Hard disk (ad esempio: C:\Programmi\gpsbabel)
  2. scomapttare  il file   gpsbabel-1.3.6.zip nella cartella appena creata; alla fine di questa operazione il contenuto di quest'ultima dovrebbe presentarsi così:  
    From gps
  3. gpsbabel ha due eseguibili: GPSBabelGUI.exe  che è un'interaccia grafica ; gpsbabel.exe che è il vero e proprio motore, eseguibile all'interno della finestra MS-DOS di Windows ;
  4. gpsbabel , pur essendo installato, può essere - al momento - eseguito solo dalla cartella ove sono collocati i files (nell'esempio: C:\Programmi\gpsbabel). Per rendere disponibili le sue funzionalità all'esterno di questa cartella e per tutti  gli altri programmi che lo utilizzano (ad es.  qgis), occorre  inserire il percorso degli eseguibili di gpsbabel nella variabile di sistema path;  per far ciò (in Windows XP):
    Impostazioni->pannello dai controllo->sistema->avanzate->variabili di ambiente:  nella casella  Variabili di sistema selezionare  Path e modificarlo, aggiungendo alla linea il percorso ove sono stati installati gli eseguibili di gpsbabel (nell'esempio: C:\Programmi\gpsbabel; ricordarsi di separare i diversi percorsi con ";") 
  5. a questo punto gpsbabel dovrebbe funzionare...
    • attraverso la riga di comando della finestra MS-DOS di Windows

    • attraverso l'interfaccia grafica di gpsbabel  (doppio click su GPSBabelGUI.exe nella cartella C:\Programmi\gpsbabel oppure creare un link sul desktop) 
      From gps
    • attraverso il plugin GPS di qgis




sabato 19 dicembre 2009

GRASS, gpsd, l'MTK e il bluetooth

Nel precedente post si è visto come connetere l'MTK 747 Transystem a  gpsd via bluetooth in Ubuntu 9.10.
Ecco un semplice memo di come utilizzare la connessione in real-time in GRASS.

Lanciare GRASS in un terminale ed entrare in una location con sistema di riferimento cartografico noto; dal terminale di GRASS lanciare gpsd e quindi gpspipe:
From Grass, gis e dintorni

Ora, con semplici comandi della shell è possibile ottenere in real-time diverse elaborazioni; per esempio con:
gpspipe -w -n 30 | grep 'GGA' | awk 'BEGIN{FS=" "}{print $5,$4}' |\
 m.proj -i --q | awk '{sx+=$1;sy+=$2}END{ print " posizione media: ", sx/NR, sy/NR " samples: ", NR }'

si riproiettano (m.proj) le coordinate (Latitudine e Longitudine) inviate dal ricevitore nel sistema cartografico di riferimento della location corrente (nell'esempio: epsg: 32632),  viene calcolata  la posizione media del ricevitore collegato al pc e si stampa a video il risultato.

Complicando ulteriormente l'esempio, è possibile editare in real-time, catturando il punto in un vettoriale:
From Grass, gis e dintorni

Per le frasi  NMEA  è possibile consultare questo link.




MTK, GPSD via bluetooth

Brevi note per connettere l'MTK 747 Transystem a  gpsd via bluetooth in Ubuntu 9.10.

Ho una chiavetta USB SURECOM per le connessioni bluetooth che viene vista dal sistema così:
lsusb
[...]
Bus 003 Device 003: ID 0a12:0001 Cambridge Silicon Radio, Ltd Bluetooth Dongle (HCI mode)
[...]

Accendere l'MTK 747 Transystem in posizione  NAV ; provare  se la perifirica Bluetooth viene riconosciuta:
hcitool scan
Scanning ...
    00:1C:88:00:CC:9C    iBT-GPS

La periferica viene riconosciuta ed associata al codice 00:1C:88:00:CC:9C.

Configuro la  connessione:
paolo@laptop:~$ sudo gedit /etc/bluetooth/rfcomm.conf

inserisco questi settaggi:
rfcomm0 {
#    # Automatically bind the device at startup
     bind yes;
#    # Bluetooth address of the device
     device 00:1C:88:00:CC:9C;
#    # RFCOMM channel for the connection
    channel    1;
#    # Description of the connection
    comment "iBT-GPS";
}

Salvo e chiudo. N.B.: non è detto che debba per forza essere rfcomm0; se si hanno altri dispositivi bluetooth associare il dispositivo ad un'altra connessione (rfcomm1, rfcomm2,  ecc.)

Quindi:
paolo@laptop:~$ rfcomm connect 0

oppure:
rfcomm connect 0 00:1C:88:00:CC:9C

se tutto è ok si avrà la seguente risposta:
Connected /dev/rfcomm0 to 00:1C:88:00:CC:9C on channel 1
Press CTRL-C for hangup

Aprire un altro terminale e digitare:
paolo@laptop:~$ gpsd /dev/rfcomm0
paolo@laptop:~$ xgps

From Grass, gis e dintorni

La disconnessione può essere effettuata dal gestore Bluetooth di Gnome:

From Grass, gis e dintorni

oppure premendo CONTROL-C nel primo terminale.
Buona cosa, inoltre, uccidere il processo gpsd
gpsd killall



 

venerdì 18 dicembre 2009

GoogleEarth e files kml

Google Earth in Linux: problema di visualizzazione non corretta dei file kml importati.

Dopo aver scaricato ed installato Google Earth in Ubuntu ho notato che i files kml importati non venivano posizionati correttamente; ho risolto il problema grazie a questo post; in pratica è sufficiente modificare il file googleearth nella cartella /home/nomeutente/google-earth (nel caso di installazione in locale, nella home utente) come segue:

LD_LIBRARY_PATH=.:${GOOGLEEARTH_DATA_PATH}:${LD_LIBRARY_PATH}
export LD_LIBRARY_PATH
export LC_ALL=us_US.UTF-8
export LANG=it_IT.UTF-8

Ecco un esempio, con una traccia denominata track_15_12_2009.gpx_trk.gpx:

step 1: conversione da gpx a kml con gpsbabel:
paolo@laptop:~$ gpsbabel -t -i gpx -f track_15_12_2009_trk.gpx -o kml -F track.kml

step 2: aprire Google Earth e caricare la traccia:

From Grass, gis e dintorni

From Grass, gis e dintorni






mercoledì 9 dicembre 2009

Codici epsg per l'Italia Nord-Occidentale.

Codici epsg: memo.

Ecco i codici di più frequente uso nell'Italia Nord-Occidentale (Piemonte, Liguria, Valle d'Aosta).

Codice
Descrizione
wkt
4806 Roma1940, in gradi.
wkt
26591
Sistema nazionale Gauss-Boaga, fuso Ovest. Datum Roma1940.  Falsa origine delle ascisse: 1.500 km. Deprecato (usare 3003).
wkt
3003
Idem c.s.
wkt
23032
European 1950 datum (ED50). Grigliato UTM, fuso 32 Nord.
wkt
4326
WGS84 datum. Lat/Long.
wkt
32632
WGS84 datum. Grigliato UTM, fuso 32 Nord. wkt

I parametri per le trasformazioni possono essere ottenuti con il comando di  GRASS g.proj.

lunedì 9 novembre 2009

Garmin GPSmap 60 Cx e GPSBabel

Usare il Garmin 60Cx con GPSBabel. Seguire le guide di OpenStreetMap e di GPSBabel.
In pratica, creare il file 51-garmin.rules:

paolo@laptop:~$ cat /etc/udev/rules.d/51-garmin.rules
# # http://wiki.openstreetmap.org/wiki/USB_Garmin_on_GNU/Linux
# # gps garmin con porta usb senza permessi
ATTRS{idVendor}=="091e", ATTRS{idProduct}=="0003", MODE="0660", GROUP="plugdev"
Inserire il modulo garmin_gps in blacklist:
paolo@laptop:~$ cat /etc/modprobe.d/blacklist
#prevent garmin gops loaded so generic USB can be used instead
blacklist garmin_gps
A questo punto il ricevitore GPS viene riconosciuto come periferica usb:.
Per testare il tutto caricando una traccia dal pc al ricevitore:

  1. cancellare ogni traccia dal ricevitore;
  2. collegare il ricevitore;
  3. caricare una traccia dal pc al ricevitore con:
    gpsbabel -t -i gpx -f miatraccia.gpx -o garmin -F usb:
    

  4. la traccia viene (dovrebbe!) essere caricata, ma il ricevitore
    non la nomina;
  5. nel menù tracce salvare la traccia dandole un nome

sabato 25 aprile 2009

Point in polygon (terza parte)

Joins spaziali e SQL

In questa sezione vengono affrontati gli aspetti principali del point in polygon maggiormente legati all'SQL.

Preparazione dei dati e prime verifiche

I dati necessari per l'esercitazione sono quelli ottenuti nell'esercitazione precedente; riassiumiamo:
  1. la location simpleloc con il mapset pointInPoly;
  2. un vettoriale poligonale polygons con attributi cat (la categoria dell'area), label (l'etichetta identificativa di ciascuna area) e GRASSRGB (il colore dell'area secondo la codifica RGB);
  3. un vettoriale di punti mypt collegato, sul layer 1, alla tabella points_code con i campi cat e code (quest'ultimo codifica la posizione del punto rispetto agli elementi dei poligoni).

Merita porre l'attenzione sull'ultimo punto: un vettoriale è stato collegato ad una tabella con un nome totalmente diverso. Spesso si è indotti a pensare che ciò non sia possibile (a causa forse dell'abitudine ad usare formati di file quali lo shapefile) ma tutti i moduli di GRASS  usati in questo esercizio funzionano benissimo anche in questo caso. Vediamo un esempio:

# # # vector and table with a different name
# # # mypt connected to points_table
table="points_code"
column="polyCat"
# # # if column not exists will be created
if ! checkColumn "$table" "$column" ; then
echo "ALTER TABLE $table ADD COLUMN $column integer" | db.execute
fi
# # upload polygons cat to table
v.what.vect vect=mypt col=polyCat \
qvector=polygons qcolumn=cat
# # # how many rows in table?
db.select sql="SELECT count(cat) FROM $table"
# # # how many points overlap polygons?
db.select sql="SELECT count(cat) FROM $table \
WHERE $column NOT NULL"
# # # how many points don't overlap polygons?
db.select sql="SELECT count(cat) FROM $table \
WHERE $column IS NULL"
# # # how many points overlaps boundaries, nodes, vertices, centroids, ecc.
db.select sql="SELECT code,count(cat) FROM $table \
WHERE $column NOT NULL GROUP BY code"

Al solito vengono lasciati i commenti.

Una nota di spiegazione meritano le funzioni checkColumn e checkTable (usata dopo). Tutto il codice riportato è contenuto in una script che viene testato ogni qualvolta si apportano modifiche; il tentativo di creare una colonna (campo) già esistente fornisce un errore e lo stesso dicasi per il tentativo di creare una tabella già presente nel database. Ripetendo quindi l'esecuzione dello script senza un preventivo controllo sull'esistenza o meno di tabelle e colonne che si volgiono creare fornirebbe numerosi errori. Per evitare questo è possibile usare queste funzioni:

checkTable()
{
# true=0
# false=1
table=$1
db.tables -p | grep '\<'$table'\>' 2>/dev/null
if [ $? -eq 0 ]; then
exists=0
else
g.message -w message="Table $table doesn't exist"
exists=1
fi
return $exists
}
checkColumn()
{
# true=0
# false=1
local table=$1
local column=$2
if checkTable $table; then
db.describe -ct $table | grep '\<'$column'\>' 2>/dev/null
if [ $? -eq 0 ]; then
exists=0
else
g.message -w message="Column $column doesn't exist"
exists=1
fi
return $exists
else
exit 1
fi
}

È chiaro che chi non ha dimestichezza con lo scripting di shell può lo stesso fare questi esercizi, senza dover per forza capire il meccanismo di funzionamento delle funzioni sopra riportate! Nell'esempio di prima basta omettere le righe che richiamano checkColumn:

# # # if column not exists will be created
echo "ALTER TABLE points_code ADD COLUMN polyCat integer" | db.execute

...ricordandosi di cancellare la colonna polyCat con db.dropcol prima di ripetere l'esercizio, per evitare errori del tipo:

GRASS 6.5.svn (simpleloc):~/script_grass > echo "ALTER TABLE points_code ADD COLUMN polyCat integer" | db.execute
DBMI-SQLite driver error:
Error in sqlite3_prepare():
duplicate column name: polyCat

ERROR: Error while executing: 'ALTER TABLE points_code ADD COLUMN polyCat
integer
'

Riprendendo l'esercitazione, se proprio si vuole avere vettoriale e tabella con lo stesso nome si può procedere semplicemente così, copiando il vettoriale con un nuovo nome:

 # # vector and table with the same name
g.copy vect=mypt,mypt2 --overwrite
db.describe -tc mypt2

oppure esportando e reimportando il vettoriale con v.out.ogr e v.in.ogr:

# # # vector and table with the same name
v.out.ogr in=mypt layer=1 type=point dsn=. \
olayer=tmp_pippo format=ESRI_Shapefile -c
v.in.ogr dsn=. output=mypt2 layer=tmp_pippo --o

Si noti che con quest'ultimo metodo si ottiene una tabella leggermente diersa rispetto all'originaria (usare db.describe... viene inserito un campo cat_; il campo polyCat viene trasformato da integer a double precision).

Per esercitarsi appieno con SQL, si crei una tabella che memorizzi la desrizione della codifica dei punti così come definita nell'esercitazione precedente

# # # create a table for points descriptions
table="points_descr"
# # # if table exists will be cancelled
if checkTable $table ; then
db.droptable -f $table
fi
echo "CREATE TABLE $table (cat integer, txtDescr varchar(25) )" | db.execute
# # insert data in table with loop and here document
while read data
do
fcat=$(echo $data | cut -d"|" -f1)
ftxtDescr=$(echo $data | cut -d"|" -f2)
echo "INSERT INTO $table (cat,txtDescr) \
values ("$fcat",'"$ftxtDescr"')" | db.execute
done<<EOF
1|out (completly out!)
2|in (but not on boundaries nor centroids)
3|on_centroids
4|on_nodes (but not on vertices)
5|on_vertices (but not on nodes)
6|on_boundaries (but not on vertices nor nodes)
EOF
db.select $table

v.distance semplificato: v.what.vect

Si è visto nell'esercitazione precedente che v.distance permetto il trasferimento della categoria o degli attributi dai poligoni ai punti; v.what.vect non è altro che uno script per semplificare v.distance in questo genere di operazioni.

Per esempio si può aggiornare il campo polyCat dei punti con la categoria dei poligoni o trasferire un attributo (in questo caso il campo GRASSRGB):

# # # upload polygons cat to points attribute
v.what.vect vect=mypt2 col=polyCat qvector=polygons qcolumn=cat
# # # coloro i punti con lo stesso colore del poligono
# # # nel quale sono contenuti
if ! checkColumn "mypt2" "GRASSRGB" ; then
v.db.addcol mypt2 col="GRASSRGB VARCHAR(11)"
fi
# # # upload polygons colors to points
v.what.vect vect=mypt2 col=GRASSRGB qvector=polygons qcolumn=GRASSRGB

Queste operazioni sono del tutto analoghe a quelle già viste in precedenza con v.distance.

Punti, poligoni e SQL

Si propone qui una semplice carrellata di esempi per dare un'idea generale sull'uso dell'SQL collegato con lo spatial join e il point in polygon; si rimanda agli ottimi tutorial presenti in rete per gli approfondimenti riguardanti l'SQL.

Verifiche sul numero dei punti

# # # how many points?
db.select sql="SELECT count(cat) FROM mypt2"
# # # how many points with polyCat not empty?
# # # how many points overlap polygons?
db.select sql="SELECT count(cat) FROM mypt2 \
WHERE polyCat NOT NULL"
# # # as above-mentioned but group by code
db.select sql="SELECT code,count(cat) FROM mypt2 \
WHERE polyCat NOT NULL GROUP BY code"
# # # how many points polyCat=NULL? (=points outside areas)
db.select sql="SELECT count(cat) FROM mypt2 \
WHERE polyCat IS NULL"

Inner Join (equal join)

Nell'Inner join vengono estratti solo i record di entrambe le tabelle nei quali il campo chiave coincide.

Nel caso di mypt2 e polygons il campo chiave è rispettivamente polyCat e cat.

db.select sql="SELECT mypt2.cat AS pointCat, \
polygons.cat, polygons.label AS polyLabel, \
polygons.GRASSRGB AS polyColor \
FROM mypt2 INNER JOIN polygons \
ON mypt2.polyCat=polygons.cat"

L'espressione con i campi chiave mypt2.polyCat=polygons.cat viene valutata true solo quando i campi, provenienti dalle due tabelle, sono uguali.

È facile verificare che è indifferente porre la clausula WHERE mypt.polyCat NOT NULL perchè nell'inner join vengono estratti solo i record di entrambe le tabelle nei quali il campo chiave coincide:

db.select sql="SELECT mypt2.cat AS pointCat, \
polygons.cat, polygons.label AS polyLabel, \
polygons.GRASSRGB AS polyColor \
FROM mypt2 INNER JOIN polygons \
ON mypt2.polyCat=polygons.cat \
WHERE mypt2.cat NOT NULL"

Si può usare anche un'altra notazione, che porta a risultati identici:

db.select sql="SELECT a.cat AS pointCat, \
b.cat, b.label AS polyLabel, \
b.GRASSRGB AS polyColor \
FROM mypt2 a, polygons b \
WHERE a.polyCat=b.cat"

È possibile ottenere i punti raggruppati per codice con la descrizione contenuta nella tabella points_descr:

db.select sql="SELECT mypt2.code AS code, count(mypt2.code) AS NumberOfPoints, \
points_descr.txtDescr AS description \
FROM mypt2 INNER JOIN points_descr \
ON mypt2.code=points_descr.cat GROUP BY mypt2.code"

...o l'elenco completo dettagliato di tutti i punti

db.select sql="SELECT mypt2.cat, mypt2.code AS code, \
points_descr.txtDescr AS description \
FROM mypt2 INNER JOIN points_descr \
ON mypt2.code=points_descr.cat ORDER BY code"

Un problema abbastanza comune è il calcolo del numero di punti contenuti in ciascun poligono:

db.select sql="SELECT COUNT(mypt2.cat) AS totalPoint, \
polygons.label AS polyLabel FROM mypt2 \
INNER JOIN polygons ON \
mypt2.polyCat=polygons.cat GROUP BY polyCat"

Questo risutato può essere memorizzato in una tabella autonoma:

table="countPointsPerPoly"
# # # if table exists will be cancelled
if checkTable $table ; then
db.droptable -f $table
fi
echo "CREATE TABLE $table AS SELECT count(mypt2.cat) AS numpoint, \
polygons.cat AS refPoly FROM mypt2 \
INNER JOIN polygons \
ON mypt2.polyCat=polygons.cat \
GROUP BY polygons.cat" | db.execute

...oppure memorizzato in un campo nella tabella attributi dei poligoni:

if ! checkColumn "polygons" "numpoints" ; then
v.db.addcol polygons col="numpoints integer"
fi
#
for i in $(v.db.select polygons col=cat -c); do
echo "UPDATE polygons SET numpoints= \
(SELECT count(mypt2.cat) FROM mypt2 \
WHERE mypt2.polyCat=$i \
GROUP BY polyCat) WHERE polygons.cat=$i" | db.execute
done

L'aggiornamento del campo numpoints è realizzato attraverso un loop; per ogni cat della tabella dei poligoni viene eseguita la subquery rappresentata dalle istruzioni all'interno delle parentesi tonde che restituisce un elenco di valori. La successiva clausula WHERE polygons.cat=$i filtra i valori ottenuti dalla subquery imponendo la condizione polygons.cat=$i ottenendo un solo valore per ogni ciclo che viene inserito nel campo numpoints.

È possibile concatenare l'inner join; nell'esempio si ottiene, per ogni poligono il punto contenuto con il suo codice e la descrizione:

# # # double inner join
db.select sql="SELECT polygons.label AS polyLabel, mypt2.cat AS pointCat,\
mypt2.code as pointCode, points_descr.txtDescr AS pointDescr \
FROM polygons \
INNER JOIN mypt2 ON polygons.cat=mypt2.polyCat \
INNER JOIN points_descr ON mypt2.code=points_descr.cat"

Left join (left outer join)

Il risultato del left join include tutti i records della prima tabella e solo i records della seconda che rendono vera l'espressione con i campi chiave:

db.select sql="SELECT mypt2.cat AS pointCat, \
polygons.cat, polygons.label AS polyLabel, \
polygons.GRASSRGB AS polyColor \
FROM mypt2 LEFT JOIN polygons \
ON mypt2.polyCat=polygons.cat"

Attualmente (aprile 2009) il right join (right outer join) non è supportato dal driver Sqlite di Grass.

Analogamente a quanto fatto in precedenza, è possibile estrarre i dati dalle tre tabelle usando sia l'inner join che il left join:

db.select sql="select mypt2.cat AS pointCat, mypt2.code AS pointCode, \
points_descr.txtDescr AS pointDescr, polygons.label AS polyLabel \
FROM mypt2 \
INNER JOIN points_descr ON mypt2.code=points_descr.cat \
LEFT JOIN polygons ON mypt2.polyCat=polygons.cat"

Casi particolari

Si voglia per esempio estrarre, per ciascun poligono avente area maggiore ad un certo valore, il numero di punti in esso contenuti:

if ! checkColumn "polygons" "area"; then
v.db.addcol polygons col="area double precision"
fi
v.to.db polygons opt=area col="area"
db.select sql="SELECT count(mypt2.polyCat) AS numPoint, \
polygons.label as polyLabel, \
polygons.area AS polyArea FROM mypt2 \
INNER JOIN polygons ON mypt2.polyCat=polygons.cat \
GROUP BY polygons.cat HAVING polyArea>=7 ORDER BY polyArea DESC"

Oppure è possibile calcolare il rapporto di densità tra punti contenuti ed area del poligono contenitore:

db.select sql="SELECT count(mypt2.polyCat) AS numPoint, \
polygons.label as polyLabel, \
polygons.area AS polyArea, round(numpoints/area,2) AS density FROM mypt2 \
INNER JOIN polygons ON mypt2.polyCat=polygons.cat \
GROUP BY polygons.cat ORDER BY density DESC"

Conclusioni

In queste tre esercitazioni si è visto a grandi linee come GRASS  affronta il problema del point in polygon. I tools base sono v.select e v.distance; v.what.vect è uno script che semplifica l'uso di v.distance; l'SQL è stato usato per l'estrazione dei dati e per creare i joins. Merita ricordare che quando l'esigenza è di unire (merge) due tabelle, c'è a disposizione v.db.join: il suo funzionamento interno, al di lá dei diversi controlli , ricorda le operazioni usate per l'aggiornamento del campo numpoints attraverso un loop.

La conoscenza dell'SQL pare uno strumento indispensabile per l'utente GIS, in particolare quando si lavora con il modello vettoriale.

Tutto il codice riportato in questo esercizio può essere scaricato da qui


giovedì 16 aprile 2009

Tre anni di Gis Open Source

Per una serie di motivi e coincidenze, nel marzo 2006 decisi di passare dal Gis proprietario in Windows al Gis Open Source in Linux.

Dopo tre anni circa mi sento di poter affermare la grande soddisfazione che questa scelta ha portato.

Poche cose mi mancano veramente del Gis-Windows:

  1. l'ambiente layout di ArcView (dalla vers. 8.2 in poi); ps.map ha grandi potenzialità ma allestimenti raffinati richiedono un ambiente dedicato WYSIWYG o WYSIWYM (What You See Is What You Mean): in ogni caso allestendo una mappa di un certo livello qualcosa (e non solo il codice) si deve vedere in tempo reale man mano che si lavora e non solo dopo compilazione. gvSIG offre un ambiente del tipo ArcView  vecchia maniera (vers. 3.2) ma come sarebbe bello avere tutto in GRASS!
  2. l'editor vettoriale di Terranova ShArc: in 3 anni GRASS  ha fatto passi da gigante dal punto di vista dell'editing vettoriale, gvSIG offre un ottimo ambiente ma non topologico; purtroppo funzionalità del tipo inseguimento linee del raster o rimozione interattiva dei dangles sono ancora lontane; nell'editing massivo di grandi quantità di elementi geometrici contano molto l'interfaccia, gli acceleratori da tastiera, la velocità di ridisegno, forma e dimensione del puntatore, ecc.
  3. la gestione delle etichette (labels) di ArcView (dalla vers. 8.2); tra le nuove caratteristiche incluse in gvSIG (vers. 1.9) vi è un motore per il labeling avanzato; speriamo in bene...


giovedì 26 marzo 2009

Point in polygon (seconda parte)

Classificazione di punti

Prima di affrontare gli aspetti del point in polygon maggiormente legati all'SQL è bene riassumere ed approfondire l'uso di v.select e v.distance.

L'esercizio che ci si propone di risolvere è il seguente: dati una copertura poligonale e un vettoriale di punti, attribuire i seguenti codici numerici:

  • 1, ai punti completamente esterni alle aree;
  • 2, ai punti completamente all'interno delle aree ma non sui boundaries né sui centroidi;
  • 3, ai punti posti sui centroidi;
  • 4, ai punti posti sui nodi ma non sui vertici dei boundaries;
  • 5, ai punti posti sui vertici ma non sui nodi;
  • 6, ai punti posti sui boundaries (ma non sui vertici né sui nodi).

In pratica si vuole una classificazione dei punti in funzione della posizione rispetto agli elementi dei poligoni (centroidi, boundaries, nodi, vertici).

La soluzione che si propone non è, probabilmente, ottimale ma... chi scrive non è riuscito a trovarne di migliori. Ben vengano i suggerimenti e le correzioni!

In generale la procedura prevede:

  1. l'attribuzione di una categoria corrispondente al codice voluto agli elementi dei poligoni (centroidi, nodi, vertici, boundaries),
  2. il trasferimento della categoria dagli elementi dei polgioni agli attributi dei punti mediante v.distance,
  3. l'estrazione dei punti con il codice desiderato (v.extract) in un nuovo vettoriale;
  4. la cancellazione dei punti nel vettoriale originario: vengono cancellati solo i punti di volta in volta già estratti con v.extract

Preparazione dei dati

Al solito, viene usata la location simpleloc. Entrare nel mapset pointInPoly creato in precedenza.
db.connect driver=sqlite database='$GISDBASE/$LOCATION_NAME/$MAPSET/pointInpoly.db'
g.copy vect=points,mypt,polyTopology,polygons --overwrite
# # # add table to mypt
v.db.addtable mypt
v.db.addcol mypt col="code integer"
g.copy vect=mypt,mypt_tmp --overwrite

Tutto il codice riportato in questo esercizio può essere scaricato da qui

Punti (solo) sui centroidi

Si attribuisce con v.category una categoria ai centroidi corrispondente al codice voluto (3); si trasferisce la categoria dai centroidi al campo code di mypt_tmp; si estraggono da mypt_tmp tutti i punti con il codice pari a 3 creando il vettoriale on_centroids; si cancellano in mypt_tmp tutti i punti con codice 3.

v.category in=polygons out=polycode \
type=centroid layer=2 cat=3 step=0 --overwrite
v.distance from=mypt_tmp to=polycode \
to_type=centroid to_layer=2 upload=cat \
col=code dmax=0
v.extract in=mypt_tmp out=on_centroids where="code=3" --overwrite
v.edit map=mypt_tmp tool=delete where="code=3"
g.remove vect=polycode

Punti sui nodi (ma non sui vertici)

Si assume la definizione di nodo fornita in questa esercitazione. Occorre estrarre i nodi con v.to.points; in questo modo però, il punto finale di un boundaries può coincidere con il punto iniziale del boundaries successivo. Si rimuovo allora i nodi duplicati con v.clean. Si attribuisce poi la categoria 4 ai nodi estratti e, così come fatto per i centroidi, si trasferisce la categoria dai nodi al campo code di mypt_tmp. Si estraggono i punti in un nuovo file on_nodes e si rimuovono nel vettoriale my_tmp.

v.to.points -n in=polygons out=multiNode \
type=boundary --overwrite
# # # remove duplicate nodes
v.clean in=multiNode out=nodeSimple \
type=point tool=rmdupl --overwrite
# # # get a unique cat for each node
v.category nodeSimple out=node \
layer=3 opt=add cat=4 step=0 --overwrite
# # # upload distance fron nearest node
v.distance from=mypt_tmp to=node to_layer=3 \
to_type=point upload=cat \
column=code dmax=0
v.extract in=mypt_tmp out=on_nodes where="code=4" --overwrite
v.edit map=mypt_tmp tool=delete where="code=4"
g.remove vect=multiNode,nodeSimple,node

Punti (solo) sui centroidi Punti sui nodi

Punti sui vertici (ma non sui nodi)

Le operazioni sono del tutto simili a quelle già svolte per i nodi.

v.to.points -v in=polygons out=multiVertex \
type=boundary --overwrite
# # # remove duplicate nodes
v.clean in=multiVertex out=vertexSimple \
type=point tool=rmdupl --overwrite
# # # set unique cat for each vertx
v.category vertexSimple out=vertex \
layer=3 opt=add cat=5 step=0 --overwrite
# # # upload distance fron nearest node
v.distance from=mypt_tmp to=vertex to_layer=3 \
to_type=point upload=cat \
column=code dmax=0
v.extract in=mypt_tmp out=on_vertices where="code=5" --overwrite
v.edit map=mypt_tmp tool=delete where="code=5"
g.remove vect=multiVertex,vertexSimple,vertex

Punti sui boundaries (ma non sui nodi né sui vertici)

A questo punto è doveroso spiegare il perchè si sono cancellati i punti da mypt_tmp ogni volta che si sono estratti in un nuovo file: se non lo si facesse sarebbe impossibile riuscire ad estrarre i punti posti sui boundaries (o meglio, lungo i boundaries) senza estrarre anche i punti posti sui vertici e i sui nodi dei boundaries medesimi (provare per credere). Come si è detto questa procedura è lungi dall'essere ottimale (sig!)...

v.category in=polygons type=boundary \
out=arcs layer=2 opt=add cat=6 step=0 --overwrite
v.distance from=mypt_tmp to=arcs to_layer=2 \
to_type=boundary to_layer=2 upload=cat \
col=code dmax=0
v.extract in=mypt_tmp out=on_boundaries where="code=6" --overwrite
v.edit map=mypt_tmp tool=delete where="code=6"
g.remove vect=arcs

Punti (solo) sui vertici Punti sui boundaries

Punti esterni alle aree

Con v.select è, a questo punto, banale estrarre i punti esterni alle aree (codice 1):
v.select -r ainput=mypt_tmp \
binput=polygons \
op=overlap out=out_area --overwrite
v.db.update out_area col=code value=1

Punti nelle aree (ma non sui centroidi né sui boundaries)

Si usa v.select con una certa sicurezza: avendo cancellato via via i punti da mypt_tmp non si dovrebbe (ancora una volta: la procedura qui proposta è lungi dall'essere ottimale) correre il rischio di estrarre punti non corretti.

v.select ainput=mypt_tmp \
binput=polygons \
op=overlap out=in_area --overwrite
v.db.update in_area col=code value=2

Punti esterni alle aree Punti nelle aree (ma non sui centroidi, nodi,vertici, boundaries)

Creazione della tabella dei punti

Verrebbe da usare, a questo punto, v.patch:

v.patch -e in=in_area,out_area,on_boundaries,\
on_nodes,on_vertices,on_centroids \
out=points_code_vpatch --overwrite
v.db.select points_code_vpatch

ma ecco il risultato nella tabella attributi:

cat|code
12|2
[...]
18|2
36|2
38|1
[...]
47|1
81|6
[...]
87|6
110|4
[...]
152|5
171|3
[...]
174|3

Le categorie originarie di mypt vengono perse; se si vogliono preservare le categorie, occorre creare una tabella unica partendo dagli attributi dei diversi files di volta in volta estratti (on_centroids, on_nodes, ecc.) e connetterla (v.db.connect) sul layer 1 sovrascrivendo la connessione esistente mypt:

# # # test if table exists 
# # # from v.db.addtable script
table="points_code"
database=`db.connect -p | grep '^database' | cut -d':' -f2`
driver=`db.connect -p | grep '^driver' | cut -d':' -f2`
db.tables database="$database" driver="$driver" 2> /dev/null | grep "^$table$"
# # # if table exists it will be deleted
if [ $? -eq 0 ]; then
db.droptable -f $table
fi

echo "CREATE TABLE $table AS \
SELECT * FROM on_boundaries \
UNION \
SELECT * FROM on_nodes \
UNION \
SELECT * FROM on_vertices \
UNION \
SELECT * FROM in_area \
UNION \
SELECT * FROM out_area \
UNION \
SELECT * FROM on_centroids" | db.execute
v.db.connect map=mypt table=points_code key=cat layer=1 -o
v.db.select mypt layer=1

In questo modo la tabella points_code viene connessa al vettoriale mypt sul layer 1 (non è detto che il nome del vettoriale e il nome della tabella ad esso collegata debbano proprio essere uguali!).

La prima parte del codice serve solamente per testare se la tabella points_code esiste già; in caso positivo viene sovrascritta con la nuova creata dal codice SQL riportato.

Come si vede dal risultato di v.db.select la tabella sul layer 2 di mypt contiene due campi: la cat che corrisponde ovviamente a quella del vettoriale originario (points) e il code contenente la codifica basata sulla localizzazione del punto rispetto agli elementi della copertura poligonale.

Un'ultima semplice e grossolana verifica sul numero dei punti nei diversi files intermedi rispetto al vettoriale contenente tutti i punti (mypt):

or vectorpoint in "on_vertices" \
"on_boundaries" \
"on_nodes" \
"on_centroids" \
"out_area" \
"in_area" \
"mypt"; do
echo "$vectorpoint $(v.info -t $vectorpoint | grep points | cut -d"=" -f2)"
done

on_vertices 7
on_boundaries 5
on_nodes 5
on_centroids 4
out_area 10
in_area 8
mypt 39

La somma dei punti nei diversi vettoriali (on_centroids, on_nodes, ecc.) corrisponde al numero dei punti contenuti in mypt.


domenica 1 marzo 2009

Point in polygon (prima parte)

Obbiettivi

Trovare se una determinata area contiene/non contiene un dato punto è tra i più classici dei problemi GIS.

GRASS  risolve il problema con alcuni tools:

v.select, v.distance e v.what.vect.

Questa esercitazione introduttiva affronta il problema del point in polygon, ed il relativo passaggio degli attributi dagli elementi poligonali ai punti contenuti (o a quelii più vicini, sovrapposti, ecc): lo spatial join.

Set di dati per l'esercitazione

Al solito l'esercitazione sarà condotta con la location simpleloc, la location creata apposta per fini didattici.

È necessario creare un nuovo mapset pointInPoly ed accedervi.

Per poter usufruire appieno dell'SQL (per la seconda parte di questa esercitazione) occorre connettere il mapset ad un database SQLite e copiare i due vettoriali dal mapset PERMANENT al mapset pointInPoly

db.connect driver=sqlite \
database='$GISDBASE/$LOCATION_NAME/$MAPSET/pointInpoly.db'
g.copy vect=points,mypt,polyTopology,polygons --overwrite
v.db.addtable mypt

Le tabelle dei due vettoriali sono ora memorizzate nel database SQLite denominato pointInpoly.db.

Il vettoriale di punti usato è stato creato apposta per consentire di avere tutte le combinazioni possibili (punto in area, fuori area, su taluni centroidi, su alcuni nodi, su alcuni vertici, sui boundaries, molto vicino ai centroidi, molto vicino ai boundaries, ecc.).

Selezione semplice di punti (in,out,on...)

Dati due vettoriali, v.select consente di selezionare ed estrarre in un nuovo vettoriale gli elementi del primo (ainput) sovrapposti al secondo (binput).

Ecco in sequenza le diverse operazioni; vengono lasciati i commenti che servono da spiegazione:

# # # select all points into polygons
v.select ainput=mypt \
binput=polygons btype=area \
op=overlap out=inarea --overwrite

# # # select all points outside polygons
v.select -r ainput=mypt \
binput=polygons \
op=overlap out=outarea --overwrite

# # # select all points on centroid
v.select ainput=mypt \
binput=polygons btype=centroid\
op=overlap out=oncentroids --overwrite

# # #select all points on boundary
v.select ainput=mypt \
binput=polygons btype=boundary \
op=overlap out=onboundaries --overwrite

Al momento overlap è l'unico operatore di v.select ma nella roabmap di GRASS si legge anche di operatori di contenimento e adiacenza.

Le mappe che seguono illustrano i risultati ottenuti con v.select e sono ottenute in modo automatico tramite lo script che può essere scaricato da questo link

Punti all'interno delle aree Punti sui boundaries
Punti sui centroidi Punti non sovrapposti alle aree

Selezione di punti con passaggio di attributi

Il modulo v.distance spaventa un po': la sua sintassi è davvero complicata (almeno per me!) ma è ampiamente giustificata dall'elevato numero di analisi possibili. Il manuale riporta che v.distance trova, per ogni elemento del vettoriale from, il più vicino elemento nel vettoriale to: decisamente diminutivo! In realtà v.distance è un formidabile tool per realizzare joins spaziali.

Una nota generale: Si impara ad usare v.distance usandolo..., come per moltissimi altri comandi di GRASS  che a prima vista spaventano per l'elevato numero di opzioni.

Upload della cat del più vicino poligono.

Il primo e più semplice dei problemi è il seguente: dato un vettoriale costituito da punti e un altro vettoriale costituito da poligoni, trovare per ciascun punto il più vicino poligono e caricare la categoria di quest'ultimo nella tabella degli attributi dei punti.

Si crea quindi una sorta di legame spaziale (=join) tra ciascun punto e il poligono più vicino.

# # # upload category of the NEAREST polygons
v.db.addcol mypt col="polyCatNearest integer"
# # # dmax=-1 (default, no limit): all points
v.distance from=mypt to=polygons \
upload=cat col=polyCatNearest dmax=-1

Settando il parametro dmax al valore di -1 non si pongono limiti alla ricerca del più vicino poligono; pertanto, anche gli attributi dei punti esterni (non sovrapposti) alle aree vengono aggiornati con la cat del poligono più vicino; se si setta il valore pari a 0, saranno trovati solo i punti sovrapposti ai poligoni.

v.db.addcol mypt col="polyCatOverlap integer"
# # # dmax=0: only points that OVERLAP polygons
v.distance from=mypt to=polygons \
upload=cat col=polyCatOverlap dmax=0

Si osservi ora la tabella degli attributi dei punti:

v.db.select mypt col=cat,polyCatNearest,polyCatOverlap

I punti 6 e 7 non sembrano essere stati presi in considerazione: la cat del poligono più vicino, per questi due punti, è nulla. I punti 6 e 7 sono in un'isola (area senza centroide e con categoria nulla): in questo caso il valore caricato da v.distance sugli attributi dei punti è nullo.

Se si vuole realizzare un upload di tutte le cat dei centroidi dei poligoni sugli attributi dei punti occorre settare il parametro to_type a centroid:

v.db.addcol mypt col="centroidCat integer"
v.distance from=mypt to=polygons to_type=centroid \
upload=cat col=centroidCat dmax=-1 \
out=distFromCentroid --overwrite

Si osservi che ora anche ai punti 6 e 7 è stata caricata la cat del centroide più vicino:

v.db.select mypt col=cat,polyCatNearest,polyCatOverlap,centroidCat \
where="cat=6 OR cat=7"

Upload della cat del centroide posto ad una distanza minima data.

Sino ad ora si è usato v.distance per trovare gli oggetti (poligoni o centroidi) più vicini; è pure possibile aggiornare la tabella attributi dei punti con la cat dei centroidi trovati ad una distanza minima data:

v.db.addcol mypt col="centroidCatDmin integer"
v.distance from=mypt to=polygons to_type=centroid \
upload=cat col=centroidCatDmin \
dmin=3 out=distFormCentroidDmin --overwrite

In questo caso v.distance cerca i centroidi distanti almeno 3 unità di mappa da ciascun punto.

Distanza dai centroidi Centroidi distanti almeno 3 unità di mappa dai punti

Upload dell'area dei poligoni sugli attributi dei punti sovrapposti.

Per creare una relazione spaziale tra punti e poligoni, è di per sè sufficiente caricare la cat dei poligoni sugli attributi dei punti con i metodi appena visti; come si vedrà in seguito, in tale maniera è possibile interrogare la tabella dei poligoni in funzione della posizione dei punti. Non sempre però è questo il risultato desiderato: si potrebbe ad esempio richiedere un passaggio diretto di taluni attributi dei poligono sulla tabella dei punti.

Quale esempio di upload di attributi si può calcolare l'area dei poligoni e caricarla nella tabella dei punti. In questo esempio verranno presi in considerazione solo i punti sovrapposti ai poligoni:

# # # upload all polygons (areas) attributes
# # # required fields (columns) in "from" vector
# # # must be present before run v.distance
v.db.addcol mypt col="polyArea double precision"
# # # add column for area polygons
v.db.addcol polygons col="area double precision"
# # # calculate polygons area size
v.to.db polygons opt=area col=area
# # # upload polygons attributes to OVERLAP points
v.distance from=mypt to=polygons \
upload=to_attr col=polyArea \
to_col=area dmax=0
v.db.select mypt

Upload della distanza fra punti e centroidi

v.distance è un calcolatore di distanze fra oggetti; ecco due esempi dell'upload delle distanze tra:

  • i punti e i centroidi dei poligoni;
  • i punti e i boundaries dei poligoni
# # # upload distance from the nearest centroid
v.db.addcol mypt col="centroidDist double precision"
v.distance from=mypt to=polygons to_type=centroid \
upload=dist col=centroidDist \
dmax=-1
# # # upload distance from the nearest boundary (ALL points)
v.db.addcol mypt col="boundaryDist double precision"
v.distance from=mypt to=polygons to_type=boundary \
upload=dist col=boundaryDist \
dmax=-1 out=distFromBondary --overwrite

v.db.select mypt col=cat,centroidDist,boundaryDist

Come si sarà certamente notato è sempre possibile materializzare la distanza fra gli oggetti, settando il parametro opzionale out con il nome del file vettoriale desiderato.

Distanze tra punti e boundaries Distanze tra punti e nodi

Casi particolari

Gli esercizi precedenti hanno consentito di apprendere le operazioni basilari con v.distance.

Ecco alcuni casi meno usuali.

Upload della distanza, della cat univoca e delle coordinate del più vicino nodo

Si assume la definizione di nodo fornita in questa esercitazione. Data quindi una copertura poligonale, estrarre i nodi rimuovendo ogni sovrapposizione, identificarli univocamente e applicare v.distance in modo da ottenere l'upload, sulla tabella attributi dei punti, di:

  1. distanza fra nodo e più vicino punto;
  2. cat univoca del nodo;
  3. coordinate del nodo.
v.db.addcol mypt col="nodeDist double precision, \
nodeCat integer, \
node_x double precision, \
node_y double precision"
# # # extract nodes from polygons
v.to.points -n in=polyTopology out=multiNode \
type=boundary --overwrite
# # # remove duplicate nodes
v.clean in=multiNode out=nodeSimple \
type=point tool=rmdupl --overwrite
# # # get a unique cat for each node
v.category nodeSimple out=node \
layer=3 opt=add --overwrite
# # # upload distance fron nearest node
v.distance from=mypt to=node to_layer=3 \
to_type=point upload=dist \
col=nodeDist dmax=-1
# # # upload cat from nearest node
v.distance from=mypt to=node to_layer=3 \
to_type=point upload=cat \
col=nodeCat dmax=-1
# # # upload coordinates from nearest node
v.distance from=mypt to=node to_layer=3 \
to_type=point upload=to_x,to_y \
col=node_x,node_y \
dmax=-1 out=distFromNode --overwrite
v.db.select mypt col=cat,nodeDist,nodeCat,node_x,node_y

Upload della cat univoca dei boundaries sugli attributi dei punti posti sui boundaries medesimi

Nel modello vettoriale implementato in GRASS  i boundaries costituiscono i contorni delle aree e, normalmente, non sono identificati univocamente con una category.

v.db.addcol mypt col="boundCatOverlap"
# # # get unique cat (on layer 2)
v.category in=polyTopology type=boundary \
out=arcs layer=2 opt=add --overwrite
v.distance from=mypt to=arcs to_layer=2 \
to_type=boundary upload=cat \
col=boundCatOverlap dmax=0
v.db.select mypt col=cat,boundCatOverlap

Upload della cat univoca dei boundaries sugli attributi dei punti posti entro due valori di distanza dati

È una veriante di quanto appena visto: i punti sui quali caricare la cat dei boundaries sono posti entro una distanza minima e massima dai boundaries medesimi.

v.db.addcol mypt col="boundCatZone double precision"
v.distance from=mypt to=arcs to_layer=2 \
to_type=boundary upload=cat \
col=boundCatZone out=distFromZone \
dmin=0.2 dmax=2 --overwrite
v.db.select mypt col=cat,boundCatOverlap,boundCatZone

Punti compresi entro due valori di distanza (0.2 e 2) dai boundaries


lunedì 9 febbraio 2009

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.