How to create a Delaunay triangulation with SciPy

Delaunay triangulation splits a coordinate set into non-overlapping triangles that follow local point neighborhoods. In SciPy, it helps Python spatial workflows prepare meshes, inspect geometry, or locate which triangle contains a test coordinate before interpolation or plotting.

scipy.spatial.Delaunay sends a floating-point coordinate array to Qhull and returns indexes back into that original array. In 2-D, every row of triangulation.simplices contains the three point indexes that make one triangle, while find_simplex() reports which triangle contains each query point.

A small point set makes the triangulation easy to check before applying the same pattern to larger data. Guard the find_simplex() value before indexing simplices because an outside point returns -1, and Python treats -1 as the last row if it is used as an index.

1. Create a Python script with a small 2-D point set.

   delaunay_demo.py

   import numpy as np
   from scipy.spatial import Delaunay
    
    
   points = np.array(
       [
           [0.0, 0.0],
           [2.0, 0.0],
           [2.0, 1.2],
           [0.0, 1.0],
           [0.9, 0.4],
           [1.3, 0.8],
       ]
   )
    
   triangulation = Delaunay(points)
    
   print(f"input points: {len(points)}")
   print(f"triangles: {len(triangulation.simplices)}")
   print("simplices:")
   for simplex_index, simplex in enumerate(triangulation.simplices):
       vertices = ", ".join(str(int(vertex)) for vertex in simplex)
       print(f"  {simplex_index}: [{vertices}]")
    
   first_simplex = triangulation.simplices[0]
   first_coordinates = [
       [float(f"{points[vertex, 0]:.1f}"), float(f"{points[vertex, 1]:.1f}")]
       for vertex in first_simplex
   ]
   print("first triangle coordinates:")
   for coordinate in first_coordinates:
       print(f"  {coordinate}")
    
   query_points = np.array(
       [
           [0.2, 0.2],
           [1.1, 0.6],
           [2.3, 0.5],
       ]
   )
    
   print("point lookup:")
   for point, simplex_index in zip(query_points, triangulation.find_simplex(query_points)):
       x, y = point
       if simplex_index == -1:
           print(f"  ({x:.1f}, {y:.1f}) -> outside triangulation")
           continue
       vertices = triangulation.simplices[simplex_index].tolist()
       print(f"  ({x:.1f}, {y:.1f}) -> simplex {int(simplex_index)} vertices {vertices}")

   The point array is shaped as rows of x, y coordinates. Each row in triangulation.simplices stores indexes back into that array, not copied coordinate values.

2. Run the script to build the triangulation and print the simplex rows.

   $ python delaunay_demo.py
   input points: 6
   triangles: 6
   simplices:
     0: [1, 4, 0]
     1: [4, 3, 0]
     2: [2, 5, 1]
     3: [5, 4, 1]
     4: [3, 5, 2]
     5: [5, 3, 4]
   first triangle coordinates:
     [2.0, 0.0]
     [0.9, 0.4]
     [0.0, 0.0]
   point lookup:
     (0.2, 0.2) -> simplex 1 vertices [4, 3, 0]
     (1.1, 0.6) -> simplex 5 vertices [5, 3, 4]
     (2.3, 0.5) -> outside triangulation

   The simplex row order is the order returned by SciPy for this input. Use the row indexes reported by find_simplex() against the same triangulation.simplices array instead of assuming a geometric ordering.

3. Match simplex rows back to coordinates before passing triangles to later code.

   Simplex 0 uses point indexes [1, 4, 0], which map to [2.0, 0.0], [0.9, 0.4], and [0.0, 0.0] in the original points array.

4. Guard outside lookup results before indexing simplices.

   find_simplex() returns -1 for a point outside the triangulation. Check for -1 first, because triangulation.simplices[-1] is a valid Python lookup for the last simplex row.

5. Check for omitted input points when the data contains duplicates or geometric degeneracy.

   Qhull can omit points from the triangulation in degenerate cases unless a resolving option such as QJ is chosen deliberately. Inspect triangulation.coplanar when duplicate or nearly collinear points may be present.

6. Remove the demo script when the check is complete.

   $ rm delaunay_demo.py