4

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)

2
  • Can you add the WKT of line 1, line 2 and points ? At least, the CRS of those ? Dec 7, 2021 at 21:19
  • 1
    WKT added to the original post.
    – Cliff
    Dec 15, 2021 at 15:27

2 Answers 2

4

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
;
3
  • Thanks! This almost works. For some reasons some of the cross section lines are in the correct location, but some go to a vertex instead of the closest point along the line.
    – Cliff
    Dec 3, 2021 at 14:40
  • You can update your question with this specific case. Dec 3, 2021 at 15:14
  • 1
    Done. See the update above.
    – Cliff
    Dec 7, 2021 at 16:22
2

So, this design is tested on the published data from your question and it works as shown in your first picture. Since you have additional conditions in your example, you have to look for workarounds, and although they are a bit complicated, they still exist.

Maybe some of the code will help you, maybe not :-)...

Run the SQL code:

WITH tbla AS (SELECT ST_Makeline(pts.geom, ST_ClosestPoint(line2.geom,pts.geom))::geometry('linestring',26917) as geom FROM line2, pts),
    tblb AS (SELECT ST_OffsetCurve(geom, -0.15, 'quad_segs=4 join=bevel') geom FROM line1),
    tblc AS (SELECT a.geom, ST_Distance(ST_StartPoint(a.geom), ST_EndPoint(a.geom)) as distance_m FROM tbla a JOIN tblb b ON ST_Disjoint(a.geom, b.geom)),
    tbld AS (SELECT a.geom geom FROM tbla a JOIN tblb b ON ST_Intersects(a.geom, b.geom)),
    tble AS (SELECT ST_MakeLine(pnt1, pnt2) geom FROM (SELECT ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) pnt1,
             ST_PointN(geom, generate_series(2, ST_NPoints(geom))) pnt2 FROM line2) AS geom),
    tblf AS (SELECT a.geom geom FROM tble a JOIN tbld b ON ST_Intersects(a.geom, b.geom)),
    tblg AS (SELECT ST_Difference(a.geom, b.geom) geom FROM tble a JOIN tblf b ON true),
    tblh AS (SELECT ST_Length(a.geom) lenght, ST_MakeLine(ST_StartPoint(a.geom), st_closestpoint(b.geom,a.geom)) AS geom FROM tbld a JOIN tblg b ON true),
    tbli AS (SELECT geom, ST_Distance(ST_StartPoint(geom), ST_EndPoint(geom)) as distance_m, MIN(lenght) min_length FROM tblh GROUP BY geom ORDER BY min_length DESC LIMIT 1)
    SELECT geom, distance_m FROM tblc UNION SELECT geom, distance_m FROM tbli

Check the result, it should match my figure below.

enter image description here

Figure

The original spatial solutions...

10
  • Thanks! This is 99.9% there. The original problem point (#5 from the right) does not have a line. Otherwise, this fixes all the issues on the western side. Any idea why there is no line there?
    – Cliff
    Jan 26 at 16:37
  • I wrote the SQL-code so that it works as shown in your and my drawings...i.e. I have one line for each point... Jan 26 at 17:15
  • If one line doesn't appear when using this query, any ideas why? If you're curious, reach out to me at cliffpatterson@gmail.com and I'll share some creds with you.
    – Cliff
    Jan 26 at 17:38
  • I did. Same exact data and I still end up without one line. It's strange.
    – Cliff
    Jan 26 at 17:41

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.