from qgis.core import QGis , QgsRaster, QgsMapLayerRegistry, QgsFields, QgsField, QgsGeometry, QgsVectorFileWriter, QgsCoordinateReferenceSystem, QgsFeature, QgsCoordinateTransform, QgsVectorLayer, QgsRasterLayer, QgsDistanceArea 
import qgis.analysis
import qgis.utils
from PyQt4.QtCore import QVariant
import os
import processing

#errors

QgsGeometryaddPart_errors = {0: "success", 
1: "not a multipolygon", 
2: "ring is not a valid geometry", 
3: "new polygon ring not disjoint with existing polygons of the feature"}

#paths

def getPath(source, name):
    source, suffix = os.path.splitext(source) # get suffix and path to file
    return source + name

# vectors

def createEmptyVectorLayer(path, formatDriver, fields=None, srcCrs = QgsCoordinateReferenceSystem(4326, QgsCoordinateReferenceSystem.EpsgCrsId), shape=QGis.WKBMultiPolygon, encoding="UTF-8"):
    """
    Create empty vector layer.
    """
    drivers = QgsVectorFileWriter.ogrDriverList().values()
    if formatDriver not in drivers:
        print "Not applicable formatDriver" + str(drivers)
        return False
    writer = QgsVectorFileWriter(path, encoding, fields, shape, srcCrs, formatDriver)
    if writer.hasError() != QgsVectorFileWriter.NoError:
        print "Error when creating shapefile: ", writer.hasError()
        del writer
        return False
    return writer

def filterRingsToHoles(to, rings, minArea, maxArea, holesfw, *holesFeatureAttributes):
    codes = []
    for ring in rings:
        ringFeature = QgsFeature()
        ring_polygon = QgsGeometry.fromPolygon([ring])
        area = ring_polygon.area()
        if area <= maxArea and area >= minArea: # add to holes
            ringFeature.setGeometry(ring_polygon)
            newfieldsvalues = [value for value in holesFeatureAttributes] # add tuple of attributes from parameter to array
            ringFeature.setAttributes(newfieldsvalues)
            holesfw.addFeature(ringFeature)
        else: # not in treshold interval, write ring to polygon (valid holes)
            code = to.addRing(ring)
            codes.append(code)
    return to, codes
# layers

def addVectorLayer(path, name, dataProvider="ogr"):
    vec = QgsVectorLayer(path, name, dataProvider)
    if vec.isValid():
        QgsMapLayerRegistry.instance().addMapLayer(vec)
        return vec
    else:
        return False

def layersIterator():
    layers = QgsMapLayerRegistry.instance().mapLayers()
    for layer in layers:
        yield layer, layers

def activeLayersIterator():
    canvas = qgis.utils.iface.mapCanvas()
    layers = canvas.layers()
    for layer in layers:
        yield layer.name(), layer

def transformGeometry(geometry, srcCrs, destCrs=QgsCoordinateReferenceSystem(4326, QgsCoordinateReferenceSystem.EpsgCrsId)):
    """
    @param geometry QgsGeometry
    """
    xform = QgsCoordinateTransform(srcCrs, destCrs)
    transformed = geometry.transform(xform)
    return transformed

def calculateArea(geom, layerCrs):
    """
    Calculate area in square km for Geographic and Projected systems in M2.
    Dont recalculate units into meters (if feet)
    """
    area = 0
    if layerCrs.geographicFlag(): # if geographic system, need to calculate in meters
        calculator = QgsDistanceArea()
        calculator.setEllipsoid(layerCrs.ellipsoidAcronym())
        calculator.setEllipsoidalMode(True)
        calculator.computeAreaInit()
        if geom.isMultipart() is False: # if only simple polygon, calculate only for this
            polyg = geom.asPolygon() # transform to list of points
            if len(polyg) > 0:
                area = calculator.measurePolygon(polyg[0])
        else: # is multipart
            multi = geom.asMultiPolygon()
            for polyg in multi:
                area = area + calculator.measurePolygon(polyg[0])
    else:
        area = geom.area()
    return round(area/1e6, 2)
    
def getFirstOccurence():
    raster = None
    vector = None
    for layername, layer in activeLayersIterator():
        if isinstance(layer, QgsRasterLayer) and raster is None:
            raster = layer
        elif isinstance(layer, QgsVectorLayer) and vector is None:
            vector = layer
        if vector is not None and raster is not None:
            break
    return vector, raster

# PUBLIC
def addLayerFromLayerExtent(path, name="BBoxes", formatDriver="ESRI Shapefile"):
    """
    Create layer with geometry of active layers extents.
    @param path path of new .shp
    """
    fieldsarr = [QgsField("name", QVariant.String), QgsField("xmax" ,QVariant.Double), QgsField("xmin" ,QVariant.Double), QgsField("ymax" ,QVariant.Double), QgsField("ymin" ,QVariant.Double), QgsField("area", QVariant.Double, prec=2)] 
    fields = QgsFields()
    for fild in fieldsarr:
        fields.append(fild)
    newvector = createEmptyVectorLayer(path, formatDriver, fields)
    if newvector:
        for layername, layer in activeLayersIterator():
            refcrs = layer.crs()
            extent = layer.extent()
            extentpolyg = QgsGeometry.fromWkt(extent.asWktPolygon())
            area = calculateArea(extentpolyg, refcrs)
            print layername + " - " + str(area)
            fet = QgsFeature()
            transformed = transformGeometry(extentpolyg, refcrs)
            if transformed == 0:
                fet.setGeometry(extentpolyg)
            else:
                print "its not possible to transform extent of layer to EPSG:4326"
                continue
            fet.setAttributes([layername ,extent.xMaximum(), extent.xMinimum(), extent.yMaximum(), extent.yMinimum(), area])
            newvector.addFeature(fet)
        added = addVectorLayer(path, name)
        if added:
            "Vector Layer added."
        else:
            "Its not possible to add vector layer to map canvas"
    else:
        print "IS not possible to creaate new vector (exists?)"
    del newvector
    
#    maptools.addLayerFromLayerExtent("C:\\marcel\\KML\\surveysbboxes.kml")

def transformActiveVectorLayers(destCrs, name="TRANSFORMED", dataProvider="ESRI Shapefile"):
    """
    copy all active layers and add spcifed name to path of layer. Not for memory layers.
    @param destCrs {String} WKT, PROJ4, EPSG: ... 
    """
    for layername, layer in activeLayersIterator():
        provider= layer.dataProvider()
        if layer.type() < 1 and provider.name() != "memory" and layer.hasGeometryType(): #if vector
            source, suffix = os.path.splitext(layer.source()) # get suffix and path to file
            path = source + name
            print path
            destCrsS = QgsCoordinateReferenceSystem(destCrs)
            newvector = createEmptyVectorLayer(path, dataProvider, layer.pendingFields(), destCrsS, layer.wkbType(), layer.dataProvider().encoding())
            if newvector:
                srcCrs = layer.crs()
                for feature in layer.getFeatures(): # iter over features in active layer
                    fet = QgsFeature(feature)
                    transform = transformGeometry(fet.geometry(), srcCrs, destCrsS)
                    if transform > 0:
                        print "not possible to transform geometry of feature" + feature.id() + "in " + layername
                        continue
                    newvector.addFeature(fet)
#                added = addVectorLayer(path, name + layername)
#                if added:
                print "Vector layer " + layername + " transformed"
#                else:
#                    print "Vector " + path + "not found."
            else:
                print "IS not possible to creaate new vector (exists?)"
            del newvector
        else:
            print layername + "is not valid Vector"
            
def drapeAndAddZAttribute():
    invector, inraster = getFirstOccurence()
    if invector is None or inraster is None:
        raise Exception("In map canvas arent vector or raster")
    processing.runalg("grass:v.drape", invector, inraster, 0, None, None, -9999999, None, None, None, 1, "C:\\marcel\\draped.shp")

def extractZFromRaster(name="EXTRACTED", attributeFieldName="TWT_SEIS", fieldType = QVariant.Double, band=1, precision=6, invert=True, dataProvider="ESRI Shapefile"):
    invector, inraster = getFirstOccurence()
    path = getPath(invector.source(), name)
    newfields = QgsFields(invector.pendingFields())
    newfields.append(QgsField(attributeFieldName, fieldType))
    if invector.type() < 1 and invector.dataProvider().name() != "memory" and invector.hasGeometryType(): #if vector
         layerWkbType = invector.wkbType()
         if layerWkbType == QGis.WKBPoint or layerWkbType == QGis.WKBPoint25D:
            newvector = createEmptyVectorLayer(path, dataProvider, newfields, invector.crs(), invector.wkbType(), invector.dataProvider().encoding())
            if newvector:
                features = invector.getFeatures()
                for feature in features:
                    oldfeatureattr = feature.attributes()
                    newfeature = QgsFeature(feature)
                    oldgeom = feature.geometry().asPoint()
#                newfeature.setGeometry(QgsGeometry.fromPoint(oldgeom))
                    newfeature.setFields(newfields)
                    ident = inraster.dataProvider().identify(oldgeom, QgsRaster.IdentifyFormatValue)
                    result = ident.results()
                    value = result[band]
                    if value is not None:
                        if fieldType == QVariant.Double:
                            value = round(value, precision)
                        if invert:
                            value = value * -1
                        oldfeatureattr.append(value)
                    newfeature.setAttributes(oldfeatureattr)
                    newvector.addFeature(newfeature)
            else:
                print "IS not possible to creaate new vector "+ path + "(exists?)"
         else:
                print invector.name() + " is not point Vector"
    else:
        print invector.name() + "is not valid vector"
                    
                

def fillHoles(minArea, maxArea,output = "FILLEDHOLES" , holes="HOLES", dataProvider="ESRI Shapefile"):
    fieldsarr = [QgsField("polygonId", QVariant.Int)] 
    fieldsholes = QgsFields()
    for fild in fieldsarr:
        fieldsholes.append(fild)
    for layername, layer in activeLayersIterator():
        provider = layer.dataProvider()
        if layer.type() < 1 and provider.name() != "memory" and layer.hasGeometryType(): #if vector
            layerWkbType = layer.wkbType()
            if layerWkbType == QGis.WKBPolygon or layerWkbType == QGis.WKBMultiPolygon or layerWkbType ==  QGis.WKBPolygon25D or layerWkbType == QGis.WKBMultiPolygon25D:
                destCrs = layer.crs()
                features = layer.getFeatures()
                lsource = layer.source()
                filledpath = getPath(lsource, output + "_" + str(minArea) + "_" + str(maxArea))
                holespath = getPath(lsource, holes  + "_" + str(minArea) + "_" + str(maxArea))
                holes = createEmptyVectorLayer(holespath, dataProvider, fieldsholes, destCrs, layer.wkbType(), provider.encoding())
                filled = createEmptyVectorLayer(filledpath, dataProvider, layer.pendingFields(), destCrs, layer.wkbType(), provider.encoding())
                if filled and holes:
                    for feature in features:
                        filledFeature = QgsFeature(feature)
                        geom = feature.geometry()
                        if not geom.isMultipart():
                            print str(feature.id()) + " is not multipart"
                            polygon = geom.asPolygon()
                            mainring = QgsGeometry.fromPolygon([polygon[0]]) # add first main ring polygon
                            filteredpolygon, ringscodes = filterRingsToHoles(mainring, polygon[1:], minArea, maxArea, holes, filledFeature.id())
                            filledFeature.setGeometry(filteredpolygon)
                        else:
                            multipolyg = geom.asMultiPolygon()
                            print str(feature.id()) + " is multipart"
                            filledGeometry = filledFeature.geometry() # multipolygon
                            for polygon in multipolyg:
                                deleteParterror = filledGeometry.deletePart(0) # delete actual first part
                                if deleteParterror:
                                    part = QgsGeometry.fromPolygon([polygon[0]])
                                    filteredpart, ringscodes = filterRingsToHoles(part, polygon[1:], minArea, maxArea, holes, filledFeature.id())
                                    partcode = filledGeometry.addPartGeometry(filteredpart) # toggle actual part
                                    if partcode > 0:
                                        print "Part of polygon with index " + partindex + " is not valid added, error: " + QgsGeometryaddPart_errors[partcode]
                                else:
                                    print "Unable to delete part " + partindex + " from multipolygon with id " + filledFeature.id() 
                        filled.addFeature(filledFeature)
                    print "Created layer with filled holes " + filledpath
                    print "Created layer with holes " + holespath
                else:
                    print "IS not possible to creaate new vectors "+ holespath + "and" + filledpath +  "(exists?)"
                del holespath
                del filledpath
            else:
                print layername + " is not polygonal valid Vector"
        else:
            layername + "is not valid vector"

def statsFromRasterToPolygon():
    invector, inraster = getFirstOccurence()
    zonalstats = qgis.analysis.QgsZonalStatistics(invector, inraster.source(), stats=qgis.analysis.QgsZonalStatistics.All)
    code = zonalstats.calculateStatistics(None)
    if code > 0:
        print "By calculating occured error " + code
    return code
    
# # selecting # #
def selectByGeometry(selectFunction):
    pass
    
# # addaptive intersect # #
## TODO:USING OGR TO WRITE NEW LAYER
def addaptiveIntersect(inputLayer, intersectLayer, addaptive=50, name = "INTERSECTED", dataProvider = "ESRI Shapefile", lines=True):
    outputpath = getPath(inputLayer.source(), name)
    print inputLayer.wkbType()
    output = createEmptyVectorLayer(outputpath, dataProvider, inputLayer.pendingFields(), inputLayer.crs(), inputLayer.wkbType(), inputLayer.dataProvider().encoding())
    if output:
        layerWkbType = inputLayer.wkbType()
        if inputLayer.type() < 1 and inputLayer.name() != "memory" and inputLayer.hasGeometryType() and intersectLayer.type() < 1 and intersectLayer.name() != "memory" and intersectLayer.hasGeometryType(): #if vector
            for feat in inputLayer.getFeatures():
                ingeom = feat.geometry()
                newfeature = QgsFeature(feat)
                for interf in intersectLayer.getFeatures():
                    intersectgeom = interf.geometry()
                    diffgeom = ingeom.difference(intersectgeom)
                    intersection = ingeom.intersection(intersectgeom)
                    if layerWkbType == QGis.WKBLineString or layerWkbType == QGis.WKBMultiLineString or layerWkbType ==  QGis.WKBLineString25D or layerWkbType == QGis.WKBMultiLineString25D:
                        if intersection.asPolyline() != [] or intersection.asMultiPolyline() != []:
                            if not intersection.isEmpty() and not diffgeom.isEmpty(): # if is not empty after difference
#                                newfeature.setGeometry(intersection)
                                ## here filter for length of intersected
                                if intersection.isMultipart() and diffgeom.isMultipart(): # only for multiparts
                                    intersectpoly = intersection.asMultiPolyline()
                                    diffpoly = diffgeom.asMultiPolyline()
                                    for index, diff in enumerate(diffpoly):# for every line from difference
                                        diffpart = QgsGeometry.fromPolyline(diff)
                                        len = diffpart.length() # check length
                                        if len < addaptive and intersection.touches(diffpart):
                                            if index == 0:
                                                intersectpoly.insert(index, diff) # insert to multi polyline
                                            else:
                                                intersectpoly.insert(index + 1, diff)
                                    newfeature.setGeometry(QgsGeometry.fromMultiPolyline(intersectpoly))
                                else:
                                    twodintersection = intersection.asPolyline()
                                    len = diffgeom.length()
                                    if len < addaptive and intersection.touches(diffgeom):
                                        if diffgeom.asPolyline() != []:
                                            diff2d = QgsGeometry.fromPolyline(diffgeom.asPolyline())
                                        else:
                                            diff2d = QgsGeometry.fromMultiPolyline(diffgeom.asMultiPolyline())
                                        union = QgsGeometry.unaryUnion([QgsGeometry.fromPolyline(twodintersection), diff2d])
                                        newfeature.setGeometry(union)
                                    else:
                                        newfeature.setGeometry(QgsGeometry.fromPolyline(twodintersection))
                                output.addFeature(newfeature)
#                                    if interf.id() == 57:
#                                        print interf.id()
#                                        print newfeature.geometry().exportToWkt()
        else:
            inputLayer.name() + " is not valid Vector"
        print outputpath + " Created "
    else:
        print "IS not possible to creaate new vectors "+ outputpath + "and(exists?)"
    del output