[Gfoss] Qgis e pyQgis: funzione contains()

Salve a tutti,
sto dando un occhio alle librerie di pyQgis e vorrei emulare in python la funzione “Contenuto” presente nel plugin Spatial Query di Qgis.

La mia necessità è quella di inserirlo nel plugin pyarchinit e i dati iniziali sono:

partendo da 2 layer, uno di poligoni e l’altro di punti, sapere per ogni poligono quanti punti ricadono in esso; che poi è quello che fa egregiamente il plugin di cui sopra.

Fino a caricare le singole geometrie dal layer chiamandole da postgis ci arrivo, però, dando un occhio al cookbook di Qgis [1], vedo che la funzione contains() sta all’interno della classe QgsGeometry, e viene chiamata a partire da poligoni, punti o linee, create ad hoc.

Partendo dai miei due layer caricati mediante la classe QgsVectorLayer e relativi metodi, come faccio a passarli a QgsGeometry? Devo ricavare i nodi di ogni poligono, e passarli a QgsGeometri e poi usare contains(), o c’è un altra strada?

Qualcuno mi sa dare una mano? Vado direttamente alla lista degli sviluppatori di Qgis?

Ciao a tutti

Luca

[1] http://www.qgis.it/pyqgis-cookbook/geometry.html

Ciao Luca,

2011/2/23 Luca Mandolesi <mandoluca@gmail.com>

partendo da 2 layer, uno di poligoni e l’altro di punti, sapere per ogni poligono quanti punti ricadono in esso; che poi è quello che fa egregiamente il plugin di cui sopra.

Fino a caricare le singole geometrie dal layer chiamandole da postgis ci arrivo, però, dando un occhio al cookbook di Qgis [1], vedo che la funzione contains() sta all’interno della classe QgsGeometry, e viene chiamata a partire da poligoni, punti o linee, create ad hoc.

Partendo dai miei due layer caricati mediante la classe QgsVectorLayer e relativi metodi, come faccio a passarli a QgsGeometry? Devo ricavare i nodi di ogni poligono, e passarli a QgsGeometri e poi usare contains(), o c’è un altra strada?

se hai i due layer, siano vlPolygons and vlPoints

vlPolygons.select( ) # recuperi tutte le geometrie senza attributi
featPoly = QgsFeature() # crei una feature vuota per il poligono
while vlPolygons.nextFeature( featPoly ): # cicli sulle feature recuperate, featPoly conterrà la feature poligonale attuale
vlPoints.select( , featPoly.geometry().boundingBox() ) # recupera i punti nel bbox del poligono
featPoint = QgsFeature() # crei una feature vuota per il punto
while vlPoints.nextFeature( featPoint ): # cicli sulle feature recuperate, featPoint conterrà la feature puntale attuale
if featPoly.geometry().contains( featPoint.geometry() ): # adesso con la contains() verifichi che effettivamente sia contenuto nel poligono

TODO: qui lavori sul punto

ad es. aggiungi il punto (l’id della feature puntale) alla lista o dizionario o quel che sia

Non l’ho provato, l’ho scritto on-the-fly, quindi potrebbe non essere perfettamente funzionante,
anche se sono abbastanza convinto di averlo scritto bene :wink:

Oppure più semplicemente potresti recuperare tutti sia tutti i poligoni che tutti i punti, ma in
quel caso saresti costretto a controllare tutti i punti per ogni poligono… no buono.

Qualcuno mi sa dare una mano? Vado direttamente alla lista degli sviluppatori di Qgis?

E’ uguale, forse però ti avrebbero risposto prima :slight_smile:
Piuttosto, guardati il plugin ClosestFeatureFinder (repo Faunalia) che usa anche l’indice spaziale
per recuperare la feature più vicina. Sicuramente ti darà qualche spunto in più.

Ciao a tutti

Saluti.

Luca

[1] http://www.qgis.it/pyqgis-cookbook/geometry.html


Iscriviti all’associazione GFOSS.it: http://www.gfoss.it/drupal/iscrizione
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e’ una lista di discussione pubblica aperta a tutti.
Non inviate messaggi commerciali.
I messaggi di questa lista non rispecchiano necessariamente
le posizioni dell’Associazione GFOSS.it.
502 iscritti all’11.2.2011


Giuseppe Sucameli

Fantastico Giuseppe ormai sei il mio best italian pyQgis guru! 8 )
Appena posso lo testo così vedo se ho capito bene il meccanismo.
Grazie mille!
ciao

Luca

Ciao Giuseppe,
ho provato il plugin e gira correttamente. Tuttavia sto sbagliando qualcosa nella creazione del dizionario.

Ho modificato lo script (vedi sotto) per le mie esigenze chiamando tutti i poligoni (strati archeologici) e le quote relative di uno scavo archeologico.
Ho creato un dizionario_strati e una lista_quote. a questo punto vorrei ottenere un dizionario tipo questo:

{“strato_q”:[“quota_1”, “quota_2”, “quota_3”], … }

tuttavia, modificando in questo modo lo script:

if featPoly.geometry().contains( featPoint.geometry() ): # adesso con la contains() verifichi che effettivamente sia contenuto nel poligono
lista_quote.append(featPoint.id())
dizionario_strati[featPoly.id()] = lista_quote

e stampando la lunghezza del dizionario, noto che il while passa in rassegna tutte le geometrie poligonali (i miei strati), ma inserisce per ogni strato in lista quote tutti i punti che ho richiamato, e non solo quelli che cadono nella singola geometria. Ho fatto una verifica manuale su una geometria

https://picasaweb.google.com/mandoluca/Qgis#5577383282759381138

qui dovrebbero comparire nella lista legata alla geometria bordata in rosso, solo tre quote, invece le ho richiamate tutte…ovviamente l’errore è in:

lista_quote.append(featPoint.id())

Come faccio a passare alla lista solo i punti che cadono nella mia geometria?

Grazie ancora

Luca

import sys
sys.path.insert(0, ‘/Applications/QGIS.app/Contents/Resources/python/’)
from qgis.core import *
from qgis.gui import *

supply path to where is your qgis installed

QgsApplication.setPrefixPath(“/Applications/QGIS.app/Contents/MacOS/”, True)

load providers

QgsApplication.initQgis()

uri = QgsDataSourceURI()

set host name, port, database name, username and password

uri.setConnection(“localhost”, “5432”, “pyarchinit”, “postgres”, “mypass”)

set database schema, table name, geometry column and optionaly subset (WHERE clause)

uri.setDataSource(“public”, “pyunitastratigrafiche”, “the_geom”, “scavo_s = ‘Palazzo Ghetti, Rimini’ AND area_s = 1”)
vlPolygons = QgsVectorLayer(uri.uri(), “US”, “postgres”)
uri.setDataSource(“public”, “pyarchinit_quote”, “the_geom”, “sito_q = ‘Palazzo Ghetti, Rimini’ AND area_q = 1”)
vlPoints = QgsVectorLayer(uri.uri(), “Quote”, “postgres”)

if vlPolygons.isValid():
print “Layer polygons loaded!”

if vlPoints.isValid():
print “Layer points loaded!”

vlPolygons.select( ) # recuperi tutte le geometrie senza attributi
featPoly = QgsFeature() # crei una feature vuota per il poligono

dizionario_strati = {}
lista_quote =

while vlPolygons.nextFeature( featPoly ): # cicli sulle feature recuperate, featPoly conterra la feature poligonale attuale
vlPoints.select( , featPoly.geometry().boundingBox() ) # recupera i punti nel bbox del poligono
featPoint = QgsFeature() # crei una feature vuota per il punto
while vlPoints.nextFeature( featPoint ): # cicli sulle feature recuperate, featPoint conterra la feature puntale attuale
if featPoly.geometry().contains( featPoint.geometry() ): # adesso con la contains() verifichi che effettivamente sia contenuto nel poligono
lista_quote.append(featPoint.id())
dizionario_strati[featPoly.id()] = lista_quote

print diz

print "Numero di poligoni controllati: ", len(diz)

Scusami, sarà l’ora tarda…non svuotavo la lista dopo averla passata al dizionario…
ciao

Luca

Ciao Luca,

2011/2/24 Luca Mandolesi <mandoluca@gmail.com>

qui dovrebbero comparire nella lista legata alla geometria bordata in rosso, solo tre quote, invece le ho richiamate tutte…ovviamente l’errore è in:

lista_quote.append(featPoint.id())

Non è esattamente quella istruzione il problema. Tuttavia fino a questo momento non avevo
capito dalle tue frasi se nella lista c’erano tutti i punti o tutti quelli dentro il bbox anche se
effettivamente stavano fuori.

Ma a vedere dal codice sotto, il problema è il primo:
il dizionario da segno che tutti i poligono contengono tutti i punti… :expressionless:
Impressionante! :slight_smile:

Ecco la soluzione:

import sys
sys.path.insert(0, ‘/Applications/QGIS.app/Contents/Resources/python/’)
from qgis.core import *
from qgis.gui import *

supply path to where is your qgis installed

QgsApplication.setPrefixPath(“/Applications/QGIS.app/Contents/MacOS/”, True)

load providers

QgsApplication.initQgis()

uri = QgsDataSourceURI()

set host name, port, database name, username and password

uri.setConnection(“localhost”, “5432”, “pyarchinit”, “postgres”, “mypass”)

set database schema, table name, geometry column and optionaly subset (WHERE clause)

uri.setDataSource(“public”, “pyunitastratigrafiche”, “the_geom”, “scavo_s = ‘Palazzo Ghetti, Rimini’ AND area_s = 1”)
vlPolygons = QgsVectorLayer(uri.uri(), “US”, “postgres”)
uri.setDataSource(“public”, “pyarchinit_quote”, “the_geom”, “sito_q = ‘Palazzo Ghetti, Rimini’ AND area_q = 1”)
vlPoints = QgsVectorLayer(uri.uri(), “Quote”, “postgres”)

if vlPolygons.isValid():
print “Layer polygons loaded!”

if vlPoints.isValid():
print “Layer points loaded!”

vlPolygons.select( ) # recuperi tutte le geometrie senza attributi
featPoly = QgsFeature() # crei una feature vuota per il poligono

dizionario_strati = {}

fino qui è ok.

lista_quote =

è questa linea che crea tutti i problemi:
tu crei una lista fuori dai 2 cicli, quindi la lista sarà sempre la stessa ad ogni ciclo e quindi
aggiungerai sempre la medesima lista a tutte le voci del dizionario. La lista alla fine conterrà
tutti i punti (o quasi tutti, dipende dall’unione dei bbox dei poligoni), quindi ogni voce del
dizionario avrà sempre come valore una lista (la medesima) con tutti i punti.

while vlPolygons.nextFeature( featPoly ): # cicli sulle feature recuperate, featPoly conterra la feature poligonale attuale

vlPoints.select( , featPoly.geometry().boundingBox() ) # recupera i punti nel bbox del poligono
featPoint = QgsFeature() # crei una feature vuota per il punto

Inserisci invece la riga di creazione della lista qui e tutto si dovrebbe risolvere :wink:

while vlPoints.nextFeature( featPoint ): # cicli sulle feature recuperate, featPoint conterra la feature puntale attuale

if featPoly.geometry().contains( featPoint.geometry() ): # adesso con la contains() verifichi che effettivamente sia contenuto nel poligono

lista_quote.append(featPoint.id())
dizionario_strati[featPoly.id()] = lista_quote

print diz

print "Numero di poligoni controllati: ", len(diz)

Saluti.


Giuseppe Sucameli

Infatti :slight_smile:
Devo dire che forse la mia risposta sarà stata un pò troppo prolissa.
Saluti.

2011/2/25 Luca Mandolesi <mandoluca@gmail.com>

Scusami, sarà l’ora tarda…non svuotavo la lista dopo averla passata al dizionario…
ciao

Luca


Giuseppe Sucameli