[postgis-users] Wrong result in spatial analysis methods like intersects due to limited coordinate precision

[Sommer Gerhard]

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#dimensionalCollapse 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