[postgis-users] Wrong result in spatial analysis methodslikeintersects due to limited coordinate precision

Wed May 24 02:17:25 PDT 2006

Hi Gerhard,

What about applying snapToGrid(geom, size) to your data to reduce
their precision, adapting it to the real precision your data were
produced with.

Keeping so many digits to data expressed in meters does not make sense
for a geographic point of view. 8 digit is close to nanometric
precision.

If data are modified or inserted automatically, a trigger may be
applied to call snapToGrip()  on the result, I think.

Just my 2 cents ;-)

Nicolas

On 5/24/06, Sommer Gerhard <som at adv.magwien.gv.at> wrote:
> Hi Martin,
>
> Thank You for Your response.
>
> The problem I wrote about has nothing to do with bad data. I send You a sketch (attachment sketch.png), which, I hope, will show You the problem better:
>
> I have double precision coordinates, which means approximately 15 significant digits. Therefore, if the coordinates have a maximum of 7 digits before the decimal point, there can be a maximum of 8 digits after the decimal point. So the theoretical grid, to which all coordinate values, that can be stored in the PostGIS database, are snapped is 0.00000001 meters wide. This grid is shown on the sketch.
>
> Point A and point B (of polygon I) were snapped to the line segment connecting point 1 and point 2 (of polygon O) and polygon I was then stored in the database. Both points A and B are not exactly on the line connecting 1 and 2, as I we can not store exact coordinates in any database. But they are as near as possible to the line and for practical purposes the line A-B touches the line 1-2 and polygon I is within polygon O.
>
> So if I call the PostGIS operation contains(O, I) or within(I, O) (or was it the other way round?), I expect true as the result. The same goes for the test if point A or point B are positioned on the line 1-2. But PostGIS would return false, because the binary predicates always use the full precision. And the calculation of the very small distance of A or B to the line 1-2 can return a value much less than 0.00000001. If You use the expression distance < epsilon (where epsilon is the smallest absolute double precision value in Your programming language) for determining, if point A is on the line 1-2, then You will get false. You would have to use the expression distance < 0.00000001 to get the correct result.
>
> I hope, that I was able to describe the problem we have. It seems to me, that this is the same problem, that Todd Jellet described in the geos-devel mailing list (Todd: please correct me, if I am wrong). I don't know the algorithms, that are used in GEOS (I am a C# developer and did not look into the C code), so I do not know, if it is possible to solve the problem at all. But if there is a solution, I am hoping, that it will be implemented in JTS/NTS/GEOS/PostGIS.
>
> Best Regards,
> Gerhard
> --
> Gerhard Sommer
> Magistrat Wien (MA 14-ADV), Rathausstraße 1, 1082 Wien, Austria
> E-Mail: som at adv.magwien.gv.at
> Tel: +43 (1) 4000 91326
> Fax: +43 (1) 4000 99 91326
>
>
> > -----Original Message-----
> > From: postgis-users-bounces at postgis.refractions.net
> > [mailto:postgis-users-bounces at postgis.refractions.net] On
> > Behalf Of Martin Davis *EXTERN*
> > Sent: Tuesday, May 23, 2006 7:30 PM
> > To: PostGIS Users Discussion
> > Subject: RE: [postgis-users] Wrong result in spatial analysis
> > methodslikeintersects due to limited coordinate precision
> >
> > > JTS seems to have solved the problem, if I got the discussion right.
> >
> > Nope - JTS has exactly the same behaviour as you're seeing in
> > PostGIS (deliberately, since PostGIS uses GEOS, which is a
> > port of JTS).
> >
> > There's no immediate plans to enhance JTS to support a
> > tolerance value for predicates, but it's an interesting idea
> > and perhaps one that should be implemented.
> >
> > I'm just not sure that using a tolerance value is a panacea
> > for bad data.  For instance, what happens if you have data
> > which has very small "zingers" (spikes folded back on the
> > main polygon)?  These would test fine using a tolerance
> > value, but would obviously not be correctly aligned.  This
> > seems to be a complex problem which isn't going to be fixed
> > simply by allowing a tolerance in predicates.
> >
> > Martin Davis, Senior Technical Architect
> > Vivid Solutions Inc.      www.vividsolutions.com
> > Suite #1A-2328 Government Street Victoria, B.C. V8T 5G5
> > Phone: (250) 385 6040 - Local 308 Fax: (250) 385 6046
> >
> >
> > > -----Original Message-----
> > > From: postgis-users-bounces at postgis.refractions.net
> > > [mailto:postgis-users-bounces at postgis.refractions.net] On
> > > Behalf Of Sommer Gerhard
> > > Sent: May 22, 2006 10:06 AM
> > > To: postgis-users at postgis.refractions.net
> > > Subject: [postgis-users] Wrong result in spatial analysis
> > > methods likeintersects due to limited coordinate precision
> > >
> > >
> > > Hi,
> > >
> > > I am using PostGIS 1.1.2 with PostgreSQL 8.1 and I have
> > > detected the following inconsistency, which can be
> > > reconstructed with the following SQL statements:
> > >
> > > select
> > > astext(line_interpolate_point(geomfromtext('LINESTRING(200
> > > 340000.00, 600 340100.00)'), 0.7345678965533442)); --Result
> > > -> "POINT(493.827158621338 340073.456789655)" select
> > > distance(geomfromtext('LINESTRING(200 340000.00, 600
> > > 340100.00)'), geomfromtext('POINT(493.827158621338
> > > 340073.456789655)')); --Result -> 0.000000000296534992623179
> > > select intersects(geomfromtext('LINESTRING(200 340000.00, 600
> > > 340100.00)'), geomfromtext('POINT(493.827158621338
> > > 340073.456789655)')); --Result -> f
> > >
> > > The last statemant shows that the point, which I have before
> > > interpolated on the line, is not on the line according to the
> > > "intersects" operation. The reason for this is clear. It is
> > > because the precision of any stored ccordinate values ist
> > > restricted by the number of bytes available for storage (in
> > > double precision this is approximately 15 digits). The
> > > problem is also described in
> > > http://www.vividsolutions.com/jts/discussion.htm#dimensionalCo
> > llapse and JTS seems to have solved the problem, if I got the
> > discussion right.
> >
> > In my opinion the intersects operation and the other spatial
> > anlysis methods should and could always return the correct
> > result. This could be done by defining a fuzzy tolerance
> > (either explicitly as a parameter of the operations or
> > implicitly depending on the precision of data storage). Then,
> > if the point in the above example was within this tolerance,
> > the function would return "true", as it should in this case.
> > In the above example the tolerance could be i.e. 0.000000001.
> >
> > It would be great, if this could be considered as an
> > enhancement to PostGIS in a future version.
> >
> > The reason, why this is a problem at our organization is,
> > that when editing polygons, the edges of the polygon are
> > snapped to another polygon. Later we need to test, which
> > polygons are "not disjoint" to the other polygon. But this
> > query often delivers the wrong results, as the snapped points
> > are partly outside and partly inside the other polygon, which
> > results in an intersection, that would not be there, if we
> > had unlimited precision available or if the limited precision
> > would be considered in the PostGIS function.
> >
> > Best Regards
> > --
> > Gerhard Sommer
> > Magistrat Wien (MA 14-ADV), Rathausstraße 1, 1082 Wien, Austria
> > E-Mail: som at adv.magwien.gv.at
> > Tel: +43 (1) 4000 91326
> > Fax: +43 (1) 4000 99 91326
> >
> > _______________________________________________
> > postgis-users mailing list postgis-users at postgis.refractions.net
> > http://postgis.refractions.net/mailman/listinfo/postgis-users
> > _______________________________________________
> > postgis-users mailing list
> > postgis-users at postgis.refractions.net
> > http://postgis.refractions.net/mailman/listinfo/postgis-users
> >
>
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>
>
>
>