Split line across multiple polygons preserving order

Question

I have polygons stored in the database as (Multi)Polygons using PostgreSQL (13.1) and PostGIS (3.0.2) with a srid of 4326 (aka states). A user generated line (LineString) (aka route) will be used to find all intersections of those polygons along that line.

I am able to get all of the intersection points using the following query (using this route / linestring):

SELECT states.*
  , (ST_Dump(ST_Intersection(ST_Boundary(states.boundary_polygon::geometry), :route::geometry))).geom AS point
FROM states
WHERE ST_Intersects(:route::geometry, states.boundary_polygon::geometry)

The issue is, I can't seem to find a way to keep the intersection points in order along the route. I have tried to use order by distance from the start point, but the issue is I need both points per each state, thus causing the order to possibly be wrong.

For example, I have a route that splits states A, B, C. I need the start point to border crossing of A. I then need the border crossing of B to C. Then C to the end point.

The reason I need this is to calculate distance, which I currently handle in a separate query as I'm not sure how to combine them into one query yet. I'm totally new to PostGIS and my main role isn't a DB admin.

I have tried:

ORDER BY ST_LineLocatePoint(ST_GeomFromText(:route_str, 4326)::geometry, (ST_Dump(ST_Intersection(ST_Boundary(states.boundary_polygon::geometry),ST_GeomFromText(:route_str, 4326)::geometry))).geom)

and

ORDER BY :start_point <-> point
  , :start_point <-> states_with_breakpoints.boundary_polygon::geometry
  , :end_point <-> point
  , :end_point <-> states_with_breakpoints.boundary_polygon::geometry

but still don't get the correct order.

Expected order: start point, A border, B border, C border, end point Output order: start point, B border, A border, C border, end point

I'm very much a beginner in this and have been searching the webs and using PostGIS docs to get as far as I have now.

Env Info: PostgreSQL 13.1, POSTGIS="3.0.2 2fb2a18" [EXTENSION] PGSQL="130" GEOS="3.8.1-CAPI-1.13.3" PROJ="7.1.1" LIBXML="2.9.10" LIBJSON="0.15", [RoR] Ruby 2.5.3, Rails 5.2.4.2 (going to be upgrading to 5.2.4.4 very soon), rgeo gem 2.1.1

---

Just so you have it, this is the query I'm using for distance:

SELECT ST_Length(
  ST_LineSubstring(
    route,
    LEAST(ST_LineLocatePoint(route, start_point), ST_LineLocatePoint(route, end_point)),
    GREATEST(ST_LineLocatePoint(route, end_point), ST_LineLocatePoint(route, start_point))
  )::geography
) AS distance
FROM (
  SELECT :route AS route
    , :start_point AS start_point
    , :end_point AS end_point
) AS distance_sub_query

---

My end goal is to have a query that splits the route in each state and gets the distance of each route segment in that state. There could be multiple segments in each state. What i'll need is the border points for each state (where the route crosses) and the length of each segment.

I'm open to other ideas or suggestions if the way I'm going about this is wrong or impossible. Would like the query to be under 1 second if possible.

---

Updated query:

SELECT s.id, s.abv, dmp.geom
FROM   states AS s,
       LATERAL ST_Dump(ST_Intersection(ST_Boundary(s.geom), :route::GEOMETRY)) AS dmp
WHERE  ST_Intersects(s.geom, :route::GEOMETRY)
ORDER BY dmp.path

Results:

Order should be: OH -> WV, WV -> PA, PA -> WV, WV -> MD, MD -> PA, PA -> NJ, NJ -> NY

---

SELECT s.id, s.abv, dmp.geom
FROM   states AS s,
       LATERAL ST_Dump(ST_Intersection(s.boundary_polygon::geometry, :route::geometry)) AS dmp
WHERE  ST_Intersects(s.boundary_polygon::geometry, :route::geometry)
ORDER BY dmp.path

This is now not returning the border points and is still out of order.

---

Update 2023-11-16

ENV: POSTGIS="3.1.4 ded6c34" [EXTENSION] PGSQL="140" GEOS="3.9.0-CAPI-1.16.2" PROJ="7.2.1" LIBXML="2.9.10" LIBJSON="0.15" LIBPROTOBUF="1.3.3" WAGYU="0.5.0 (Internal)"

Testing route line here

Answer 1

It turns out that in recent versions of PostGIS the ST_Intersection function has been enhanced to:

• preserve order of sections along a line
• allow intersecting with a invalid MultiPolygon with adjacent elements

So the following query will: create the route sections (in order); link them back to the parent states (via intersecting the section midpoint; and extract the start and end points.

WITH states(id, geom) AS (VALUES
  (1, 'POLYGON ((-120 45, -115 45, -115 40, -120 40, -120 45))'::geometry),
  (2, 'POLYGON ((-110 45, -110 40, -115 40, -115 45, -110 45))'::geometry),
  (3, 'POLYGON ((-110 40, -110 45, -105 42, -103 45, -100 40, -110 40))'::geometry),
  (4, 'POLYGON ((-110 45, -103 45, -105 42, -110 45))'::geometry)
),
route(geom) AS (VALUES
  ('LINESTRING (-114 41, -112 43, -109 42, -106 44, -102 42)'::geometry)
),
route_states AS (SELECT ST_Collect(s.geom) AS geom
  FROM route r
  JOIN states s ON ST_Intersects(r.geom, s.geom)
),
route_segs AS (SELECT (dmp).path[1] AS seq, (dmp).geom AS seg
  FROM (
    SELECT ST_Dump(ST_Intersection(r.geom, rs.geom)) AS dmp
    FROM route r
    JOIN route_states rs ON true) AS t
)
SELECT rs.seq, s.id,
      ST_StartPoint(rs.seg) AS start_pt, ST_EndPoint(rs.seg) AS end_pt
  FROM route_segs rs
  JOIN states s ON ST_Intersects(s.geom, ST_LineInterpolatePoint(rs.seg, 0.5));

Answer 2

ST_Dump returns a geometry_dump composite type, having

• path INT[]
• geom GEOMETRY

fields; you should be able to retain order of border crossing Points via the path field even across MultiPolygons:

SELECT dmp.geom
FROM   states AS s,
       LATERAL ST_Dump(ST_Intersection(ST_Boundary(s.geom), :route::GEOMETRY)) AS dmp
WHERE  ST_Intersects(s.geom, :route::GEOMETRY)
ORDER BY
       dmp.path
;

---

To add both the ST_StartPoint and ST_EndPoint of the :route, simply UNION them in a set of query results:

SELECT ST_Startpoint(:route) AS geom
UNION ALL
<above_query>
UNION ALL
SELECT ST_EndPoint(:route) AS geom
;

Answer 3

If you want the length of the route in each state, the query is simple. Compute the length of the intersection of the route line with each state polygon that it crosses:

WITH states(id, geom) AS (VALUES
  (1, 'POLYGON ((-120 45, -115 45, -115 40, -120 40, -120 45))'::geometry),
  (2, 'POLYGON ((-110 45, -110 40, -115 40, -115 45, -110 45))'::geometry), 
  (3, 'POLYGON ((-110 40, -110 45, -105 42, -103 45, -100 40, -110 40))'::geometry), 
  (4, 'POLYGON ((-110 45, -103 45, -105 42, -110 45))'::geometry)
),
route(geom) AS (VALUES
  ('LINESTRING (-114 41, -112 43, -109 42, -106 44, -102 42)'::geometry)
)
SELECT s.id, ST_Length( ST_Intersection(r.geom, s.geom))
  FROM route r 
  JOIN states s ON ST_Intersects(r.geom, s.geom);

Result:

 id |     st_length      
----+--------------------
  2 |  4.936612231525109
  3 | 5.8468593073171355
  4 |  3.284920476535893