ST_HexagonGrid produces overlapping hexagons

[dbaston]

SELECT * FROM  (SELECT (ST_HexagonGrid(5, '01030000000100000005000000fdffffffff7f66c0fcffffffff7f56c0fdffffffff7f66c07cc0e71a95e8544000000000008066407cc0e71a95e854400000000000806640fcffffffff7f56c0fdffffffff7f66c0fcffffffff7f56c0')).*) sq  WHERE i = 5 AND j IN (0, 1);

[Algunenano]

An example that could be used as test based on that:

WITH geoms AS
(
    SELECT * FROM  (SELECT (ST_HexagonGrid(5, '01030000000100000005000000fdffffffff7f66c0fcffffffff7f56c0fdffffffff7f66c07cc0e71a95e8544000000000008066407cc0e71a95e854400000000000806640fcffffffff7f56c0fdffffffff7f66c0fcffffffff7f56c0')).*) sq  WHERE i = 5 AND j IN (0, 1)
),
j0 AS
(
    SELECT * FROM geoms WHERE i = 5 AND j = 0
),
j1 AS 
(
    SELECT * FROM geoms WHERE i = 5 AND j = 1
)
SELECT
    ST_Intersects(j0.geom, j1.geom),
    ST_AsText(ST_Intersection(j0.geom, j1.geom))
FROM j0,j1;

The current output is:

 st_intersects |                                                   st_astext                                                    
---------------+----------------------------------------------------------------------------------------------------------------
 t             | POLYGON((35 8.66025403784439,40 8.66025403784439,40 8.66025403784439,35 8.66025403784439,35 8.66025403784439))

When transformed into text, the latitude is the same for all points but probably due to some rounding the binary value will be different, thus the intersection is a polygon and not a point.

I'm not sure how hard this would be to fix.

[Algunenano]

I have an improved version in ​https://github.com/postgis/postgis/pull/544 but I'm certain it is architecture dependent.

I think that if it doesn't break Windows we could push it to improve the precision between tiles.

[pramsey]

Ug, well an improvement is nice. I fear the only real full-on solution is the one Komzpa put in the old PR: feeding in any pre-existing coordinates to the generation function so they can be re-used… so the HexagonGrid function would have to carry a bunch of prior state, and the Hexagon function would have to be able to make use of said state. And that would still leave cases where independently generated hexagons from the same tiling space didn't have perfectly matching vertex values.

[Algunenano]

I've modified the PR to change how the coordinates are calculated to use something like:

pt.y = origin_u + tile_height * n

N will hopefully more precise value (an integer or an integer.5), but what's more important when changing to a different/adjacent tile it will do the same calculations.

For example cell_j_10 + 0.5 (10.5) should be the same as cell_j_11 - 0.5 (10.5). Technically there could still be some small differences (likely with floating point numbers), which could be solved by using ints (everything multiplied by 2, then at the very end transform it to float and divide it by 2) but I think it's an unlikely scenario.

[Raúl Marín <git@…>]

In 95bff56/git:

ST_HexagonGrid: Improve precision

Change the order of the operations so adjacent tiles have the exact
same coordinates

Closes #4627
Closes ​https://github.com/postgis/postgis/pull/544