Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Division de polígonos sin errores topológicos en QGIS usando Python
# -*- coding: utf-8 -*-
def calcula_lado_linea(x0, x1, x2, y0, y1, y2):
valor = (x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0)
if valor > 0 and valor >= (1 * 10**(-10)):
aux_valor = 1
elif valor > 0 and valor < (1 * 10**(-10)):
aux_valor = 0
elif valor < 0 and abs(valor) < (1 * 10**(-10)):
aux_valor = 0
elif valor == 0:
aux_valor = 0
else:
aux_valor = -1
return aux_valor
def intersec_segment_4pts(segm1_p1x, segm1_p1y, segm1_p2x, segm1_p2y, segm2_p1x, segm2_p1y, segm2_p2x, segm2_p2y):
deltaOne_x = segm1_p2x - segm1_p1x
deltaOne_y = segm1_p2y - segm1_p1y
deltaTwo_x = segm2_p2x - segm2_p1x
deltaTwo_y = segm2_p2y - segm2_p1y
s = (-deltaOne_y * (segm1_p1x - segm2_p1x) + deltaOne_x * (segm1_p1y - segm2_p1y)) / (-deltaTwo_x * deltaOne_y + deltaOne_x * deltaTwo_y)
t = ( deltaTwo_x * (segm1_p1y - segm2_p1y) - deltaTwo_y * (segm1_p1x - segm2_p1x)) / (-deltaTwo_x * deltaOne_y + deltaOne_x * deltaTwo_y)
if (s >= 0 and s <= 1) and (t >= 0 and t <= 1):
intersec_x = segm1_p1x + (t * deltaOne_x)
intersec_y = segm1_p1y + (t * deltaOne_y)
return True
else:
return False
def check_misma_linea(s1_x1, s1_y1, s1_x2, s1_y2, s2_x1, s2_y1, s2_x2, s2_y2):
s1x1, s1x2, s2x1, s2x2 = s1_x1, s1_x2, s2_x1, s2_x2
s1y1, s1y2, s2y1, s2y2 = s1_y1, s1_y2, s2_y1, s2_y2
if s1_x1 > s1_x2:
s1x1 = s1_x2
s1y1 = s1_y2
s1x2 = s1_x1
s1y2 = s1_y1
if s2_x1 > s2_x2:
s2x1 = s2_x2
s2y1 = s2_y2
s2x2 = s2_x1
s2y2 = s2_y1
pendiente_pol = (s1y2 - s1y1) / (s1x2 - s1x1)
pendiente_lin = (s2y2 - s2y1) / (s2x2 - s2x1)
if pendiente_pol == pendiente_lin:
return True
else:
return False
def check_vertice_repetido(listaPoligono, listaPoligonoDos, pointsLinea):
for vertexPol in range(len(listaPoligono) - 1):
s1_x1 = listaPoligono[vertexPol][0]
s1_y1 = listaPoligono[vertexPol][1]
s1_x2 = listaPoligono[vertexPol + 1][0]
s1_y2 = listaPoligono[vertexPol + 1][1]
for vertexLin in range(len(pointsLinea) - 1):
boolean_intersec = False
s2_x1 = pointsLinea[vertexLin][0]
s2_y1 = pointsLinea[vertexLin][1]
s2_x2 = pointsLinea[vertexLin + 1][0]
s2_y2 = pointsLinea[vertexLin + 1][1]
bool_pendiente = check_misma_linea(s1_x1, s1_y1, s1_x2, s1_y2, s2_x1, s2_y1, s2_x2, s2_y2)
if bool_pendiente == False:
aux_booleanIntersec = intersec_segment_4pts(s1_x1, s1_y1, s1_x2, s1_y2, s2_x1, s2_y1, s2_x2, s2_y2)
if aux_booleanIntersec == True:
lado_valorOne = calcula_lado_linea(s2_x1, s2_x2, s1_x1, s2_y1, s2_y2, s1_y1)
lado_valorTwo = calcula_lado_linea(s2_x1, s2_x2, s1_x2, s2_y1, s2_y2, s1_y2)
if lado_valorOne != lado_valorTwo:
if lado_valorOne != 0 and lado_valorTwo != 0:
vertexOne = [s1_x1, s1_y1]
vertexTwo = [s1_x2, s1_y2]
for vertexPolTwo in listaPoligonoDos:
if vertexOne == vertexPolTwo:
return vertexPol
elif vertexTwo == vertexPolTwo:
return vertexPol + 1
return -1
nombre_lineaCorte = "Linea_Corte"
capa_linea = QgsMapLayerRegistry.instance().mapLayersByName(nombre_lineaCorte)[0]
capa_lineaCorte_DP = QgsMapLayerRegistry.instance().mapLayersByName(nombre_lineaCorte)[0].dataProvider()
capa_manzana = QgsMapLayerRegistry.instance().mapLayersByName("Poligono")[0]
iface.setActiveLayer(capa_manzana)
capa_manzana_DP = QgsMapLayerRegistry.instance().mapLayersByName("Poligono")[0].dataProvider()
campos = [ ]
aux_campos = None
atributos = None
resultl = QgsVectorLayer("Polygon?crs=epsg:4326", "virtual_one", "memory")
resultpr = resultl.dataProvider()
resultl_dos = QgsVectorLayer("Polygon?crs=epsg:4326", "virtual_two", "memory")
resultpr_dos = resultl_dos.dataProvider()
puntos_poligono = [ ]
auxFeature_one = QgsFeature()
auxFeature_two = QgsFeature()
for feature in capa_manzana.selectedFeatures():
auxFeature_one = feature
auxFeature_two = feature
atributos = feature.attributes()
aux_campos = feature.fields()
for polygonFeat in feature.geometry().asMultiPolygon():
for lineaFeat in polygonFeat:
for pointFeat in lineaFeat:
puntos_poligono.append(pointFeat)
auxGeometry = QgsGeometry.fromPolygon([puntos_poligono])
auxFeature_one.setGeometry(auxGeometry)
auxFeature_two.setGeometry(auxGeometry)
resultpr.addFeatures([auxFeature_one])
for feature in resultpr.getFeatures():
geometria_original = QgsGeometry.fromPolygon(feature.geometry().asPolygon())
for line in capa_lineaCorte_DP.getFeatures():
intersect_boolean = feature.geometry().reshapeGeometry(line.geometry().asPolyline())
if (intersect_boolean == 0):
diferencia = QgsFeature()
diferencia.setGeometry(geometria_original.difference(feature.geometry()))
resultpr_dos.addFeatures([feature])
resultpr_dos.addFeatures([diferencia])
points_line = [ ]
for line in capa_linea.getFeatures():
for point in line.geometry().asPolyline():
points_line.append([point.x(), point.y()])
polygon_counter = 0
points_polygon_one = [ ]
points_polygon_two = [ ]
for feature in resultpr_dos.getFeatures():
polygon_counter += 1
for vertices in feature.geometry().asPolygon():
for point in vertices:
if polygon_counter == 1:
points_polygon_one.append([point.x(), point.y()])
else:
points_polygon_two.append([point.x(), point.y()])
vertice_repetido = None
index_polygon_one = -1
index_polygon_two = -1
index_polygon_one = check_vertice_repetido(points_polygon_one, points_polygon_two, points_line)
if index_polygon_one == -1:
index_polygon_two = check_vertice_repetido(points_polygon_two, points_polygon_one, points_line)
if index_polygon_one == -1 and index_polygon_two == -1:
pass
else:
if index_polygon_one == -1:
if index_polygon_two == 0:
points_polygon_two.pop()
points_polygon_two.pop(0)
aux_vertex = points_polygon_two[0]
points_polygon_two.append(aux_vertex)
else:
points_polygon_two.pop(index_polygon_two)
else:
if index_polygon_one == 0:
points_polygon_one.pop()
points_polygon_one.pop(0)
aux_vertex = points_polygon_one[0]
points_polygon_one.append(aux_vertex)
else:
points_polygon_one.pop(index_polygon_one)
list_points_one = [ ]
for i in points_polygon_one:
point = QgsPoint(i[0], i[1])
list_points_one.append(point)
list_points_two = [ ]
for i in points_polygon_two:
point = QgsPoint(i[0], i[1])
list_points_two.append(point)
resultl = QgsVectorLayer("Polygon?crs=epsg:4326", "virtual_one", "memory")
resultpr = resultl.dataProvider()
auxFeature_one = QgsFeature()
auxFeature_two = QgsFeature()
auxFeature_one.setGeometry(QgsGeometry.fromPolygon([list_points_one]))
auxFeature_two.setGeometry(QgsGeometry.fromPolygon([list_points_two]))
resultpr.addFeatures([auxFeature_one])
resultpr.addFeatures([auxFeature_two])
QgsMapLayerRegistry.instance().addMapLayer(resultl)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.