PostGIS Distance – How to Determine Shortest Distance Between Points and Line 2 Without Intersecting Line 1

postgis

I am generating cross sections between each points and the closest point along Line 2 (red line). These cross sections are generated with the following statement in PostGIS:

        SELECT point.id,
st_makeline(point.geom,st_closestpoint(line.geom,point.geom))::geometry('linestring',26917) as geom,
    st_distance(point.geom,st_closestpoint(line.geom,point.geom)) as distance_m
    FROM line2, point

The problem is that while this works for >90% of the cross sections (i.e. most of the black lines below), there are a few cross sections that intersect with Line 1, which I need to avoid.

How do I find the shortest distance from each point to Line 2 (red line) that DOES NOT intersect with Line 1 (green line)?

Note that the points are all on Line 1, so they already intersect with that geometry once.

enter image description here

UPDATE

I ran the following query and the results are better, but still a few issues to resolve.

WITH segments AS (
    SELECT
        ROW_NUMBER() OVER(PARTITION BY id) AS id,
        ST_MakeLine(lag((pt).geom, 1, NULL) OVER (PARTITION BY id ORDER BY id, (pt).path), (pt).geom) AS geom
  FROM (SELECT id, ST_DumpPoints(geom) AS pt FROM line2 WHERE id = 1) as dumps
)

SELECT
    point.id,
    ST_ShortestLine(point.geom, closest_segment.geom) AS geom
FROM
    point
    LEFT JOIN LATERAL (
        SELECT
          segments.id,
          segments.geom,
          ST_Distance(point.geom, segments.geom) as dist
          FROM segments
            WHERE NOT ST_CROSSES(ST_ShortestLine(point.geom, segments.geom), (SELECT geom FROM line WHERE id = 1))
          ORDER BY segments.geom <-> point.geom
          LIMIT 1
    ) AS closest_segment ON TRUE
;

There are 3 issues:

  1. When segmenting line1, the result produces a NULL value for the geom.
  2. There are two points that do not have a cross-section.
  3. A few cross sections meet at the vertex instead of the closest location along line2.

enter image description here

WKT for the lines and points (all in EPSG:26917)

Points
POINT (517828.640365819 4739956.45175504)
POINT (517770.320250362 4739922.84355291)
POINT (517715.328384183 4739891.15332494)
POINT (517664.845481695 4739871.94301691)
POINT (517602.563855499 4739856.91274439)
POINT (517562.453236079 4739855.04713418)
POINT (517564.090909638 4739895.05315968)
POINT (517565.603003744 4739931.99145857)
POINT (517567.103531322 4739968.64720368)
POINT (517533.29866912 4739990.56644695)
POINT (517485.227046466 4739996.3494993)
POINT (517465.657916174 4739998.70368039)
POINT (517349.273067239 4739967.20679058)
POINT (517343.498718221 4740015.32636574)
POINT (517336.633726275 4740056.00105978)
POINT (517326.864327411 4740090.19395581)
POINT (517319.674977505 4740115.35668048)
POINT (517441.859575148 4739992.26320526)
POINT (517408.695055905 4739983.28798883)

Line 1

LINESTRING (517828.640365819 4739956.45175504,517715.328384183 4739891.15332494,517628.519882998 4739858.12000148,517562.453236079 4739855.04713418,517567.830753852 4739986.4122112,517465.657916174 4739998.70368039,517349.273067239 4739967.20679058,517340.054465343 4740044.02847304,517315.471526955 4740130.06875741)

Line 2

LINESTRING (517859.369038805 4739864.26573608,517732.99737115 4739774.38436759,517558.612151956 4739745.96034508,517484.479228377 4739779.76188537,517537.87029769 4739973.35252518,517479.485819017 4739971.0478747,517423.405990818 4739933.4052503,517366.173837381 4739819.70916025,517283.20642032 4739875.02077162,517255.550614632 4739990.25329532,517235.576977191 4740105.48581902)

Best Answer

It's a little bit more complex, haven't test it on all cases but it works :

-- explode the line id = 2 into segments
WITH segments AS (
    SELECT
        ROW_NUMBER() OVER(PARTITION BY gid) AS id,
        ST_MakeLine(lag((pt).geom, 1, NULL) OVER (PARTITION BY gid ORDER BY gid, (pt).path), (pt).geom) AS geom
  FROM (SELECT gid, ST_DumpPoints(geom) AS pt FROM line WHERE gid = 2) as dumps
)

SELECT
    point.id,
    ST_ShortestLine(point.geom, closest_segment.geom) AS geom
FROM
    point
    -- compute all distances between points and segments,
    -- keep the closest that doesn't cross the line id = 1
    LEFT JOIN LATERAL (
        SELECT
          segments.id,
          segments.geom,
          ST_Distance(point.geom, segments.geom) as dist
          FROM segments
            WHERE NOT ST_CROSSES(ST_ShortestLine(point.geom, segments.geom), (SELECT geom FROM line WHERE gid = 1))
          ORDER BY segments.geom <-> point.geom
          LIMIT 1
    ) AS closest_segment ON TRUE
;
Related Question