En posts anteriores (1, 2) se consideró la intersección de rasgos (features) adyacentes de polígonos seleccionados en PyQGIS para crear un plugin de QGIS 2 que permitiera generar las fronteras comunes como multipolilíneas. En éste se adaptó el código allí propuesto para que pudiese funcionar en QGIS 3 y generar un QDialog, con sus elementos Qt equivalentes, en una clase personalizada que se ejecuta desde la Python Console.
Se descubrió, además, que el referido código no funcionaba adecuadamente con los polígonos de los estados Arizona, Colorado, Nuevo México y Utah de USA, porque contiene un punto de unión cuádruple que genera dos puntos (no polilíneas) para la selección de dos rásgos adyacentes. Aunque se corrigió la falta de funcionalidad no se añadió como excepción.El código producido se incluye a continuación:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 | from PyQt5.QtCore import Qt import itertools class Dlg(QDialog): def __init__(self): QDialog.__init__(self) self.layout = QGridLayout(self) self.label1 = QLabel('Select polygon layer: ') self.wcb = QgsMapLayerComboBox() self.wcb.setFixedWidth(200) self.btn1 = QPushButton('OK', self) self.btn1.setFixedWidth(50) self.layout.addWidget(self.label1, 0, 0) self.layout.addWidget(self.wcb, 0, 1) self.layout.addWidget(self.btn1, 1, 1) self.btn1.clicked.connect(self.borders_features) def borders_features(self): layer = self.wcb.currentLayer() try: features = layer.selectedFeatures() list = range(len(features)) types = [ features[i].geometry().intersection(features[j].geometry()).type() for i,j in itertools.combinations(list, 2) if features[i].geometry().intersects(features[j].geometry()) ] intersections = [ features[i].geometry().intersection(features[j].geometry()).asWkt() for i,j in itertools.combinations(list, 2) if features[i].geometry().intersects(features[j].geometry()) ] n = len(intersections) if n == 0: iface.messageBar().pushMessage("Warning: ", "There is not adjacent selected features", Qgis.Warning, 5) return #creating a memory layer for multilinestring crs = layer.crs() epsg = crs.postgisSrid() if crs.isGeographic() is False: string = "length_m" else: string = "length_km" uri = "MultiLineString?crs=epsg:" + str(epsg) + "&field=id:integer&field=" + string + ":double""&index=yes" mem_layer = QgsVectorLayer(uri, "multipolyline", "memory") QgsProject.instance().addMapLayer(mem_layer) mem_layer.startEditing() #Set features feature = [QgsFeature() for i in range(n)] k = 0 dao = QgsDistanceArea() dao.setEllipsoid('WGS84') for i in range(n): if types[i] == 1: #set geometry geom = QgsGeometry.fromWkt(intersections[i]) if crs.isGeographic() is False: l = geom.length() else: l = dao.measureLength(geom)/1000 feature[i].setGeometry(geom) #set attributes values feature[i].setAttributes([k, l]) mem_layer.addFeature(feature[i]) k += 1 #stop editing and save changes mem_layer.commitChanges() except AttributeError: self.iface.messageBar().pushMessage("Warning: ", "There is not adjacent selected features", Qgis.Warning, 5) return w = Dlg() w.setWindowTitle('Borders Features') w.setWindowFlags(Qt.WindowStaysOnTopHint) w.move(300,200) w.show() |
Se probó con el shapefile ya mencionado seleccionando los cuatro estados; tal como se aprecia en la imagen siguiente:
Después de ejecutado el código anterior en la Python Console de QGIS, se puede observar en la imagen a continuación las cuatro líneas producidas. La tabla de atributos incluye adicionalmente la longitud, en kilómetros, de éstas.
No hay comentarios:
Publicar un comentario