2

I have two polygons called 1, 2. I have a linestring that crosses the two polygons but does not remain only inside them.

enter image description here

What I'm trying to do:

Find the linestring that is from A to D in the case of polygon 1 and Find the line that is from E to H in the case of polygon number 2

So from that big linstring I only need these two linstrings and their size in KM

Now as I tried to do:

def geospatial_intersection(pg_conn, geom_A, geom_B):
   
    sql1 = """SELECT st_intersection('%s' , '%s') as rez""" % (geom_A, geom_B)
    #print(sql1)
    c1 = pg_conn.cursor()
    c1.execute(sql1)
    row1 = c1.fetchall() 
    
    return row1[0][0]

It's just that it only returns these to me: MULTILINESTRING with segments A-B and C-D in the case of polygon 1 and MULTILINESTRING with segments E-F and G-H in the case of polygon 2

I miss polygon 1 segment B-C and polygon 2 segment F-G ...

    SELECT
  ply.ogc_fid AS pid, 
  lns.ogc_fid AS lid,
  ST_LineSubstring(lns.wkb_geometry, frc.st, frc.ed) AS geom
FROM
  limite_pug_bacau AS ply
  JOIN
  traseu AS lns ON ST_Intersects(ply.wkb_geometry, lns.wkb_geometry),
  LATERAL (
    SELECT
      ST_Intersection(ply.wkb_geometry, lns.wkb_geometry) AS geom
  ) AS its,
  LATERAL (
    SELECT
      ST_MakeLine(dmp.geom ORDER BY dmp.path) AS geom
    FROM
      ST_DumpPoints(its.geom) AS dmp
  ) AS ln,
  LATERAL (
    SELECT
      ST_LineLocatePoint(lns.wkb_geometry, ST_StartPoint(ln.geom)) AS st,
      ST_LineLocatePoint(lns.wkb_geometry, ST_EndPoint(ln.geom)) AS ed
  ) AS frc
;

I tried as geozelot said but it returns this error to me:

ERROR:  line_locate_point: 1st arg isn't a line
   SQL state: XX000

The two tables in which the geometries are located are not linked with each other

4
  • In SQL terms, the intersection of a line and a polygon is the parts of the line within the polygon, so your result is correct. If your want the parts of the line which are outside the polygon, you need a different query (using ST_Difference to erase the parts of the line within the polygons).
    – Vince
    Jun 15, 2022 at 19:40
  • If you split the linestring at intersections the desired result contains the parts that are inside any polygon, plus the parts with both start point and end point touching the boundary of the same polygon. Would that be correct?
    – user30184
    Jun 15, 2022 at 19:57
  • Yes. This is corect! @user30184 Jun 15, 2022 at 20:14
  • Then for finding the missing pieces you could use postgis.net/docs/ST_ExteriorRing.html, postgis.net/docs/ST_StartPoint.html, postgis.net/docs/ST_EndPoint.html, and maybe postgis.net/docs/ST_Touches.html but postgis.net/docs/ST_Intersects.html will do as well. There will be a practical problem because it is difficult to make artificially created points to touch exactly a line when used in another function so it may be necessary to use ST_Buffer for start/endpoints, or to use postgis.net/docs/ST_DWithin.html instead ST_Touches.
    – user30184
    Jun 15, 2022 at 20:48

1 Answer 1

1

In theory, simple, but a little wordy in SQL - still fast:

Edited to use OPs relation/column names.

SELECT
  ply.ogc_fid AS pid, 
  lns.ogc_fid AS lid,
  ST_LineSubstring(lns.wkb_geometry, frc.st, frc.ed) AS geom
FROM
  limite_pug_bacau AS ply
  JOIN
  traseu AS lns ON ST_Intersects(ply.wkb_geometry, lns.wkb_geometry),
  LATERAL (
    SELECT
      ST_Intersection(ply.wkb_geometry, lns.wkb_geometry) AS geom
  ) AS its,
  LATERAL (
    SELECT
      ST_MakeLine(dmp.geom ORDER BY dmp.path) AS geom
    FROM
      ST_DumpPoints(its.geom) AS dmp
  ) AS ln,
  LATERAL (
    SELECT
      ST_LineLocatePoint(ST_LineMerge(lns.wkb_geometry), ST_StartPoint(ln.geom)) AS st,
      ST_LineLocatePoint(ST_LineMerge(lns.wkb_geometry), ST_EndPoint(ln.geom)) AS ed
  ) AS frc
;
8
  • I edited the question and asked the problem I encountered with your version @geozelot. Can we transcribe a query because we don't need to link the tables together? Jun 16, 2022 at 6:14
  • 1
    @CoderCoder42 the FROM clause specifies tables (or set rerurning expressions), nor columns I assume wkb_geometry is the column name in both tables? See my update.
    – geozelot
    Jun 16, 2022 at 8:13
  • Yes its the geometry column name its the both tables Jun 16, 2022 at 11:38
  • 1
    Mixed value MULTILINESTRING and LINESTRING uplied this ST_LineMerge and resolv.... Jun 16, 2022 at 12:26
  • 1
    Thank you! @geozelot Jun 16, 2022 at 12:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.