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


Nessun commento:

Posta un commento