[postgis-users] ST_Intersection returning offset point on two lines which intersect

Fri Apr 24 07:54:27 PDT 2020

Thanks for getting back to me.

I see you have the problem solved as true is returned for both, but I have
no idea how to implement this theory into my current problem.

I boiled down what was happening to the smallest scale I could and posted
the query that replicates it. In actuality, I have a table of records with
thousands of linestring geometries. I have extended certain lines using
ST_Extend and populated a new field called extend_geom in the same table.
Now I am trying to create a new line in my table connecting the original
geom to the locations where the extend geom intersects another geom using:

SELECT DISTINCT (CASE WHEN ST_Intersects(a.extend_geom, b.geom) THEN
ST_Shortestline(a.geom, ST_Intersection(a.extend_geom, b.geom)) ELSE NULL
END) geom
FROM linestring_table a JOIN linestring_table b ON
ST_DWithin(a.extend_geom, b.geom, 50)
WHERE a.id <> b.id AND a.extend_geom IS NOT NULL

The returned geometry in many cases was offset from where the
st_intersection actually was.

How do I mold what you (Paul) have solved into my case?

David

On Thu, Apr 23, 2020 at 1:02 PM Paul Ramsey <pramsey at cleverelephant.ca>
wrote:

> You don't. If you introduce a vertex into the original line at the
> intersection point, then the vextexes can match up, but the new point is
> not necessarily representable as being on the mathematical intersection of
> the two inputs, it's at the closest representable point.
>
> This works, for example.
>
> WITH lines AS (
> SELECT
> 'LINESTRING(470874.945140537 6000126.53834916,470825.026548551
> 6000129.39039651,470813.970641131 6000131.67770864,470798.339625488
> 6000130.91512638,470770.509325729 6000128.62757298,470716.373201037
> 6000122.90909733,470592.358279565 6000112.18733429,470584.082725689
> 6000111.4719529,470471.616772575 6000098.31957132,470442.261380633
> 6000108.99375658,470402.993405271 6000132.6307069,470339.326451592
> 6000140.06487686,470261.553382752 6000143.87732586,470226.86058284
> 6000141.9710548,470158.999734862 6000139.4932543,470075.507824541
> 6000123.86262835,470042.340010929 6000112.04367615,470005.359698272
> 6000095.84140767,469987.441344206 6000088.9789069,469981.323382904
> 6000088.34605401,469976.38532627 6000087.83546738,469928.191247079
> 6000074.51890886)'::geometry AS l1,
> 'LINESTRING(470835.701330449 6000121.00308247,470835.612988186
> 6000125.16357655,470834.757333107 6000165.42524755,470834.681958091
> 6000168.96974537,470841.821713877 6000270.84345442,470838.892222296
> 6000668.7623786,470838.348652514 6000835.0334808)'::geometry AS l2
> ),
> outputs AS (
>         SELECT
>                 ST_Snap(l1, ST_Intersection(l1, l2), 0.00001) AS l1,
>                 ST_Snap(l2, ST_Intersection(l1, l2), 0.00001) AS l2,
>                 ST_Intersection(l1, l2) AS p
>         FROM lines
> )
> SELECT ST_Intersects(p, l1) AS p_l1, ST_Intersects(p, l2) as p_l2
> FROM outputs;
>
>
>
>
> > On Apr 23, 2020, at 11:42 AM, David Campbell <david.campbell at forcorp.com>
> wrote:
> >
> > Hi,
> >
> > I have two lines which intersect but ST_Intersection returns an offset
> point. Offset is around 0.000026m.
> >
> > WITH line1 AS (SELECT ST_GeomFromText('LINESTRING(470874.945140537
> 6000126.53834916,470825.026548551 6000129.39039651,470813.970641131
> 6000131.67770864,470798.339625488 6000130.91512638,470770.509325729
> 6000128.62757298,470716.373201037 6000122.90909733,470592.358279565
> 6000112.18733429,470584.082725689 6000111.4719529,470471.616772575
> 6000098.31957132,470442.261380633 6000108.99375658,470402.993405271
> 6000132.6307069,470339.326451592 6000140.06487686,470261.553382752
> 6000143.87732586,470226.86058284 6000141.9710548,470158.999734862
> 6000139.4932543,470075.507824541 6000123.86262835,470042.340010929
> 6000112.04367615,470005.359698272 6000095.84140767,469987.441344206
> 6000088.9789069,469981.323382904 6000088.34605401,469976.38532627
> 6000087.83546738,469928.191247079 6000074.51890886)') as the_geom),
> >
> >      line2
> > AS (SELECT ST_GeomFromText('LINESTRING(470835.701330449
> 6000121.00308247,470835.612988186 6000125.16357655,470834.757333107
> 6000165.42524755,470834.681958091 6000168.96974537,470841.821713877
> 6000270.84345442,470838.892222296 6000668.7623786,470838.348652514
> 6000835.0334808)') as the_geom)
> > SELECT st_intersects(ST_Intersection(line1.the_geom,
> line2.the_geom),line1.the_geom)
> >
> >
> > FROM
> >  line1
> >
> > JOIN line2 ON true
> >
> >
> https://gis.stackexchange.com/questions/359210/line-intersection-resulting-in-offset-point
>
> >
> > I have just tried to introduce ST_SnapToGrid, but all attempts at syntax
> still return false. in example:
> >
> > WITH line1 AS (SELECT ST_SnapToGrid(the_geom, 0.001) the_geom FROM
> (SELECT ST_GeomFromText('LINESTRING(470874.945140537
> 6000126.53834916,470825.026548551 6000129.39039651,470813.970641131
> 6000131.67770864,470798.339625488 6000130.91512638,470770.509325729
> 6000128.62757298,470716.373201037 6000122.90909733,470592.358279565
> 6000112.18733429,470584.082725689 6000111.4719529,470471.616772575
> 6000098.31957132,470442.261380633 6000108.99375658,470402.993405271
> 6000132.6307069,470339.326451592 6000140.06487686,470261.553382752
> 6000143.87732586,470226.86058284 6000141.9710548,470158.999734862
> 6000139.4932543,470075.507824541 6000123.86262835,470042.340010929
> 6000112.04367615,470005.359698272 6000095.84140767,469987.441344206
> 6000088.9789069,469981.323382904 6000088.34605401,469976.38532627
> 6000087.83546738,469928.191247079 6000074.51890886)') as the_geom) a),
> >
> >      line2
> > AS (SELECT ST_SnapToGrid(the_geom, 0.001) the_geom FROM (SELECT
> ST_GeomFromText('LINESTRING(470835.701330449
> 6000121.00308247,470835.612988186 6000125.16357655,470834.757333107
> 6000165.42524755,470834.681958091 6000168.96974537,470841.821713877
> 6000270.84345442,470838.892222296 6000668.7623786,470838.348652514
> 6000835.0334808)') as the_geom) a)
> > SELECT st_intersects(ST_SnapToGrid(ST_Intersection(line1.the_geom,
> line2.the_geom), 0.001), line1.the_geom)
> >
> >
> > FROM
> >  line1
> >
> > JOIN line2 ON true
> >
> > how do I get this query to return true?
> >
> >
> > David
> > _______________________________________________
> > postgis-users mailing list
> > postgis-users at lists.osgeo.org
> > https://lists.osgeo.org/mailman/listinfo/postgis-users
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20200424/a99edb25/attachment.html>