Triangulation produces wrong results

[GeorgySk]

For some reason shapely.ops.triangulate doesn't like some geometries and loses their parts.

Simple example:

from shapely.geometry import Polygon
from shapely.ops import triangulate

p = Polygon([(0, 0), (1, 0), (-1, 0.05)])
triangulate(p)
# []

And another one:

p = Polygon([(0, 0), (0, 486), (1, 486), (1, 22), (2, 22), (2, 0)])
for x in triangulate(p):
    print(x.wkt)
# POLYGON ((0 486, 1 22, 1 486, 0 486))
# POLYGON ((1 486, 1 22, 2 22, 1 486))
# POLYGON ((0 0, 2 0, 1 22, 0 0))
# POLYGON ((1 22, 2 0, 2 22, 1 22))

There is a missing triangle here.
But if I replace 486 by, for example, 200, then it will work fine:

p = Polygon([(0, 0), (0, 200), (1, 200), (1, 22), (2, 22), (2, 0)])
for x in triangulate(p):
    print(x.wkt)
# POLYGON ((0 200, 0 0, 1 22, 0 200))
# POLYGON ((0 200, 1 22, 1 200, 0 200))
# POLYGON ((1 200, 1 22, 2 22, 1 200))
# POLYGON ((0 0, 2 0, 1 22, 0 0))
# POLYGON ((1 22, 2 0, 2 22, 1 22))

---

Shapely version: 1.6.4.post1, installed from conda.

[drnextgis]

In the first example we have just a simple triangle so getting an empty list is expected behavior, isn't it?

Second and third examples also give us a valid Delaunay triangulation:

[GeorgySk]

@drnextgis shapely.ops.triangulate returns Delaunay triangulation of vertices of a given geometry. And one of the properties of Delaunay triangulation is that the union of all resulting triangles equals the convex hull of given points. So, when giving this function a triangle polygon, I think I'm right to expect the triangle itself as a result.

And even if it would be a default behavior to return empty lists then it is definitely not consistent:

from shapely.geometry import Polygon
from shapely.ops import triangulate

# this is broken:
p = Polygon([(0, 0), (1, 0), (-1, 0.05)])
triangulate(p)
# []

# but this works fine:
p = Polygon([(0, 0), (1, 0), (-1, 0.1)])
triangulate(p)[0].wkt
# 'POLYGON ((-1 0.1, 0 0, 1 0, -1 0.1))'

---

Regarding the second and the third example. In your plot there is a missing line, hence missing triangle:

Taking in mind the aforementioned property with the convex hull, we can also perform a simple check of resulting areas:

# there is a bug here:
p = Polygon([(0, 0), (0, 486), (1, 486), (1, 22), (2, 22), (2, 0)])
p.convex_hull.area, sum(triangle.area for triangle in triangulate(p))
# (740.0, 497.0)

# this works fine though:
p = Polygon([(0, 0), (0, 200), (1, 200), (1, 22), (2, 22), (2, 0)])
p.convex_hull.area, sum(triangle.area for triangle in triangulate(p))
# (311.0, 311.0)

[dbaston]

Confirmed for both MULTIPOINT ((0 0), (1 0), (-1 0.05), (0 0)) and MULTIPOINT ((0 0), (0 486), (1 486), (1 22), (2 22), (2 0)) in JTS master @dr-jts

[dr-jts]

This appears to be fixed in latest GEOS (by libgeos/geos#728). Can Shapely confirm?

[jorisvandenbossche]

Yes, I can confirm this is fixed with latest released GEOS and shapely. Thanks for the note!

First example no longer gives an empty list:

>>> p = Polygon([(0, 0), (1, 0), (-1, 0.05)])
>>> triangulate(p)
[<POLYGON ((-1 0.05, 0 0, 1 0, -1 0.05))>]

And the second example also now gives all expected triangles:

>>> p = Polygon([(0, 0), (0, 486), (1, 486), (1, 22), (2, 22), (2, 0)])
>>> for x in triangulate(p):
...     print(x.wkt) 
POLYGON ((0 486, 0 0, 1 22, 0 486))
POLYGON ((0 486, 1 22, 1 486, 0 486))
POLYGON ((1 486, 1 22, 2 22, 1 486))
POLYGON ((0 0, 2 0, 1 22, 0 0))
POLYGON ((1 22, 2 0, 2 22, 1 22))

and confirming the above result visually:

[dr-jts]

Should be fixed by libgeos/geos#953