equality of SquareGrid polygons

[David Anderson]

I am working out a process to convert a set of irregular features into a set of regular features.  The irregular features are a layer of vegetation cover.
There are the multi colored areas.

I run the SquareGrid function against the vegetation layer.  I use the default 0,0 for the origin.
It correctly generates 7 squares that cover the vegetation poygons in the area shown below.

I only want 1 regular/square polygon.

My current approach is to group by the geometry field, then select the one with the square that has the maximum area of a vegetation polygon.
The problem is that the groupby returns 3 square polygons for this area.
Using the ST_EQUALS function the 3 polygons are not equal.
Yet the areas are the same out to 10 decimal places, the centroids return the same X and Y values.  So I am wondering why the squares are not equal.

I am also wondering why 4 of the squares are equal to at least one of the three squares returned.  It seems like they should all be equal or not should be equal.  





It might be floating point problem.  Here is the WKT of two the almost equal polygons:
POLYGON((530154.569 5256880.132,530296.8159999999 5256880.132,530296.8159999999 5257022.379,      530154.569 5257022.379,      530154.569 5256880.132))
POLYGON((530154.569 5256880.132,530296.8159999999 5256880.132,530296.8159999999 5257022.379000001,530154.569 5257022.379000001,530154.569 5256880.132))

They differ only by 0.000000001 for two points.  Which appears to be enough to be not equal.

Perhaps they could be some rounding in the SquareGrid function?

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

On Fri, 21 Jul 2017 11:44:27 -0700 (PDT), David Anderson wrote:
> I am working out a process to convert a set of irregular features
> into
> a set of regular features.  The irregular features are a layer of
> vegetation cover.
> There are the multi colored areas.
>
> I run the SquareGrid function against the vegetation layer.  I use
> the default 0,0 for the origin.
> It correctly generates 7 squares that cover the vegetation poygons in
> the area shown below.
>
> I only want 1 regular/square polygon.
>
> My current approach is to group by the geometry field, then select
> the
> one with the square that has the maximum area of a vegetation
> polygon.
> The problem is that the groupby returns 3 square polygons for this
> area.
> Using the ST_EQUALS function the 3 polygons are not equal.
> Yet the areas are the same out to 10 decimal places, the centroids
> return the same X and Y values.  So I am wondering why the squares
> are not equal.
>

Hi David,

I not sure to correctly understand what are you doing so to
replicate a testcase.
the attached figure and WKT samples aren't enough: a detailed
specification of all SQL statements you're using will surely be
more useful.

> I am also wondering why 4 of the squares are equal to at least one of
> the three squares returned.  It seems like they should all be equal
> or not should be equal.  
>

> It might be floating point problem.  Here is the WKT of two the
> almost equal polygons:
> POLYGON((530154.569 5256880.132,530296.8159999999
> 5256880.132,530296.8159999999 5257022.379,      530154.569
> 5257022.379,      530154.569 5256880.132))
> POLYGON((530154.569 5256880.132,530296.8159999999
> 5256880.132,530296.8159999999 5257022.379000001,530154.569
> 5257022.379000001,530154.569 5256880.132))
>
> They differ only by 0.000000001 for two points.  Which appears to be
> enough to be not equal.
>
> Perhaps they could be some rounding in the SquareGrid function?
>

I strongly doubt that SquareGrid could be specifically prone to
rounding errors, but for sure all SQL functions involving some
GEOS operation are, most notably on 32 bit x86 platforms, as it
recently emerged in a discussion with GEOS developers.

a little known technical fact is that all modern x86 CPUs
(starting since year 2000, Pentium III generation) include
in their HW _TWO_ different floating point processors:
- the old 387, supporting 80 bit internal registers
- the new SSE, supporting 64 bit internal registers
(and the SSE fp processor is someway faster than the 387)

all C/C++ compilers (gcc, msvc and so on) always use by default
the 387 fp processor when compiling 32 bit code (mainly to
preserve historical compatibility), whilst they automatically
switch to SSE for 64 bit code.

the problem with 387 is in that its odd 80 bit internal
registers doesn't fit standard DOUBLE program variables
that always are 64 bit.
so every time that a 387 fp result is transferred from
CPU to memory a conversion from 80 to 64 bit happens, and
this can easily cause weird truncation/rounding effects.

the next version of GEOS will probably fix this issue
(SSE should be always switched on, also when compiling
32 bit code).

bye Sandro

[David Anderson]

Sandro,
The work flow is rather cryptic.
I am adding everything that I have done to get to the point of the original message.

The data is here:

https://www.cloudvault.usda.gov/index.php/s/Lq8M9gdPFuc296z

Here is the SQL that I used

create table square_step1
as
select objectid,ht_grp,species,density,size_class,squaregrid(shape,142.247) as square_geom,shape as orig_geom
from veg4simpplle

select recovergeometrycolumn('square_step1','square_geom',
(select srid(square_geom) from square_step1 limit 1),
(select geometrytype(square_geom) from square_step1 limit 1))

select elementarygeometries('square_step1','square_geom','square_step2','new_oid','old_oid')

create table square_step3
as
select square_geom,max(st_area(st_intersection(square_geom,orig_geom))) as maxarea
from square_step2
group by square_geom

select recovergeometrycolumn('square_step3','square_geom',
(select srid(square_geom) from square_step3 limit 1),
(select geometrytype(square_geom) from square_step3 limit 1))

After getting more squares than I thought should be in the resultant data set,  I did some looking at the data in GIS where I found some polygons that should be the same and combined in the groupby.
The SQL to produce the WKT values is

select a.gid,b.gid,

aswkt(a.square_geom),

aswkt(b.square_geom)

from square_step3 a,square_step3 b

where a.gid=48769

and b.gid=48770

I hope this helps in any error tracking.

I am running a 32 bit SQLite database front end.  I'd like to switch to a 64 bit one but I can't find one that I like.

David

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

Hi David,

yes, I can confirm that the reported issue is real
and indifferently affects both 32 and 64 bit.

just a casual finding: using a slightly different
cell size value (142.25 instead of 142.247) the
issue magically vanishes and everything works
exactly as expected (thus indirectly confirming
that the SpatiaLite's code itself is sound and
safe).

after an extenuating debugging session I was
finally able to identify the real cause, and
I was very surprised ... it's SQLite that
sometimes badly converts SQL floating point
constants into C double variables.

the following is a simple test that you can
easily replicate by yourself:

SELECT SetDecimalPrecision(20);
SELECT 142.5, 142.25, 142.247, 142.244, 142.255;

note: this is a pure SQL expression just involving
few float constants; spatialite hasn't any role in
this, except for printing all returned values
with an extended precision (20 decimal digits).

and here are the results:
142.500000000000000000
142.250000000000000000
142.247000000000014097
142.243999999999999773
142.254999999999995453

as you can easily check, in many cases SQLite
has arbitrarily introduced an infinitesimal
difference (sometimes in excess, other times
in defect).

the difference is so small that only the
less significant bit is affected, but this
is enough to introduce unexpected shifts
in all subsequent calculations based on
that constant value.

possible workarounds:
=======================================
A) build a single grid covering all your polygons;
something like this:

CREATE TABLE fusion AS
SELECT ST_Union(shape) AS geom
FROM veg4simpplle;
SELECT RecoverGeometryColumn('fusion', 'geom', 26912,
'MULTIPOLYGON', 'XY')
CREATE TABLE square_grid AS
SELECT SquareGrid(geom, 142.247) AS geom
FROM fusion;
SELECT RecoverGeometryColumn('square_grid', 'geom',
26912, 'MULTIPOLYGON', 'XY');
SELECT ElementaryGeometries('square_grid', 'geom',
'square_grid_elem', 'new_oid', 'old_oid');

B) use a cell size value immune from nasty rounding
effects (as e.g. 142.25).

bye Sandro

[David Anderson]

Sandro, Thanks for looking into this questions.  I apologize for not responding earlier but I got caught up in implementing your solution.

David