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.


Nessun commento:

Posta un commento