# -*- coding: latin1 -*-
#---------------------------------------------------------------------
# 
# FreeFrame - A QGIS plugin to execute planimetry transformations (LV03 -> LV95)
#                                      for vector and raster data.
#
# Copyright (C) 2010 Stefan Ziegler
#
# EMAIL: stefan.ziegler (at) bd.so.ch
# WEB  : www.catais.org
#
#---------------------------------------------------------------------
# 
# licensed under the terms of GNU GPL 2
# 
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
# 
#---------------------------------------------------------------------

## Import the PyQt and the QGIS libraries
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *

## Import own stuff.
from freeframegui import FreeFrameGui
import utils

## Import python stuff.
import os, sys, math,  copy

## Initialize Qt resources from file resources.py.
import resources

## Our main class for the plugin.
class FreeFrame:

    def __init__(self, iface):
        self.iface = iface
        self.canvas = self.iface.mapCanvas()
          
        self.settings = QSettings("CatAIS","freeframe")      
        
    
    def initGui(self):
        self.index = None
      
        self.action = QAction(QIcon(":/plugins/freeframe/freeframe.png"), "Free Frame", self.iface.mainWindow())
    
        QObject.connect(self.action, SIGNAL("triggered()"), self.run)
                   
        self.iface.addToolBarIcon(self.action)
        self.iface.addPluginToMenu("&FreeFrame", self.action)
        
        
    def run(self):
        self.ctrl = FreeFrameGui(self.iface.mainWindow())
        self.ctrl.initGui()
        self.ctrl.show()        
        
        QObject.connect(self.ctrl, SIGNAL("okClicked()"), self.transform)                       


    def unload(self):
        #QObject.disconnect(self.ctrl, SIGNAL("okClicked()"), self.transform)        
        self.iface.removePluginMenu("&FreeFrame",self.action)
        self.iface.removeToolBarIcon(self.action)      

    def transform(self):
        print "********************** begin"
        self.ctrl.progressBar.setValue(0)
        
        # projectEPSG-Parameter will be used for checking
        # in which reference system we want to transform.
        if self.ctrl.rBtnLV03LV95.isChecked():
            self.projectEPSG = 21781
        elif self.ctrl.rBtnLV95LV03.isChecked():
            self.projectEPSG = 2056
        else:
            self.projectEPSG = 21781
        
        layerName = self.ctrl.inLayer.currentText()
        
        ## Raster transformation is not yet supported.
        type = utils.getLayerTypeByName(layerName)
        if type == 1:
            # Prüfen ob GDAL installiert ist....
            self.transformRasterLayer(layerName)
            
            return 
        else:
            # Prüfen ob spatialite support vorhanden.
            print "Vector Transformation"
            self.transformVectorLayer(layerName)
        
    
    def transformRasterLayer(self,  layerName):
        rl = utils.getRasterLayerByName(layerName)        
        dir = self.settings.value("gui/dirpath")     
        self.ctrl.progressBar.setEnabled(False)
        
        if self.projectEPSG == 21781:
            src_srs = "+proj=somerc +lat_0=46.952405555555555N +lon_0=7.439583333333333E +ellps=bessel +x_0=600000 +y_0=200000 +towgs84=674.374,15.056,405.346  +units=m +nadgrids=" + QDir.convertSeparators(QDir.cleanPath(QgsApplication.qgisSettingsDirPath() + '/python/plugins/freeframe/chenyx06/chenyx06a.gsb')) + " +k_0=1 +wktext"
            dst_srs = "+proj=somerc +lat_0=46.9524055555555N +lon_0=7.439583333333333E +ellps=bessel +x_0=2600000 +y_0=1200000 +towgs84=674.374,15.056,405.346  +units=m +nadgrids=@null +k_0=1 +wktext"
            dst_epsg = "2056"
        else:
            dst_srs = "\"+proj=somerc +lat_0=46.952405555555555N +lon_0=7.439583333333333E +ellps=bessel +x_0=600000 +y_0=200000 +towgs84=674.374,15.056,405.346  +units=m +nadgrids=@null +k_0=1 +wktext\""
            srs_srs = "\"+proj=somerc +lat_0=46.9524055555555N +lon_0=7.439583333333333E +ellps=bessel +x_0=2600000 +y_0=1200000 +towgs84=674.374,15.056,405.346  +units=m +nadgrids=" + QDir.convertSeparators(QDir.cleanPath(QgsApplication.qgisSettingsDirPath() + '/python/plugins/freeframe/chenyx06/chenyx06a.gsb')) + " +k_0=1 +wktext\""
            dst_epsg = "21781"
            
        #self.proc = QProcess()
        #self.proc.start("cp ",  str(rl.source()) + " /home/stefan/tmp/bubu.tif")
        #self.proc.waitForFinished()           
       
#        cmd = "gdalwarp -of GTiff -s_srs " + str(src_srs ) + " -t_srs " + str(dst_srs)  + " "  + str(rl.source()) + " " + str(dir.toString())+os.sep+str(layerName) + "_tmp.tif"
#        print cmd
#        return
            
        arguments = QStringList()
        arguments << "-of"
        arguments << "GTiff"
        arguments << "-s_srs"
        arguments << src_srs
        arguments << "-t_srs"
        arguments << dst_srs
        arguments << str(rl.source())
        arguments << str(dir.toString())+os.sep+str(layerName) + "_tmp.tif"
        try:
            os.remove(str(dir.toString())+os.sep+str(layerName) + "_tmp.tif")        
        except OSError:
            print "tmp file not deleted"
        QApplication.setOverrideCursor(Qt.WaitCursor)
        try:
            # Vielleicht hier noch mit Timeout (-> nicht statische Methode verwenden).    
            QProcess().execute("gdalwarp",  arguments)
        except:
            QApplication.restoreOverrideCursor()
            
        arguments = QStringList()
        arguments << "-of"
        arguments << "GTiff"
        arguments << "-a_srs"
        arguments << "epsg:"+dst_epsg
        arguments << str(dir.toString())+os.sep+str(layerName) + "_tmp.tif"
        arguments << str(dir.toString())+os.sep+str(layerName) + "_ntv2.tif"
        try:
            os.remove(str(dir.toString())+os.sep+str(layerName) + "_ntv2.tif")        
        except OSError:
            print "ntv2 file not deleted"
        # Vielleicht hier noch mit Timeout (-> nicht statische Methode verwenden).    
        try:
            QProcess().execute("gdal_translate",  arguments)
        except:
            QApplication.restoreOverrideCursor()
        QApplication.restoreOverrideCursor()

        try:
            os.remove(str(dir.toString())+os.sep+str(layerName) + "_tmp.tif")
        except OSError:
            print "tmp file not deleted"
    
        self.ctrl.progressBar.setEnabled(True)
        QMessageBox.information(None, 'Information', "Raster dataset transformed.")   


    def transformVectorLayer(self,  layerName):
        ## Layer and project must have same crs.
        vl = utils.getVectorLayerByName(layerName)
        if vl.featureCount() == 0:
            QMessageBox.information(None, 'Information', "No features found.")   
            return
        
#        layerEPSG = vl.srs().epsg()
#        if layerEPSG != self.projectEPSG:
#            QMessageBox.information(None, 'Information', str("Layer (EPSG:" + str(layerEPSG) + str(") and project (EPSG:") + str(self.projectEPSG) +") do not have same crs."))   
#            return 

        if self.projectEPSG == 21781:
            frame_src= "lv03"
            frame_dst = "lv95"
            srs_src = 21781
            srs_dst = 2056
        else:
            frame_src = "lv95"
            frame_dst = "lv03"
            srs_src = 2056
            srs_dst = 21781
            
        ## Load the spatialite db and
        ## build a spatial index.
        uri = QgsDataSourceURI()
        uri.setDatabase(QDir.convertSeparators(QDir.cleanPath(QgsApplication.qgisSettingsDirPath() + '/python/plugins/freeframe/chenyx06/chenyx06.sqlite')))
        
        uri.setDataSource('','chenyx06', 'the_geom_'+str(frame_src))
        self.vl_src = QgsVectorLayer(uri.uri(), 'chenyx06_'+str(frame_src), 'spatialite')  
        self.vl_src.select([], QgsRectangle())
        
        uri.setDataSource('','chenyx06', 'the_geom_'+str(frame_dst))
        self.vl_dst = QgsVectorLayer(uri.uri(), 'chenyx06_'+str(frame_dst), 'spatialite')          
        self.vl_dst.select([], QgsRectangle())
        
        if not (self.vl_src.isValid() or self.vl_dst.isValid()):
            print "Layer failed to load!"
            QMessageBox.information(None, 'Information', str("Error while loading transformation file."))   
            return
        #QgsMapLayerRegistry.instance().addMapLayer(self.vl_dst)
                
        #if self.index == None:
        QApplication.setOverrideCursor(Qt.WaitCursor)
        try:
            if 1 == 1:            
                self.index = QgsSpatialIndex()
                feat = QgsFeature()
                for feat in self.vl_src:
                    self.index.insertFeature(feat)
                self.ctrl.progressBar.setValue(10.0)
            else:
                self.ctrl.progressBar.setValue(10.0)
        except:
            QApplication.restoreOverrideCursor()
        QApplication.restoreOverrideCursor()


        ## Transform the features
        provider = vl.dataProvider()
        allAttrs = provider.attributeIndexes()
        provider.select(allAttrs)

        vtype = vl.wkbType()

        ## Create a vector writer with destination srs (21781 or 2056)
        ## TODO: encoding.... use a qgsfiledialog....
        dir = self.settings.value("gui/dirpath").toString()
        time = QDateTime.currentDateTime()
        shpSuffix = str(time.toString(Qt.ISODate)).replace(":", "").replace("-", "")
        shpFilePath = os.path.join(str(dir), str(layerName) + str("_") + shpSuffix + str(".shp"))
        print shpFilePath
        writer = QgsVectorFileWriter(shpFilePath, "CP1250", provider.fields(), vtype, QgsCoordinateReferenceSystem(srs_dst, QgsCoordinateReferenceSystem().EpsgCrsId))  
        
        if writer.hasError() != QgsVectorFileWriter.NoError:
            print "Error when creating shapefile: ", writer.hasError()
            QMessageBox.information(None, 'Information', str("Error when creating shapefile."))   
            return

        if self.ctrl.useSelected.isChecked():
            start = 10.00
            add = 90.00 / len(vl.selectedFeatures())            
            
            selectedFeatures = vl.selectedFeatures()
            for f in selectedFeatures:
                g = f.geometry()            
                fet = QgsFeature()
                fet.setGeometry(self.transformGeometry(g,  vtype))
                fet.setAttributeMap(f.attributeMap())
                writer.addFeature(fet)
            
                start = start + add
                self.ctrl.progressBar.setValue(start)            
                
                if self.ctrl.progressBar.value() == 99:
                    self.ctrl.progressBar.setValue(100)                  

        else:
            start = 10.00
            add = 90.00 / provider.featureCount()
            QApplication.setOverrideCursor(Qt.WaitCursor)
            try:
                f = QgsFeature()
                while provider.nextFeature(f):  
                    QApplication.processEvents()  
                    g = f.geometry()   
                    fet = QgsFeature()
                    tgeom = self.transformGeometry(g, vtype)
                    if tgeom == None:
                        continue
                    fet.setGeometry(tgeom)
                    fet.setAttributeMap(f.attributeMap())
                    writer.addFeature(fet)
                    
                    start = start + add
                    self.ctrl.progressBar.setValue(start)
                    
                if self.ctrl.progressBar.value() == 99:
                    self.ctrl.progressBar.setValue(100)  
                    
            except:
                QApplication.restoreOverrideCursor()
        QApplication.restoreOverrideCursor()
            
        del writer
        QMessageBox.information(None, 'Information', "Vector dataset transformed.")   


    def transformGeometry(self,  g,  vtype):
        if vtype == QGis.WKBPoint or vtype == QGis.WKBPoint25D:
            p = self.transformPoint(g.asPoint())
            if p == None:
                print "none"
                return None
            return QgsGeometry().fromPoint(p)

        elif vtype == QGis.WKBLineString or vtype == QGis.WKBLineString25D:
            coords = g.asPolyline()
            coords_transformed = []
            for i in coords:
                p = self.transformPoint(i)
                if p == None:
                    return None
                coords_transformed.append(p)
            return QgsGeometry().fromPolyline(coords_transformed)
        
        elif vtype == QGis.WKBPolygon or vtype == QGis.WKBPolygon25D:
            coords = g.asPolygon()
            coords_transformed = []
            ring = []
            for i in coords:
                for k in i: 
                    p = self.transformPoint(k)
                    if p == None:
                        return None
                    ring.append(p)
                coords_transformed .append(ring)
                ring = []
            return QgsGeometry().fromPolygon(coords_transformed )
                
        elif vtype == QGis.WKBMultiPoint or vtype == QGis.WKBMultiPoint25D:
            coords = g.asMultiPoint()
            coords_transformed = []
            for i in coords:
                p = self.transformPoint(i)
                if p == None:
                    return None                
                coords_transformed.append(p)
            return QgsGeometry().fromMultiPoint(coords_transformed)
            
        elif vtype == QGis.WKBMultiLineString or vtype == QGis.WKBMultiLineString25D:
            coords = g.asMultiPolyline()
            coords_transformed = []
            singleline = [] 
            for i in coords:
                for j in i:
                    p = self.transformPoint(j)
                    if p == None:
                        return None
                    singleline.append(p)
                coords_transformed.append(singleline)
                singleline = []
            return QgsGeometry().fromMultiPolyline(coords_transformed)
                
        elif vtype == QGis.WKBMultiPolygon or vtype == QGis.WKBMultiPolygon25D:
            coords = g.asMultiPolygon()
            coords_transformed = []
            ring = []
            for i in coords:
                for j in i:
                    for k in j:
                        p = self.transformPoint(k)
                        if p == None:
                            return None
                        ring.append(p)
                    coords_transformed.append(ring)
                    ring = []
            return QgsGeometry().fromMultiPolygon([coords_transformed])
            
        else:
            QMessageBox.information(None, 'Information', str("Vector type is not yet supported."))   
            return None
    

    def transformPoint(self,  point):
        intersect = self.index.intersects(QgsRectangle(point.x(), point.y(), point.x(), point.y()))        
        for id in intersect:
            feat_src = QgsFeature()
            self.vl_src.featureAtId(int(id),  feat_src,  True,  True)
            if feat_src.geometry().distance(QgsGeometry.fromPoint(point)) == 0:
                
                feat_dst = QgsFeature()
                self.vl_dst.featureAtId(int(id),  feat_dst,  True,  True)
                
                poly_src = feat_src.geometry().asPolygon()
                poly_dst = feat_dst.geometry().asPolygon()
            
                return self.fineltra(point,  poly_src,  poly_dst)
                    
        return None


    def fineltra(self,  point,  poly_src,  poly_dst):
        x1_src = poly_src[0][0].x()
        y1_src = poly_src[0][0].y()
        x2_src = poly_src[0][1].x()
        y2_src = poly_src[0][1].y()
        x3_src = poly_src[0][2].x()
        y3_src = poly_src[0][2].y()      
   
#        print "src point"
#        print "%12.5f" % x1_src
#        print "%12.5f" % y1_src
#        print "%12.5f" % x2_src
#        print "%12.5f" % y2_src
#        print "%12.5f" % x3_src
#        print "%12.5f" % y3_src

        x1_dst = poly_dst[0][0].x()
        y1_dst = poly_dst[0][0].y()
        x2_dst = poly_dst[0][1].x()
        y2_dst = poly_dst[0][1].y()
        x3_dst = poly_dst[0][2].x()
        y3_dst = poly_dst[0][2].y()                     

#        print "dst point"
#        print "%12.5f" % x1_dst 
#        print "%12.5f" % y1_dst 
#        print "%12.5f" % x2_dst 
#        print "%12.5f" % y2_dst 
#        print "%12.5f" % x3_dst 
#        print "%12.5f"  % y3_dst 

        ## Siehe Fineltra Manual Seiten 4-5
        x0_src = point.x()
        y0_src = point.y()

        P1 = math.fabs( 0.5 * ( x0_src*(y2_src-y3_src) + x2_src*(y3_src-y0_src) + x3_src*(y0_src - y2_src) ) )
        P2 = math.fabs( 0.5 * ( x0_src*(y1_src-y3_src) + x1_src*(y3_src-y0_src) + x3_src*(y0_src - y1_src) ) )                
        P3 = math.fabs( 0.5 * ( x0_src*(y1_src-y2_src) + x1_src*(y2_src-y0_src) + x2_src*(y0_src - y1_src) ) )

#        print "P1: " + str(P1*1000)
#        print "P2: " + str(P2*1000)
#        print "P3: " + str(P3*1000)

        ## vxi und vyi berechnen.
        vx1 = x1_dst - x1_src
        vy1 = y1_dst - y1_src
        vx2 = x2_dst - x2_src
        vy2 = y2_dst - y2_src
        vx3 = x3_dst - x3_src
        vy3 = y3_dst - y3_src        
        
#        print "vx1: " + str(vx1)
#        print "vy1: " + str(vy1)
#        print "vx2: " + str(vx2)
#        print "vy2: " + str(vy2)
#        print "vx3: " + str(vx3)
#        print "vy3: " + str(vy3)

        ## Interpolationsverbesserungen.
        DX = (vx1*P1 + vx2*P2 + vx3*P3) / ( P1 + P2 + P3 )
        DY = (vy1*P1 + vy2*P2 + vy3*P3) / ( P1 + P2 + P3 )
                        
        ## Definitive Koordinaten.
        x0_dst = x0_src + DX
        y0_dst = y0_src + DY
   
        return QgsPoint(x0_dst,  y0_dst)       

  
    
 
