[QGIS-it-user] pyQgis: come verificare una gemetria id=X nel corso di un ciclo?

Salve,

su spatialite ho un layer di poligoni e dato un certo valore contenuto in un campo che mi segnala se quel poligono si sovrappone ad un altro nel medesimo layer, vorrei verificare con intersect se le due geometrie rispettano quanto segnalato nel campo.

Facendo un ciclo tra le feature, quando trovo il valore di sovrapposizione, come si richiama una geometria all’interno dello stesso layer senza uscire dal ciclo?

#print a[0].geometry().intersects(b[0].geometry())

Sotto metto un po’ di codice con la parte di pseudocodice.

Scusate, partito l’invio, il print di prima era un errore.

Io pensavo ad una cosa tipo [0], ma ricevo poi errore perchè a e b non accedono a geometry:

[0]
request = QgsFeatureRequest().setFilterExpression( ‘“id_us” = 2’ )
a = vlayer.getFeatures( request )
request = QgsFeatureRequest().setFilterExpression( ‘“id_us” = 3’ )
b = vlayer.getFeatures( request )

a.geometry().intersects(b.geometry())

Grazie per i suggerimenti.

Luca

Ciao,
non ho capito se vuoi sapere soltanto quali siano le features che si intersecano tra loro (ad es. una lista) oppure se devi fare un ulteriore controllo su di esse una volta verificata l’intersezione (ad esempio: se la feature “a” interseca altre 3 features all’interno del layer, verificare che a[“id_us”] sia pari a 3).
A prescindere dall’obiettivo, vuoi solo stampare qualcosa o modificare gli attributi?
Marco

···

Il giorno 5 novembre 2016 09:53, Luca Mandolesi <mandoluca@gmail.com> ha scritto:

Scusate, partito l’invio, il print di prima era un errore.

Io pensavo ad una cosa tipo [0], ma ricevo poi errore perchè a e b non accedono a geometry:

[0]
request = QgsFeatureRequest().setFilterExpression( ‘“id_us” = 2’ )
a = vlayer.getFeatures( request )
request = QgsFeatureRequest().setFilterExpression( ‘“id_us” = 3’ )
b = vlayer.getFeatures( request )

a.geometry().intersects(b.geometry())

Grazie per i suggerimenti.

Luca


QGIS-it-user mailing list
QGIS-it-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/qgis-it-user

Ciao Marco.

per ora devo solo mandare fuori l’id delle 2 geometrie e la verifica della presenza di intersezione. Dal cookbook si arriva bene qua:

iter = vlayer.getFeatures()

Faccio un loop sulle singole geometrie e se becco il rapporto copre

ricavare da iter la feature id=x senza fare un loop.

Tipo con una chiamata:

#singFeat = vlayer.featureAtId(1)

Ma mi dice che vlayer non ha featureAtId come metodo

E’ possibile?

Grazie

Luca

···

Il giorno 5 novembre 2016 17:11, Marco Grisolia <marco.grisolia5@gmail.com> ha scritto:

Ciao,
non ho capito se vuoi sapere soltanto quali siano le features che si intersecano tra loro (ad es. una lista) oppure se devi fare un ulteriore controllo su di esse una volta verificata l’intersezione (ad esempio: se la feature “a” interseca altre 3 features all’interno del layer, verificare che a[“id_us”] sia pari a 3).
A prescindere dall’obiettivo, vuoi solo stampare qualcosa o modificare gli attributi?
Marco

Il giorno 5 novembre 2016 09:53, Luca Mandolesi <mandoluca@gmail.com> ha scritto:

Scusate, partito l’invio, il print di prima era un errore.

Io pensavo ad una cosa tipo [0], ma ricevo poi errore perchè a e b non accedono a geometry:

[0]
request = QgsFeatureRequest().setFilterExpression( ‘“id_us” = 2’ )
a = vlayer.getFeatures( request )
request = QgsFeatureRequest().setFilterExpression( ‘“id_us” = 3’ )
b = vlayer.getFeatures( request )

a.geometry().intersects(b.geometry())

Grazie per i suggerimenti.

Luca


QGIS-it-user mailing list
QGIS-it-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/qgis-it-user

On 11/6/16, Luca Mandolesi <mandoluca@gmail.com> wrote:

Ciao Marco.

ciao Luca,

per ora devo solo mandare fuori l'id delle 2 geometrie e la verifica della
presenza di intersezione. Dal cookbook si arriva bene qua:
.......

credo tu abbia due possibilità: o confronti le geometrie e quando
trovi intersezioni verifichi la bontà dell'attributo interessato o,
come mi sembrava indicasse il codice che avevi postato, prendi le
feature con attributo positivo e verifichi l'effettiva intersezione;
non so se ho capito bene fino a questo punto, ma confesso di non aver
capito il tuo problema da quì in avanti:

ciao,
giuliano

Il giorno 6 novembre 2016 13:33, Giuliano Curti <giulianc51@gmail.com> ha
scritto:

On 11/6/16, Luca Mandolesi <mandoluca@gmail.com> wrote:
> Ciao Marco.
il tuo problema da quì in avanti:

Ciao Giuliano,
il mio unico problema è:

mentre sono in un ciclo for e analizzo il campo di un record, come faccio a
chiamare un altro record e caricarne la geometria. La parte del perchè devo
farlo è inerente all'archeologia e non è importante [0].

In pseudo codice:

per ogni feature in tutte le features:
      se il feature.campo == "controlla":
            id_altra_feature = feature.campo_con_id_da_verificare
            altra_geometria = features.ricava_record(id_altra_feature) <-
questo è il passaggio che non riesco a ricavare
            return feature.geometria.interseca(altra_geometria.geometria)

Ciao
Luca

[0]
Per approfondire ma non necessario:
- nella tabella alfanumerica sono indicati i rapporti tra i record dentro
ad una lista di liste id=1, campo = [[copre, 2], [copre, 3], [coperto da,
4], [...]]

- devo verificare solo certi rapporti

- non posso sapere a priori se è buono il disegno oppure l'attributo
inserito (il data entry crea errori senza possibilità di verifica runtime)

- ad un record alfanumerico corrispondono molte geometrie

Quindi quello che si vuole ottenere, come per il plugin di verifica della
topologia è un bel listato con il tipo di rapporto e l'assenza di
corrispondenza da un punto di vista topografico...ma per fare questo ho
tutto l'occorrente.

On 11/6/16, Luca Mandolesi <mandoluca@gmail.com> wrote:

Il giorno 6 novembre 2016 13:33, Giuliano Curti <giulianc51@gmail.com> ha
scritto:

On 11/6/16, Luca Mandolesi <mandoluca@gmail.com> wrote:

mentre sono in un ciclo for e analizzo il campo di un record, come faccio a
chiamare un altro record e caricarne la geometria.

rimango sempre nel dubbio di aver capito bene :slight_smile: questo spezzone di
codice cmq ti dà le due feature f1 e f2, le rispettive geometrie, le
confronta (scartando ovviamente il controllo
con sè stessa) e puoi accedere ad ogni altro loro attributo;

  for f1 in layer.getFeatures():
    ......
    for f2 in layer.getFeatures():
      ........
      if f2.id() <> f1.id():
        .........
        if f1.geometry().intersects(f2.geometry()):
          (controllo se l'attributo delle due feature è corretto)
          .......

NB: per il tuo caso forse è ridondante perchè confronta A con B e B
con A, mentre a te potrebbe bastare confrontare A con B, ad es.
modificando il test in
      if f2.id() > f1.id():
o qualche modo migliore che sapranno indicarti altri;

Ciao
Luca

ciao,
giuliano

Sperando di essere ancora in tempo, ti parlo della mia idea.
Premetto di non avere capito del tutto come siano strutturati i tuoi dati, per cui la soluzione che ti propongo potrebbe essere leggermente modificata in funzione di come essi siano presentati.
Ho creato questa situazione, sperando che sia rappresentativa del tuo caso:


Lo script che ti incollo nel seguito esegue queste operazioni:

  1. cerco le features che hanno, tra gli attributi, un determinato valore ( o “rapporto”, come lo chiami tu);
  2. salvo queste features in un dizionario;
  3. estrapolo dalla colonna del rapporto (in questo caso, quella con indice [1]) i valori numerici che fanno riferimento alle geometrie da confrontare;
  4. verifico che la feature in esame intersechi le altre features di confronto;
  5. stampo in ordine: il codice della feature, il rapporto dichiarato, la veridicità di quanto riportato nel rapporto (VERO o FALSO).

Per questo esempio, lo script stampa a video il risultato:
Feature No. 3, 3 copre 2, FALSO
Feature No. 5, 5 copre 3 e 4, VERO
Feature No. 4, 4 copre 5, VERO
Feature No. 6, 6 coperto da 7, VERO
Feature No. 7, 7 copre 6, VERO

Come puoi notare, ho scritto volutamente un rapporto sbagliato per la feature 3 (“3 copre 2”) per validare lo script.

Sperando di non essere andato troppo fuori strada (o di aver commesso altri errori formali), ti incollo il codice che ho utilizzato io, ciao!

Marco

···

Il giorno 6 novembre 2016 15:06, Giuliano Curti <giulianc51@gmail.com> ha scritto:

On 11/6/16, Luca Mandolesi <mandoluca@gmail.com> wrote:

Il giorno 6 novembre 2016 13:33, Giuliano Curti <giulianc51@gmail.com> ha
scritto:

On 11/6/16, Luca Mandolesi <mandoluca@gmail.com> wrote:

mentre sono in un ciclo for e analizzo il campo di un record, come faccio a
chiamare un altro record e caricarne la geometria.

rimango sempre nel dubbio di aver capito bene :slight_smile: questo spezzone di
codice cmq ti dà le due feature f1 e f2, le rispettive geometrie, le
confronta (scartando ovviamente il controllo
con sè stessa) e puoi accedere ad ogni altro loro attributo;

for f1 in layer.getFeatures():

for f2 in layer.getFeatures():

if f2.id() <> f1.id():

if f1.geometry().intersects(f2.geometry()):
(controllo se l’attributo delle due feature è corretto)

NB: per il tuo caso forse è ridondante perchè confronta A con B e B
con A, mentre a te potrebbe bastare confrontare A con B, ad es.
modificando il test in
if f2.id() > f1.id():
o qualche modo migliore che sapranno indicarti altri;

Ciao
Luca

ciao,
giuliano


QGIS-it-user mailing list
QGIS-it-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/qgis-it-user

Grazie Marco,

credo che hai risolto il mio problema, poi lo riadatto e vedo se mi funzia.

Dal codice mi pare che fai quello che voglio io ma partendo da una strutturazione diversa del codice.

Io iteravo sulle feature e trovando la condizione copre == vero, cercavo di richiamare la feature tramite il suo id per confrontarla dentro al ciclo.

Nel tuo caso se capisco bene prima si itera su tutte e si ricava la condizione copre == vero, poi si reitera su tutte le feature che hanno solo quella condizione e salti il confronto con se stessa con if feat2 != feat1.

Invece io volevo evitare questo passaggio: for feat2 in vlayer.getFeatures():

e mentre ero nel primmo ciclo, trovando che 1 copre 2, dirgli di chiamare da vlayer solo la feat == 2 e verificae l’intersect.

Ora non mi resta che provare.

Grazie mille

Luca

Ho cambiato il codice di prima, omettendo il ciclo for che hai
segnalato: mi è bastato cambiare "l'ordine degli addendi" e dare una
diversa definizione del dizionario "allfeatures". Trovi il codice
commentato qui sotto, ma aggiungo una nota: in presenza di alcune
intersezioni 'vere' e di altre intersezioni 'false', il vecchio codice
generava falsi positivi assegnado il valore 'vero' (in sostanza,
bastava che ci fosse almeno un'intersezione tra le features nel
rapporto perché il risultato fosse 'VERO'). Spero che ora sia a posto!
Ciao,
Marco

#####################

# se vuoi verificare l'intersezione, 'copre' e 'coperto' significano
la stessa cosa
test_string1 = 'copre'
test_string2 = 'coperto'

# Crea un dizionario con tutte le features
allfeatures={}
for item in vlayer.getFeatures():
    allfeatures[item[0]] = QgsFeature(item)

for feat1 in allfeatures.values():
    feat = 'Feature No. %s,' %(feat1[0])
    string = '%s %s,' %(feat1[0], feat1[1])
    inGeom1 = feat1.geometry()
    if test_string1 in item[1] or test_string2 in item[1]: # faccio il
controllo sul rapporto
        test = feat1[1]
        numb_test1=[int(s) for s in test.split() if s.isdigit()] #
estraggo i valori numerici dal rapporto
        if numb_test1: # vado avanti solo se numb_test1 e' una lista
non nulla (cosi' evito di portarmi dietro i rapporti con formattazione
sbagliata o non richiesta
            numbs = # qui inserisco tutti i risultati delle
intersezioni (vedi dopo)
            for k in numb_test1: # verifico l'intersezione tra feat1 e
tutti i valori numerici presenti nel rapporto
                inGeom2 = allfeatures[k].geometry()
                if inGeom1.intersects(inGeom2):
                    check = 'VERO'
                    numbs.append(check)
                else:
                    check = 'FALSO'
                    numbs.append(check)
            if numbs.count(numbs[0]) != len(numbs): # se i valori di
numbs non sono tutti uguali (tutti 'VERO' o tutti 'FALSO', allora il
rapporto e' 'FALSO'
                check = 'FALSO'
            print feat, string, check

Il 7 novembre 2016 20:40, Luca Mandolesi <mandoluca@gmail.com> ha scritto:

Grazie Marco,
credo che hai risolto il mio problema, poi lo riadatto e vedo se mi funzia.
Dal codice mi pare che fai quello che voglio io ma partendo da una
strutturazione diversa del codice.
Io iteravo sulle feature e trovando la condizione copre == vero, cercavo di
richiamare la feature tramite il suo id per confrontarla dentro al ciclo.

Nel tuo caso se capisco bene prima si itera su tutte e si ricava la
condizione copre == vero, poi si reitera su tutte le feature che hanno solo
quella condizione e salti il confronto con se stessa con if feat2 != feat1.

Invece io volevo evitare questo passaggio: for feat2 in
vlayer.getFeatures():
e mentre ero nel primmo ciclo, trovando che 1 copre 2, dirgli di chiamare da
vlayer solo la feat == 2 e verificae l'intersect.

Ora non mi resta che provare.

Grazie mille
Luca

Luca,

... e meno male che sei un "archeologo" !!!!
E se eri informatico che facevi?? .... letteralmente spaccavi !!

Bravo non ce che dire !!

Salutoni
Nino

Il 7 novembre 2016 20:40, Luca Mandolesi <mandoluca@gmail.com> ha scritto:

Grazie Marco,
credo che hai risolto il mio problema, poi lo riadatto e vedo se mi funzia.
Dal codice mi pare che fai quello che voglio io ma partendo da una
strutturazione diversa del codice.
Io iteravo sulle feature e trovando la condizione copre == vero, cercavo di
richiamare la feature tramite il suo id per confrontarla dentro al ciclo.

Nel tuo caso se capisco bene prima si itera su tutte e si ricava la
condizione copre == vero, poi si reitera su tutte le feature che hanno solo
quella condizione e salti il confronto con se stessa con if feat2 != feat1.

Invece io volevo evitare questo passaggio: for feat2 in
vlayer.getFeatures():
e mentre ero nel primmo ciclo, trovando che 1 copre 2, dirgli di chiamare da
vlayer solo la feat == 2 e verificae l'intersect.

Ora non mi resta che provare.

Grazie mille
Luca

_______________________________________________
QGIS-it-user mailing list
QGIS-it-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/qgis-it-user