[postgis-devel] Fixed precision for ST_Intersection, ST_Union, etc.

[Daniel Baston]

The relevant code is in GeometryPrecisionReducer.cpp, around line 72.  It
can be paraphrased as:

geom_reduced = reduce_coordinate_precision(geom)
if (!geom_reduced.isValid())
   geom_reduced = geom_reduced.buffer(0)
return geom_reduced

The reason I chose to bypass everything after the first line is that (a) an
isValid() check is always performed, and (b) if that fails, our geometry
gets the bazooka of .buffer(0).  There's no algorithm available to
surgically remove only degeneracies caused by the precision reduction.
It's very easy to imagine a case (and I've seen them) where:

- user has a large polygon at full double-precision
- coordinates are snapped or truncated to 1e-6, producing a tiny bow-tie
effect along an edge somewhere
- in an attempt to fix the problem, user performs a .buffer(0), causing
99.9% of the original polygon to disappear, and keeping only the tiny part
of the bowtie

This is not the kind of thing I'd want to have silently and automatically
in every overlay function.  I really don't see the problem with having
overlay functions fail on invalid inputs, even if the input is invalid only
because the user supplied an incorrect precision argument.

The only change I can see making to GEOS is in the C API, renaming the
existing GEOSGeom_setPrecision to GEOSGeom_reducePrecision, and adding a
new GEOSGeom_setPrecision that does nothing more than set the
PrecisionModel.  The PostGIS overlays could then use the new
GEOSGeom_setPrecision, while a PostGIS ST_ReducePrecision could use
GEOSGeom_reducePrecision.

Dan


On Tue, Mar 15, 2016 at 7:22 AM, Sandro Santilli <strk at keybit.net> wrote:

> On Tue, Mar 15, 2016 at 05:52:10AM -0400, Daniel Baston wrote:
>
> > The cost of the checking is the same as ST_IsValid,
>
> Are you sure ? Do you have numbers to show ?
> Had you considered improving GEOS to make it irrelevant if
> none of the coordinates were really changed ?
> (as that would be your use case).
>
> > Here's a query that I think illustrates the problem:
> >
> > -- No preicsion model
> > postgres=# SELECT ST_AsText(ST_Intersection('POINT (0 0)', 'POLYGON ((0
> 0,
> > 1 0, 0 1, 1 1, 0 0))'));
> > ERROR:  Error performing intersection: TopologyException: Input geom 1 is
> > invalid: Self-intersection at or near point 0.5 0.5 at 0.5 0.5
> > CONTEXT:  SQL function "st_intersection" statement 1
> >
> > -- With precision model, GEOSGeom_setPrecision flags GEOS_PREC_NO_TOPO |
> > GEOS_PREC_KEEP_COLLAPSED
> > postgres=# SELECT ST_AsText(ST_Intersection('POINT (0 0)', 'POLYGON ((0
> 0,
> > 1 0, 0 1, 1 1, 0 0))', op_precision := 1e-12));
> > ERROR:  Error performing intersection: TopologyException: Input geom 1 is
> > invalid: Self-intersection at or near point 0.5 0.5 at 0.5 0.5
> >
> > -- WIth precision model, no flags to GEOSGeom_setPrecision
> > postgres=# SELECT ST_AsText(ST_Intersection('POINT (0 0)', 'POLYGON ((0
> 0,
> > 1 0, 0 1, 1 1, 0 0))', op_precision := 1e-12));
> >         st_astext
> > --------------------------
> >  GEOMETRYCOLLECTION EMPTY
>
> The final result is unexpected, how was the polygon changed ?
> Reduction should have not moved any vertex, right ?
>
> I guess I just don't like how the handling is performed by
> GEOS. Ideally, it should only fix what was broken by precision
> reduction, not what was broken in the input.
>
> Remember the GEOS side of this is recent, so might need to
> still be tweaked.
>
> --strk;
> _______________________________________________
> postgis-devel mailing list
> postgis-devel at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/postgis-devel
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-devel/attachments/20160315/13e5ee67/attachment.html>