# -*- coding: utf-8 -*-

"""
Gnss_NMEA2layer.py
***************************************************************************
*                                                                         *
*   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.                                   *
*                                                                         *
***************************************************************************
"""
__author__ = 'Leandro França'
__date__ = '2022-02-27'
__copyright__ = '(C) 2022, Leandro França'

from qgis.PyQt.QtCore import QMetaType
from qgis.core import *
from lftools.geocapt.imgs import Imgs
from lftools.translations.translate import translate
from lftools.geocapt.cartography import raioMedioGauss
import numpy as np
import datetime
import codecs
import os
from qgis.PyQt.QtGui import QIcon, QColor


class NMEA2layer(QgsProcessingAlgorithm):

    LOC = QgsApplication.locale()[:2]

    def tr(self, *string):
        return translate(string, self.LOC)

    def createInstance(self):
        return NMEA2layer()

    def name(self):
        return 'nmea2layer'

    def displayName(self):
        return self.tr('NMEA to layer', 'NMEA para camada')

    def group(self):
        return self.tr('GNSS', 'GNSS')

    def groupId(self):
        return 'gnss'

    def tags(self):
        return 'GeoOne,gps,navigation,satellites,surveying,glonass,beidou,compass,galileu,track,kinematic,rtk,ntrip,static'.split(',')

    def icon(self):
        return QIcon(os.path.join(os.path.dirname(os.path.dirname(__file__)), 'images/satellite.png'))

    txt_en = '''Loads a NMEA file (protocol 0183) from GNSS receivers as a point layer.
Modes:
◼️ Kinematic - generates all tracked points with their accuracies (PDOP, HDOP and VDOP) and number of satellites.
◼️ Static - calculates the mean and standard deviation of the observed points, for all points or only for fixed solution points (best result).'''
    txt_pt = '''Carrega um arquivo NMEA de rastreio GNSS (protocolo 0183) como uma camada do tipo ponto.
Modos:
◼️ Cinemático - gera todos os pontos rastreados com suas precisões (PDOP, HDOP e VDOP) e número de satélites.
◼️ Estático - calcula a média e desvio-padrão dos pontos observados, para todos os pontos ou somente para os pontos de solução fixa (melhor resultado).'''

    figure = 'images/tutorial/gnss_nmea.jpg'

    def shortHelpString(self):
        social_BW = Imgs().social_BW
        footer = '''<div align="center">
                      <img src="'''+ os.path.join(os.path.dirname(os.path.dirname(__file__)), self.figure) +'''">
                      </div>
                      <div align="left">
                      <p>
                      <b><a href="'''+ self.tr('https://portal.geoone.com.br/m/lessons/pvgnss?classId=3639') + '''" target="_blank">'''+ self.tr('Click here to watch a full class on this tool',
                                    'Clique aqui para assistir uma aula completa sobre esta ferramenta') +'''</a></b>
                      </p>
                      <p>
                      <b><a href="'''+ self.tr('https://geoone.com.br/pvgnss/') + '''" target="_blank">'''+ self.tr('Sign up for the GNSS with RTKLib and QGIS course',
                                    'Inscreva-se no curso de GNSS com RTKLib e QGIS') +'</a><br><br>'+ self.tr('Author: Leandro Franca', 'Autor: Leandro França')+'''</b>
                      </p>'''+ social_BW + '''
                    </div>'''
        return self.tr(self.txt_en, self.txt_pt) + footer


    FILE = 'FILE'
    TYPE = 'TYPE'
    HEIGHT = 'HEIGHT'
    OUTPUT = 'OUTPUT'
    CRS = 'CRS'

    def initAlgorithm(self, config=None):

        self.addParameter(
            QgsProcessingParameterFile(
                self.FILE,
                self.tr('NMEA file .nmea', 'Arquivo NMEA .nmea'),
                behavior = Qgis.ProcessingFileParameterBehavior.File,
                fileFilter = 'NMEA (*.nmea)'
            )
        )

        tipo = [self.tr('Kinematic','Cinemático'),
                self.tr('Static, fixed solution (best result)','Estático, solução fixa  (melhor resultado)'),
                self.tr('Static, all observations','Estático, todas as observações')
               ]

        self.addParameter(
            QgsProcessingParameterEnum(
                self.TYPE,
                self.tr('Attributes', 'Atributos'),
				options = tipo,
                defaultValue= 0
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.HEIGHT,
                self.tr('Antenna height', 'Altura da antena'),
                type = QgsProcessingParameterNumber.Type.Double,
                defaultValue = 0.0,
                minValue = 0
                )
            )

        self.addParameter(
            QgsProcessingParameterCrs(
                self.CRS,
                self.tr('CRS', 'SRC'),
                eval(self.tr("QgsCoordinateReferenceSystem('EPSG:4326')", "QgsCoordinateReferenceSystem('EPSG:4674')"))
                )
            )

        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.OUTPUT,
                self.tr('Point layer', 'Camada de ponto(s)')
            )
        )


    def processAlgorithm(self, parameters, context, feedback):

        caminho = self.parameterAsFile(
            parameters,
            self.FILE,
            context
        )
        if not caminho:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.FILE))

        path, file = os.path.split(caminho)
        nome = file[:-5]

        aa = self.parameterAsDouble(
            parameters,
            self.HEIGHT,
            context
        )

        crs = self.parameterAsCrs(
            parameters,
            self.CRS,
            context
        )
        if not crs.isGeographic():
            raise QgsProcessingException(self.tr('Choose a geographic CRS!', 'Escolha um SRC geográfico!'))

        arq = codecs.open(caminho, 'r', encoding='utf-8', errors='ignore')

        tipo = self.parameterAsEnum(
            parameters,
            self.TYPE,
            context
        )

        lista = []
        PDOP, HDOP, VDOP = -1, -1, -1
        nmea_quality = {0:'0: Invalid', 1:'1: Standalone', 2:'2: DGPS', 3:'3: n/a', 4:'4: RTK fixed', 5:'5: RTK float'}

        for line in arq.readlines():
            try:
                if line[3:6] == 'GSA': # GPS DOP and active satellites
                    partes = line.split(',')
                    PDOP = float(partes[-4])
                    HDOP = float(partes[-3])
                    VDOP = float(partes[-2])

                if line[3:6] == 'GGA': # global position system fix data
                    partes = line.split(',')
                    hora = int(partes[1][0:2])
                    min = int(partes[1][2:4])
                    seg = float(partes[1][4:10])
                    lat = (-1 if partes[3] == 'S' else 1)*( float(partes[2][0:2]) + float(partes[2][2:-1])/60)
                    lon = (-1 if partes[5] == 'W' else 1)*( float(partes[4][0:3]) + float(partes[4][3:-1])/60)
                    quality = nmea_quality[int(partes[6])]
                    num_sat = int(partes[7])
                    HDOP = float(partes[8])
                    H = float(partes[9]) - aa
                    N = float(partes[11])
                    h = N + H
                    if quality != nmea_quality[0]: # Only parse if is valid
                        lista += [[lat, lon, h, H, N, HDOP, VDOP, PDOP, hora, min, seg, quality, num_sat]]

                if line[3:6] == 'ZDA': # Time & Date – UTC, Day, Month, Year and Local Time Zone
                    partes = line.split(',')
                    dia = int(partes[2])
                    mes = int(partes[3])
                    ano = int(partes[4])

                if line[3:6] == 'RMC': # Position, Time & Date – UTC, Day, Month, Year -- (Get only date)
                    partes = line.split(',')
                    if partes[2] == 'A':            # Only parse if is valid
                        dia = int(partes[9][0:2])
                        mes = int(partes[9][2:4])
                        ano = int(partes[9][4:6])+2000
            except:
                pass

        arq.close()
        if len(lista) == 0:
            raise QgsProcessingException(self.tr('No valid observations found!', 'Nenhuma observação válida encontrada!'))
        else:
            valores = np.array(lista)
        # Campos
        if tipo == 0:
            itens  = {"lat": QMetaType.Type.Double,
                      "lon": QMetaType.Type.Double,
                      "h_ellip": QMetaType.Type.Double,
                      "H_orto": QMetaType.Type.Double,
                      "N": QMetaType.Type.Double,
                      "datetime": QMetaType.Type.QString,
                      "HDOP": QMetaType.Type.Double,
                      "VDOP": QMetaType.Type.Double,
                      "PDOP": QMetaType.Type.Double,
                      "quality": QMetaType.Type.QString,
                      "num_sat": QMetaType.Type.Int
                 }
        else:
            itens  = {"lat": QMetaType.Type.Double,
                      "lon": QMetaType.Type.Double,
                      "h_ellip": QMetaType.Type.Double,
                      "H_orto": QMetaType.Type.Double,
                      "N": QMetaType.Type.Double,
                      "start_datetime": QMetaType.Type.QString,
                      "end_datetime": QMetaType.Type.QString,
                      "sigma_x": QMetaType.Type.Double,
                      "sigma_y": QMetaType.Type.Double,
                      "sigma_z": QMetaType.Type.Double,
                      "num_obs": QMetaType.Type.Int
                 }
        Fields = QgsFields()
        for item in itens:
            Fields.append(QgsField(item, itens[item]))

        (sink, dest_id) = self.parameterAsSink( parameters, self.OUTPUT, context, Fields, Qgis.WkbType.PointZ, crs)
        if sink is None:
            raise QgsProcessingException(self.invalidSinkError(parameters, self.OUTPUT))

        if tipo == 0: # Cinemático
            feat = QgsFeature()
            for pnt in lista:
                feat.setGeometry(QgsPoint(pnt[1],pnt[0], pnt[2]))
                try:
                    data_hora = str(datetime.datetime(ano, mes, dia, pnt[-5], pnt[-4], int(pnt[-3])))
                except:
                    data_hora = None
                feat.setAttributes(pnt[0:5] + [data_hora] + pnt[5:8] + pnt[-2:]  ) # lat, lon, h, H, N, HDOP, VDOP, PDOP, hora, min, seg, quality, num_sat
                sink.addFeature(feat, QgsFeatureSink.Flag.FastInsert)

        elif tipo == 1: # Estático com solução fixa
            valores = valores[valores[:,-2] == nmea_quality[4] ] # eliminando observações de baixa qualidade
            if len(valores) == 0:
                raise QgsProcessingException(self.tr('There is no observation with RTK correction.', 'Não existe observação com correção RTK.'))
            # calculo de valores médios
            lat = valores[:,0].astype('float').mean()
            s_lat = valores[:,0].astype('float').std()
            lon = valores[:,1].astype('float').mean()
            s_lon = valores[:,1].astype('float').std()
            h = valores[:,2].astype('float').mean()
            s_h = valores[:,2].astype('float').std()
            H = valores[:,3].astype('float').mean()
            N = valores[:,4].astype('float').mean()
            try:
                data_hora_ini = str(datetime.datetime(ano, mes, dia, int(valores[0, -5]), int(valores[0, -4]), int(float(valores[0, -3]))))
                data_hora_fim = str(datetime.datetime(ano, mes, dia, int(valores[-1, -5]), int(valores[-1, -4]), int(float(valores[-1, -3]))))
            except:
                data_hora_ini, data_hora_fim = None, None
            # Raio Médio de Gauss
            R = raioMedioGauss(lat, self.tr('EPSG:4326','EPSG:4674'))
            sigma_x = (R+h)*np.radians(s_lon)
            sigma_y = (R+h)*np.radians(s_lat)

            feat = QgsFeature()
            feat.setGeometry(QgsPoint(lon, lat, h))
            feat.setAttributes([float(lat),
                                float(lon),
                                float(h),
                                float(H),
                                float(N),
                                data_hora_ini,
                                data_hora_fim,
                                float(sigma_x),
                                float(sigma_y),
                                float(s_h),
                                len(valores)])
            sink.addFeature(feat, QgsFeatureSink.Flag.FastInsert)

        elif tipo == 2: # Estático com todos os valores
            # calculo de valores médios
            lat = valores[:,0].astype('float').mean()
            s_lat = valores[:,0].astype('float').std()
            lon = valores[:,1].astype('float').mean()
            s_lon = valores[:,1].astype('float').std()
            h = valores[:,2].astype('float').mean()
            s_h = valores[:,2].astype('float').std()
            H = valores[:,3].astype('float').mean()
            N = valores[:,4].astype('float').mean()
            try:
                data_hora_ini = str(datetime.datetime(ano, mes, dia, int(valores[0, -5]), int(valores[0, -4]), int(float(valores[0, -3]))))
                data_hora_fim = str(datetime.datetime(ano, mes, dia, int(valores[-1, -5]), int(valores[-1, -4]), int(float(valores[-1, -3]))))
            except:
                data_hora_ini, data_hora_fim = None, None
            # Raio Médio de Gauss
            R = raioMedioGauss(lat, self.tr('EPSG:4326','EPSG:4674'))
            sigma_x = (R+h)*np.radians(s_lon)
            sigma_y = (R+h)*np.radians(s_lat)

            feat = QgsFeature()
            feat.setGeometry(QgsPoint(lon, lat, h))
            feat.setAttributes([float(lat),
                                float(lon),
                                float(h),
                                float(H),
                                float(N),
                                data_hora_ini,
                                data_hora_fim,
                                float(sigma_x),
                                float(sigma_y),
                                float(s_h),
                                len(valores)])
            sink.addFeature(feat, QgsFeatureSink.Flag.FastInsert)

        feedback.pushInfo(self.tr('Operation completed successfully!', 'Operação finalizada com sucesso!'))
        feedback.pushInfo(self.tr('Leandro Franca - Cartographic Engineer', 'Leandro França - Eng Cart'))

        global renamer
        renamer = Renamer(nome)
        context.layerToLoadOnCompletionDetails(dest_id).setPostProcessor(renamer)
        self.SAIDA = dest_id

        return {self.OUTPUT: dest_id}
    
    
    def postProcessAlgorithm(self, context, feedback):
        layer = QgsProcessingUtils.mapLayerFromString(self.SAIDA, context)
        if layer is not None:

           # Criar a simbologia baseada em regras
            root_rule = QgsRuleBasedRenderer.Rule(None)

            # Helper: cria regra com filtro, label e cor
            def add_rule(filter_expr, label, color_hex, size=2.1):
                rule = QgsRuleBasedRenderer.Rule(QgsSymbol.defaultSymbol(layer.geometryType()))
                rule.setFilterExpression(filter_expr)
                rule.setLabel(label)
                rule.symbol().setColor(QColor(color_hex))
                rule.symbol().setSize(size)
                root_rule.appendChild(rule)

            # FIX (verde)
            add_rule('"quality" = \'4: RTK fixed\'',   "4: RTK fixed",    "#33a02c")  # verde

            # FLOAT (laranja)
            add_rule('"quality" = \'5: RTK float\'', "5: RTK float",  "#FF8C00")  # laranja

            # SINGLE (vermelho)
            add_rule('"quality" = \'1: Standalone\'',"1: Standalone", "#FF0000")  # vermelho

            # OUTROS (cinza) - pega tudo que não caiu nas regras acima
            add_rule('ELSE', self.tr('others', 'outras'), "#b0b0b0")      # cinza

            # Aplicar a simbologia baseada em regras na camada
            renderer = QgsRuleBasedRenderer(root_rule)
            layer.setRenderer(renderer)
            layer.triggerRepaint()

        return {}


class Renamer (QgsProcessingLayerPostProcessorInterface):
    def __init__(self, layer_name):
        self.name = layer_name
        super().__init__()

    def postProcessLayer(self, layer, context, feedback):
        layer.setName(self.name)
