Dividiendo polígonos irregulares en dos partes con diferentes proporciones y con diferente orientación de la línea de corte con PyQGIS

Un requerimiento común en las opciones de geoproceso es la división de polígonos irregulares en diferentes proporciones y con diferente orientación de la línea de corte. Aunque parece difícil, la opción más sencilla para realizar ésto es mediante la bounding box del polígono rotado por el ángulo correspondiente a la línea de corte para la proporción deseada.

Supongamos que queremos dividir un polígono en dos partes; una de ellas representa el 35 % del área total y la otra el 65 % restante. Además, hay que considerar la orientación de la línea de corte, por ejemplo, 45 grados en sentido horario.

El algoritmo a emplear considera entonces la bounding box del polígono rotado con la orientación de -45 grados (los ángulos negativos corresponden a la orientación horaria en PyQGIS) y realiza el corte de la porción occidental del polígono con la proporción requerida con base en el área total de éste. Las diferencias se van corrigiendo moviendo la línea de corte de oeste a este y de este a oeste hasta lograr convergencia con base en una tolerancia < 1x10^-5 % con base en el área total. Un procedimiento similar es realizado con la porción oriental y las dos partes resultantes se combinan en una memory layer con los dos rasgos correspondientes. Es lo que refleja el código que se presenta 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 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140

import math registry = QgsProject.instance() #layer = registry.mapLayersByName('polygon1') layer = registry.mapLayersByName('polygon2') angle = -45 feats = [ feat for feat in layer[0].getFeatures() ] pol_geom = feats[0].geometry() area = pol_geom.area() area2 = pol_geom.area() pt = pol_geom.centroid().asPoint() pol_geom.rotate(angle, pt) bbox = pol_geom.boundingBox() pol = bbox.asWktPolygon() xmin = bbox.xMinimum() xmax = bbox.xMaximum() ymin = bbox.yMinimum() ymax = bbox.yMaximum() xmean = (xmin + xmax)/2. xmean2 = (xmin + xmax)/2. width = bbox.width() increment2 = width/20 xold2 = xmean2 new_area2 = 0 old_area2 = 0 width = bbox.width() increment = width/20 xold = xmean n = 1 new_area = 0 old_area = 0 int1 = 0 int2 = 0 while(n > 0): left_half = QgsRectangle(xmin, ymin, xmean, ymax).asWktPolygon() left_half_geom = QgsGeometry.fromWkt(left_half) half_intersect = left_half_geom.intersection(pol_geom) half_intersect_area = half_intersect.area() right_half = QgsRectangle(xmean2, ymin, xmax, ymax).asWktPolygon() right_half_geom = QgsGeometry.fromWkt(right_half) half_intersect2 = right_half_geom.intersection(pol_geom) half_intersect_area2 = half_intersect2.area() if n == 500: break if int1 == 1 and int2 == 1: break fact = 0.35 fact2 = 1 - fact if area*fact > half_intersect_area and int1 == 0: xold = xmean xmean += increment old_area = half_intersect_area var = new_area-old_area if area*fact < half_intersect_area and int1 == 0: xmean = xold increment = increment/2 new_area = half_intersect_area var = new_area-old_area if new_area-old_area < 1e-5: int1 = 1 if int1 == 1: print("error1", new_area-old_area) if n == 1 and var > 0 and int1 == 0: xmean = xmin + increment if area2*fact2 > half_intersect_area2 and int2 == 0: xold2 = xmean2 xmean2 -= increment2 old_area2 = half_intersect_area2 var2 = new_area2-old_area2 if area2*fact2 < half_intersect_area2 and int2 == 0: xmean2 = xold2 increment2 = increment2/2 new_area2 = half_intersect_area2 var2 = new_area2-old_area2 if new_area2-old_area2 < 1e-5: int2 = 1 if int2 == 1: print("error2", new_area2-old_area2) if n == 1 and var2 > 0 and int2 == 0: xmean2 = xmax - increment2 n += 1 geoms = [] new_geom = QgsGeometry.fromWkt(half_intersect.asWkt()) new_geom.rotate(-angle, pt) geoms.append(new_geom) new_geom2 = QgsGeometry.fromWkt(half_intersect2.asWkt()) new_geom2.rotate(-angle, pt) geoms.append(new_geom2) epsg = layer[0].crs().postgisSrid() uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes" mem_layer = QgsVectorLayer(uri, 'divided_polygon', 'memory') prov = mem_layer.dataProvider() feats = [ QgsFeature() for i in range(len(geoms)) ] for i, feat in enumerate(feats): feat.setAttributes([i]) feat.setGeometry(geoms[i]) prov.addFeatures(feats) registry.addMapLayer(mem_layer)

Después de ejecutado el código anterior en la Python Console de QGIS el polígono original muestra el resultado a continuación. Ha sido dividido en dos partes. La occidental representa el 35 % del área total y la oriental el 65 %. La orientación de la línea de corte corresponde a 45 grados; tal como se esperaba.