[GIS] adding points to a line shapefile at specified distance in python

gdalpythonshapefile

I want to add create nodes\points at an interval on a line shapefile let say at every 100 meters and save those points as shapefile.

I came across a plugin for Qgis called Qchainage, this exactly performs the task i want to do, but I am not sure how to adapt this plugin into my python script.

Best Answer

In PyQGIS you can use interpolate method from QgsGeometry to do that. Generated points can be stored as a memory layer. Complete example code is:

layer = iface.activeLayer()

feat = layer.getFeatures().next()

geom = feat.geometry()

length = geom.length()

distance = 100

points = []

iter = distance

while iter <= length:

    pt = feat.geometry().interpolate(iter).exportToWkt()

    points.append(pt)

    iter += distance

epsg = layer.crs().authid()

uri = "Point?crs=" + epsg + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'points',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(points)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromWkt(points[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

I tried it out with line shapefile of next image; where it is also observed memory point layer produced.

Editing Note:

A shorter code could be generated if you have installed in your system fiona and shapely python modules.

import fiona
from shapely.geometry import shape, mapping

line = fiona.open('pyqgis_data/new_line.shp')

crs = line.crs

line = line.next()

geom = shape(line['geometry'])

# length of the LineString
length = geom.length

# creation of the resulting shapefile
schema = {'geometry': 'Point','properties': {'id': 'int'}}

with fiona.open('pyqgis_data/new_shape.shp', 'w', 'ESRI Shapefile', schema, crs=crs) as output:
    # create points every 100 meters along the line
    for i, distance in enumerate(range(0, int(length), 100)):
         point = geom.interpolate(distance)   
         output.write({'geometry':mapping(point),'properties': {'id':i}})

It produces same result layer (but it's saved in disk as shapefile; not as memory layer).

Related Solutions

python – Getting Lines and Nodes Around the Confluence Point in a Network System (Line Shapefile)

You work here with a Graph (Graph Theory) and your confluence point is simply a node with a degree (number of edges adjacent to that node) > 2.

As you mentioned a Python tag, a solution in Python with Fiona (for reading the shapefile) and NetworkX (Graph Theory)

The Polylines

Open the shapefile and convert all the segments of the Polylines in egdes of a Graph

import fiona
import networkx as nx
import itertools # for extracting all the line segments
# creation of an undirected Graph
G = nx.Graph()
for line in fiona.open("lines.shp"):
   # convert the line segments to Graph edges
   for seg_start, seg_end in itertools.izip(line['geometry']['coordinates'],line['geometry']['coordinates'][1:]):
       G.add_edge(seg_start, seg_end)

The Graph with Nodes(x) and edges

Now, in NetworkX, the degrees of a node are obtained by degree (I use Shapely here but you can save the resulting shapefiles with Fiona and even Networkx)

from shapely.geometry import Point, LineString
for node in G.nodes_iter():
    if G.degree(node) > 2:
        print Point(node) #the nodes
        for edge in G.edges(node):
            print LineString(edge) # the edges (>2)

Result

You can translate easily the Fiona solution in PyQGIS.

[GIS] Splitting line shapefile into segments of equal length using Python

Shapely can do what you ask. Here's an example of splitting a line string into 4 equal length parts.

from shapely.geometry import LineString, MultiPoint
from shapely.ops import split

line = LineString([(0, 0), (10, 10)])
splitter = MultiPoint([line.interpolate((i/4), normalized=True) for i in range(1, 4)])
split(line, splitter).wkt
# 'GEOMETRYCOLLECTION (LINESTRING (0 0, 2.5 2.5), LINESTRING (2.5 2.5, 5 5), LINESTRING (5 5, 7.5 7.5), LINESTRING (7.5 7.5, 10 10))'

It uses Shapely's interpolate() to find the points at which to split. Shapely's split() function is new in version 1.6.0. This method is general and doesn't depend on the number of points that define your line or whether it bends.

In practice, split() is sensitive to the points being precisely on the geometry you want to split. Floating point issues might confound you in some cases.

Best Answer

Related Solutions

python – Getting Lines and Nodes Around the Confluence Point in a Network System (Line Shapefile)

[GIS] Splitting line shapefile into segments of equal length using Python

Related Question