Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have this polygon

POLYGON((1.2E-15 1.1E-15,1.31E-15 1.21E-15,-180 -90,180 90,1.2E-15 1.1E-15)) 

It uses an extreme coordinates like (1.2E-15 1.1E-15) and (180 90) When I call ST_isValidReason, it says that polygon is invalid:

SELECT ST_isValidReason( ST_GeomFromText(
    'POLYGON((1.2E-15 1.1E-15,1.31E-15 1.21E-15,-180 -90,180 90,1.2E-15 1.1E-15))') );
st_isvalidreason       
Self-intersection[-180 -90] 

But when I reduce number precision form E-15 to E-12 it works ok:

SELECT ST_isValidReason( ST_GeomFromText(
    'POLYGON((1.2E-12 1.1E-12,1.31E-12 1.21E-12,-180 -90,180 90,1.2E-12 1.1E-12))') );
st_isvalidreason 
Valid Geometry

The question is: is it a PostGis bug or I shouldn't use number like 1.1E-15 for coordinates?

share|improve this question

1 Answer 1

up vote 3 down vote accepted

There are two general reasons to avoid tiny coordinates (or tiny differences of coordinates). Because I cannot (readily) find any PostGIS documentation about its internal storage format for coordinates, I will give them both.

Integer representation

Many GISes convert coordinates into indexes into an extremely fine grid of points covering a study region or "area of interest." This enables them to do all geometric calculations using integral coordinates (just row and column indexes). Textbook algorithms can be used almost unchanged when integral coordinates are involved; storage requirements may be halved (an integer takes four bytes while a double takes eight bytes); and possibly the calculations will be faster on some processors.

Using four bytes (=32 bits) for integer values limits these grids to 2^32 rows and columns. For geographic coordinates, values up to 180 in size may have to be stored. So that the grid can cover such an extent, the interval between neighboring grid points cannot be any less than 180 / 2^32, which is about 4.2E-8. Values between -2.1E-8 and +2.1E-8 therefore would all coincide with each other and with 0.

Whence, if there is any possibility that your data might be processed by a GIS using integers for its calculations, you cannot rely on it distinguishing numbers like 1.2E-15 and 1.31E-15 from each other.

Double-precision floating point issues

I suspect this is not the problem, though. Although the error message is strange, one interpretation is that PostGIS detected a self-intersection during sequential processing of the polygon coordinates and the detection occurred when it was processing (-180, -90). This means it thinks (somehow) that there is a problem with the first two edges.

The equations to compute a self-intersection come down to inverting a matrix. The matrix coefficients in this case are

-180.     -90.
-1.1E-16   -1.1E-16

This matrix is extraordinarily ill-conditioned: its eigenvalues of -180. and -5.5E-17 cover 18.5 orders of magnitude, giving a condition number of 3.3E+18. This enormous condition number means that numerical procedures to invert the matrix will be unstable, even with double precision arithmetic. It is likely this instability led to a spurious intersection of the first two edges of the polygon, resulting in the error.

18.5 orders of magnitude is greater than the 16 to 17 significant figures of precision available with double precision floats. In effect, all bets are off concerning what solution the software will find.

When E-15 is changed to E-12, the condition number drops three orders of magnitude. This time the software can find a solution--but probably only to one significant figure. Perhaps by accident, that solution does not indicate a self-intersection. But you're not out of the woods yet: one significant figure is terrible precision for GIS work (or almost any other numerical work). There will be concern about these precision issues until the E-15 is increased to around E-6 or so. The latter value represents about 0.11 meter (at the most): see http://gis.stackexchange.com/a/8674 for how this distance was computed. For almost any application, getting a precision of 0.1 meters when mapping at the scale of the entire earth is more than good enough.

share|improve this answer
    
great explanation, thanks! –  GetUsername Jul 5 '12 at 7:11

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.