# -*- coding: utf-8 -*-
"""Various utility functions.

Copyright (c) 2020, 3D geoinformation group, Delft University of Technology

The MIT License (MIT)

Permission is hereby granted, free of charge, to any person obtaining a copy of 
this software and associated documentation files (the "Software"), to deal in 
the Software without restriction, including without limitation the rights to 
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies 
of the Software, and to permit persons to whom the Software is furnished to do 
so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all 
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 
SOFTWARE.
"""
import json
import math
from statistics import mean
from typing import Iterable, Tuple, Mapping, TextIO, Union
import logging
import zipfile, gzip
from platform import platform
from pathlib import Path

log = logging.getLogger(__name__)

def create_rectangle_grid(bbox: Iterable[float], hspacing: float,
                          vspacing: float) -> Iterable:
    """

    :param bbox: (xmin, ymin, xmax, ymax)
    :param hspacing:
    :param vspacing:
    :return: A MultiPolygon of rectangular polygons (the grid) as Simple Feature
    """
    xmin, ymin, xmax, ymax = bbox
    width = math.ceil(xmax - xmin)
    height = math.ceil(ymax - ymin)
    cols = math.ceil(width / hspacing)
    rows = math.ceil(height / vspacing)
    multipolygon = list()

    for col in range(cols):
        x1 = float(xmin) + (col * hspacing)
        x2 = x1 + hspacing

        for row in range(rows):
            y1 = float(ymax) - (row * vspacing)
            y2 = y1 - vspacing

            # A polygon with a single (outer) ring
            polygon = [[(x1, y1), (x1, y2), (x2, y2), (x2, y1), (x1, y1)]]
            multipolygon.append(polygon)

    return multipolygon

def create_rectangle_grid_morton(bbox: Iterable[float], hspacing: float,
                                 vspacing: float) -> Mapping:
    """Creates a grid of rectangular polygons and computes their Morton code.

    If the width or height of the ``bbox`` is not divisible by 4 without a
    remainder, then the extent of the grid is expanded until it is. The reason
    for doing so is that the output grid is meant to be used as the leafs of a
    quadtree.

    :param bbox: (xmin, ymin, xmax, ymax)
    :param hspacing: Width of a cell
    :param vspacing: Height of a cell
    :return: A dictionary of {morton code: Polygon}. Polygon is represented as
        Simple Feature.
    """
    xmin, ymin, xmax, ymax = bbox
    width = math.ceil(xmax - xmin)
    height = math.ceil(ymax - ymin)
    _c = math.ceil(width / hspacing)
    cols = _c + (_c % 4)
    _r = math.ceil(height / vspacing)
    rows = _r + (_r % 4)
    if rows < cols:
        rows += cols - rows
    elif cols < rows:
        cols += rows - cols
    # Expand extent until we get enough cells for a full quadtree
    exponent = math.ceil(math.log(rows*cols, 4))
    full_cells = 4**exponent
    rows = int(math.sqrt(full_cells))
    cols = int(math.sqrt(full_cells))
    grid = dict()


    for col in range(cols):
        x1 = float(xmin) + (col * hspacing)
        x2 = x1 + hspacing

        for row in range(rows):
            y1 = float(ymax) - (row * vspacing)
            y2 = y1 - vspacing

            ring = [(x1, y1), (x1, y2), (x2, y2), (x2, y1), (x1, y1)]
            # A polygon with a single (outer) ring
            polygon = [ring,]

            centroid = mean_coordinate(ring)
            morton_key = morton_code(*centroid)

            grid[morton_key] = polygon

    return dict((k, grid[k]) for k in sorted(grid))


def index_quadtree(grid):
    """Create indices for the leafs of the quadtree.

    Based on AHN's tile indexing.

    :param grid: A rectangular grid of polygons which has 4**x cells. The cells
        must be sorted in Morton-order.
    """
    quadtree = dict()
    nr_cells = len(grid)
    if not math.log(nr_cells, 4).is_integer():
        raise ValueError(f"There are {nr_cells} in the grid. The grid must "
                         f"contain 4**x cells to form a full quadtree. ")
    # Nr. levels in the quadtree
    nr_lvls = int(math.log(nr_cells, 4))
    log.debug(f"Nr. levels={nr_lvls}, cells={nr_cells}")

    id_map = {
        0 : ('1', '2', '3', '4'),
        1 : ('e', 'f', 'g', 'i'),
        2 : ('a', 'b', 'c', 'd'),
        3 : ('1', '2', '3', '4'),
        4 : ('1', '2', '3', '4')
    }

    # Extend the ID map to have as many levels as there are in the quadtree
    if nr_lvls > len(id_map):
        diff = nr_lvls - len(id_map)
        for i in range(diff):
            id_map[5+i] = id_map[i]

    # Compose the cell IDs per level
    for i, mcode in enumerate(grid):
        cell_id = ""
        for j in range(nr_lvls, 0, -1):
            lvl_id = id_map[j-1]
            lvl_idx = int((i % 4**j) / 4**(j-1))
            cell_id += lvl_id[lvl_idx]
        if cell_id in quadtree:
            raise IndexError(f"ID {cell_id} already exists in the quadtree")
        else:
            quadtree[cell_id] = mcode

    return quadtree


def bbox(polygon: Iterable) -> Tuple[float, float, float, float]:
    """Compute the Bounding Box of a polygon.

    :param polygon: A Simple Feature Polygon, defined as [[[x1, y1], ...], ...]
    """
    x,y = 0,1
    vtx = polygon[0][0]
    minx, miny, maxx, maxy = vtx[x], vtx[y], vtx[x], vtx[y]
    for ring in polygon:
        for vtx in ring:
            if vtx[x] < minx:
                minx = vtx[x]
            elif vtx[y] < miny:
                miny = vtx[y]
            elif vtx[x] > maxx:
                maxx = vtx[x]
            elif vtx[y] > maxy:
                maxy = vtx[y]
    return minx, miny, maxx, maxy


def distance(a,b) -> float:
    """Distance between point a and point b"""
    x,y = 0,1
    return math.sqrt((a[x] - b[x])**2 + (a[y] - b[y])**2)


def is_between(a,c,b) -> bool:
    """Return True if point c is on the segment ab

    Ref.: https://stackoverflow.com/a/328193
    """
    return math.isclose(distance(a,c) + distance(c,b), distance(a,b))


def in_bbox(point: Tuple[float, float], bbox: Tuple) -> bool:
    """Evaluates if a point is in the provided bounding box.

    A poin is in the BBOX if it is either completely within
    the BBOX, or overlaps with the South (lower) or West (left) boundaries
    of the BBOX.

    :param point: A point defined as a tuple of cooridnates of (x,y)
    :param bbox: Bounding Box as (minx, miny, maxx, maxy)
    """
    if not bbox or not point:
        return False
    x,y = 0,1
    minx, miny, maxx, maxy = bbox
    within = ((minx < point[x] < maxx) and
              (miny < point[y] < maxy))
    on_south_bdry = is_between((minx, miny), point, (maxx, miny))
    on_west_bdry = is_between((minx, miny), point, (minx, maxy))
    return any((within, on_south_bdry, on_west_bdry))


def mean_coordinate(points: Iterable[Tuple]) -> Tuple[float, float]:
    """Compute the mean x- and y-coordinate from a list of points.

    :param points: An iterable of coordinate tuples where the first two elements
        of the tuple are the x- and y-coordinate respectively.
    :returns: A tuple of (mean x, mean y) coordinates
    """
    mean_x = mean(pt[0] for pt in points)
    mean_y = mean(pt[1] for pt in points)
    return mean_x, mean_y

# Computing Morton-code. Reference: https://github.com/trevorprater/pymorton ---

def __part1by1_64(n):
    """64-bit mask"""
    n &= 0x00000000ffffffff                  # binary: 11111111111111111111111111111111,                                len: 32
    n = (n | (n << 16)) & 0x0000FFFF0000FFFF # binary: 1111111111111111000000001111111111111111,                        len: 40
    n = (n | (n << 8))  & 0x00FF00FF00FF00FF # binary: 11111111000000001111111100000000111111110000000011111111,        len: 56
    n = (n | (n << 4))  & 0x0F0F0F0F0F0F0F0F # binary: 111100001111000011110000111100001111000011110000111100001111,    len: 60
    n = (n | (n << 2))  & 0x3333333333333333 # binary: 11001100110011001100110011001100110011001100110011001100110011,  len: 62
    n = (n | (n << 1))  & 0x5555555555555555 # binary: 101010101010101010101010101010101010101010101010101010101010101, len: 63

    return n


def __unpart1by1_64(n):
    n &= 0x5555555555555555                  # binary: 101010101010101010101010101010101010101010101010101010101010101, len: 63
    n = (n ^ (n >> 1))  & 0x3333333333333333 # binary: 11001100110011001100110011001100110011001100110011001100110011,  len: 62
    n = (n ^ (n >> 2))  & 0x0f0f0f0f0f0f0f0f # binary: 111100001111000011110000111100001111000011110000111100001111,    len: 60
    n = (n ^ (n >> 4))  & 0x00ff00ff00ff00ff # binary: 11111111000000001111111100000000111111110000000011111111,        len: 56
    n = (n ^ (n >> 8))  & 0x0000ffff0000ffff # binary: 1111111111111111000000001111111111111111,                        len: 40
    n = (n ^ (n >> 16)) & 0x00000000ffffffff # binary: 11111111111111111111111111111111,                                len: 32
    return n


def interleave(*args):
    """Interleave two integers to create a Morton key."""
    if len(args) != 2:
        raise ValueError('Usage: interleave2(x, y)')
    for arg in args:
        if not isinstance(arg, int):
            print('Usage: interleave2(x, y)')
            raise ValueError("Supplied arguments contain a non-integer!")

    return __part1by1_64(args[0]) | (__part1by1_64(args[1]) << 1)


def deinterleave(morton_key):
    """Deinterleave a Morton key to get the original coordinates."""
    if not isinstance(morton_key, int):
        print('Usage: deinterleave2(n)')
        raise ValueError("Supplied arguments contain a non-integer!")

    return __unpart1by1_64(morton_key), __unpart1by1_64(morton_key >> 1)


def morton_code(x: float, y: float):
    """Takes an (x,y) coordinate tuple and computes their Morton-key.

    Casts float to integers by multiplying them with 100 (millimeter precision).
    """
    return interleave(int(x * 100), int(y * 100))


def rev_morton_code(morton_key: int) -> Tuple[float, float]:
    """Get the coordinates from a Morton-key"""
    x,y = deinterleave(morton_key)
    return float(x)/100.0, float(y)/100.0


def read_geojson_polygon(fo: TextIO) -> Iterable:
    """Reads a single polygon from a GeoJSON file.
    :returns: A Simple Feature representation of the polygon
    """
    polygon = list()
    # Only Polygon is allowed (no Multi-)
    gjson = json.load(fo)
    if gjson['features'][0]['geometry']['type'] != 'Polygon':
        raise ValueError(f"The first Feature in GeoJSON is "
                         f"{gjson['features'][0]['geometry']['type']}. Only "
                         f"Polygon is allowed.")
    else:
        polygon = gjson['features'][0]['geometry']['coordinates']
    return polygon


def polygon_to_ewkt(polygon, srid) -> str:
    """Creates a EWKT representation of a Simple Feature polygon.
    :returns: The EWKT string of ``polygon``
    """
    ring = [" ".join(map(str, i)) for i in polygon[0]]
    ewkt = f'SRID={srid};POLYGON(({",".join(ring)}))'
    return ewkt


def polyline_to_ewkt(polyline, srid) -> str:
    """Creates a EWKT representation of a Simple Feature polyline.
    :returns: The EWKT string of ``polyline``
    """
    linestring = [" ".join(map(str, i)) for i in polyline]
    ewkt = f'SRID={srid};LINESTRING({",".join(linestring)})'
    return ewkt


def rectangle_sw_boundary(rectangle):
    """Extracts the South-West edges of a rectangle polygon into a polyline."""
    # rectangle[0] is the outer ring of the polygon
    bbox = [*rectangle[0][0], *rectangle[0][1]]
    for v in rectangle[0]:
        if v[0] < bbox[0]:
            bbox[0] = v[0]
        elif v[0] > bbox[2]:
            bbox[2] = v[0]
        if v[1] < bbox[1]:
            bbox[1] = v[1]
        elif v[1] > bbox[3]:
            bbox[3] = v[1]
    # polyline in CCW of the polygon boundary
    return [(bbox[2], bbox[1]), (bbox[0], bbox[1]), (bbox[0], bbox[3])]

def lod_to_string(lod: Union[int, float]) -> Union[str, None]:
    """Convert and LoD integer or float to string.
    """
    if lod is None:
        return None
    elif isinstance(lod, str):
        return lod
    elif isinstance(lod, int):
        return str(lod)
    elif isinstance(lod, float):
        return str(round(lod, 1))
    else:
        raise ValueError(f"Type {type(lod)} is not allowed as input")

def parse_lod_value(lod_key: str) -> str:
    """Extract the LoD value from an LoD parameter key (eg. lod13).

    For example 'lod13' -> '1.3'
    """
    pos = lod_key.lower().find('lod')
    if pos != 0:
        raise ValueError(f"The key {lod_key} does not begin with 'lod'")
    value = lod_key[3:]
    if len(value) == 1:
        return value
    elif len(value) == 2:
        return f"{value[0]}.{value[1]}"
    else:
        raise ValueError(f"Invalid LoD value '{value}' in key {lod_key}")


def write_zip(data: bytes, filename: str, outdir: Path):
    """Write out a citymodel to a zip file.

    On Linux and MacOS it uses Gzip, on Windows it uses Zip.

    :param data: Data to compress into a file
    :param filename: Filename to write
    :param outdir: Output directory
    """
    outfile = outdir / filename
    if "windows" in platform().lower():
        outzip = outfile.with_suffix(".zip")
        with zipfile.ZipFile(file=outzip, mode="w") as zout:
            zout.writestr(zinfo_or_arcname=filename,
                          data=data)
    else:
        outzip = outfile.with_suffix(".json.gz")
        with gzip.open(outzip, "w") as zout:
            zout.write(data)
    return outzip