[postgis-users] Trouble with ST_ShortestLine: the returned line DOES NOT start in g1 and end in g2

Wed Feb 16 15:45:14 PST 2022

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 at 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 at 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 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
>
> 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 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/20220216/6e5b383b/attachment.html>