5 votos

Seleccionar y cargar polígonos que intersecan a otro con PyQGIS (polígonos almacenados en PostGIS)

Todavía soy un novato en Python y estoy teniendo problemas para seleccionar y cargar polígonos que se cruzan con otros con PyQGIS.

He conseguido utilizar PyQGIS para lograr esa operación con shapefiles almacenados en mi portátil pero no he encontrado una solución para proceder al mismo tratamiento en shapefiles almacenados en PostGIS.

Esto es lo que he encontrado hasta ahora:

from qgis.core import *
from PyQt4.QtCore import *
from PyQt4.QtCore import QSettings
from PyQt4.QtCore import QSettings
from qgis.core import QgsVectorLayer, QgsDataSourceURI
import processing

uri = QgsDataSourceURI()
uri.setConnection("localhost", "5432", "name of my database", "my user name", "my password")
uri.setDataSource("database schema", "table name", "geom")
couche = "PM"
layer1 = QgsVectorLayer(uri.uri(), couche, "postgres")
if not layer1.isValid:
    print ('Chargement de la couche: echec !')
else:
    QgsMapLayerRegistry.instance().addMapLayer(layer1)

uri.setDataSource("database schema", "table name2", "geom")
couche = "batiments"
layer2 = QgsVectorLayer(uri.uri(), couche, "postgres" )

subset_layer2 = processing.runalg('qgis:extractbylocation', layer2,layer1,['intersects'],0.0,None)

uri.setDataSource("rbal", "parcelles", "geom")
couche = "parcelles"
layer3 = QgsVectorLayer(uri.uri(), couche, "postgres" )

subset_layer3 = processing.runalg('qgis:extractbylocation', layer3,layer1,['intersects'],0.0,None)

Cuando ejecuto el código no me da ningún error pero el tratamiento se carga eternamente sin devolver ningún resultado...


EDITAR :

Con ese script he conseguido seleccionar todos los polígonos contenidos en la capa2 que intersectan a la capa1 pero tengo que cargar la capa2 para hacerlo. ¿Hay alguna manera de seleccionar y cargar esos polígonos que se cruzan sin cargar toda la capa en primer lugar?

uri = QgsDataSourceURI()
uri.setConnection("localhost", "5432", "name of my database", "my user name", "my password")
uri.setDataSource("database schema", "table name", "geom")
couche = "PM"
layer1 = QgsVectorLayer(uri.uri(), couche, "postgres")
if not layer1.isValid:
    print ('Chargement de la couche: echec !')
else:
    QgsMapLayerRegistry.instance().addMapLayer(layer1)

uri.setDataSource("database schema", "table name2", "geom")
couche = "batiments"
layer2 = QgsVectorLayer(uri.uri(), couche, "postgres" )
if not layer2.isValid:
    print ('Chargement de la couche: echec !')
else:
    QgsMapLayerRegistry.instance().addMapLayer(layer2)

feats_layer1 = layer1.getFeatures()
for feat in feats_layer1:
    geom_layer1 = feat.geometry()
feats_layer2 = layer2.getFeatures()

for feature in feats_layer2:
    if feature.geometry().intersects(geom_layer1):
        layer2.select(feature.id())

4voto

wanderzen Puntos 87

Gracias a J.Monticolo he encontrado una solución existente en stackoverflow:

PyQGIS - Cargar una capa PostgreSQL/PostGIS desde una consulta SQL

QgsDataSourceuRI y left outer join

EDITAR :

Aunque los enlaces anteriores proporcionan (en mi ejemplo) la mejor solución, hay una manera de hacer lo mismo sólo con PyQGIS:

from qgis.core import *
from PyQt4.QtCore import *
from PyQt4.QtCore import QSettings
from PyQt4.QtCore import QSettings
from qgis.core import QgsVectorLayer, QgsDataSourceURI
from qgis.gui import QgsMapCanvas
import processing

uri = QgsDataSourceURI()
uri.setConnection("localhost", "5432", "name of my database", "my user name", "my password")
uri.setDataSource("database schema", "table name", "geom")
couche = "PM"
layer1 = QgsVectorLayer(uri.uri(), couche, "postgres")
if not layer1.isValid:
    print ('Chargement de la couche: echec !')
else:
    QgsMapLayerRegistry.instance().addMapLayer(layer1)

uri.setDataSource("database schema", "table name2", "geom")
couche = "batiments"
layer2 = QgsVectorLayer(uri.uri(), couche, "postgres" )
if not layer2.isValid:
    print ('Chargement de la couche: echec !')
else:
    QgsMapLayerRegistry.instance().addMapLayer(layer2)

inter_feats = []

feats_layer1 = layer1.getFeatures()
for feat in feats_layer1:
    geom_layer1 = feat.geometry()
feats_layer2 = layer2.getFeatures()

for feature in feats_layer2:
    if feature.geometry().intersects(geom_layer1):
        inter_feats.append(feature.id())

layer2.select(inter_feats)

intersecting_features = QgsVectorFileWriter.writeAsVectorFormat(layer2, 'C:/Users/Users/Documents/test.shp', "utf-8",  None, "ESRI Shapefile", onlySelected=True)

bati = QgsVectorLayer('C:/Users/Users/Documents/test.shp', "test", "ogr")
if not bat.isValid():
    print("Impossible de charger la couche !")

# add layer to the registry
QgsMapLayerRegistry.instance().addMapLayer(bati)

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X