1

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:

query version 3 result set

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.

query version 4 result set

3
  • If you just want the length of the route in each state, you can just compute the length of the intersection of the route with each state polygon.
    – dr_jts
    Mar 19 at 4:41
  • @dr_jts how would I go about doing this?
    – AJ_
    Mar 21 at 13:12
  • See answer below.
    – dr_jts
    Mar 21 at 19:19

3 Answers 3

1

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
;
2
  • ST_Intersection doesn't preserve order of intersections along a line. (Perhaps it should, but it doesn't currently). This can probably be done by using ST_LineLocatePoint to compute the distance fraction and order by that.
    – dr_jts
    Mar 19 at 4:33
  • Actually my last comment is incorrect. In recent versions of PostGIS ST_Intersection does preserve order along a line. So this technique can be part of the solution. Still need to associate line fragments with the state they cross, however (probably by joining on the midpoint of the fragment).
    – dr_jts
    Mar 21 at 21:02
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));
3
  • What version of PostgreSQL and PostGIS are you using?
    – AJ_
    Mar 22 at 11:28
  • 1
    POSTGIS="3.1.4 v2.4.3-44-g7ef659a" [EXTENSION] PGSQL="140" GEOS="3.9.1-CAPI-1.14.2" PROJ="8.1.1" LIBXML="2.9.11" LIBJSON="0.15" LIBPROTOBUF="1.4.0" WAGYU="0.5.0 (Internal)"
    – dr_jts
    Mar 22 at 18:07
  • Thanks! I have added this to our internal issue for when we update everything to the latest version (since our current setup is older than yours). I'll then test out the above code and if it works as expected I'll mark this as the answer. Thanks again!
    – AJ_
    Mar 23 at 12:15
0

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
1
  • Thanks for this! This does answer part of my question, but how would I keep these points in proper order (like in my example above)?
    – AJ_
    Mar 21 at 20:00

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.