[prev in list] [next in list] [prev in thread] [next in thread]
List: postgis-users
Subject: Re: [postgis-users] ST_Intersection returning offset point on two lines which intersect
From: Martin Davis <mtnclimb () gmail ! com>
Date: 2020-04-26 4:18:59
Message-ID: CAK2ens25oai2NxCsL+voseTo-DAU8VwH7SprwmeKsi37x++30w () mail ! gmail ! com
[Download RAW message or body]
[Attachment #2 (multipart/alternative)]
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@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@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@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@lists.osgeo.org
>> > https://lists.osgeo.org/mailman/listinfo/postgis-users
>>
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users@lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/postgis-users
>
> _______________________________________________
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
[Attachment #5 (text/html)]
<div dir="ltr">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.</div><br><div class="gmail_quote"><div dir="ltr" \
class="gmail_attr">On Fri, Apr 24, 2020 at 7:55 AM David Campbell <<a \
href="mailto:david.campbell@forcorp.com">david.campbell@forcorp.com</a>> \
wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px \
0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Thanks \
for getting back to me.<div><br></div><div>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. </div><div><br></div><div>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:</div><div><br></div><div><div \
style="color:rgb(0,0,0);background-color:rgb(255,255,254);font-family:SFMono-Medium,"SF \
Mono","Segoe UI Mono","Roboto Mono","Ubuntu \
Mono",Menlo,monospace;font-size:13px;line-height:18px;white-space:pre-wrap"><div><span \
style="color:rgb(9,101,69)">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</span></div><div><span style="color:rgb(9,101,69)">FROM \
linestring_table a JOIN linestring_table b ON ST_DWithin(a.extend_geom, b.geom, \
50)</span></div><div><span style="color:rgb(9,101,69)">WHERE <a href="http://a.id" \
target="_blank">a.id</a> <> <a href="http://b.id" target="_blank">b.id</a> AND \
a.extend_geom IS NOT NULL</span></div></div></div><div><br></div><div>The returned \
geometry in many cases was offset from where the st_intersection actually was. \
</div><div><br></div><div>How do I mold what you (Paul) have solved into my case? \
</div><div><br></div><div><br></div><div><div><div dir="ltr"><div dir="ltr">David \
</div></div></div></div></div><br><div class="gmail_quote"><div dir="ltr" \
class="gmail_attr">On Thu, Apr 23, 2020 at 1:02 PM Paul Ramsey <<a \
href="mailto:pramsey@cleverelephant.ca" \
target="_blank">pramsey@cleverelephant.ca</a>> wrote:<br></div><blockquote \
class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid \
rgb(204,204,204);padding-left:1ex">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.<br> <br>
This works, for example.<br>
<br>
WITH lines AS (<br>
SELECT <br>
'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,<br> '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<br> ),<br>
outputs AS (<br>
SELECT <br>
ST_Snap(l1, ST_Intersection(l1, l2), 0.00001) AS l1,<br>
ST_Snap(l2, ST_Intersection(l1, l2), 0.00001) AS l2,<br>
ST_Intersection(l1, l2) AS p<br>
FROM lines<br>
)<br>
SELECT ST_Intersects(p, l1) AS p_l1, ST_Intersects(p, l2) as p_l2<br>
FROM outputs;<br>
<br>
<br>
<br>
<br>
> On Apr 23, 2020, at 11:42 AM, David Campbell <<a \
href="mailto:david.campbell@forcorp.com" \
target="_blank">david.campbell@forcorp.com</a>> wrote:<br> > <br>
> Hi,<br>
> <br>
> I have two lines which intersect but ST_Intersection returns an offset point. \
Offset is around 0.000026m.<br> > <br>
> 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),<br> > \
<br> > line2 <br>
> 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)<br> > SELECT \
st_intersects(ST_Intersection(line1.the_geom, line2.the_geom),line1.the_geom)<br> \
> <br> > <br>
> FROM<br>
> line1<br>
> <br>
> JOIN line2 ON true<br>
> <br>
> <a href="https://gis.stackexchange.com/questions/359210/line-intersection-resulting-in-offset-point" \
rel="noreferrer" target="_blank">https://gis.stackexchange.com/questions/359210/line-intersection-resulting-in-offset-point</a> \
<br> > <br>
> I have just tried to introduce ST_SnapToGrid, but all attempts at syntax still \
return false. in example:<br> > <br>
> 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),<br> > <br>
> line2 <br>
> 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)<br> > \
SELECT st_intersects(ST_SnapToGrid(ST_Intersection(line1.the_geom, line2.the_geom), \
0.001), line1.the_geom)<br> > <br>
> <br>
> FROM<br>
> line1<br>
> <br>
> JOIN line2 ON true<br>
> <br>
> how do I get this query to return true?<br>
> <br>
> <br>
> David<br>
> _______________________________________________<br>
> postgis-users mailing list<br>
> <a href="mailto:postgis-users@lists.osgeo.org" \
target="_blank">postgis-users@lists.osgeo.org</a><br> > <a \
href="https://lists.osgeo.org/mailman/listinfo/postgis-users" rel="noreferrer" \
target="_blank">https://lists.osgeo.org/mailman/listinfo/postgis-users</a><br> <br>
_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@lists.osgeo.org" \
target="_blank">postgis-users@lists.osgeo.org</a><br> <a \
href="https://lists.osgeo.org/mailman/listinfo/postgis-users" rel="noreferrer" \
target="_blank">https://lists.osgeo.org/mailman/listinfo/postgis-users</a></blockquote></div>
_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@lists.osgeo.org" \
target="_blank">postgis-users@lists.osgeo.org</a><br> <a \
href="https://lists.osgeo.org/mailman/listinfo/postgis-users" rel="noreferrer" \
target="_blank">https://lists.osgeo.org/mailman/listinfo/postgis-users</a></blockquote></div>
[Attachment #6 (text/plain)]
_______________________________________________
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic