Oriented minimum bounding box con PyQGIS 3

Las oriented minimum bounding box ya se encuentran disponibles como oferta en QGIS 3 a través de la ToolBox de Processing. Sin embargo, los inconvenientes a la hora de acceder a los archivos temporales o en memoria pueden hacernos pensar porqué no tener nuestra propia rutina de cálculo. El fundamento del algoritmo puede ser visualizado aquí.

Con base en lo anterior se desarrolló el siguiente código:

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

registry = QgsProject.instance() polygon = registry.mapLayersByName('polygon1_convex_hull') feats_polygon = [ feat for feat in polygon[0].getFeatures() ] orig_pol = [ feat.geometry().asMultiPolygon() for feat in feats_polygon ] points_pol = [ item for item in orig_pol[0][0][0] ] points_pol = [ QgsPoint(point) for point in points_pol ] angles = [ 90 - points_pol[i].azimuth(points_pol[i+1]) for i in range(len(points_pol)-2) ] epsg = polygon[0].crs().postgisSrid() uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes" areas = [] for j, item in enumerate(angles): mem_layer = QgsVectorLayer(uri, 'rotated', 'memory') prov = mem_layer.dataProvider() feats = [ QgsFeature() for i in range(len(feats_polygon)) ] for i, feat in enumerate(feats): feat.setAttributes([i]) geom = feats_polygon[i].geometry() pt = geom.centroid().asPoint() geom.rotate(item, pt) feat.setGeometry(geom) prov.addFeatures(feats) # registry.addMapLayer(mem_layer) feats_mem = [ feat for feat in mem_layer.getFeatures() ] rotated_pol = [ feat.geometry().asMultiPolygon() for feat in feats_mem ] points = [ item for item in rotated_pol[0][0][0] ] x = [ point.x() for point in points ] y = [ point.y() for point in points ] #creating bounding box p1 = QgsPointXY(min(x), max(y)) p2 = QgsPointXY(max(x), min(y)) rect= QgsRectangle(p1, p2) new_geom = rect.asWktPolygon() areas.append([j, rect.area(), p1, p2]) #find minimum area by second column idx, min_area, p1, p2 = min(areas, key = lambda i : i[1]) new_mem_layer = QgsVectorLayer(uri, 'minimum bounding box', 'memory') prov = new_mem_layer.dataProvider() feats = [ QgsFeature() for i in range(len(feats_polygon)) ] rect= QgsRectangle(p1, p2) new_geom = rect.asWktPolygon() for i, feat in enumerate(feats): feat.setAttributes([i]) geom = feat.geometry().fromWkt(new_geom) geom.rotate(-angles[idx], pt) feat.setGeometry(geom) prov.addFeatures(feats) registry.addMapLayer(new_mem_layer)

La ejecución del código anterior en la Python Console de QGIS 3 produce la deseada bounding box.

Si desmarcamos la instrucción que permite visualizar todas las memory layers usadas para determinar las bounding box del vectorial rotado se obtiene lo siguiente:

De las nueve rotaciones posibles sólo la que ha rotado 16.676180725954737 grados es la que produce la bounding box visualizada (área mínima).