Reorder and Orient MULTILINESTRING Components in PostGIS

README.md

Reorder and Orient MULTILINESTRING Components in PostGIS

This PostgreSQL/PostGIS function, reorder_multilinestring, takes a MULTILINESTRING geometry and returns a new MULTILINESTRING with components reordered and oriented based on geometric criteria. It is particularly useful for GIS applications where the order and orientation of line segments matter, such as routing, network analysis, or pipeline management.

Features

1. Component extraction: Splits the input multilinestring into individual lines using ST_Dump.
2. Sequence ordering: Orders line components by minimum distance between consecutive endpoints, creating a coherent path.
3. Orientation of components: Ensures each line is oriented correctly to connect smoothly to the previous line.
4. Merging of close components: Optionally merges consecutive components into a single LINESTRING if the distance between them is below a configurable threshold.
5. Global orientation check: Optionally accepts two reference points. If the first point appears after the second along the assembled line, the entire multilinestring is reversed.
6. Full inversion: When reversing, both the order of components and the direction of each line are inverted to maintain logical traversal.
7. Flexible: Works for any number of line components and supports Z-dimension geometries.

Usage

SELECT reorder_multilinestring(
    geom,
    ST_SetSRID(ST_MakePoint(x1, y1), 25832),
    ST_SetSRID(ST_MakePoint(x2, y2), 25832),
    5  -- merge_threshold in map units
)
FROM route_table
WHERE road_code = 'SP42';

• geom: input MULTILINESTRING geometry.
• reference_point_start / reference_point_end: optional points used to check and enforce global orientation.
• merge_threshold: optional numeric distance (in the same SRID units as geom). If > 0, consecutive components closer than this threshold will be merged into a single linestring.

Benefits

• Produces a clean, ordered, and oriented multilinestring suitable for downstream GIS analysis.
• Handles Z-dimension geometries correctly.
• Can simplify geometry topology by merging adjacent segments within tolerance.
• Maintains internal orientation while allowing full inversion when required.
• Ready for integration in PostGIS pipelines, spatial ETL scripts, or GIS applications.

Changelog

• v1.0 – Initial implementation:

  • Component extraction, sequence ordering, orientation handling, and optional global inversion with reference points.

• v1.1 – New functionality:

  • Added merge_threshold parameter.
  • If consecutive components are closer than the given threshold, they are automatically fused into a single linestring.
  • Improved flexibility for handling fragmented networks, survey data, or digitized features with small gaps.

reorder_multilinestring.sql

CREATE OR REPLACE FUNCTION progressive.reorder_multilinestring(
    input_multilinestring geometry,
    reference_point_start geometry DEFAULT NULL,
    reference_point_end geometry DEFAULT NULL,
    merge_threshold double precision DEFAULT 0
)
RETURNS geometry
LANGUAGE plpgsql AS
$$
DECLARE
    line_components geometry[];
    ordered_sequence geometry[];
    merged_sequence geometry[];       -- sequence after optional merging
    num_components int;
    component_index int;
    comparison_index int;
    max_distance double precision := -1;
    first_line_index int;
    second_line_index int;
    first_line_starts boolean;
    second_line_starts boolean;
    current_line geometry;
    next_line geometry;
    next_line_index int;
    min_distance double precision;
    result_multilinestring geometry;
    concatenated_line geometry;
    loc_start double precision;
    loc_end double precision;
    component_used boolean[];
    start_points geometry[];
    end_points geometry[];
BEGIN
    -- Ensure the input geometry is a MULTILINESTRING
    IF GeometryType(input_multilinestring) NOT IN ('MULTILINESTRING', 'MULTILINESTRINGZ') THEN
        RAISE EXCEPTION 'Input must be a MULTILINESTRING';
    END IF;

    -- Extract individual line components
    SELECT array_agg((d).geom)
    INTO line_components
    FROM ST_Dump(input_multilinestring) AS d;

    num_components := array_length(line_components,1);
    IF num_components IS NULL OR num_components = 0 THEN
        RETURN NULL;
    END IF;

    start_points := ARRAY(SELECT ST_StartPoint(component) FROM unnest(line_components) component);
    end_points   := ARRAY(SELECT ST_EndPoint(component) FROM unnest(line_components) component);

    component_used := ARRAY(SELECT false FROM generate_series(1,num_components));
    ordered_sequence := ARRAY[]::geometry[];

    -- Step 1: find farthest pair
    FOR component_index IN 1..num_components LOOP
        FOR comparison_index IN 1..num_components LOOP
            IF component_index = comparison_index THEN CONTINUE; END IF;

            IF ST_Distance(start_points[component_index], start_points[comparison_index]) > max_distance THEN
                max_distance := ST_Distance(start_points[component_index], start_points[comparison_index]);
                first_line_index := component_index;
                second_line_index := comparison_index;
                first_line_starts := true;
                second_line_starts := true;
            END IF;
            IF ST_Distance(start_points[component_index], end_points[comparison_index]) > max_distance THEN
                max_distance := ST_Distance(start_points[component_index], end_points[comparison_index]);
                first_line_index := component_index;
                second_line_index := comparison_index;
                first_line_starts := true;
                second_line_starts := false;
            END IF;
            IF ST_Distance(end_points[component_index], start_points[comparison_index]) > max_distance THEN
                max_distance := ST_Distance(end_points[component_index], start_points[comparison_index]);
                first_line_index := component_index;
                second_line_index := comparison_index;
                first_line_starts := false;
                second_line_starts := true;
            END IF;
            IF ST_Distance(end_points[component_index], end_points[comparison_index]) > max_distance THEN
                max_distance := ST_Distance(end_points[component_index], end_points[comparison_index]);
                first_line_index := component_index;
                second_line_index := comparison_index;
                first_line_starts := false;
                second_line_starts := false;
            END IF;
        END LOOP;
    END LOOP;

    -- Step 2: initialize sequence
    current_line := line_components[first_line_index];
    IF first_line_starts = false THEN
        current_line := ST_Reverse(current_line);
    END IF;
    ordered_sequence := ARRAY[current_line];
    component_used[first_line_index] := true;

    -- Step 3: greedy sequence building
    WHILE array_length(ordered_sequence,1) < num_components LOOP
        min_distance := 1e20;
        next_line_index := NULL;

        FOR component_index IN 1..num_components LOOP
            IF component_used[component_index] THEN CONTINUE; END IF;

            IF ST_Distance(ST_EndPoint(current_line), ST_StartPoint(line_components[component_index])) < min_distance THEN
                min_distance := ST_Distance(ST_EndPoint(current_line), ST_StartPoint(line_components[component_index]));
                next_line_index := component_index;
                next_line := line_components[component_index];
            END IF;
            IF ST_Distance(ST_EndPoint(current_line), ST_EndPoint(line_components[component_index])) < min_distance THEN
                min_distance := ST_Distance(ST_EndPoint(current_line), ST_EndPoint(line_components[component_index]));
                next_line_index := component_index;
                next_line := ST_Reverse(line_components[component_index]);
            END IF;
        END LOOP;

        ordered_sequence := ordered_sequence || next_line;
        component_used[next_line_index] := true;
        current_line := next_line;
    END LOOP;

    -- Step 4: merge components if closer than threshold
    merged_sequence := ARRAY[]::geometry[];
    current_line := ordered_sequence[1];

    FOR component_index IN 2..array_length(ordered_sequence,1) LOOP
        next_line := ordered_sequence[component_index];
        IF merge_threshold > 0 AND ST_Distance(ST_EndPoint(current_line), ST_StartPoint(next_line)) <= merge_threshold THEN
            -- fuse with current
            current_line := ST_MakeLine(current_line, next_line);
        ELSE
            -- push current into merged_sequence
            merged_sequence := merged_sequence || current_line;
            current_line := next_line;
        END IF;
    END LOOP;

    -- push last one
    merged_sequence := merged_sequence || current_line;

    result_multilinestring := ST_Multi(ST_Collect(merged_sequence));

    -- Step 5: orientation check with reference points
    IF reference_point_start IS NOT NULL AND reference_point_end IS NOT NULL THEN
        concatenated_line := ST_MakeLine(
            ARRAY(
                SELECT (vertex).geom
                FROM unnest(merged_sequence) AS line_component,
                     LATERAL ST_DumpPoints(line_component) AS vertex
            )
        );

        loc_start := ST_LineLocatePoint(concatenated_line, reference_point_start);
        loc_end   := ST_LineLocatePoint(concatenated_line, reference_point_end);

        IF loc_start > loc_end THEN
            merged_sequence := ARRAY(
                SELECT ST_Reverse(merged_sequence[array_idx])
                FROM generate_series(array_length(merged_sequence,1),1,-1) AS array_idx
            );
            result_multilinestring := ST_Multi(ST_Collect(merged_sequence));
        END IF;
    END IF;

    RETURN result_multilinestring;
END;
$$;