Quadkey UDFs

redshift_quadkey_udf.sql

create or replace function latlng_to_quadkey (lat float, lng float, zoom int)
  returns varchar
STABLE
as $$

    import math
    if lng is None or lat is None:
        return None
    lat = float(lat)
    lng = float(lng)
    zoom = int(zoom)
    lat = math.radians(lat)
    n = 2.0 ** zoom
    xtile = (lng + 180.0) / 360.0 * n

    try:
        ytile = (
            (1.0 - math.log(math.tan(lat) + (1.0 / math.cos(lat))) / math.pi) / 2.0 * n
        )
    except:
        return None

    xtile = int(math.floor(xtile))
    ytile = int(math.floor(ytile))

    qk = []
    for z in range(zoom, 0, -1):
        digit = 0
        mask = 1 << (z - 1)
        if xtile & mask:
            digit += 1
        if ytile & mask:
            digit += 2
        qk.append(str(digit))

    return "".join(qk)

$$ language plpythonu;




create or replace function quadkey_to_latlng (qk varchar)
  returns varchar
STABLE
as $$

    import math
    if qk is None:
        return None

    if len(qk) == 0:
        return 0.0
    xtile, ytile = 0, 0
    for i, digit in enumerate(reversed(qk)):
        mask = 1 << i
        if digit == "1":
            xtile = xtile | mask
        elif digit == "2":
            ytile = ytile | mask
        elif digit == "3":
            xtile = xtile | mask
            ytile = ytile | mask
        elif digit != "0":
            None
    zoom = i+1

    n = 2.0 ** zoom
    ul_lon_deg = xtile / n * 360.0 - 180.0
    ul_lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
    ul_lat_deg = math.degrees(ul_lat_rad)

    xtile +=1
    ytile +=1
    br_lon_deg = xtile / n * 360.0 - 180.0
    br_lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
    br_lat_deg = math.degrees(br_lat_rad)

    lat_deg = (br_lat_deg + ul_lat_deg) /2
    lon_deg = (ul_lon_deg + br_lon_deg) /2

    return ",".join([str(lat_deg), str(lon_deg)])

$$ language plpythonu;