[prev in list] [next in list] [prev in thread] [next in thread]
List: postgis-users
Subject: Re: [postgis-users] Trouble with ST_ShortestLine: the returned line DOES NOT start in g1 and end in
From: Martin Davis <mtnclimb () gmail ! com>
Date: 2022-02-16 23:45:14
Message-ID: CAK2ens3uuqxo9Skc-6FtXFPHTXhQZE2Oivs=To1ut0BY9aRMbA () mail ! gmail ! com
[Download RAW message or body]
[Attachment #2 (multipart/alternative)]
Paul's comment about constructed points not necessarily lying on lines due
to numerical rounding is 100% correct. However, in this particular case
it turns out that the fact that ST_Intersects reports false is a bug in
GEOS. See:
https://github.com/libgeos/geos/issues/565
The good news is that a fix has been identified and will be implemented in
GEOS soon.
Handling robustness problems in computation geometry algorithms is a
long-standing research problem. So any suggested solution should be made
in that perspective.
My thinking is that a way to deal with this is to support a distance
tolerance on spatial predicates. This will require a new algorithm for
evaluating predicates, which hopefully can be worked on in the near/medium
term.
On Wed, Feb 9, 2022 at 4:51 PM Jorge Gustavo Rocha <jgr@di.uminho.pt> wrote:
> Hi Paul,
>
> Thanks for the feedback.
>
> I would expect, as a assertion:
> ST_Intersects( ST_ClosestPoint(line, point), line) = true
>
> In our current implementation, the ST_ClosestPoint(geometry g1, geometry
> g2) can return a point that is not on g1.
>
> If ST_ClosestPoint is returning a point that is not exactly on g1, due
> to the FP limitations, the ST_Intersects will have the same limitations
> and could agree that point is on the line.
>
> St_distance is reporting zero as the distance between the
> ST_ClosestPoint and the line, as expected. It seems like a exact zero.
>
> My point is, as a user, even knowing that the fp types are inexact, why
> st_distance is returning zero and st_intesects, st_touches, st_crosses,
> st_disjoint, are returning unexpected values?
>
> Maybe this can not be solved using the current fp representation, but
> maybe we can improve the logic to overcome some fp limitations.
>
> I'll push this issue to my Postgis stack of "todo".
>
> Regards,
>
> Jorge
>
> test code
>
> --8<-----------------
>
> --test code
>
> with
> p as (
> select st_geometryfromtext('SRID=3763;POINT (-29802.23305474065
> 138518.7938969792)') as geom
> ),
> l as (
> select st_geometryfromtext('SRID=3763;LINESTRING
> (-29796.696826656284
> 138522.76848210802, -29804.3911369969 138519.3504205817)') as geom
> ),
> short as (
> select ST_ShortestLine(l.geom, p.geom) as geom,
> ST_ClosestPoint(l.geom,
> p.geom) as closestpoint
> from l, p
> )
> select short.geom, short.closestpoint, ST_Intersects(short.closestpoint,
> l.geom), ST_StartPoint(short.geom),
> st_distance(short.closestpoint, l.geom) * 1E+24,
> ST_Disjoint(short.geom, l.geom),
> ST_Intersects(short.geom, l.geom),
> ST_Touches(short.geom, l.geom),
> ST_Crosses(short.geom, l.geom)
> from short, l;
>
>
> On 09/02/22 19:58, Paul Ramsey wrote:
> > If the shortest line runs vertext-to-vertext I would expect it to touch
> both lines. If on the other hand it runs vertex-to-midsegment I would
> expect it is possible it would not touch. Check the distance of the
> shortest line to the two parent lines. It should be zero or very very very
> small. That would be expected. It's not possible to place constructed
> points exactly on the lines they are constructed from, except in some cases
> and on vertices, thanks to the granularity of floating point
> representation. (Discrete math, not continuous math)
> >
> > P
> >
> > > On Feb 9, 2022, at 11:31 AM, Jorge Gustavo Rocha <jgr@di.uminho.pt>
> wrote:
> > >
> > > Hi,
> > >
> > > Help needed :-) I'm having trouble with lines returned by
> ST_ShortestLine function.
> > >
> > > The computed lines *should* start and end on the geometries [1], but it
> is not happening.
> > >
> > > In fact, the computed geometry sometimes crosses the original
> geometries or does not touch the geometries at all.
> > >
> > > EXAMPLE 1: (of a computed geometry between a point and a line, that
> does not touches the line)
> > >
> > > with
> > > p as (
> > > select st_geometryfromtext('SRID=3763;POINT (-29802.23305474065
> 138518.7938969792)') as geom
> > > ),
> > > l as (
> > > select st_geometryfromtext('SRID=3763;LINESTRING
> (-29796.696826656284 138522.76848210802, -29804.3911369969
> 138519.3504205817)') as geom
> > > ),
> > > short as (
> > > select ST_ShortestLine(l.geom, p.geom) as geom
> > > from l, p
> > > )
> > > select short.geom,
> > > ST_Disjoint(short.geom, l.geom),
> > > ST_Intersects(short.geom, l.geom),
> > > ST_Touches(short.geom, l.geom),
> > > ST_Crosses(short.geom, l.geom)
> > > from short, l;
> > >
> > > > geom
> > st_disjoint|st_intersects|st_touches|st_crosses|
> > >
> > -----------------------------------------------------------------------------------------|-----------|-------------|----------|----------|
> >
> > > > LINESTRING (-29802.795222153436 138520.05937757515, -29802.23305474065
> 138518.7938969792)|true |false |false |false |
> > >
> > > EXAMPLE 2: (of a computed geometry between a point and a line. The
> computed line crosses the original line.)
> > >
> > > l as (
> > > select st_geometryfromtext('SRID=3763;LINESTRING
> (-29790.144327955728 138508.0183209069, -29798.28270998299
> 138504.40298837385)') as geom
> > > ),
> > > short as (
> > > select ST_ShortestLine(l.geom, p.geom) as geom
> > > from l, p
> > > )
> > > select short.geom,
> > > ST_Disjoint(short.geom, l.geom),
> > > ST_Intersects(short.geom, l.geom),
> > > ST_Touches(short.geom, l.geom),
> > > ST_Crosses(short.geom, l.geom)
> > > from short, l;
> > >
> > > > geom
> > st_disjoint|st_intersects|st_touches|st_crosses|
> > >
> > -----------------------------------------------------------------------------------------|-----------|-------------|----------|----------|
> >
> > > > LINESTRING (-29795.33635541747 138505.7118543719, -29794.774188004685
> 138504.44637377595)|false |true |false |true |
> > >
> > > In both cases, I was expecting that the computed line touch the
> original line.
> > >
> > > Am I understanding the ST_ShortestLine incorrectly or should I report a
> possible error? Can this be an (math approximation) error?
> > >
> > > I'm using Postgis 3.3 dev, GEOS 3.11
> > >
> > > POSTGIS="3.3.0dev 3.1.0rc1-1150-g05d3150c9" [EXTENSION] PGSQL="140"
> GEOS="3.11.0dev-CAPI-1.16.0" PROJ="8.2.0" LIBXML="2.9.12"
> LIBPROTOBUF="1.3.3" WAGYU="0.5.0 (Internal)"
> > >
> > >
> > > Thanks in advance,
> > >
> > > Jorge Gustavo
> > >
> > > [1] https://postgis.net/docs/ST_ShortestLine.html
> > > --
> > > Jorge Gustavo Rocha
> > > Departamento de Informática
> > > Universidade do Minho
> > > 4710-057 Braga
> > > Gabinete 3.29 (Piso 3)
> > > Tel: +351 253604480
> > > Fax: +351 253604471
> > > Móvel: +351 910333888
> > > skype: nabocudnosor
> > > _______________________________________________
> > > 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
>
> J. Gustavo
> --
> Jorge Gustavo Rocha
> Departamento de Informática
> Universidade do Minho
> 4710-057 Braga
> Gabinete 3.29 (Piso 3)
> Tel: +351 253604480
> Fax: +351 253604471
> Móvel: +351 910333888
> skype: nabocudnosor
> _______________________________________________
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
>
[Attachment #5 (text/html)]
<div dir="ltr">Paul's comment about constructed points not necessarily lying on \
lines due to numerical rounding is 100% correct. However, in this particular case \
it turns out that the fact that ST_Intersects reports false is a bug in GEOS. \
See:<div><br></div><div><a \
href="https://github.com/libgeos/geos/issues/565">https://github.com/libgeos/geos/issues/565</a><br></div><div><br></div><div>The \
good news is that a fix has been identified and will be implemented in GEOS soon. \
</div><div><br></div><div>Handling robustness problems in computation geometry \
algorithms is a long-standing research problem. So any suggested solution should be \
made in that perspective.</div><div><br></div><div>My thinking is that a way to deal \
with this is to support a distance tolerance on spatial predicates. This will \
require a new algorithm for evaluating predicates, which hopefully can be worked on \
in the near/medium term. </div><div> </div></div><br><div \
class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Feb 9, 2022 at 4:51 PM \
Jorge Gustavo Rocha <<a href="mailto:jgr@di.uminho.pt">jgr@di.uminho.pt</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">Hi Paul,<br> <br>
Thanks for the feedback.<br>
<br>
I would expect, as a assertion:<br>
ST_Intersects( ST_ClosestPoint(line, point), line) = true<br>
<br>
In our current implementation, the ST_ClosestPoint(geometry g1, geometry <br>
g2) can return a point that is not on g1.<br>
<br>
If ST_ClosestPoint is returning a point that is not exactly on g1, due <br>
to the FP limitations, the ST_Intersects will have the same limitations <br>
and could agree that point is on the line.<br>
<br>
St_distance is reporting zero as the distance between the <br>
ST_ClosestPoint and the line, as expected. It seems like a exact zero.<br>
<br>
My point is, as a user, even knowing that the fp types are inexact, why <br>
st_distance is returning zero and st_intesects, st_touches, st_crosses, <br>
st_disjoint, are returning unexpected values?<br>
<br>
Maybe this can not be solved using the current fp representation, but <br>
maybe we can improve the logic to overcome some fp limitations.<br>
<br>
I'll push this issue to my Postgis stack of "todo".<br>
<br>
Regards,<br>
<br>
Jorge<br>
<br>
test code<br>
<br>
--8<-----------------<br>
<br>
--test code<br>
<br>
with<br>
p as (<br>
select st_geometryfromtext('SRID=3763;POINT (-29802.23305474065 <br>
138518.7938969792)') as geom<br>
),<br>
l as (<br>
select st_geometryfromtext('SRID=3763;LINESTRING (-29796.696826656284 \
<br> 138522.76848210802, -29804.3911369969 138519.3504205817)') as geom<br>
),<br>
short as (<br>
select ST_ShortestLine(l.geom, p.geom) as geom, ST_ClosestPoint(l.geom, \
<br> p.geom) as closestpoint<br>
from l, p<br>
)<br>
select short.geom, short.closestpoint, ST_Intersects(short.closestpoint, <br>
l.geom), ST_StartPoint(short.geom),<br>
st_distance(short.closestpoint, l.geom) * 1E+24,<br>
ST_Disjoint(short.geom, l.geom),<br>
ST_Intersects(short.geom, l.geom),<br>
ST_Touches(short.geom, l.geom),<br>
ST_Crosses(short.geom, l.geom)<br>
from short, l;<br>
<br>
<br>
On 09/02/22 19:58, Paul Ramsey wrote:<br>
> If the shortest line runs vertext-to-vertext I would expect it to touch both \
lines. If on the other hand it runs vertex-to-midsegment I would expect it is \
possible it would not touch. Check the distance of the shortest line to the two \
parent lines. It should be zero or very very very small. That would be expected. \
It's not possible to place constructed points exactly on the lines they are \
constructed from, except in some cases and on vertices, thanks to the granularity of \
floating point representation. (Discrete math, not continuous math)<br> > <br>
> P<br>
> <br>
>> On Feb 9, 2022, at 11:31 AM, Jorge Gustavo Rocha <<a \
href="mailto:jgr@di.uminho.pt" target="_blank">jgr@di.uminho.pt</a>> wrote:<br> \
>><br> >> Hi,<br>
>><br>
>> Help needed :-) I'm having trouble with lines returned by \
ST_ShortestLine function.<br> >><br>
>> The computed lines *should* start and end on the geometries [1], but it is \
not happening.<br> >><br>
>> In fact, the computed geometry sometimes crosses the original geometries or \
does not touch the geometries at all.<br> >><br>
>> EXAMPLE 1: (of a computed geometry between a point and a line, that does not \
touches the line)<br> >><br>
>> with<br>
>> p as (<br>
>> select st_geometryfromtext('SRID=3763;POINT (-29802.23305474065 \
138518.7938969792)') as geom<br> >> ),<br>
>> l as (<br>
>> select st_geometryfromtext('SRID=3763;LINESTRING \
(-29796.696826656284 138522.76848210802, -29804.3911369969 138519.3504205817)') \
as geom<br> >> ),<br>
>> short as (<br>
>> select ST_ShortestLine(l.geom, p.geom) as geom<br>
>> from l, p<br>
>> )<br>
>> select short.geom,<br>
>> ST_Disjoint(short.geom, l.geom),<br>
>> ST_Intersects(short.geom, l.geom),<br>
>> ST_Touches(short.geom, l.geom),<br>
>> ST_Crosses(short.geom, l.geom)<br>
>> from short, l;<br>
>><br>
>> |geom \
|st_disjoint|st_intersects|st_touches|st_crosses|<br> >> \
|-----------------------------------------------------------------------------------------|-----------|-------------|----------|----------|<br>
>> |LINESTRING (-29802.795222153436 138520.05937757515, -29802.23305474065 \
138518.7938969792)|true |false |false |false |<br> \
>><br> >> EXAMPLE 2: (of a computed geometry between a point and a line. \
The computed line crosses the original line.)<br> >><br>
>> l as (<br>
>> select st_geometryfromtext('SRID=3763;LINESTRING \
(-29790.144327955728 138508.0183209069, -29798.28270998299 138504.40298837385)') \
as geom<br> >> ),<br>
>> short as (<br>
>> select ST_ShortestLine(l.geom, p.geom) as geom<br>
>> from l, p<br>
>> )<br>
>> select short.geom,<br>
>> ST_Disjoint(short.geom, l.geom),<br>
>> ST_Intersects(short.geom, l.geom),<br>
>> ST_Touches(short.geom, l.geom),<br>
>> ST_Crosses(short.geom, l.geom)<br>
>> from short, l;<br>
>><br>
>> |geom \
|st_disjoint|st_intersects|st_touches|st_crosses|<br> >> \
|-----------------------------------------------------------------------------------------|-----------|-------------|----------|----------|<br>
>> |LINESTRING (-29795.33635541747 138505.7118543719, -29794.774188004685 \
138504.44637377595)|false |true |false |true \
|<br> >><br>
>> In both cases, I was expecting that the computed line touch the original \
line.<br> >><br>
>> Am I understanding the ST_ShortestLine incorrectly or should I report a \
possible error? Can this be an (math approximation) error?<br> >><br>
>> I'm using Postgis 3.3 dev, GEOS 3.11<br>
>><br>
>> POSTGIS="3.3.0dev 3.1.0rc1-1150-g05d3150c9" [EXTENSION] \
PGSQL="140" GEOS="3.11.0dev-CAPI-1.16.0" PROJ="8.2.0" \
LIBXML="2.9.12" LIBPROTOBUF="1.3.3" WAGYU="0.5.0 \
(Internal)"<br> >><br>
>><br>
>> Thanks in advance,<br>
>><br>
>> Jorge Gustavo<br>
>><br>
>> [1] <a href="https://postgis.net/docs/ST_ShortestLine.html" rel="noreferrer" \
target="_blank">https://postgis.net/docs/ST_ShortestLine.html</a><br> >> -- \
<br> >> Jorge Gustavo Rocha<br>
>> Departamento de Informática<br>
>> Universidade do Minho<br>
>> 4710-057 Braga<br>
>> Gabinete 3.29 (Piso 3)<br>
>> Tel: +351 253604480<br>
>> Fax: +351 253604471<br>
>> Móvel: +351 910333888<br>
>> skype: nabocudnosor<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><br> <br>
J. Gustavo<br>
-- <br>
Jorge Gustavo Rocha<br>
Departamento de Informática<br>
Universidade do Minho<br>
4710-057 Braga<br>
Gabinete 3.29 (Piso 3)<br>
Tel: +351 253604480<br>
Fax: +351 253604471<br>
Móvel: +351 910333888<br>
skype: nabocudnosor<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> \
</blockquote></div>
_______________________________________________
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