Espesor mínimo de un polígono no convexo con agujeros mediante PyQGIS

En este hilo de gis.stackexchange.com se consideran, aparentemente sin éxito o con un tiempo de cómputo excesivamente largo, diferentes procedimientos para determinar el espesor mínimo de un polígono no convexo con agujeros.

Me interesé por el tema porque desde hace tiempo quería desarrollar un algoritmo que, una vez generados puntos aleatorios dentro de la superficie de un polígono, determinará el punto más cercano a uno de los lados o vértice del polígono y luego "creciera" en sentido contrario hasta intersectar con el lado opuesto.

Como puede desprenderse de lo anterior, la primera intersección que se alcance será la que, para los puntos generados, le corresponda el espesor mínimo. Es obvio que los puntos aleatorios generados tienen que tener una buena distribución en todas las áreas del polígono; algo difícil de obtener con los métodos por defecto en los SIGs pero la imagen a continuación sólo refleja el intento de probar la validez del algoritmo.

El código a continuación contiene el intento de determinar el segmento más cercano por contexto a la periferia del polígono (que puede ser también un agujero) y desde el punto comenzar a "crecer" en el sentido opuesto basado en los cosenos directores. El crecimiento se hace con base a una distancia fija en un loop hasta que alguno de esos segmentos intersecte nuevamente con el polígono. Es de hacer notar que la distancia de crecimiento no puede ser muy elevada (algo que se nos puede ocurrír para disminuir el tiempo de cómputo) porque podría dar falsos positivos.

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

import math def azimuth(point1, point2): return point1.azimuth(point2) #in degrees def cosdir_azim(azim): azim = math.radians(azim) cosa = math.sin(azim) cosb = math.cos(azim) return cosa,cosb registry = QgsMapLayerRegistry.instance() polygon = registry.mapLayersByName('polygon_with_holes') point_layer = registry.mapLayersByName('Random_points') points = [ feat.geometry().asPoint() for feat in point_layer[0].getFeatures() ] feat_polygon = polygon[0].getFeatures().next() #producing rings polygons rings_polygon = feat_polygon.geometry().asPolygon() segments = [] epsg = point_layer[0].crs().authid() uri = "LineString?crs=" + epsg + "&field=id:integer""&index=yes" mem_layer = QgsVectorLayer(uri, 'increasing_segments', 'memory') prov = mem_layer.dataProvider() length = 10 pt2 = 0 k = 0 while pt2 == 0: for i, point in enumerate(points): #determining closest distance to vertex or side polygon dist1 = feat_polygon.geometry().closestSegmentWithContext(point)[0] #determining point with closest distance to vertex or side polygon pt = feat_polygon.geometry().closestSegmentWithContext(point)[1] cosa, cosb = cosdir_azim(azimuth(pt, point)) #extending segment in opposite direction based in director cosine and length op_pt = QgsPoint(point.x() + (length*cosa), point.y() + (length*cosb)) segments.append([pt,op_pt]) geom = QgsGeometry.fromPolyline([point,op_pt]) for ring in rings_polygon: geom_ring = QgsGeometry.fromPolyline(ring) if geom.intersects(geom_ring): pt3 = geom.intersection(geom_ring) pt2 = pt3.distance(QgsGeometry.fromPoint(point)) ms = [pt3.asPoint(), pt] length += 100 k += 1 new_segments = segments[len(segments) -1 - len(segments)/k: len(segments) - 1] feats = [ QgsFeature() for i in range(len(new_segments)) ] for i,feat in enumerate(feats): feat.setAttributes([i]) geom = QgsGeometry.fromPolyline(new_segments[i]) feat.setGeometry(geom) prov.addFeatures(feats) QgsMapLayerRegistry.instance().addMapLayer(mem_layer) minimum_segment = QgsGeometry().fromPolyline(ms).exportToWkt() print minimum_segment, k

Después de ejecutado el código anterior se obtiene el resultado de la imagen siguiente. Puede observarse que el primer segmento creciente alcanza el lado opuesto del polígono en el área esperada (área delgada central).

Un zoom a la zona referida y la utilización del plugin QuickWKT permiten corroborar lo anterior. No obstante, se puede observar una ligera inclinación del segmento (no es perpendicular). Ello se debe a que el "segmento más cercano con contexto" se obtuvo para un vértice; no para un lado del polígono. Esto puede arreglarse con una excepción o generando más puntos.