[Gfoss] Trovare extent del layer da console di python

Ciao a tutti,

sto provando a creare uno script nella console di qgis che come prima cosa estragga da ciascuna cartella lo shape chiamato “particella.shp” e poi esegua un controllo sull’estensione del layer per verificare se è all’interno della provincia che mi interessa oppure no, nel primo caso vorrei che il layer venisse aggiunto nella toc altrimenti scartato.
Ammesso che mi sia riuscita a spiegare vi allego parte dello script

import sys,os,string
import processing

folder = ‘C:\Dati_geografici\Particelle catastali\Catastale’
lista_folder = (os.listdir(folder))

for sub_folder in lista_folder:
new_folder = folder+‘\’+sub_folder
lista_shp = (os.listdir(new_folder))

for shp in lista_shp:
if string.find(shp,‘particella.shp’) >=0:
MANCA FUNZIONE PER CONTROLLARE L’ESTENSIONE DEL LAYER
qgis.utils.iface.addVectorLayer(shp,shp,“ogr”)

cercando tra le API di qgis ho visto che esiste questa classe QgsVectorLayer che tra le varie cose calcola anche l’extent ma non so come richiamarla e soprattutto quale operatore di contronto usare per verificare che il mio layer cada dentro un determinato extent.
Altro problema il caricamento del layer nella toc mediante addVectorLayer mi restituisce questo errore “Il layer non è valido: Il layer D197__particella.shp non è valido e non può essere aggiunto alla mappa” che nonc apisco da cosa dipenda.

Grazie in anticipo!

senza scomodare gdal e rimanendo in pyqgis

1) carichi il layer in memoria (senza visualizzarlo)

vl = QgsVectorLayer(shp, "mioshape", "ogr")

2) prendi l'extent del vettore

extent = vl.extent()

come vedrai extent e' un QgsRectangle =>

3) ho bisogno di un QgsGeometry per usarne gli operatori spaziali

geomExtent = QgsGeometry.fromRect(extent)

4) alla stessa cosa per generarti la geometria dell'extent che vuoi contollare

geomExtentDiControllo = QgsGeometry.fromRect(...un qgis rectangle... )

ti consiglio di vederti: http://qgis.org/api/classQgsGeometry.html
per costruirti il QgsREctangle di controllo

5) usi l'operatore di QgsGeometry che meglio ti si addice

isInside = geomExtent.within(geomExtentDiControllo)

oppure

isInside = geomExtent.intersects(geomExtentDiControllo)

if isInside:
    QgsMapLayerRegistri.instance().addMapLayer(vl)

ciao Luigi Pirelli

2015-03-06 10:36 GMT+01:00 Romina Di Paolo <romi.dipaolo@gmail.com>:

Ciao a tutti,

sto provando a creare uno script nella console di qgis che come prima cosa
estragga da ciascuna cartella lo shape chiamato "particella.shp" e poi
esegua un controllo sull'estensione del layer per verificare se è
all'interno della provincia che mi interessa oppure no, nel primo caso
vorrei che il layer venisse aggiunto nella toc altrimenti scartato.
Ammesso che mi sia riuscita a spiegare vi allego parte dello script

import sys,os,string
import processing

folder = 'C:\\Dati_geografici\\Particelle catastali\\Catastale'
lista_folder = (os.listdir(folder))

for sub_folder in lista_folder:
    new_folder = folder+'\\'+sub_folder
    lista_shp = (os.listdir(new_folder))

    for shp in lista_shp:
        if string.find(shp,'particella.shp') >=0:
       MANCA FUNZIONE PER CONTROLLARE L'ESTENSIONE DEL LAYER
            qgis.utils.iface.addVectorLayer(shp,shp,"ogr")

cercando tra le API di qgis ho visto che esiste questa classe QgsVectorLayer
che tra le varie cose calcola anche l'extent ma non so come richiamarla e
soprattutto quale operatore di contronto usare per verificare che il mio
layer cada dentro un determinato extent.
Altro problema il caricamento del layer nella toc mediante addVectorLayer mi
restituisce questo errore "Il layer non è valido: Il layer
D197__particella.shp non è valido e non può essere aggiunto alla mappa" che
nonc apisco da cosa dipenda.

Grazie in anticipo!

_______________________________________________
Gfoss@lists.gfoss.it
http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss
Questa e' una lista di discussione pubblica aperta a tutti.
I messaggi di questa lista non hanno relazione diretta con le posizioni
dell'Associazione GFOSS.it.
666+40 iscritti al 5.6.2014

folder = ‘S:\Dati_geografici_RP\I.6_Particelle catastali\Catastale’
lista_folder = (os.listdir(folder))
geomExtentDiControllo = QgsGeometry.fromRect(QgsRectangle (xmin =442882.70,ymin = 5015200.28,xmax =487902.38,ymax =5087429.92))

for sub_folder in lista_folder:
new_folder = folder+‘\’+sub_folder
lista_shp = (os.listdir(new_folder))

for shp in lista_shp:
if string.find(shp,‘PARTICELLA.shp’) >=0:
vl = QgsVectorLayer(shp, “mioshape”, “ogr”)
extent = vl.extent()
geomExtent = QgsGeometry.fromRect(extent)
isInside = geomExtent.within(geomExtentDiControllo)
if isInside:
QgsMapLayerRegistri.instance().addMapLayer(vl)

In grassetto la parte che mi è ancora lacunosa infatti il programma finisce senza restituire nulla!

Sto usando come riferimento questo perchè sono ancora alle prime armi come si può capire!
http://www.padido.eu/gfoss/qgis/workshop/workshop/esplorare/python_in_qgis_tutorial2.html

···

Grazie Luigi per l’esempio fornito,
Quindi da quello che hai riportato nell’esempio da una parte ho la mia lista di shp (extent = vl.extent) che deve essere confrontata con un extent di controllo ovvero geomExtentDiControllo. Il mio problema è che non so come accedere alle proprietà geometriche di quest’ultimo.

Ho provato così

2015-03-06 11:10 GMT+01:00 Luigi Pirelli <luipir@gmail.com>:

senza scomodare gdal e rimanendo in pyqgis

1) carichi il layer in memoria (senza visualizzarlo)

vl = QgsVectorLayer(shp, “mioshape”, “ogr”)

2) prendi l’extent del vettore

extent = vl.extent()

come vedrai extent e’ un QgsRectangle =>

3) ho bisogno di un QgsGeometry per usarne gli operatori spaziali

geomExtent = QgsGeometry.fromRect(extent)

4) alla stessa cosa per generarti la geometria dell’extent che vuoi contollare

geomExtentDiControllo = QgsGeometry.fromRect(…un qgis rectangle… )

ti consiglio di vederti: http://qgis.org/api/classQgsGeometry.html
per costruirti il QgsREctangle di controllo

5) usi l’operatore di QgsGeometry che meglio ti si addice

isInside = geomExtent.within(geomExtentDiControllo)

oppure

isInside = geomExtent.intersects(geomExtentDiControllo)

if isInside:
QgsMapLayerRegistri.instance().addMapLayer(vl)

ciao Luigi Pirelli

non avendo i dati sottomano, su due piedi farei il debug della nonna
1) controllerei che le x e y di geomExtentDiControllo siano lon e lat
2) controllerei che gli shape siano nella stessa proiezione con cui
hai espresso il bounding box di geomExtentDiControllo
3) controllerei che gli extent degli shape (fanne un print) per
verificare a mano se within debba darti true o false
4) verificherei anche con l'operatore intersects di QgsGeometry

ciao Luigi Pirelli

2015-03-09 12:44 GMT+01:00 Romina Di Paolo <romi.dipaolo@gmail.com>:

Grazie Luigi per l'esempio fornito,
Quindi da quello che hai riportato nell'esempio da una parte ho la mia lista
di shp (extent = vl.extent) che deve essere confrontata con un extent di
controllo ovvero geomExtentDiControllo. Il mio problema è che non so come
accedere alle proprietà geometriche di quest'ultimo.

Ho provato così

folder = 'S:\\Dati_geografici_RP\\I.6_Particelle catastali\\Catastale'
lista_folder = (os.listdir(folder))
geomExtentDiControllo = QgsGeometry.fromRect(QgsRectangle (xmin
=442882.70,ymin = 5015200.28,xmax =487902.38,ymax =5087429.92))

for sub_folder in lista_folder:
    new_folder = folder+'\\'+sub_folder
    lista_shp = (os.listdir(new_folder))

    for shp in lista_shp:
        if string.find(shp,'PARTICELLA.shp') >=0:
            vl = QgsVectorLayer(shp, "mioshape", "ogr")
            extent = vl.extent()
            geomExtent = QgsGeometry.fromRect(extent)
            isInside = geomExtent.within(geomExtentDiControllo)
            if isInside:
                QgsMapLayerRegistri.instance().addMapLayer(vl)

In grassetto la parte che mi è ancora lacunosa infatti il programma finisce
senza restituire nulla!

Sto usando come riferimento questo perchè sono ancora alle prime armi come
si può capire!
http://www.padido.eu/gfoss/qgis/workshop/workshop/esplorare/python_in_qgis_tutorial2.html

2015-03-06 11:10 GMT+01:00 Luigi Pirelli <luipir@gmail.com>:

senza scomodare gdal e rimanendo in pyqgis

1) carichi il layer in memoria (senza visualizzarlo)

vl = QgsVectorLayer(shp, "mioshape", "ogr")

2) prendi l'extent del vettore

extent = vl.extent()

come vedrai extent e' un QgsRectangle =>

3) ho bisogno di un QgsGeometry per usarne gli operatori spaziali

geomExtent = QgsGeometry.fromRect(extent)

4) alla stessa cosa per generarti la geometria dell'extent che vuoi
contollare

geomExtentDiControllo = QgsGeometry.fromRect(...un qgis rectangle... )

ti consiglio di vederti: http://qgis.org/api/classQgsGeometry.html
per costruirti il QgsREctangle di controllo

5) usi l'operatore di QgsGeometry che meglio ti si addice

isInside = geomExtent.within(geomExtentDiControllo)

oppure

isInside = geomExtent.intersects(geomExtentDiControllo)

if isInside:
    QgsMapLayerRegistri.instance().addMapLayer(vl)

ciao Luigi Pirelli