Error in the implementation of tile box boundaries

[ysig]

Hi,

thank you for this really useful library.

I noticed an unexpected behavior, when I run:

from pyquadkey2 import quadkey
def in_boundaries(lat, lon, res):
    qk = quadkey.from_geo((lat, lon), res)
    min_lat, min_lon = qk.to_geo(quadkey.TileAnchor.ANCHOR_SW)
    max_lat, max_lon = qk.to_geo(quadkey.TileAnchor.ANCHOR_NE)
    assert min_lat <= lat <= max_lat                   # error

for

in_boundaries(-27.052395, 152.97702, 6)

I get an assertion returns AssertionError is being raised, which is a behavior that clearly shouldn't happen as each point should be included in its associated rectangular box (or I'm missing something).

Thank you,

[ymorin-orange]

assert min_lat <= lat <= max_lat

~This looks fishy. One of the comparison will return True, and that will be compared to the last element.~

~I.e. if you add parentheses, it is the same as (min_lat <= lat) <= max_lat, which would be evaluated as True <= max_lat, which then can't evaluate to True.~

~Have you tried printing the three values, and splitting the assert in two?~

Forget that, that was a lack of coffee speaking...

[muety]

Sorry, didn't receive any e-mail notifications on this repo in the past time and thus totally missed all recent issues. Will have a look into this later today!

[muety]

Looks like we're indeed picking up the wrong tile here. Your code yields quadkey 311213, while testing with jquad returns 311211 (one tile above). Perhaps a floating point precision / rounding issue? Definitely something we'll need to look into! Thanks for bringing this up!

[muety]

This is what it looks like on the map: blue (1) is your query point, green (2 and 3) are SW and NE reported by pyquadkey2 and red (4 and 5) are SW and NE reported by jquad (which are probably the correct ones).

fid,lat,lon,color
1,-27.052395,152.97702,1
2,-31.952162238025,151.875,2
3,-27.059125784374,157.5,2
4,-27.05912578437406,151.875,3
5,-21.94304553343818,157.5,3

[muety]

Actually, the key returned by pyquadkey2 is correct, at least Microsoft's reference implementation yields the same result:

using Microsoft.MapPoint;

double latitude = -27.052395;
double longitude = 152.97702;
int levelOfDetail = 6;

string quadKey = GetQuadKey(latitude, longitude, levelOfDetail);

Console.WriteLine("QuadKey: " + quadKey);  // returns '311213'

static string GetQuadKey(double latitude, double longitude, int levelOfDetail)
{
    int pixelX, pixelY;
    TileSystem.LatLongToPixelXY(latitude, longitude, levelOfDetail, out pixelX, out pixelY);

    int tileX, tileY;
    TileSystem.PixelXYToTileXY(pixelX, pixelY, out tileX, out tileY);

    return TileSystem.TileXYToQuadKey(tileX, tileY, levelOfDetail);
}

But this would mean both ethlo/jquad and CartoDB/python-quadkey are wrong? :thinking:

So, the bug must be somewhere in the tile -> bbox calculation part.

[muety]

I went through every line of the code, but couldn't find the issue. I honestly have no clue what's wrong here. Would love to get some help on this!

[cheginit]

I think the issue is here:

https://github.com/muety/pyquadkey2/blob/8222ca5b65e2435ad1e35916323018519c583e23/src/pyquadkey2/quadkey/tilesystem/tilesystem.pyx#L63

You can check this out in the code below:

import math

def get_quadkey(lat: float, lng: float, zoom: int)->str:
    x = int((lng + 180) / 360 * (1 << zoom))
    y = int((1 - math.log(math.tan(math.radians(lat)) + 1 / math.cos(math.radians(lat))) / math.pi) / 2 * (1 << zoom))
    quadkey = ''
    for i in range(zoom, 0, -1):
        digit = 0
        mask = 1 << (i - 1)
        if (x & mask) != 0:
            digit += 1
        if (y & mask) != 0:
            digit += 2
        quadkey += str(digit)
    return quadkey

get_quadkey(-27.052395, 152.97702, 6)

If you replace int with round you'll get 311213 instead of 311211. In python, int truncates, but round either rounds up or down.

[muety]

Thanks a lot for your help on this!