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.
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.
WORLD = Quadtile(-180, 180, 90, -90) WORLD
Quadtile(xFrom=-180, xTo=180, yFrom=90, yTo=-90)
Operations on Quads
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.
def quad_center(quad): x = (quad.xFrom + quad.xTo) / 2 y = (quad.yFrom + quad.yTo) / 2 return (x, y)
quad_center(WORLD)
(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
def quad_split_x(quad, right): mid = (quad.xFrom + quad.xTo) / 2 if right: xf = mid xt = quad.xTo else: xf = quad.xFrom xt = mid return Quadtile(xf, xt, quad.yFrom, quad.yTo) quad_center(quad_split_x(WORLD, False)) quad_center(quad_split_x(WORLD, True))
(-90.0, 0.0)
(90.0, 0.0)
Now let’s do it in the other axis.
def quad_split_y(quad, bottom): mid = (quad.yFrom + quad.yTo) / 2 if not bottom: yf = mid yt = quad.yTo else: yf = quad.yFrom yt = mid return Quadtile(quad.xFrom, quad.xTo, yf, yt) quad_center(quad_split_y(WORLD, False)) quad_center(quad_split_y(WORLD, True))
(0.0, -45.0)
(0.0, 45.0)
A simple encoding
blah blah blah
def decode_bits(bits): is_x = True quad = WORLD for bit in bits: bit = int(bit) if is_x: quad = quad_split_x(quad, bit) else: quad = quad_split_y(quad, bit) is_x = not is_x return quad_center(quad) decode_bits("01001") decode_bits("01100111")
(-112.5, 22.5)
(-56.25, 39.375)
Blah blah blah
def encode_location(coord, N): bits = "" quad = WORLD x, y = coord is_x = True for i in range(N): mid = quad_center(quad) if is_x: if x > mid[0]: bit = 1 else: bit = 0 quad = quad_split_x(quad, bit) else: if y > mid[1]: bit = 1 else: bit = 0 quad = quad_split_y(quad, bit) 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
HOUSE = (-8.577507, 52.664838)
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()
plot_path(encode_location(HOUSE, 8), 0.55)
Path length: 8
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.
plot_path(encode_location(HOUSE, 16), 0.01)
Path length: 16
16 bits (2 bytes) gets us in the county, really close to the actual city. Let’s keep going.
plot_path(encode_location(HOUSE, 24), 0.008) plot_path(encode_location(HOUSE, 24), 0.0005)
Path length: 24
Path length: 24
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.
path = encode_location(HOUSE, 32) plot_path(path, 0.000019) dist = haversine(decode_bits(path), HOUSE) f"{dist} meters from target"
Path length: 32
'55.883351336259565 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.
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')
'https://openstreetmap.org/go/esb8PMRe'
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.
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]
7776
['abacus', 'abdomen', 'abdominal', 'abide', 'abiding']
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
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)
'01111010110001101111110000111100'
'grooving.familiar.clasp'
Let’s turn it back.
def to_coords(text): bits = "" for word in text.split('.'): bits += word2bin[word] return bits loc to_coords('grooving.familiar.clasp')
'01111010110001101111110000111100'
'01111010110001101111110000111100'
We can even test this with random coordinates.
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
flinch.errand.hula => -123.11,18.66 136 m freebie.energetic.grape => -144.44,62.74 141 m mobster.napkin.enable => 168.58,11.03 183 m kite.entomb.error => 159.67,-4.47 68 m jargon.monotype.jellied => 94.01,-39.69 66 m majority.favored.fled => 63.31,70.27 103 m name.mandate.headboard => 173.76,77.19 75 m matrimony.laborious.maker => 132.66,42.53 62 m handcuff.luxury.gestation => 21.32,-63.75 124 m hardiness.fancied.jurist => 69.06,-76.43 26 m fraying.flannels.landmass => -153.37,56.57 148 m lapdog.luridness.mystified => 22.49,37.08 160 m hut.duty.lanky => 56.52,-5.30 86 m gravy.factual.ferris => -48.44,69.56 60 m junkyard.flattop.edgy => 143.47,-32.95 159 m moonbeam.kerchief.liking => 97.86,55.61 122 m graves.hush.magazine => -59.12,71.17 131 m gainfully.effort.garbage => -112.29,63.18 94 m legwork.esquire.fade => 85.49,7.16 152 m foil.embody.flame => -121.81,22.74 146 m facelift.mooned.everyday => -85.71,-1.22 146 m landing.eggshell.galley => 39.93,12.47 136 m encrypt.handful.ground => -150.91,-8.10 123 m gulf.musky.narrow => -17.58,78.66 83 m exert.mounted.erasable => -33.86,-46.41 136 m
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)"
'Average error: 99 meters (50000 samples)'