Precision of coordinates in spatial index

[Roman]

Hello,

I am fighting with some spatial operations that don't work as expected. I have found that this is (probably) due to incorrect coordinates in spatial index. To explain it quickly, take a look at attached example. Red rectangle is an extent from spatial index of a yellow line. Extent coordinates are not exactly min a max coordinates of a line but they are somehow rounded. In this particular example, for a linestring

LINESTRING(-651396.04 -1069001.71, -651403.23 -1069028.47)

there is a record in spatial index with this extent

xmin=-651403.250000, xmax=-651396.062500, ymin=-1069028.500000, ymax=-1069001.750000

This causes that this linestring is ignored when finding relationships with other geometries. 

Why coordinates in spatial index are different from min, max coordinates of geometries?

Why is the extent shifted? Is it causes by negative values of coordinates?

I tried to build spatial index in Spatialite GUI 1.4.0 and 1.4.1 (spatialite 2.4.0 RC4 and RC5) with same results.

Roman

[a.furieri]

Hi Roman,

this is a well known problem; please see my previous post here:
http://groups.google.com/group/spatialite-users/browse_thread/thread/282f27d38995e96d

very shortly resumed:
- geometry coordinates are stored as 64 bit floating point (double
precision)
- R*Tree coordinates are stored as 32 bit floating point (single
precision)
obviously, some slight precision loss occurs.

Anyway, caring for this is absolutely simple: you are simply required
to "enlarge" a little bit the BBOX used to query the R*Tree.

starting since v.3.0.0-alpha you can safely use the VirtualIndex
interface
and/or "geometry callback" methods (RTreeIntersects, RTreeWithin and
RTreeContain), because all them implicitly take care to add such BBOX
"extra margins" so to recover for any adverse rounding/truncation
side effect.

bye Sandro

[Roman]

Sandro,

64bit vs. 32bit floating point values could explain a difference after seventh significant digit (in my described example it is 4cm) but I have to enlarge BBOX by about 6 meters to get the related line. So it has to be some other problem. On the other hand, even with not enlarged BBOX, spatial query returns geometries whose BBOX are several meters far (disjoint) from requested BBOX. I will investigate it further and try to prepare testing case with very few geometries.

Roman

P.S. VirtualSpatialIndex is great simplification, I am looking forward to version 3.0.0.

[Roman]

My first attempt with enlarged BBOX was false. Now it's working OK.

Thank you for suggestions.

Roman

[a.furieri]

Hi Roman,

this elementary C snippet exactly corresponds
to the LINESTRING in your first report:
-------------------------------------------
#include <stdio.h>
void main(void)
{
double dbl_a = -651396.04;
double dbl_b = -1069001.71;
double dbl_c = -651403.23;
double dbl_d = -1069028.47;
float flt_a = dbl_a;
float flt_b = dbl_b;
float flt_c = dbl_c;
float flt_d = dbl_d;
printf("DBL: %1.4f %1.4f %1.4f %1.4f\n",
dbl_a, dbl_b, dbl_c, dbl_d);
printf("FLT: %1.4f %1.4f %1.4f %1.4f\n",
flt_a, flt_b, flt_c, flt_d);
}
-------------------------------------------

here is the generated output:
-------------------------------------------
DBL: -651396.0400 -1069001.7100 -651403.2300 -1069028.4700
FLT: -651396.0625 -1069001.7500 -651403.2500 -1069028.5000

As you can easily check, FLT strictly corresponds to values
you are reporting for xmin, ymin, xmax and ymax (R*Tree).
And this is the final and definitive answer: a slight
rounding (truncation) actually occurs when converting
64 bit doubles into 32 bit floats.
The average difference is very narrow: about 2-4 cm;
anyway, that's enough to cause some inconsistency.

"enlarging" the query bounding box by 10 cm (four sides) surely
will definitively eradicate any possible R*Tree related issue.

bye Sandro

[Roman]

"enlarging" the query bounding box by 10 cm (four sides) surely
will definitively eradicate any possible R*Tree related issue.

Sandro,

I came to the same conclusion and added exactly 0.1 (10cm). Question is, how to solve it generally for queries to VirtualIndex where extra margin is added implicitly. Margin should be large enough to solve the problem with rounding to 32bit floating point but not too large to select unnecessary geometries. Adding 0.1 is suitable for my particular dataset in meters but not for datasets in other coordinate systems, e.g.. SRID 4326 where appropriate value should be 0.00001 (eighth significant digit - counting 3 possible digits before decimal point). Possible solution: 1. find the extent of a spatial table 2. take largest value of extent and find at what position is eighth significant digit 3. add appropriate margin in VirtualIndex query. 

Roman

[a.fu...@lqt.it]

The latest v.3.3.0 *automatically* adds
an appropriate extra-margin to the BBOX,
so to compensate for any rounding/truncation.

This is silently applied by VirtualSpatialIndex
and by RTreeWithin(), RTreeContains() and
RTreeIntersects() [Geometry CallBacks].

for the most curious about implementation details,
please see the following code snippet:

-------------------
/* adjusting the MBR so to compensate for DOUBLE/FLOAT truncations */
float fminx = xmin;
float fminy = ymin;
float fmaxx = xmax;
float fmaxy = ymax;
double tic = fabs (xmin - fminx);
double tic2 = fabs (ymin - fminy);
if (tic2 > tic)
tic = tic2;
tic2 = fabs (xmax - fmaxx);
if (tic2 > tic)
tic = tic2;
tic2 = fabs (ymax - fmaxy);
if (tic2 > tic)
tic = tic2;
tic *= 2.0;
mbr->minx = xmin - tic;
mbr->miny = ymin - tic;
mbr->maxx = xmax + tic;
mbr->maxy = ymax + tic;
-------------------

the above method doesn't depends on scale:
it works using UTM coordinates (e.g. 1234567.23);
and it works as well using long-lat coordinates
(e.g. 12.3456789).

bye Sandro