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

Sat Apr 25 21:18:59 PDT 2020

You cannot rely on a point computed by ST_Intersection returning
ST_Intersects = true.  Instead you need to test using a small
distance tolerance, using ST_DWithin.

On Fri, Apr 24, 2020 at 7:55 AM David Campbell <david.campbell at forcorp.com>
wrote:

> 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
>
> _______________________________________________
> 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/20200425/eeffbf0f/attachment.html>