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

Wed May 24 04:45:50 PDT 2006

Hi Nicolas,

If I would snap the data to the real precision (in our case 0.005 meters), I would have the same problem. What I would need is a storage with infinite precision, so that I could store the exact position of any calculated point on a line segment. It does not matter how wide the grid is. The problem is, that there is a grid at all and that any calculated point that I store might be moved from the actual position to the next grid point, if it is not already by coincidence on a grid point.

But if the binary predicates were aware of the fact, that there is only a finite precision awailable for stored geometry, it could deliver correct results in spite of the limited precision.

If a distance calculation is used to determine, if a point is exactly on a line segment, and the distance is very small (0.0000000523457989383956, which has 15 significant digits in double precision), this distance has a much higher precision than the stored coordinates and the result is, that the point is not exactly on the line segment. But in fact for my coordinate precision, which might be 0.005 meters, the point is on the line segment, and I want to get true as a result. But this is only possible, if the software allows me to tell it, what my precision is. Or if the software automatically calculates the grid size, which depends on the number of digits before the decimal point of the biggest of my coordinate values (using 7 to 8 digits before the decimal point would be enough for any projection worldwide and would result in a grid size of 0.0000001 meters, using geographic coordinates the grid size would be 0.000000000001 degrees). 

For practical purposes the default tolerance value, if the user does not enter his own, could be derived from the SRID of the geometries, which is stored with the geometries anyway (is it a geographic coordinate system or any projected coordinate system?). This would be faster than calculating the maximum coordinate value and then calculating the "exact" precision for these geometries.

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 Nicolas Ribot *EXTERN*
> Sent: Wednesday, May 24, 2006 11:17 AM
> To: PostGIS Users Discussion
> Cc: todd.jellet at caris.com
> Subject: Re: [postgis-users] Wrong result in spatial 
> analysismethodslikeintersects due to limited coordinate precision
> 
> 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
> >
> >
> >
> >
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
> 