Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It only takes a minute to sign up.
Sign up to join this community
I want to split a line shapefile into segments of equal length, let say into lengths of x meters each. There is tool v.split.length as mentioned in this answer: How to split line into specific number of parts? Is there a way to do that in python without using QGIS?
Maybe using shapely and Fiona? as there is an interpolate method in shapely which could be used to add points at a specified distance on a line. But I couldn't find a way to divide a line into segments of a specified length.
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.
Here's a method using just ogr, it hasn't been tested extensively though. The split_line_multiple function takes an ogr LineString Geometry and returns a list of ogr LineString Geometries. Either the number of sub-strings or their individual length can be specified.
from osgeo import ogr
import math
def _distance(a, b):
""" Return the distance separating points a and b.
a and b should each be an (x, y) tuple.
Warning: This function uses the flat surface formulae, so the output may be
inaccurate for unprojected coordinates, especially over large distances.
"""
dx = abs(b[0] - a[0])
dy = abs(b[1] - a[1])
return (dx ** 2 + dy ** 2) ** 0.5
def _get_split_point(a, b, dist):
""" Returns the point that is <<dist>> length along the line a b.
a and b should each be an (x, y) tuple.
dist should be an integer or float, not longer than the line a b.
"""
dx = b[0] - a[0]
dy = b[1] - a[1]
m = dy / dx
c = a[1] - (m * a[0])
x = a[0] + (dist**2 / (1 + m**2))**0.5
y = m * x + c
# formula has two solutions, so check the value to be returned is
# on the line a b.
if not (a[0] <= x <= b[0]) and (a[1] <= y <= b[1]):
x = a[0] - (dist**2 / (1 + m**2))**0.5
y = m * x + c
return x, y
def split_line_single(line, length):
""" Returns two ogr line geometries, one which is the first length
<<length>> of <<line>>, and one one which is the remainder.
line should be a ogr LineString Geometry.
length should be an integer or float.
"""
line_points = line.GetPoints()
sub_line = ogr.Geometry(ogr.wkbLineString)
while length > 0:
d = _distance(line_points[0], line_points[1])
if d > length:
split_point = _get_split_point(line_points[0], line_points[1], length)
sub_line.AddPoint(line_points[0][0], line_points[0][1])
sub_line.AddPoint(*split_point)
line_points[0] = split_point
break
if d == length:
sub_line.AddPoint(*line_points[0])
sub_line.AddPoint(*line_points[1])
line_points.remove(line_points[0])
break
if d < length:
sub_line.AddPoint(*line_points[0])
line_points.remove(line_points[0])
length -= d
remainder = ogr.Geometry(ogr.wkbLineString)
for point in line_points:
remainder.AddPoint(*point)
return sub_line, remainder
def split_line_multiple(line, length=None, n_pieces=None):
""" Splits a ogr wkbLineString into multiple sub-strings, either of
a specified <<length>> or a specified <<n_pieces>>.
line should be an ogr LineString Geometry
Length should be a float or int.
n_pieces should be an int.
Either length or n_pieces should be specified.
Returns a list of ogr wkbLineString Geometries.
"""
if not n_pieces:
n_pieces = int(math.ceil(line.Length() / length))
if not length:
length = line.Length() / float(n_pieces)
line_segments = []
remainder = line
for i in range(n_pieces - 1):
segment, remainder = split_line_single(remainder, length)
line_segments.append(segment)
else:
line_segments.append(remainder)
return line_segments
LINESTRING (656772.081066615 -185004.101691471 0,656823.94704498 -185003.710250124 0,656861.329693538 -185004.101691471 0,656870.798486507 -184999.01659895 0) any idea why we have these zeros?
Aug 18, 2017 at 8:38
Here is a modified shapely cut function for splitting a line into segments.
from shapely.geometry import LineString, Point
def cut(line, distance, lines):
# Cuts a line in several segments at a distance from its starting point
if distance <= 0.0 or distance >= line.length:
return [LineString(line)]
coords = list(line.coords)
for i, p in enumerate(coords):
pd = line.project(Point(p))
if pd == distance:
return [
LineString(coords[:i+1]),
LineString(coords[i:])
]
if pd > distance:
cp = line.interpolate(distance)
lines.append(LineString(coords[:i] + [(cp.x, cp.y)]))
line = LineString([(cp.x, cp.y)] + coords[i:])
if line.length > distance:
cut(line, distance, lines)
else:
lines.append(LineString([(cp.x, cp.y)] + coords[i:]))
return lines
Some test
line = LineString([(6348053.180476057, 6652085.281694949), (6348206.3069751775, 6651979.8705129465)])
print(round(line.length, 4))
#185.9012
lines = cut(line, 30, list())
print([round(line.length, 4) for line in lines])
#[30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 5.9012]
You can use the shapely.ops.substring() function. This cuts out one substring at an arbitrary location given by distance along the line. I use np.arange to compute multiples of the distance. Then I apply Python's built-in map to get the substring of the line for each pair of distances (combined by slicing appropriately). Since the line is always the same, I need itertools.repeat to feed it to map.
from itertools import repeat
import numpy as np
from shapely import LineString
from shapely.ops import substring
def cut(line, distance):
distances = np.arange(0, line.length + distance, distance)
return list(map(substring, repeat(line), distances[:-1], distances[1:]))
Test slightly adapted from Grigory Nevsky's test:
line = LineString([(6348053.180476057, 6652085.281694949), (6348206.3069751775, 6651979.8705129465)])
print(round(line.length, 4))
#185.9012
lines = cut(line, 30)
print([round(line.length, 4) for line in lines])
#[30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 5.9012]
