Proyectar la longitud de polígonos sobre segmentos de rectas con PyQGIS 3

Para proyectar la longitud de un polígono sobre una recta es necesario encontrar el ángulo que ésta mantiene con el eje X (será negativo si la pendiente es negativa), rotar cada rasgo ese ángulo con centro de rotación en el centroide, calcular la bounding box del rasgo rotado y, finalmente para comprobar los resultados, rotar la bounding box en sentido contrario por el mismo valor del ángulo (ángulo negativo) y con centro de rotación en el centroide del feature (no de la bounding box).

La situación que se va a probar se encuentra a continuación:

Consta de un polígono con dos features que se van a proyectar alternativamente sobre los segmentos de recta my_line y my_line2. El código desarrollado para tal fin se encuentra 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

registry = QgsProject.instance() my_line = registry.mapLayersByName('my_line') polygon1 = registry.mapLayersByName('polygon1') for feature in my_line[0].getFeatures(): fgeom = feature.geometry() # check if feature geometry is multipart if fgeom.isMultipart(): print("yes multipart") new_features = [] temp_feature = QgsFeature(feature) # set the geometry of each part for part in fgeom.asGeometryCollection(): temp_feature.setGeometry(part) new_features.append(temp_feature) slines = [ feat.geometry().asPolyline() for feat in new_features ] else: print("not multipart") slines = [ feature.geometry().asPolyline() for feature in my_line[0].getFeatures() ] points = [ point for line in slines for point in line ] azimuth = points[0].azimuth(points[1]) angle = 90 - azimuth feats_polygon = [ feat for feat in polygon1[0].getFeatures() ] epsg = polygon1[0].crs().postgisSrid() uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer&field=width:double""&index=yes" mem_layer = QgsVectorLayer(uri, 'rotated bounding box', 'memory') prov = mem_layer.dataProvider() feats = [ QgsFeature() for i in range(len(feats_polygon)) ] for i, feat in enumerate(feats): geom = feats_polygon[i].geometry() pt = geom.centroid().asPoint() geom.rotate(angle, pt) new_geom = QgsGeometry.fromWkt(geom.boundingBox().asWktPolygon()) width = new_geom.boundingBox().width() feat.setAttributes([i, width]) new_geom.rotate(-angle, pt) feat.setGeometry(new_geom) prov.addFeatures(feats) registry.addMapLayer(mem_layer)

Cuando se ejecuta el código anterior con base en my_line (en rojo) se obtiene el resultado de la imagen siguiente:

Si se toma como referencia my_line2 (en azul) el resultado de ejecución se presenta de la manera siguiente:

En la tabla de atributos de las bounding box rotadas se encuentra como 'width' el valor projectado sobre la recta considerada.