Dividir la línea en varios polígonos preservando el orden

Question

Tengo polígonos almacenados en la base de datos como (Multi)Polígonos usando PostgreSQL (13.1) y PostGIS (3.0.2) con un srid de 4326 (también conocido como states). Una línea generada por el usuario (LineString) (llamada ruta) se usará para encontrar todas las intersecciones de esos polígonos a lo largo de esa línea.

Puedo obtener todos los puntos de intersección utilizando la siguiente consulta (utilizando esta ruta / línea):

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

El problema es que no puedo encontrar una forma de mantener los puntos de intersección en orden a lo largo de la ruta. He intentado ordenar por distancia desde el punto de inicio, pero el problema es que necesito ambos puntos por cada state, lo que posiblemente cause que el orden sea incorrecto.

Por ejemplo, tengo una ruta que atraviesa los estados A, B, C. Necesito el punto de inicio al cruce de la frontera de A. Luego necesito el cruce de la frontera de B a C. Luego de C al punto final.

La razón por la que necesito esto es para calcular la distancia, la cual actualmente manejo en una consulta separada ya que todavía no estoy seguro de cómo combinarlas en una sola consulta. Soy totalmente nuevo en PostGIS y mi rol principal no es como administrador de BD.

He intentado:

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)

y

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

pero aún así no obtengo el orden correcto.

Orden esperado: punto de inicio, frontera de A, frontera de B, frontera de C, punto final Orden de salida: punto de inicio, frontera de B, frontera de A, frontera de C, punto final

Soy muy principiante en esto y he estado buscando en la web y utilizando la documentación de PostGIS para llegar hasta donde estoy ahora.

Info del entorno: 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 (próximamente actualizaré a 5.2.4.4), rgeo gem 2.1.1

---

Solo para que lo tengas, esta es la consulta que estoy usando para calcular la distancia:

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 distancia
FROM (
  SELECT :route AS route
    , :start_point AS start_point
    , :end_point AS end_point
) AS subconsulta_distancia

---

Mi objetivo final es tener una consulta que divida la ruta en cada estado y obtenga la distancia de cada segmento de ruta en ese estado. Puede haber varios segmentos en cada estado. Lo que necesitaré son los puntos de la frontera para cada estado (donde la ruta se cruza) y la longitud de cada segmento.

Estoy abierto a otras ideas o sugerencias si la forma en que estoy abordando esto es incorrecta o imposible. Me gustaría que la consulta se ejecute en menos de 1 segundo si es posible.

---

Consulta actualizada:

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

Resultados:

El orden debería ser: 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

Ahora esto no está devolviendo los puntos de la frontera y sigue fuera de orden.

---

Actualización 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)"

Aquí se encuentra la línea de ruta de prueba

Answer 1

ST_Dump devuelve un geometry_dump tipo compuesto, teniendo

• path INT[]
• geom GEOMETRY

campos; deberías poder conservar el orden de cruce de Puntos de borde a través del campo path incluso a través de MultiPolígonos:

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
;

---

Para añadir tanto el ST_StartPoint como el ST_EndPoint de la :route, simplemente UNION en un conjunto de resultados de consulta:

SELECT ST_Startpoint(:route) AS geom
UNION ALL

UNION ALL
SELECT ST_EndPoint(:route) AS geom
;

Answer 2

Resulta que en las versiones recientes de PostGIS la función ST_Intersection ha sido mejorada para:

• preservar el orden de secciones a lo largo de una línea
• permitir la intersección con un MultiPolygon inválido con elementos adyacentes

Entonces, la siguiente consulta creará las secciones de ruta (en orden); las vinculará de vuelta a los estados padres (a través de la intersección del punto medio de la sección; y extraerá los puntos de inicio y fin.

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 3

Si desea la longitud de la ruta en cada estado, la consulta es simple. Calcule la longitud de la intersección de la línea de la ruta con cada polígono estatal que cruza:

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);

Resultado:

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