11

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.

0

4 Answers 4

14

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.

1
  • I´ve tried this. It's awesome. :) BUT: How can i access the data of one member of the geometrycollection? For example, i would like to have the first linestring as an array. Commented Jul 14, 2021 at 8:27
5

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
7
  • thanks for the solution, how to save the line_segments into shapefile?
    – FJ_Abbasi
    Commented Aug 17, 2017 at 11:50
  • 1
    Have a look at the answer given here: gis.stackexchange.com/questions/129745/…, or here: pcjericks.github.io/py-gdalogr-cookbook/… (although change the driver to ESRI Shapefile instead of geojson). Commented Aug 17, 2017 at 12:01
  • how to convert the list of ogr wkbLineString Geometries (line_segments) to list of wtkLineString?
    – FJ_Abbasi
    Commented Aug 17, 2017 at 20:06
  • An ogr Geometry has a ExportToWkt() method which returns a wkt string, see here: gdal.org/python/osgeo.ogr.Geometry-class.html Commented Aug 17, 2017 at 21:36
  • the segments produced by this method have 0 with the x ,y values in the Linestring like this: 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?
    – FJ_Abbasi
    Commented Aug 18, 2017 at 8:38
1

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 lines + [LineString(line)]

    coords = list(line.coords)
    for i, p in enumerate(coords):
        pd = line.project(Point(p))

        if pd == distance:
            lines.append(LineString(coords[: i + 1]))
            line = LineString(coords[i:])

            return cut(line, distance, lines)

        if pd > distance:
            cp = line.interpolate(distance)
            lines.append(LineString(coords[:i] + [(cp.x, cp.y)]))
            line = LineString([(cp.x, cp.y)] + coords[i:])

            return cut(line, distance, lines)

Some test 1

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]

Some test 2

line = LineString([(0, 0), (1, 0), (2, 0), (3, 0), (4, 0), (5, 0)])
print(round(line.length, 4))
#5.0
lines = cut(line, 2, list())
print([round(line.length, 4) for line in lines])
#[2.0, 2.0, 1.0]
0

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]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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