#,!/usr/bin/env python

import math
from functools import partial
from operator import itemgetter


import fiona
import pyproj
from ujson import dumps  # noqa
from shapely.ops import transform
from shapely.geometry import mapping, shape  # noqa

import shapely.speedups
shapely.speedups.available and shapely.speedups.enable()


def meters_per_pixel(zoom, lat):
    """
    ground resolution = cos(latitude * pi/180) * earth circumference / map width
    """

    return (math.cos(lat * math.pi/180.0) * 2 * math.pi * 6378137) / (256 * 2**zoom)


def map_scale(zoom, lat, dpi=96.0):
    """
    map scale = 1 : ground resolution * screen dpi / 0.0254 meters/inch
    """

    res = meters_per_pixel(zoom, lat)
    return (res * dpi) / 0.0254


def zoom_to_scaledenom(z1, z2=False):
    """
    Converts zoom level to mapnik's scaledenominator pair for EPSG:3857
    """
    if not z2:
        z2 = z1
    s = 279541132.014
    z1 = (s / (2 ** (z1 - 1)) + s / (2 ** (z1 - 2))) / 2
    z2 = (s / (2 ** (z2 - 1)) + s / (2 ** z2)) / 2
    # return 100000000000000, 1
    return z1, z2


def pixel_size_at_zoom(z, l=1):
    """
    Converts l pixels on tiles into length on zoom z
    """
    return int(math.ceil(l * 20037508.342789244 / 256 * 2 / (2 ** z)))


zoom = 3
lat = -28.226970038918328
mpp = meters_per_pixel(zoom, lat)
scale = map_scale(zoom, lat)

ds = fiona.open("data/regions/LGA", "r")
fs = list(ds)
gs = map(shape, filter(None, map(itemgetter("geometry"), fs)))


epsg_4283_to_epsg_3857 = partial(
    pyproj.transform,
    pyproj.Proj(init="epsg:4283"),
    pyproj.Proj(init="epsg:3857")
)

epsg_3857_to_epsg_4326 = partial(
    pyproj.transform,
    pyproj.Proj(init="epsg:3857"),
    pyproj.Proj(init="epsg:4326")
)

g1 = gs[510]  # largest geometry in LGA regions ~165KB GeoJSON
g2 = transform(epsg_4283_to_epsg_3857, g1)
#g3 = g2.simplify(pixel_size_at_zoom(zoom, 0.5))
g3 = g2.buffer(pixel_size_at_zoom(zoom, 10)).simplify(pixel_size_at_zoom(zoom, 3)).buffer(pixel_size_at_zoom(zoom, -10)
g4 = transform(epsg_3857_to_epsg_4326, g3)

print "size(g1):", len(dumps(mapping(g1)))
print "size(g4):", len(dumps(mapping(g4)))

#Map resolution = 156543.04 meters/pixel * cos(latitude) / (2 ^ zoomlevel)
#Komzzzpa> well, algo is:
#          1) convert to EPSG:3857
#          2) make simplification for 0.5px with .simplify(pixel_size_at_zoom(0.5))
#          3) convert to EPSG:4326
#          4) serialize/show
#<Komzzzpa> for polygons, buffer(10px).simplify(3px).buffer(-10px) might be a lot better
#<Komzzzpa> for adjacent polygons, whole other ideas are needed :3