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

"""
/***************************************************************************
 GeoINCRA
                                 A QGIS plugin
 Georreferenciamento de Imóveis Rurais
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2022-02-13
        copyright            : (C) 2022 by Tiago Prudencio e Leandro França
        email                : contato@geoone.com.br
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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__ = 'Tiago Prudencio e Leandro França'
__date__ = '2025-05-29'
__copyright__ = '(C) 2025 by Tiago Prudencio e Leandro França'


from qgis.PyQt.QtCore import QCoreApplication
from qgis.core import *
import os, processing
from qgis.PyQt.QtGui import QIcon
import numpy as np


class DividedByRoad(QgsProcessingAlgorithm):

    PARCELA = 'PARCELA'
    VERTICE = 'VERTICE'
    ESTRADA = 'ESTRADA' # linha ou polígono
    DOMINIO = 'DOMINIO'
    SALVAR = 'SALVAR'

    def initAlgorithm(self, config):

        self.addParameter(
			QgsProcessingParameterVectorLayer(
				self.PARCELA,
				self.tr('Camada Parcela (GeoRural)'),
				[QgsProcessing.TypeVectorPolygon],
                optional = False
			)
		)

        self.addParameter(
            QgsProcessingParameterVectorLayer(
                self.VERTICE,
                self.tr('Camada Vértice (GeoRural)'),
                [QgsProcessing.TypeVectorPoint],
                optional = True
                )
        )

        self.addParameter(
			QgsProcessingParameterFeatureSource(
				self.ESTRADA,
				self.tr('Estrada (Linha ou Polígono)'),
				[QgsProcessing.TypeVectorPolygon, QgsProcessing.TypeVectorLine],
                optional = False
			)
		)

        self.addParameter(
            QgsProcessingParameterNumber(
                self.DOMINIO,
                self.tr('Faixa de domínio em metros'),
                type = QgsProcessingParameterNumber.Type.Double,
                minValue = 1,
                optional = True
                )
        )

        self.addParameter(
            QgsProcessingParameterBoolean(
                self.SALVAR,
                self.tr('Salvar modificação'),
                defaultValue = False
            )
        )

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

        parcela = self.parameterAsVectorLayer(
            parameters,
            self.PARCELA,
            context)

        vertice = self.parameterAsVectorLayer(
            parameters,
            self.VERTICE,
            context)

        estrada = self.parameterAsVectorLayer(
            parameters,
            self.ESTRADA,
            context
        )

        dominio = self.parameterAsDouble(
            parameters,
            self.DOMINIO,
            context
        )

        salvar = self.parameterAsBool(
           parameters,
           self.SALVAR,
           context
        )

        # VALIDAÇÕES

        # Camada parcela deve ter apenas uma polígono
        if parcela.featureCount() != 1:
            raise QgsProcessingException('Camada parcela deve ter apenas um único polígono!')

        # Camada estrada deve ter apenas uma feição (linha ou polígono)
        try:
            if estrada.featureCount() != 1:
                raise QgsProcessingException('Camada estrada deve ter apenas uma feição!')
        except:
            raise QgsProcessingException('Camada estrada deve ter apenas uma feição!')

        # Camada da estrada não pode ser igual a camada da parcela
        if parcela.source() == estrada.source():
            raise QgsProcessingException('Camada da estrada não pode ser igual a camada da parcela!')

        # Verificar se camada estrada for do tipo Linha, se a medida da faixa de domínio foi preenchida
        eh_linha = False
        if estrada.wkbType() in (QgsWkbTypes.LineString, QgsWkbTypes.LineStringZ, QgsWkbTypes.MultiLineString, QgsWkbTypes.MultiLineStringZ):
            eh_linha = True
            if not dominio:
                raise QgsProcessingException('Insira a medida da faixa de domínio!')

        # Camada estrada deve ter o mesmo SRC da camada parcela
        reprojetar = False
        if parcela.crs() != estrada.crs():
            crsSrc = estrada.crs()
            crsDest = parcela.crs()
            coordinateTransformer = QgsCoordinateTransform(crsSrc, crsDest, QgsProject.instance())
            reprojetar = True

        # Pegar geometria do polígono da parcela
        for feat in parcela.getFeatures():
            pol_parcela = feat.geometry()
            feat1 = feat

        # Verificar o tipo de geometria da estrada obter geometria da estrada
        feedback.pushInfo('Verificando o tipo de geometria da estrada...')
        for feat in estrada.getFeatures():
            if eh_linha:
                lin_estrada = feat.geometry()
                if estrada.crs().isGeographic():
                    # Transformar a medida do faixa de domínio de metros para graus usando o raio médio de Gauss
                    centroide = lin_estrada.centroid()
                    lat = centroide.asPoint().y()
                    dist = self.metrosParaGraus(dominio, lat)
                    # Reprojetar camada para SIRGAS2000
                    params = {
                        'INPUT': estrada,  # pode ser QgsVectorLayer ou string
                        'TARGET_CRS': QgsCoordinateReferenceSystem('EPSG:4674'),
                        'OUTPUT': 'memory:reprojetado'
                    }
                    resultado = processing.run('native:reprojectlayer', params)
                    camada = resultado['OUTPUT']
                else:
                    camada = estrada
                    dist = dominio
                # Gerar buffer a partir da linha central da estrada
                params_buffer = {
                    'INPUT': camada,
                    'DISTANCE': float(dist),
                    'SEGMENTS': 5,
                    'END_CAP_STYLE': 2,  # Square
                    'JOIN_STYLE': 1,     # Miter , pontiagudo
                    'MITER_LIMIT': 2,
                    'DISSOLVE': False,
                    'OUTPUT': 'memory:buffer_estrada'
                }
                buffer_result = processing.run('native:buffer', params_buffer, context=context, feedback=feedback)
                buffer_layer = buffer_result['OUTPUT']
                for feat_buff in buffer_layer.getFeatures():
                    pol_estrada = feat_buff.geometry()
                if reprojetar:
                    lin_estrada.transform(coordinateTransformer)
                    pol_estrada.transform(coordinateTransformer)
            else:
                pol_estrada = feat.geometry()
                if reprojetar:
                    pol_estrada.transform(coordinateTransformer)

        # Fazer a diferença do polígono pela área da faixa de domínio, o resultado deve ser um multipolígono com 2 aneis (partes)
        if pol_estrada.intersects(pol_parcela):
            diferenca = pol_parcela.difference(pol_estrada)
        else:
            raise QgsProcessingException('Estrada e área do imóvel não se cruzam!')

        if eh_linha:
            # Criar camada temporária para visualização da diferença
            camada_temp = QgsVectorLayer('Polygon?crs=' + parcela.crs().authid(), 'Faixa de domínio', 'memory')
            prov = camada_temp.dataProvider()
            # Criar e adicionar feição com a geometria da diferença
            feat_dif = QgsFeature()
            feat_dif.setGeometry(pol_estrada)
            prov.addFeature(feat_dif)
            # Adicionar camada ao projeto
            QgsProject.instance().addMapLayer(camada_temp)

        # Verificar se a diferença tem 2 partes...

        # Alterar a geometria original da camada parcela
        feedback.pushInfo('Editando geometria da camada parcela...')
        parcela.startEditing() # coloca no modo edição
        parcela.changeGeometry(feat1.id(), diferenca)

        # Preencher vértices tipo V (virtual) PA1 - Paralela na camada Vértice
        if vertice:
            feedback.pushInfo('Adicionando vértices virtuais PA1...')
            campos_vertice = vertice.fields()
            vertice.startEditing()

            # Obter geometria do polígono da parcela como uma lista de pontos
            pontos_parcela = []
            for parte in pol_parcela.asMultiPolygon():
                for anel in parte:
                    for ponto in anel:
                        pontos_parcela.append(ponto)

            # Obter geometria da diferença como uma lista de pontos
            pontos_diferenca = []
            for parte in diferenca.asMultiPolygon():
                for anel in parte:
                    for ponto in anel:
                        pontos_diferenca.append(ponto)

            # Identificar quais vértices do multipolígono da diferença não estão na geometria original da parcela
            for ponto in pontos_diferenca:
                if ponto not in pontos_parcela:
                    fet = QgsFeature(campos_vertice)
                    fet['tipo_verti'] = 'V'
                    fet['metodo_pos'] = 'PA1'
                    fet.setGeometry(QgsGeometry.fromPointXY(ponto))
                    vertice.addFeature(fet)
                    pontos_parcela.append(ponto) # Evitar vértices duplicados

        if salvar: # salva as edições na camada
            parcela.commitChanges()
            if vertice:
                vertice.commitChanges()

        return {'RESULTADO': 'Geometria alterada com sucesso'}


    def name(self):
        return 'DividedByRoad'.lower()

    def displayName(self):
        return self.tr('d. Dividir parcela por Estrada')

    def group(self):
        return self.tr('2. ⚙️ Fluxo de processamento')

    def groupId(self):
        return 'fluxo_processamento'

    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return DividedByRoad()

    def icon(self):
        return QIcon(os.path.join(os.path.dirname(os.path.dirname(__file__)), 'images/geoincra_pb.png'))
    
    def tags(self):
        return 'GeoOne,GeoRural,INCRA,Sigef,separado,estrada,rodovia,glebas,divisão,cortar,multipoligono,parcelar,parcelas,regularização,fundiária'.split(',')
    
    def shortHelpString(self):
        txt = '''Divide o polígono da Parcela de imóvel rural por uma Estrada do tipo linha ou polígono diretamente no banco de dados GeoRural, preenchendo automaticamente os vértices virtuais "V" do tipo "PA1 - Paralela" (opcional).
        A estrada pode ser do tipo linha ou polígono. Se for do tipo linha, é necessário definir a sua medida da faixa de domínio em metros para o cálculo do offset.'''

        footer = '''<div>
                      <div align="center">
                      <a target="_blank" rel="noopener noreferrer" href="https://geoone.com.br/pvgeoincra/"><img title="Inscreva-se no curso de GeoINCRA" style="width: 100%; height: auto;" src="'''+ os.path.join(os.path.dirname(os.path.dirname(__file__)), 'images/INCRA_GeoOne.jpg') +'''"></a>
                      </div>
                      <div align="right">
                      <p align="right">
                      <a href="https://geoone.com.br/pvgeoincra/"><span style="font-weight: bold;">Conheça o curso de GeoINCRA no QGIS</span></a>
                      </p>
                      <p align="right">
                      <a href="https://portal.geoone.com.br/m/lessons/geoincra?classId=5808"><span style="font-weight: bold;">Acesse a aula sobre esta ferramenta no curso de GeoINCRA no GeoOne</span></a>
                      </p>
                      <a target="_blank" rel="noopener noreferrer" href="https://geoone.com.br/"><img title="GeoOne" width="280"  src="'''+ os.path.join(os.path.dirname(os.path.dirname(__file__)), 'images/GeoOne.png') +'''"></a>
                      <p><i>"Mapeamento automatizado, fácil e direto ao ponto é na GeoOne!"</i></p>
                      </div>
                    </div>'''
        return txt + footer

    def metrosParaGraus(self, dist, lat):
        crsGeo = QgsCoordinateReferenceSystem('EPSG:4674')
        ellipsoid_id = crsGeo.ellipsoidAcronym()
        ellipsoid = QgsEllipsoidUtils.ellipsoidParameters(ellipsoid_id)
        a = ellipsoid.semiMajor
        f_inv = ellipsoid.inverseFlattening
        f = 1/f_inv
        e2 = f*(2-f)
        N = a/np.sqrt(1-e2*(np.sin(lat))**2) # Raio de curvatura 1º vertical
        M = a*(1-e2)/(1-e2*(np.sin(lat))**2)**(3/2.) # Raio de curvatura meridiana
        R = np.sqrt(M*N) # Raio médio de Gauss
        theta = dist/R
        theta = np.degrees(theta) # Radianos para graus
        return theta
