Skip to content

Instantly share code, notes, and snippets.

@BibMartin
Created May 15, 2020 17:56
Show Gist options
  • Save BibMartin/e7da03789f175efe43cd484fc789591a to your computer and use it in GitHub Desktop.
Save BibMartin/e7da03789f175efe43cd484fc789591a to your computer and use it in GitHub Desktop.
Three tiles systems : Squares, Triangles, Hexagons ; with a hash system to encode them
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment