leo.blog();

Quadtiles

Some tasks require to you to communicate or store your geographical location. In a local application, you can store it as two floating point numbers and never worry about anything else.

For this purpose, we can use Quadtiles. In our Python program, we can represent a quadtile as the four corners of a rectangular area. I’ll call these xFrom, xTo, yFrom and yTo, but it might be more accurate to use terminology like latitude and longitude.

In [2]:
Quadtile = namedtuple("Quadtile", "xFrom xTo yFrom yTo")

Quadtiles are about splitting the current Quad into 4 tiles, picking one of them and repeating. Since we want to be able to point to any location on Earth, our initial Quad will always be the whole world.

In [3]:
WORLD = Quadtile(-180, 180, 90, -90)

WORLD
Out [3]:
Quadtile(xFrom=-180, xTo=180, yFrom=90, yTo=-90)

Getting the center point

When working with locations, we often need a single point instead of a Quad, so let’s write a function to get the center point. We will calculate this by getting the midpoint value of the X and Y ranges.

In [4]:
@property
def center(quad):
    x = (quad.xFrom + quad.xTo) / 2
    y = (quad.yFrom + quad.yTo) / 2
    return (x, y)
Quadtile.center = center
In [5]:
WORLD.center
Out [5]:
(0.0, 0.0)

Makes sense, the center of the world is at (0, 0). This is the value we will get when we try to compute the location represented by the empty string, i.e. the default location.

Splitting quads

In [7]:
@property
def left(quad):
    xMid, yMid = quad.center
    
    return Quadtile(quad.xFrom, xMid, quad.yFrom, quad.yTo)
Quadtile.left = left

@property
def right(quad):
    xMid, yMid = quad.center
    
    return Quadtile(xMid, quad.xTo, quad.yFrom, quad.yTo)
Quadtile.right = right

WORLD.left.center
WORLD.right.center
Out [7]:
(-90.0, 0.0)
Out [7]:
(90.0, 0.0)

Now let’s do it in the other axis.

In [12]:
@property
def top(quad):
    xMid, yMid = quad.center
    
    return Quadtile(quad.xFrom, quad.xTo, quad.yFrom, yMid)
Quadtile.top = top

@property
def bottom(quad):
    xMid, yMid = quad.center
    
    return Quadtile(quad.xFrom, quad.xTo, yMid, quad.yTo)
Quadtile.bottom = bottom

WORLD.top.center
WORLD.bottom.center
Out [12]:
(0.0, 45.0)
Out [12]:
(0.0, -45.0)

We can use these to split the world into increasingly smaller chunks in order to represent the location we are interested in. For example;

In [22]:
WORLD.top.left.center
WORLD.top.left.bottom.right.center
WORLD.top.left.bottom.right.top.left.top.right.center
Out [22]:
(-90.0, 45.0)
Out [22]:
(-45.0, 22.5)
Out [22]:
(-56.25, 39.375)

A simple encoding

Using the Quadtile methods shown above, we will make

In [23]:
def decode_bits(bits):
    is_x = True
    
    quad = WORLD
    for bit in bits:
        bit = int(bit)
        
        if is_x:
            quad = quad.right if bit else quad.left
        else:
            quad = quad.bottom if bit else quad.top
        
        is_x = not is_x
        
    return quad.center

decode_bits("01001")
decode_bits("01100111")
Out [23]:
(-112.5, -22.5)
Out [23]:
(-56.25, -39.375)

Blah blah blah

In [15]:
def encode_location(coord, N):
    bits = ""
    quad = WORLD
    x, y = coord
    
    is_x = True
    
    for i in range(N):
        mid = quad.center
        if is_x:
            if x > mid[0]:
                bit = 1
            else:
                bit = 0
            quad = quad.right if bit else quad.left
        else:
            if y > mid[1]:
                bit = 1
            else:
                bit = 0
            quad = quad.bottom if bit else quad.top
        bits += str(bit)
        is_x = not is_x
    
    return bits

Let’s make a plot of how accurate our positions gets with more bits added. Lower numbers are more accurate.

Map visualizations

In [16]:
HOUSE = (-8.577507, 52.664838)
In [17]:
def plot_path(path, zoom):
    print(f"Path length: {len(path)}")
    fig, ax = plt.subplots(figsize=(8, 8), dpi=100)
    ax.xaxis.set_visible(False)
    ax.yaxis.set_visible(False)

    xs = [0.5]
    ys = [0.5]

    bits = ""
    for bit in path:
        bits += bit
        x, y = decode_bits(bits)
        x, y = tilemapbase.project(x, y)
        xs.append(x)
        ys.append(y)
        
    extent = tilemapbase.Extent.from_centre_lonlat(*decode_bits(bits), xsize=zoom)
    extent = extent.to_aspect(1.0)
    
    plotter = tilemapbase.Plotter(extent, t, width=400)
    plotter.plot(ax, t)
    
    _ = ax.plot(xs, ys)
    plt.show()
In [18]:
plot_path(encode_location(HOUSE, 8), 0.55)
Out:
Path length: 8
Out:
[ERROR] NameError: name 'plt' is not defined
---------------------------------------------------------------------------NameError                                 Traceback (most recent call last)<ipython-input-18-3899a8b637ef> in <module>
----> 1 plot_path(encode_location(HOUSE, 8), 0.55)
<ipython-input-17-cf8bea3f6b56> in plot_path(path, zoom)
      1 def plot_path(path, zoom):
      2     print(f"Path length: {len(path)}")
----> 3     fig, ax = plt.subplots(figsize=(8, 8), dpi=100)
      4     ax.xaxis.set_visible(False)
      5     ax.yaxis.set_visible(False)
NameError: name 'plt' is not defined

We can see that 8 bits is enough to get us next to the correct country. Let’s do the same for getting to the city.

In:
plot_path(encode_location(HOUSE, 16), 0.01)

16 bits (2 bytes) gets us in the county, really close to the actual city. Let’s keep going.

In:
plot_path(encode_location(HOUSE, 24), 0.008)
plot_path(encode_location(HOUSE, 24), 0.0005)

With another byte, 24-bits puts us in the correct suburb. Let’s see if we can pinpoint the house with a few more bits.

In:
path = encode_location(HOUSE, 32)
plot_path(path, 0.000019)
dist = haversine(decode_bits(path), HOUSE)
f"{dist} meters from target"

OpenStreetMap shortlink

It turns out OpenStreetMap short URLs use the same encoding. Therefore we can generate these links without using the website.

In:
def bitstring_to_bytes(s):
    return int(s, 2).to_bytes((len(s) + 7) // 8, byteorder='big')

loc = encode_location(HOUSE, 48)
by = bitstring_to_bytes(loc)

import base64
"https://openstreetmap.org/go/" + base64.b64encode(by).decode('ascii')

Encoding coordinates as English words

We can also have some fun and encode these coordinates as combinations of English words. While this might sound like a weird thing to do, you might end up making millions with a product like this.

In:
import requests

URL = "https://www.eff.org/files/2016/07/18/eff_large_wordlist.txt"
text = requests.get(URL).text.split("\n")

wordlist = []

for line in text:
    if not line:
        break
    number, word = line.split()
    wordlist.append(word)
wordlist_orig = list(wordlist)
len(wordlist)
wordlist[:5]
In:
wordlist = list(wordlist_orig)

word2bin = {}
bin2word = {}
num = 0

for nbits in range(1, 11 + 1):
    for n in range(2 ** nbits):
        word = wordlist.pop(0)
        bits = bin(n)[2:].rjust(nbits, '0')
        word2bin[word] = bits
        bin2word[bits] = word
In:
loc = encode_location(HOUSE, 32)

def to_words(bits):
    result = []
    
    group = ""
    for bit in bits:
        group += bit
        
        if len(group) == 11:
            result.append(bin2word[group])
            group = ""
    if group:
        result.append(bin2word[group])
    return '.'.join(result)

loc
to_words(loc)

Let’s turn it back.

In:
def to_coords(text):
    bits = ""
    for word in text.split('.'):
        bits += word2bin[word]
    return bits

loc
to_coords('grooving.familiar.clasp')

We can even test this with random coordinates.

In:
import random

for _ in range(25):
    x = random.uniform(-180, 180)
    y = random.uniform(90, -90)
    
    loc = encode_location((x, y), 33)
    words = to_words(loc)
    coords = to_coords(words)
    
    print(f"{words.ljust(35)} => {x:.2f},{y:.2f}\t{int(haversine((x, y), decode_bits(loc)))} m")
    assert loc == coords
In:
error = 0
n = 0

for _ in range(50_000):
    x = random.uniform(-180, 180)
    y = random.uniform(-90, 90)
    loc = encode_location((x, y), 33)
    words = to_words(loc)
    coords = to_coords(words)

    # Every coordinate can be represented by 3 words
    assert len(words.split(".")), 3

    assert loc == coords

    error += haversine((x, y), decode_bits(loc))
    n += 1
    
error /= n
f"Average error: {int(error)} meters ({n} samples)"

Leave a Comment