Shapely: Unexpected triangulation results for a concave shape

[kuanb]

When triangulating a concave shape, I would expect that the edges of the shape would be preserved such that the outside edges of the geometry would become vertices for 1 or more of the resulting geometries.

If this assumption is incorrect, is there a method or parameter that could force this to be true for the result? Is that easily achievable?

Example plots:

Code to create the same situation (using same below geometry):

triangles = MultiPolygon(triangulate(largest, tolerance=0.001))

Here's the WKT for this example geometry:

POLYGON ((-74.05644319847762 4.664371152795165, -74.05701264773319 4.663503533579181, -74.05770573357918 4.662810447733186, -74.05896428283818 4.662056102337443, -74.05990224838993 4.661771573597983, -74.06224145161008 4.661772473597984, -74.06317941716183 4.662057002337444, -74.06443796642083 4.662811347733187, -74.06572065226682 4.664052233579182, -74.06921901369725 4.668360960588712, -74.07141674761461 4.670691972117472, -74.07635895116509 4.673818151938486, -74.07894493390593 4.675834266094067, -74.08192435226682 4.679424333579181, -74.08891615226682 4.688383433579182, -74.08958587724112 4.689463067989243, -74.09086467349047 4.690719228823506, -74.09790275116509 4.694460551938487, -74.10034036642082 4.696114147733187, -74.10386724657296 4.698958762835078, -74.10814346936803 4.700870863662334, -74.10930545161006 4.700957773597984, -74.11043741716183 4.701320402337444, -74.11139045116509 4.701828051938487, -74.11214813390593 4.702449866094067, -74.11445885226682 4.704984433579183, -74.11521319766256 4.706242982838176, -74.11590442640203 4.70845234838992, -74.11611382363337 4.710865185701647, -74.11554750632175 4.71273208368413, -74.11467343390593 4.713910633905932, -74.11305131716183 4.714994497662556, -74.11211335161008 4.715279026402015, -74.11073328554708 4.715355222566556, -74.11041129766257 4.716545717161825, -74.1094660522668 4.718057466420818, -74.10756213390592 4.720161033905933, -74.10568445116508 4.721873348061513, -74.10420772338625 4.722584201678661, -74.10275629999998 4.7227995, -74.10130487661371 4.722584201678661, -74.09995094883489 4.721945648061514, -74.09403374773316 4.716245366420819, -74.09338219367824 4.71525808368413, -74.09288787359796 4.71387655161008, -74.09288787359796 4.711925648389919, -74.0936469519385 4.710096948834901, -74.09548963357916 4.708043347733187, -74.09763068874044 4.70677727898599, -74.09460666784895 4.704312469919607, -74.09203589291714 4.702595752799392, -74.08765371631586 4.700417706321742, -74.08459603357916 4.698508252266813, -74.08273894773316 4.696641066420818, -74.08206922275886 4.695561432010757, -74.08118604773318 4.694727366420818, -74.0717936025263 4.682825210334004, -74.06562484883489 4.678834648061513, -74.06390001213198 4.677314420684328, -74.0535710885149 4.739037649974117, -74.05520778283815 4.737937702337444, -74.05712119999998 4.7375571, -74.06045152338628 4.738003698321339, -74.06217206642083 4.738923347733187, -74.06340970632175 4.74043141631587, -74.06393542640203 4.74187004838992, -74.06400742363337 4.743335585701648, -74.063537026402 4.74698295161008, -74.06279044806149 4.748785351165098, -74.06139073691359 4.75018144712255, -74.0609743487264 4.751890490745711, -74.06128649999999 4.753612599999999, -74.06107120167866 4.755064023386272, -74.06039964806149 4.756487551165098, -74.0594142664208 4.757574752266813, -74.05806621716184 4.758372197662556, -74.05615280000001 4.7587528, -74.05377854838993 4.758482026402016, -74.0507277649826 4.757844764394359, -74.05055722640202 4.75945045161008, -74.04983964400328 4.762451751932777, -74.04873032640204 4.76925805161008, -74.04829449766255 4.770818617161826, -74.04766074258318 4.772140089033794, -74.04726809766255 4.773711517161825, -74.04651375226682 4.774970066420818, -74.04542655116509 4.775955448061513, -74.04410012338627 4.776582801678661, -74.04215861429836 4.776774023633361, -74.04029171631588 4.776207706321742, -74.03878364773318 4.774970066420818, -74.03800540233745 4.773666917161825, -74.03729177636662 4.771281485701648, -74.03727967636664 4.770004014298352, -74.03747549832134 4.769002476613728, -74.03843161386389 4.767199327773277, -74.04047532279198 4.75512284897233, -74.0369408023945 4.754256982815782, -74.03025007661371 4.753148101678661, -74.02876484883491 4.752415548061514, -74.02767764773319 4.751430166420819, -74.02692330233745 4.750171617161826, -74.02654087636664 4.748620285701648, -74.02683367636664 4.736912214298353, -74.02739999367826 4.73504531631587, -74.02863763357918 4.733537247733187, -74.03035817661373 4.732617598321339, -74.03229968570164 4.732426376366639, -74.03372301716183 4.732782902337444, -74.03534513390593 4.733866766094067, -74.03621920632175 4.73504531631587, -74.03676491551958 4.736773285680133, -74.04355019487681 4.738452884397855, -74.04767016881966 4.714185707761746, -74.04885077359798 4.706409648389919, -74.05027314410287 4.698677293675041, -74.05045487359796 4.696561148389919, -74.05077748485968 4.695365046361829, -74.05289107166563 4.682999393176747, -74.05422957359798 4.67404324838992, -74.05473458359565 4.671889824337649, -74.05543207359798 4.667273848389919, -74.05644319847762 4.664371152795165))

Shapely-1.6.1
Python 3.6

[sgillies]

@kuanb I don't have a ton of experience with the triangulation algorithm in GEOS. This result surprises me, too. From what I see in https://postgis.net/docs/ST_DelaunayTriangles.html (also based on GEOS), this might not be unexpected.

[dr-jts]

Delaunay Triangulation does not respect edges - it only takes into account points. To respect edges you need to use Constrained or Conforming Delaunay Triangulation. JTS has a Conforming implementation; not sure of the status of this in GEOS.

[sgillies]

🙏 @dr-jts!

[dr-jts]

Nice plot, by the way! JTS TestBuilder needs something like this... :)

[dr-jts]

Here's a image of the Conforming Delaunay computed by JTS.

[dr-jts]

Just to beat this one into the ground, note that:

• A Conforming Delaunay is a valid DT, with sites added along constraint edges to approximate the shape of the constraints
• A Constrained Delaunay has the constraint edges as edges of the triangulation, but is potentially NOT a valid DT

[kuanb]

Thanks so much for beating it into the ground! I appreciated the additional information. Any chance you are aware of a Python implementation of conforming? I've been using https://pypi.python.org/pypi/tri/0.2 which only performs constrained. I suppose I could also follow along with the JTS triangulation logic and attempt to implement it in Py. Thanks again.

[dr-jts]

Sorry, I don't keep up with the Python spatial world much. I'm impressed that there is a constrained DT implementation!

It's possible that the JTS CDT is exposed in GEOS, in which case you could link to that. Porting the JTS logic might be a challenge - I'm all too aware that it is tricky code to get right.

[dr-jts]

If you're wanting to triangulate just the interior of the polygon, you could look at a Polygon Triangulation algorithm (which is very different to DT, surprising though that may seem). The classic algorithm is Ear Clipping:

In general, PT-EC is not Delaunay. However, there is the possibility of doing Delaunay Refinement of a Polygon Triangulation:
(Unfortunately this code has not yet appeared in JTS).

[mhtb32]

Conforming triangulation has been put on Milestone 3.9 of GEOS. Maybe Shapely can have this after GEOS implements it?