6

This seems like it should be easy, but I'm not figuring it out: How do I convert a polygon layer into a line layer representing the polygon edges, where each edge belongs only to a single feature (or part)? When I convert from polygons to lines, each edge is represented twice (because each edge belongs to two adjacent polygons). Even if I dissolve the features, the duplicates remain as separate parts, so I cannot just dissolve, then split into singleparts. This causes problems with symbolization, for example if I want to use a dashed line style.

My thought was to split each polygon boundary into separate edge segments, as shown in these two questions:
Split polygon into lines
How do I split polygons into line segments?

…and then to remove duplicates, for example using v.clean with rmdupl. But I can't figure out how to do either of those steps.

For splitting segments, either the answer isn't in those two linked threads, or else I'm not understanding the explanation. Using v.split doesn't work, because I can only select to split by length, or by number of vertices, not by intersections with other line segments.

For removing duplicates, rmdupl isn't working; the duplicate edges always remain after I run the algorithm. It seems that the edges aren't recognized as duplicate geometries, even though they have identical points. (Running snap first didn't help; the layer appears to be topologically correct already.)

What am I missing?

6

If you are not forced to use QGIS, another Open Source GIS software OpenJUMP http://openjump.org/ has a Planer Graph tool that may be exactly what you need.

Here you can find the tool.

enter image description here

If you need only the edges you can uncheck all extra options.

enter image description here

The result contains the common edges only once. With real data the result may not be perfect because adjacent polygons do not necessarily share exactly the same vertices. In that case you must fix the source data or edit the planar graph layer.

enter image description here

If you want to create planar graph with just QGIS this answer may be useful QGIS Polygon edge style based on adjoining feature?.

  • 1
    The beauties of open source. – Alexandre Neto Apr 15 '17 at 13:38
  • Yes, that is exactly what I need! – Nathan Apr 15 '17 at 15:00
2

Here's a python solution using the Fiona and Shapely libraries.. It might be possible to implement in pyqgis without the need for those libraries...

It only works with POLYGON layers at the moment, so MULTIPOLYGON layers would need to be converted.

It breaks polygons into distinct segments. It handles:-

  • snapping to fixed number of decimal places
  • combining line segments if they happen to be identical but going in different directions
  • counting numbers of times each equivalent geometry is seen

import fiona
from shapely.geometry import LineString, Point

def explode_polygons_to_line_segments(filename, dps=5):
    """
    Explode POLYGON layer to distinct edge segments, together
    with counts. (e.g. for a country map, 1=coastline and 2=land border)
    :param filename: Filename for source
    :param dps: Number decimal places to round to (suggest min 5)
    :return: 
    """
    segs = {}
    geoms = {}
    with fiona.open(filename, "r") as source:
        for feature in source:
            # only works on exterior for now
            coords = feature["geometry"]["coordinates"][0]
            coords_rounded = []
            for x, y in coords:
                rounded_x = round(x, dps)
                rounded_y = round(y, dps)
                coords_rounded.append((rounded_x, rounded_y))
            for i in range(0, len(coords_rounded)-1):
                x1, y1 = coords_rounded[i]
                x2, y2 = coords_rounded[i+1]
                # deduplicate lines which overlap but go in different directions
                if (x1 < x2):
                    key = (x1, y1, x2, y2)
                else:
                    if (x1 == x2):
                        if (y1 < y2):
                            key = (x1, y1, x2, y2)
                        else:
                            key = (x2, y2, x1, y1)
                    else:
                        key = (x2, y2, x1, y1)
                if key not in segs:
                    segs[key] = 1
                else:
                    segs[key] += 1
                line = LineString([Point(x1,y1), Point(x2,y2)])
                geoms[key] = line.wkt
    return geoms, segs

For a quick-and-dirty export to CSV you could call it with

geoms, segs = explode_polygons_to_line_segments("/path/to/singleparts.shp", dps=5)
with open("/tmp/exploded.csv","w") as fo:
    fo.write("COUNT\tGEOM\n")
    for key in segs:
        fo.write("{}\t{}\n".format(segs[key], geoms[key]))
print("Found {} unique segments".format(len(segs)))

And here's an example... I've styled by the overlap count, so the exterior borders can be styled differently to the adjacent borders.

enter image description here

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.