Trovare extent del layer da console di python

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Trovare extent del layer da console di python

Romina Di Paolo
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!

_______________________________________________
[hidden email]
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
Reply | Threaded
Open this post in threaded view
|

Re: Trovare extent del layer da console di python

Luigi Pirelli-2
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 <[hidden email]>:

> 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!
>
> _______________________________________________
> [hidden email]
> 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
_______________________________________________
[hidden email]
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
Reply | Threaded
Open this post in threaded view
|

Re: Trovare extent del layer da console di python

Romina Di Paolo

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!




2015-03-06 11:10 GMT+01:00 Luigi Pirelli <[hidden email]>:

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


_______________________________________________
[hidden email]
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
Reply | Threaded
Open this post in threaded view
|

Re: Trovare extent del layer da console di python

Luigi Pirelli-2
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 <[hidden email]>:

>
> 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 <[hidden email]>:
>
> 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
_______________________________________________
[hidden email]
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