Division de polígonos sin errores topológicos en QGIS usando Python

divide_poligono_qgis.py

# -*- 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)