'''
Created on 15.03.2017

@author: luroth
'''

import os
import math
import numpy as np
from matplotlib import path
from osgeo import gdal

DELTA_x = 2693000
DELTA_y = 1256000

def rotation_matrix(axis, theta):
    axis = np.asarray(axis)
    axis = axis/math.sqrt(np.dot(axis, axis))
    a = math.cos(theta/2.0)
    b, c, d = -axis*math.sin(theta/2.0)
    aa, bb, cc, dd = a*a, b*b, c*c, d*d
    bc, ad, ac, ab, bd, cd = b*c, a*d, a*c, a*b, b*d, c*d
    return np.array([[aa+bb-cc-dd, 2*(bc+ad), 2*(bd-ac)],
                     [2*(bc-ad), aa+cc-bb-dd, 2*(cd+ab)],
                     [2*(bd+ac), 2*(cd-ab), aa+dd-bb-cc]])
    
def find_mins_maxs(mesh):
    minx = maxx = miny = maxy = minz = maxz = None
    for p in mesh.points:
        # p contains (x, y, z)
        if minx is None:
            minx = p[stl.Dimension.X]
            maxx = p[stl.Dimension.X]
            miny = p[stl.Dimension.Y]
            maxy = p[stl.Dimension.Y]
            minz = p[stl.Dimension.Z]
            maxz = p[stl.Dimension.Z]
        else:
            maxx = max(p[stl.Dimension.X], maxx)
            minx = min(p[stl.Dimension.X], minx)
            maxy = max(p[stl.Dimension.Y], maxy)
            miny = min(p[stl.Dimension.Y], miny)
            maxz = max(p[stl.Dimension.Z], maxz)
            minz = min(p[stl.Dimension.Z], minz)
    return minx, maxx, miny, maxy, minz, maxz
    
def get_polygon_normale(points):
    a = np.array(points[0])
    b = np.array(points[1])
    c = np.array(points[2])

    v1 = a-b
    v2 = a-c
    normale = np.cross(v1, v2)
    
    return normale

def getWKT_PRJ (epsg_code):
    try:
        file_dir = os.path.dirname(os.path.realpath(__file__))
    except:
        file_dir = '/home/luroth/git_repos/MSc_Lukas_Roth/'
    
    with open(os.path.join(file_dir, 'EPSG_files', 'EPSG' + str(epsg_code) + '.txt'), 'r') as file:
        remove_spaces = file.read()
        remove_spaces = remove_spaces.replace(' ', '')
        output = remove_spaces.replace("\n", "")
    return output

def point_in_polygon(x,y,poly,precision):        
    path_ = path.Path(poly)
    return path_.contains_point((x, y), radius=precision)

from shapely.geometry import Point
from shapely.geometry.polygon import Polygon

def point_in_polygon2(x,y,poly):
    point = Point(x, y)
    polygon = Polygon(poly)
    return polygon.contains(point)

def point_in_poly(x,y,poly):
    n = len(poly)
    inside = False

    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        if y > min(p1y,p2y):
            if y <= max(p1y,p2y):
                if x <= max(p1x,p2x):
                    if p1y != p2y:
                        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                    if p1x == p2x or x <= xints:
                        inside = not inside
        p1x,p1y = p2x,p2y

    return inside

def create_GeoTiff(file_path, array, data_type):
    
    driver = gdal.GetDriverByName('GTiff')
    file = file_path + '.tif'
    
    ds = driver.Create(file, array.shape[1], array.shape[0], 1, data_type)

    band_1 = ds.GetRasterBand(1)
    band_1.WriteArray(array)
    band_1.SetNoDataValue(0)
    
    ds = None
    
    return file_path