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