[GIS] How to create point labels from midpoints of road polylines in QGIS

centroidsdissolveqgis

I have a road polyline layer (singlepart ie has several features per road, grouped by an id field), and i want to (after merging the parts) find the mid point of the resulting lines. This is for roads geolocation.

The problem is that the QGIS' dissolve function only works on polygons.

I've tried:
a) mmqgis plugin "create label layer" (handles multiparts nicely)
b) buffer road lines to 1m then dissolve by roadid, then ftools centroids

Both of which give a 'sort of' result, but the main problem is that the point is not actually on the road line. Is there a way to either force polygon centroids to be within the buffer polygon, or alternatively, some way to actually join the polylines and find mid points directly.

Note that while there are (a small number of) 'Y' shaped lines, most of the multipart lines start and end at the same place, (maybe gridifyed to 1cm).

Best Answer

You can join the polylines by editing the layer, selecting all the lines (all of them or by attribute) and then click the 'Merge Features' tool:

To find the midpoint, calculate distance of lines, then use point along line tools to put a point at distance/2 along each line.

Related Solutions

[GIS] How to create midpoints for a set of polylines in QGIS

There is a GDAL tool Points along lines (in Processing Toolbox | GDAL | Vector geoprocessing).

Start this tool, then set the Distance from line start represented as fraction of line length parameter to 0.5.

This 0.5, being the fraction, puts a point at the middle of the line. For you most important part of this tool is that the points are placed on the line(s).

QGIS – How to Snap a Road Network to a Hexagonal Grid in QGIS

My solution involves a PyQGIS script that is faster and more effective than a workflow involving snapping (I gave it a try too). Using my algorithm I've obtained these results:

You can run the following code snippets in sequence from within QGIS (in the QGIS Python console). At the end you get a memory layer with the snapped routes loaded into QGIS.

The only prerequisite is to create a multipart road Shapefile (use Processing->Singleparts to multipart, I used the field fictitiuos as Unique ID field parameter). This will give us a roads_multipart.shp file with a single feature.

Here is the algorithm explained:

Get the nearest hexagon sides where routes cross. For each hexagon we create 6 triangles between each pair of neighbour vertices and the corresponding centroid. If any road intersects a triangle, the segment shared by the hexagon and the triangle is added to the final snapped route. This is the heavier part of the whole algorithm, it takes 35 seconds running on my machine. In the first two lines there are 2 Shapefile paths, you should adjust them to fit your own file paths.

hexgrid = QgsVectorLayer("/docs/borrar/hex_grid_question/layers/normal-hexgrid.shp", "hexgrid", "ogr")
roads = QgsVectorLayer("/docs/borrar/hex_grid_question/layers/roads_multipart.shp", "roads", "ogr")  # Must be multipart!

roadFeat = roads.getFeatures().next() # We just have 1 geometry
road = roadFeat.geometry() 
indicesHexSides = ((0,1), (1,2), (2,3), (3,4), (4,5), (5,0))

epsilon = 0.01
# Function to compare whether 2 segments are equal (even if inverted)
def isSegmentAlreadySaved(v1, v2):
    for segment in listSegments:        
        p1 = QgsPoint(segment[0][0], segment[0][1])
        p2 = QgsPoint(segment[1][0], segment[1][1])
        if v1.compare(p1, epsilon) and v2.compare(p2, epsilon) \
            or v1.compare(p2, epsilon) and v2.compare(p1, epsilon):
            return True
    return False

# Let's find the nearest sides of hexagons where routes cross
listSegments = []
for hexFeat in hexgrid.getFeatures():
    hex = hexFeat.geometry()
    if hex.intersects( road ):
        for side in indicesHexSides:
            triangle = QgsGeometry.fromPolyline([hex.centroid().asPoint(), hex.vertexAt(side[0]), hex.vertexAt(side[1])])
            if triangle.intersects( road ):
                # Only append new lines, we don't want duplicates!!!
                if not isSegmentAlreadySaved(hex.vertexAt(side[0]), hex.vertexAt(side[1])): 
                    listSegments.append( [[hex.vertexAt(side[0]).x(), hex.vertexAt(side[0]).y()], [hex.vertexAt(side[1]).x(),hex.vertexAt(side[1]).y()]] )

Get rid of disconnected (or 'open') segments by using Python lists, tuples, and dictionaries. At this point, there are some disconnected segments left, i.e., segments that have one vertex disconnected but the other one connected to at least other 2 segments (see red segments in the next figure). We need to get rid of them.

# Let's remove disconnected/open segments
lstVertices = [tuple(point) for segment in listSegments for point in segment]
dictConnectionsPerVertex = dict((tuple(x),lstVertices.count(x)-1) for x in set(lstVertices))

# A vertex is not connected and the other one is connected to 2 segments
def segmentIsOpen(segment):
    return dictConnectionsPerVertex[tuple(segment[0])] == 0 and dictConnectionsPerVertex[tuple(segment[1])] >= 2 \
        or dictConnectionsPerVertex[tuple(segment[1])] == 0 and dictConnectionsPerVertex[tuple(segment[0])] >= 2

# Remove open segments
segmentsToDelete = [segment for segment in listSegments if segmentIsOpen(segment)]        
for toBeDeleted in segmentsToDelete:
    listSegments.remove( toBeDeleted )

Now we can create a vector layer from the list of coordinates and load it to the QGIS map:

# Create a memory layer and load it to QGIS map canvas
vl = QgsVectorLayer("LineString", "Snapped Routes", "memory")
pr = vl.dataProvider()
features = []
for segment in listSegments:
    fet = QgsFeature()
    fet.setGeometry( QgsGeometry.fromPolyline( [QgsPoint(segment[0][0], segment[0][1]), QgsPoint(segment[1][0], segment[1][1])] ) )
    features.append(fet)

pr.addFeatures( features )
vl.updateExtents()
QgsMapLayerRegistry.instance().addMapLayer(vl)

Another part of the result:

Should you need attributes in the snapped routes, we could use a Spatial Index to rapidly evaluate intersections (such as in https://gis.stackexchange.com/a/130440/4972 ), but that's another story.

Hope this helps!

Best Answer

Related Solutions

[GIS] How to create midpoints for a set of polylines in QGIS

QGIS – How to Snap a Road Network to a Hexagonal Grid in QGIS

Related Question