5

I want to split a line shapefile into segments of equal length, let say into length of x meters each. There is tool v.split.length as mentioned in this answer: how to split a line into a specific no. of parts? Is there a way to do that in python without using QGIS? May be 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 specified length.

1

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
  • thanks for the solution, how to save the line_segments into shapefile? – FJ_Abbasi Aug 17 '17 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). – lpd2 Aug 17 '17 at 12:01
  • how to convert the list of ogr wkbLineString Geometries (line_segments) to list of wtkLineString? – FJ_Abbasi Aug 17 '17 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 – lpd2 Aug 17 '17 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 Aug 18 '17 at 8:38
7

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.

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.