PostGIS Splitting – Splitting Lines into Non-Overlapping Subsets Based on Points Using PostGIS

linear-referencingpostgissplitting

Given a table with line geometry, and one or more points that are snapped to this line in a separate table, I would like to split each line with one or more intersecting points at each of the locations where the line intersects a point.

For example, there is a line, L, with three intersecting points, A, B, and C in order along the line geometry. I would like to return L as four distinct geometries: from the beginning of L to A, from A to B along L, from B to C along L, and from C to the end of L.

In the past I have used shapely for this task, which is a linear referencing problem (http://sgillies.net/blog/1040/shapely-recipes/). However, this would not be practicable in this case, which has many millions of lines and points. Instead, I'm looking for a solution using PostgreSQL/PostGIS.

Note that points are constrained to be on a line. Further, a point can validly be on the start or end of a line, in which case the line need not be split (unless there are other points which are not coincident with the same line's start or end points). The subset lines need to retain their direction and their attributes, but the attributes of the point features do not matter.

Best Answer

The ST_Split PostGIS function is probably what you want.

PostGIS 2.2+ now supports Multi* geometries in ST_Split.

For older versions of PostGIS, read on:

To get a single line split by multiple points, you could use something like this multipoint wrapper plpgsql function. I've simplified it just to the "split (multi)lines with (multi)points" case below:

DROP FUNCTION IF EXISTS split_line_multipoint(input_geom geometry, blade geometry);
CREATE FUNCTION split_line_multipoint(input_geom geometry, blade geometry)
  RETURNS geometry AS
$BODY$
    -- this function is a wrapper around the function ST_Split 
    -- to allow splitting multilines with multipoints
    --
    DECLARE
        result geometry;
        simple_blade geometry;
        blade_geometry_type text := GeometryType(blade);
        geom_geometry_type text := GeometryType(input_geom);
    BEGIN
        IF blade_geometry_type NOT ILIKE 'MULTI%' THEN
            RETURN ST_Split(input_geom, blade);
        ELSIF blade_geometry_type NOT ILIKE '%POINT' THEN
            RAISE NOTICE 'Need a Point/MultiPoint blade';
            RETURN NULL;
        END IF;

        IF geom_geometry_type NOT ILIKE '%LINESTRING' THEN
            RAISE NOTICE 'Need a LineString/MultiLineString input_geom';
            RETURN NULL;
        END IF;

        result := input_geom;           
        -- Loop on all the points in the blade
        FOR simple_blade IN SELECT (ST_Dump(ST_CollectionExtract(blade, 1))).geom
        LOOP
            -- keep splitting the previous result
            result := ST_CollectionExtract(ST_Split(result, simple_blade), 2);
        END LOOP;
        RETURN result;
    END;
$BODY$
LANGUAGE plpgsql IMMUTABLE;

-- testing
SELECT ST_AsText(split_line_multipoint(geom, blade))
    FROM (
        SELECT ST_GeomFromText('Multilinestring((-3 0, 3 0),(-1 0, 1 0))') AS geom,
        ST_GeomFromText('MULTIPOINT((-0.5 0),(0.5 0))') AS blade
        --ST_GeomFromText('POINT(-0.5 0)') AS blade
    ) AS T;

Then to create a multipoint geometry to cut by, use ST_Collect and either create it manually from inputs:

SELECT ST_AsText(ST_Collect(
  ST_GeomFromText('POINT(1 2)'),
  ST_GeomFromText('POINT(-2 3)')
));

st_astext
----------
MULTIPOINT(1 2,-2 3)

Or aggregate it from a subquery:

SELECT stusps,
  ST_Multi(ST_Collect(f.the_geom)) as singlegeom
FROM (SELECT stusps, (ST_Dump(the_geom)).geom As the_geom
      FROM somestatetable ) As f
GROUP BY stusps
Related Question