Three tiles systems : Squares, Triangles, Hexagons ; with a hash system to encode them

tiles demo.ipynb

tiles.py

import numpy as np
import pandas as pd

def mercator(lat):
    """Map latitudes (in deg) to Mercator projection.
    (Between -1 and 1 for latitudes between -85.05 and +85.05).
    """
    return np.arcsinh(np.tan(lat * np.pi / 180.)) / np.pi

def mercator_inv(x):
    """Map Mercator "longitudes" to DEG-longitudes."""
    return np.arctan(np.sinh(x * np.pi)) * 180 / np.pi

porte_de_dijon = (45.8381, 3.1205)  # A GPS position to make tests

def to_hash(x, y, z, prefix=1):
    """Hash the values of x and y in an single int.
    If x = 0bXXXXXX and y = 0bYYYYYY,
    then the output's binary string will be : bin(prefix) + "XYXYXYXYXY".
    """
    assert prefix in [1, 2, 3]
    return (
        int('0'.join(bin(x)[2:])+'0', 2) + int('0'.join(bin(y)[2:]), 2) + prefix * 2**(2*z)
    )

def from_hash(n):
    if n in (1,2,3):
        return 0,0,0,n
    s = bin(n)[2:]
    if len(s)%2:
        prefix, suffix = int(s[:1], 2), s[1:]
    else:
        prefix, suffix = int(s[:2], 2), s[2:]
    return int(suffix[::2], 2), int(suffix[1::2], 2), len(suffix)//2, prefix


def tile_from_hash(n):
    x, y, z, prefix = from_hash(n)
    if prefix == 1:
        return SquareTile(x, y, z)
    if prefix == 2:
        return TriangleTile(x, y, z)
    if prefix == 3:
        return HexagonalTile(x, y, z)


def test_hash_function():
    assert bin(to_hash(0,0,5)) == '0b10000000000'
    assert bin(to_hash(31,0,5)) == '0b11010101010'
    assert bin(to_hash(0,31,5)) == '0b10101010101'
    assert bin(to_hash(31,31,5)) == '0b11111111111'

    assert bin(to_hash(0,0,5,2)) == '0b100000000000'
    assert bin(to_hash(31,0,5,2)) == '0b101010101010'
    assert bin(to_hash(0,31,5,2)) == '0b100101010101'
    assert bin(to_hash(31,31,5,2)) == '0b101111111111'

    for n in range(4, 10000):
        assert to_hash(*from_hash(n)) == n, n

    i = 0
    for z in range(8):
        for prefix in range(1,4):
            for x in range(2**z):
                for y in range(2**z):
                    i += 1
                    assert from_hash(to_hash(x, y, z, prefix)) == (x, y, z, prefix), (x, y, z, prefix)

class Tile(object):
    def __init__(self, x, y, z):
        self.x = int(x)
        self.y = int(y)
        self.z = int(z)

    def to_string(self, include_zoom=True):
        if include_zoom:
            return '{}/{}/{}'.format(self.z, self.x, self.y)
        else:
            return '{}/{}'.format(self.x, self.y)

    def to_tuple(self, include_zoom=True):
        if include_zoom:
            return (self.z, self.x, self.y)
        else:
            return (self.x, self.y)

    def vertex(self):
        raise NotImplementedError()

    def to_geojson(self):
        return { "type": "Polygon", "coordinates": [[x[::-1] for x in self.vertex()]]}

class SquareTile(Tile):
    def from_lat_lng(lat, lng, zoom):
        x, y = (lng + 180) / 360, (1 - mercator(lat)) / 2  # x=0-1, y=1-0
        return SquareTile(int(x * 2**zoom), int(y * 2**zoom), zoom)

    def central_point(self):
        x = (self.x + 0.5) / 2**self.z
        y = (self.y + 0.5) / 2**self.z
        return mercator_inv(1 - y * 2), x * 360 - 180

    def vertex(self):
        return [
            [mercator_inv(1 - (self.y + j) / 2**self.z * 2), (self.x + i) / 2**self.z * 360 - 180]
            for i, j in zip([0, 1, 1, 0], [0, 0, 1, 1])]

    def to_hash(self):
        prefix = 1
        return to_hash(self.x, self.y, self.z, prefix=prefix)

def test_square_tile():
    z = 18
    s = SquareTile.from_lat_lng(*porte_de_dijon, z)
    sc = s.central_point()
    assert SquareTile.from_lat_lng(*sc, z).central_point() == sc


class TriangleTile(Tile):
    def from_lat_lng(lat, lng, zoom):
        x, y = (lng + 180) / 360, (0.5 - mercator(lat)*0.5/np.sqrt(3))
        dx, mx = divmod(x * 2**zoom, 1)
        dy, my = divmod(y * 2**zoom, 1)
        if (dx + dy) % 2 == 0:
            return TriangleTile(int(dx + (mx>my)), int(dy), zoom)
        else:
            return TriangleTile(int(dx + (mx>1-my)), int(dy), zoom)

    def _transform(self, lat, lng):
        """transforms latitude and longitude into variables x and y,
        so that a tile will have width and heaight equal to 1."""
        return (
            (lng + 180) / 360 * 2**self.z,
            (1 - mercator(lat)/np.sqrt(3)) / 2 * 2**self.z
            )
            
    def _transform_inv(self, x, y):
        """transforms latitude and longitude into variables x and y,
        so that a tile will have width and heaight equal to 1."""
        return (
            mercator_inv((1 - 2 * y / 2**self.z) * np.sqrt(3)),
            x / 2**self.z * 360 - 180
            )

    def central_point(self):
        x = self.x
        if int(self.x + self.y) % 2 == 0:  # triangle pointing up
            y = (self.y + 2/3)
        else:
            y = (self.y + 1/3)
        return self._transform_inv(x, y)

    def vertex(self):
        if int(self.x + self.y) % 2 == 0:  # triangle pointing up
            dys = [1, 0, 1]
        else:
            dys = [0, 1, 0]
        return [
            self._transform_inv(self.x + dx, self.y + dy)
            for dx, dy in zip([-1, 0, 1], dys)]

    def to_hash(self):
        prefix = 2
        return to_hash(self.x, self.y, self.z, prefix=prefix)


def test_triangle_tile():
    z = 18
    t = TriangleTile.from_lat_lng(*porte_de_dijon, z)
    tc = t.central_point()
    assert TriangleTile.from_lat_lng(*tc, z).central_point() == tc
    assert t._transform(*t._transform_inv(t.x, t.y)) == (t.x, t.y)
    assert (
        np.round(t._transform_inv(*t._transform(*porte_de_dijon)), 10) == np.round(porte_de_dijon, 10)
        ).all()


class HexagonalTile(Tile):
    def from_lat_lng(lat, lng, zoom):
        x, y = (lng + 180) / 360, (0.5 - mercator(lat)*0.5/np.sqrt(3))
        dx, mx = divmod(x * 2**zoom, 1)
        dy, my = divmod(y * 2**zoom, 1)
        if (dx + dy) % 2 == 0:
            if mx < 2 - 3 * my:
                return HexagonalTile(int(dx), int(dy), zoom)
            else:
                return HexagonalTile(int(dx+1), int(dy+1), zoom)
        else:
            if 1 - mx < 2 - 3 * my:
                return HexagonalTile(int(dx+1), int(dy), zoom)
            else:
                return HexagonalTile(int(dx), int(dy+1), zoom)

    def _transform(self, lat, lng):
        """transforms latitude and longitude into variables x and y,
        so that a tile will have width and heaight equal to 1."""
        return (
            (lng + 180) / 360 * 2**self.z,
            (1 - mercator(lat)/np.sqrt(3)) / 2 * 2**self.z
            )
            
    def _transform_inv(self, x, y):
        """transforms latitude and longitude into variables x and y,
        so that a tile will have width and heaight equal to 1."""
        return (
            mercator_inv((1 - 2 * y / 2**self.z) * np.sqrt(3)),
            x / 2**self.z * 360 - 180
            )

    def central_point(self):
        x = self.x
        y = self.y
        #if int(self.x + self.y) % 2 == 0:  # triangle pointing up
        #    y = (self.y + 2/3)
        #else:
        #    y = (self.y + 1/3)
        return self._transform_inv(x, y)

    def vertex(self):
        dxs = np.array([0, -1, -1, 0, 1, 1])
        dys = np.array([-2, -1, 1, 2, 1, -1]) / 3
        #if int(self.x + self.y) % 2 == 0:  # triangle pointing up
        #    dys = [1, 0, 1]
        #else:
        #    dys = [0, 1, 0]
        return [
            self._transform_inv(self.x + dx, self.y + dy)
            for dx, dy in zip(dxs, dys)]

    def to_hash(self):
        prefix = 3
        return to_hash(self.x, self.y, self.z, prefix=prefix)


def test_hexagon_tile():
    z = 18
    t = HexagonalTile.from_lat_lng(*porte_de_dijon, z)
    tc = t.central_point()
    assert HexagonalTile.from_lat_lng(*tc, z).central_point() == tc
    assert t._transform(*t._transform_inv(t.x, t.y)) == (t.x, t.y)
    assert (
        np.round(t._transform_inv(*t._transform(*porte_de_dijon)), 10) == np.round(porte_de_dijon, 10)
        ).all()