2

I have a TIN model where I remove some of the triangles (see fig. 1). Most of the triangles touch on one edge, some on one vertice. I would now like to generate a polygon shapefile containing the remaining triangles, dissolved into one polygon (of multiple rings, see fig. 2).

Fig 1. - TIN with some triangles removed. Vertices of touching triangles are exactly identical, no sliver polygons.

Fig. 2 - Note, that this is only one polygon with a single line in the attribute table

Pictures were made by using the QGIS "Vector-->Geoprocessing Tools-->Dissolve" tool. Is this possible with the ogr python bindings? Here's what I've tried so far:

from osgeo import ogr

tri = ogr.Geometry(ogr.wkbLinearRing)
tri.AddPoint(0,0)
tri.AddPoint(0,1)
tri.AddPoint(1,1)
tri.CloseRings()
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(tri)

tri2 = ogr.Geometry(ogr.wkbLinearRing)
tri2.AddPoint(1,1)
tri2.AddPoint(0,1)
tri2.AddPoint(0,4)
poly2 = ogr.Geometry(ogr.wkbPolygon)
poly2.AddGeometry(tri2)

poly.Union(poly2) # returns None
poly.Intersect(poly2) # returns False
poly.Touches(poly2) # returns False (?)
poly.Disjoint(poly2) # returns False

I'm missing something like poly.Dissolve(), maybe used on a ogr.wkbGeometryCollection?

  • 1
    I guess you just forgot to close the second ring, in your case with tri2.CloseRings() , but the method you are looking for is ogr.Geometry.Union() as you tried in your code (or maybe UnionCascaded()). Anyway with second ring closed, your example work perfectly. In [3]: union_poly = poly.Union(poly2) In [4]: union_poly.ExportToWkt() Out[4]: 'POLYGON ((0 0 0,0 1 0,0 4 0,1 1 0,0 0 0)) – mgc Aug 26 '15 at 17:40
  • Yup, that's it. Thanks. Would you mind posting as an answer, so I can accept it? (For reference: with a third triangle only connected on one vertice, this will create a MULTIPOLYGON (((0 0 0,0 1 0,0 2 0,1 1 0,0 0 0)),((1 1 0,2 2 0,2 0 0,1 1 0)))) – LuWi Aug 27 '15 at 9:29
3

As asked by the OP, I put this in answer. The problem wasn't with the method used, as ogr.Geometry.Union() is the appropriate one but with the second polygon which was not closed.
It can be done with :
tri2.CloseRings()
or by repeating the first point :
tri2.AddPoint(1,1)
If there is many features to dissolve (like the entire layer) it can be done with UnionCascaded() which is supposed to provide significant time savings compared to an iterative union.

  • I can confirm the significantly lower runtime of UnionCascaded(). It also helps if you don't dissolve the polygons one-by-one, but rather in a binary tree fashion. – LuWi Sep 1 '15 at 12:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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