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