Decoding lat/long in file name

[johnmcd]

Hi,

I'm something of a novice Python programmer so please bear with me. OpenStreetMaps uses a 'shortlink' protocol to create file names from lat/long coordinates, which is described here: https://wiki.openstreetmap.org/wiki/Shortlink. They also include links to various coding implementations; however, the link to the Python implementation only includes an encoding routine, not a decoding one, and I need to be able to decode filenames. Bit twiddling in Python is a bit beyond my current skill level, so I grabbed the decode routine from the .js example and ran it through an online language converter to get Python. I ran it and the decoded lat/long coordinates were incorrect, so I tried it with the .java and .pl examples and got the exact same incorrect results, so I assume it's something to do with the formula. Here's the code I'm using (_decode was converted from .js):

ARRAY = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789_~'
import math

def _encode(lat, lon, z):
    """given a location and zoom, return a short string representing it."""
    x = int((lon + 180.0) * 2**32 / 360.0)
    y = int((lat +  90.0) * 2**32 / 180.0)
    code = _interleave(x, y)
    str = ''
    # add eight to the zoom level, which approximates an accuracy of
    # one pixel in a tile.
    for i in range(int(math.ceil((z + 8) / 3.0))):
        digit = (code >> (56 - 6 * i)) & 0x3f;
        str += ARRAY[digit]
    # append characters onto the end of the string to represent
    # partial zoom levels (characters themselves have a granularity
    # of 3 zoom levels).
    for i in range((z + 8) % 3):
        str += "-"
    return str
    
def _decode(sc):
    char_array = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789_~"
    i = 0
    x = 0
    y = 0
    z = -8
    for i in range(len(sc)):
        ch = sc[i]
        digit = char_array.index(ch)
        if digit == -1:
            break
        # distribute 6 bits into x and y
        x <<= 3
        y <<= 3
        for j in range(2, -1, -1):
            x |= (0 if (digit & (1 << (j+j+1)) == 0) else (1 << j))
            y |= (0 if (digit & (1 << (j+j)) == 0) else (1 << j))
        z += 3

    x = x * pow(2, 2-3*i) * 90 - 180
    y = y * pow(2, 2-3*i) * 45 - 90

    # adjust z
    if i < len(sc) and sc[i] == "-":
        z -= 2
        if i+1 < len(sc) and sc[i+1] == "-":
            z += 1

    return [y, x, z]

def _interleave(x, y):
    """combine 2 32 bit integers to a 64 bit integer"""
    c = 0
    for i in range(31, 0, -1):
      c = (c << 1) | ((x >> i) & 1)
      c = (c << 1) | ((y >> i) & 1)
    return c
    
encoded = _encode(50.671530961990356, 6.09715461730957, 16)
print(encoded)

decoded = _decode(encoded)
print(decoded)

Here are the results I got:

>python decode_short.py
0GAjIv8h
[1035.3722476959229, 1308.7772369384766, 16]

The encoding is correct, but the decoded values are obviously way off. Can anyone provide and hints/tips/guidance on what I should look at?

Thanks

John

[deanhystad]

The problem with your code is here:

    x = x * pow(2, 2-3*i) * 90 - 180
    y = y * pow(2, 2-3*i) * 45 - 90

i should be 8, but it is only 7, so x and y are roughly 8 times larger than they should be.

The reason for the discrepancy is that range works differently than a for loop in js or C.
Python

for i in range(8): pass
print(i)

prints 7

C++

for(i=0; i< 8; i++) ;
cout << i

prints 8

This appears to work.

def decodejs(sc):
    index_of = {x: y for y, x in enumerate("ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789_~")}
    x = y = 0
    z = -8
    i = 0
    for ch in sc:
        if (digit := index_of.get(ch)) is None:
            break

        # distribute 6 bits into x and y
        x <<= 3
        y <<= 3
        for j in range(2, -1, -1):
            if digit & 1 << (2 * j + 1):
                x |= 1 << j
            if digit & 1 << (2 * j):
                y |= 1 << j
        z += 3
        i += 1

    x = x * 2**(2 - 3 * i) * 90 - 180
    y = y * 2**(2 - 3 * i) * 45 - 90

    if i < len(sc) and sc[i] == "-":
        z -= 2
        if i + 1 < len(sc) and sc[i + 1] == "-":
            z += 1

    return z, y, x


print(decodejs('0GAjIv8h'))

[johnmcd]

Dean,

That's fantastic - thank you! I'm going to sit down and go through everything so I can learn what you did and why - my goal is to learn something new with Python every day.

John

[deanhystad]

This may be more "pythonic". Added some comments now that I understand what is going on.

def decodejs(sc):
    codes = "-ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789_~"
    x = y = z = i = 0
    for ch in sc:
        # Convert asci letter to int using codes.  '-' signals early end of bytes
        if (digit := codes.index(ch)-1) < 0:
            break

        # Each digit contains 3 bits of latitude (y) and 3 bits of longitude (x)
        # interleaved like this y2x2y1x1y0x0.  Extract y2y1y0 and x2x1x0
        # and add to the latitude and longitude.
        x <<= 3
        y <<= 3
        for xbit, ybit, bit in zip((2, 8, 32), (1, 4, 16), (1, 2, 4)):
            x |= bit if digit & xbit else 0
            y |= bit if digit & ybit else 0
        z += 3
        i += 1

    # This scales latitude (x) from 0...0xFFFFFF to -180...180 and
    # the latitude (y) from from 0...0xFFFFFF to -90...90.
    x = x * 2**(2 - 3 * i) * 90 - 180
    y = y * 2**(2 - 3 * i) * 45 - 90

    if i < len(sc) and sc[i] == "-":
        z -= 2
        if i + 1 < len(sc) and sc[i + 1] == "-":
            z += 1

    return z-8, y, x


print(decodejs('0GAjIv8h'))

[johnmcd]

Thanks again - the comments are extremely useful.