1

I am interested in using OpenStreetMaps building data in PyQGIS to generate triangulated VTK or STL files. I am able to get the features I want as QgsGeometry objects. The problem is that these geometries are often concave, and often have holes (interior rings), so delaunayTriangulation and voronoiDiagram have not produced good results.

So for example, directly calculating the area of the GEOS geometry directly gives a good value. Calculating the edge-constrained Delaunay or Voronoi fills returns an empty geometry (area = 0). Calculating the un-constrained Delaunay or Voronoi fills in the concave regions and holes.

for f in features:
    geom = f.geometry()
    geom.transform(qgis.core.QgsCoordinateTransform(
        qgis.core.QgsCoordinateReferenceSystem('EPSG:4326'),
        qgis.core.QgsCoordinateReferenceSystem('EPSG:27572'),
        qgis.core.QgsProject.instance()))
    print("True area = %.3f"%(geom.area()))

    delaunay_geom = geom.delaunayTriangulation(0.1, True)
    print("Delaunay constrained area = %.3f"%(delaunay_geom.area()))

    delaunay_geom = geom.delaunayTriangulation(0.1, False)
    print("Delaunay unconstrained area = %.3f"%(delaunay_geom.area()))

    voronoi_geom = geom.voronoiDiagram(qgis.core.QgsGeometry(), 0.1, True)
    print("Voronoi constrained area = %.3f"%(voronoi_geom.area()))

    voronoi_geom = geom.voronoiDiagram(qgis.core.QgsGeometry(), 0.1, False)
    print("Voronoi unconstrained area = %.3f"%(voronoi_geom.area()))

The above code block returns these values, which suggests that the Delaunay and Voronoi are not producing something I can use.

> True area = 6912.680
> Delaunay constrained area = 0.000
> Delaunay unconstrained area = 8819.325
> Voronoi constrained area = 0.000
> Voronoi unconstrained area = 150244.332

I have tried extracting the points and exterior/interior rings, and doing the triangulation myself.

for f in features:
    geom = f.geometry()

    point_index = 0
    rings = []
    vertices = []
    for part in geom.parts():
        for n in range(part.childCount()):
            child = part.childGeometry(n)
            rings.append([])
            for j in range(child.childCount()):
                vertices.append([child.childPoint(j).x(), child.childPoint(j).y(), 0])
                rings[-1].append(point_index)
                point_index = point_index + 1
    print(rings)

This output corresponds to the large building in the attached pictures with 4 holes. It has the exterior loop, and 4 interior loops.

[[0, 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]]

Doing this triangulation manually is a very hard problem, and I'm hoping PyQGIS has a better solution. I can see in the attached QGIS GUI rendering that the building geometries are filled in correctly -- the holes are holes, and the concave regions are respected. I can also see that calculating the area by feature.geometry().area() is giving good results. So I assume that somehow, QGIS must be splitting the geometry into well-behaved polygons to do the rendering and area calculations. Is there any way to access more well-behaved versions of the polygons, or convert the feature geometry into simpler (convex) polygons or triangles?

QGIS render view of concave building geometry with holes from OSM

0

1 Answer 1

1

I ended up finding the solution to this myself. The first trick is to densify the geometry before running delaunayTriangulation. This helps force delaunayTriangulation to respect holes and concave regions. But note that this triangulated geometry will be convex and have no holes, which is solved in the second step.

for f in features:
    geom = f.geometry()
    geom.transform(qgis.core.QgsCoordinateTransform(
        qgis.core.QgsCoordinateReferenceSystem('EPSG:4326'),
        qgis.core.QgsCoordinateReferenceSystem('EPSG:27572'),
        qgis.core.QgsProject.instance()))

    delaunay_geom = geom.densifyByDistance(2.0)
    delaunay_geom = delaunay_geom.delaunayTriangulation(1e-6, False)

The second trick is to loop over all the triangulated cells, and delete any whose centroids are not contained within the original geometry. This step recovers the concave regions and the holes. I had to delete the triangulated cells in the reverse order, because it seemed like the parts list resized whenever a cell was deleted.

    part_delete_list = []
    for part in delaunay_geom.parts():
        centroid = part.centroid()
        centroid_xy = qgis.core.QgsPointXY(centroid.x(), centroid.y())
        is_inside = geom.contains(qgis.core.QgsGeometry.fromPointXY(centroid_xy))
        part_delete_list.append(not is_inside)

    for part_index in reversed(range(len(part_delete_list))):
        if part_delete_list[part_index]:
            delaunay_geom.deletePart(part_index)

After that, you can extract the points and triangles, and generate a VTK or STL mesh using your favorite VTK or STL library.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.