[Python] Segmentation fault when using 'triangle' library ONLY when processing multiple...

I am trying to process shapely polygon geometries and decompose them into triangles with the constraint that triangle vertices can only exist on the polygon edge(s). For some strange reason the code works most of the time, but then I have a specific polygon, that results in a segmentation fault

(Process finished with exit code 139 (interrupted by signal 11:SIGSEGV): POLYGON ((545155.9709965222 5378933.039350484, 545155.8786686319 5378932.570194424, 545153.9828366871 5378934.823291995, 545154.3704756413 5378934.7733652275, 545155.9709965222 5378933.039350484))
​

and this ONLY when I run the code for the entire geodataframe containing all geometries to be triangulated. I am at a loss what the issue is, since it works if I execute the code line by line for this specific polygon, but keeps crashing if run within my loop. I am running python 3.12 and use triangle library version 20230923.

Also, if someone has a workaround in python I would be very grateful! Unfortunately, the shapely triangulation does not have any constraints as far as I understand, so I cannot use it for my purpose.

My code:

import os
import geopandas as gpd
from shapely.geometry import Point, MultiPolygon, Polygon, GeometryCollection
from shapely.validation import explain_validity
import shapely.geometry as geom
import numpy as np
import triangle as tr


def is_triangle(geometry):
if isinstance(geometry, Polygon):
# For a single Polygon, check if it has 3 vertices
return len(list(geometry.exterior.coords)) == 4 # 4 because the last point is the same as the first (closed polygon)
elif isinstance(geometry, MultiPolygon):
# For MultiPolygon, check if each polygon is a triangle
return all(len(list(poly.exterior.coords)) == 4 for poly in geometry.geoms)
else:
return False # In case of unexpected geometry types


def process_geometry(geometry):
# Check if the geometry is valid
if not geometry.is_valid:
print(f"Invalid geometry: {explain_validity(geometry)}")
geometry = geometry.buffer(0) # Try fixing the geometry by buffering with 0

if is_triangle(geometry):
return [geometry] # Return as a list for consistency
elif isinstance(geometry, Polygon):
return triangulate_polygon(geometry) # Triangulate and return the list of triangles
elif isinstance(geometry, MultiPolygon):
triangles = []
for poly in geometry.geoms:
triangles.extend(triangulate_polygon(poly)) # Triangulate each part of the MultiPolygon
return triangles
else:
raise ValueError("Geometry type not supported")

def triangulate_polygon(geometry):
# Remove near-duplicate vertices to prevent errors caused by floating-point precision issues:
cleaned_coords = remove_near_duplicate_vertices(np.array(geometry.exterior.coords))
polygon = Polygon(cleaned_coords)

coords = np.array(polygon.exterior.coords)
# coords = np.array(geometry.exterior.coords)

# Check for Degenerate or Collinear Segments
for i in range(len(coords) - 1):
dist = np.linalg.norm(coords - coords[i + 1])
if dist < 1e-8: # Check if distance is too small
print(f"Nearly degenerate segment between points {i} and {i + 1}")

# Make sure there are enough unique points for triangulation
if len(coords) < 3:
return []
segments = [[i, (i + 1) % len(coords)] for i in range(len(coords) - 1)]

# Ensure that no index is out of bounds
for segment in segments:
if any(idx >= len(coords) for idx in segment):
print(f"Invalid segment index {segment}")



# Perform triangulation using the triangle library
triangulation = tr.triangulate({'vertices': coords, 'segments': segments}, 'p')

triangles = create_triangles_from_triangulation(triangulation)

return triangles


def create_triangles_from_triangulation(triangulation):
"""Convert the triangulation output into Shapely Polygon objects."""
triangles = []
if 'triangles' in triangulation:
coords = triangulation['vertices']
for tri_indices in triangulation['triangles']:
triangle_coords = coords[tri_indices]
triangles.append(geom.Polygon(triangle_coords))
return triangles



all_triangles = []

for id, tri in non_triangles_gdf.iterrows():
print(tri)
print(id)
print(tri['geometry'])
geom = tri['geometry']
all_triangles.extend(process_geometry(geom))

Continue reading...