import rasterio
from os import getcwd, path
import sys
sys.path.append(path.join(getcwd(), "thirdparty"))
from core.Data import *
from core.Parameters import *
from core.PreProcessing2D import *
from tool_mesh import *
from scripts.only_banks import  *
from cephee_cli import *



param_init= read_param()
work_directory = '/Users/cassan/PycharmProjects/cephee/tmp/stOmer/work_network_filtered2'
param_init.work_path = work_directory
param_init.verbose = False
##Input
param_init['I']['DEM_path'] = '/Users/cassan/PycharmProjects/cephee/tmp/stOmer/data/DEM_1m_T2D_dig.tif'
param_init['C']['DEM_CRS_ID'] = 'EPSG:2154'
param_init['I']['network_filepath'] = '/Users/cassan/PycharmProjects/cephee/tmp/stOmer/data/centerline_topage_all.shp'
param_init['I']['outlet_point'] = Point(644130 ,7082661)
boundary_line = '/Users/cassan/PycharmProjects/cephee/tmp/stOmer/data/masque_BV_T2D.shp'
size_mesh = 20


#output
bank_file=os.path.join(param_init.work_path, 'banks_lines_Himposed_h0.2.shp') #/Users/cassan/PycharmProjects/cephee/tmp/vicdessos_complet/banks_lines.shp'
shift_bank_name = '/Users/cassan/PycharmProjects/cephee/tmp/stOmer/data/shift_banks_all_step1.shp'#/Users/cassan/PycharmProjects/cephee/tmp/vicdessos_complet/banks_lines.shp'#
shift_bank_name2 = '/Users/cassan/PycharmProjects/cephee/tmp/stOmer/data/shift_banks_all_step2.shp'

##type of process
max_water_extent = True
detect_subbassin = False
detect_dike = False #trouve les berge
is_mesh_dike = False
detect_and_correct_bank =True#True#True# sans calcul hydraulique
adjust_polygone = True#écarte les lignes de size mesh
export_line_for_bluekenue =True
is_density = True
find_centerline = True
is_mesh_channel = True

create_mesh = False

#constant parameters
tolerance = 2*size_mesh
dilate_distance = 1.1 * size_mesh
param_init['C']['computeGlobal'] = False

if max_water_extent:
    '''Détecte les zones potentiellement en eau pour les débits max
    '''

    param_init['XS']['creation_step'] = 5000
    param_init['XS']['width'] = 1000
    param_init['N']['npoly'] = 2
    # parametres à garder fixes
    param_init['XS']['number_of_points'] = 1000
    param_init['H']['dxlat'] = -1
    param_init['H']['himposed'] = 0.5
    param_init['H']['createBanks'] = True
    param_init['H']['createBanksMethods'] = 'Normal'
    param_init['H']['output_resolution'] = 10
    param_init['H']['mapping_method'] = 'by_polygon'
    param_init['H']['outletDischarge'] = 100
    BV1 = cephee_cli(param_init)

    water_extent_file = water_extent(BV1, param_init,size_mesh)

    print('water mask created')
dqfs
if detect_subbassin:
    '''
    Découpe le DEM en sous bassin pour avoir les lignes de contraintes hautes
    '''
    BV = ModelCatchment(qgis_mode=False)
    param_init['C']['computeGlobal'] = True
    param_init['N']['classeMax'] = 4
    param_init['N']['minAccumulativeArea'] = 1000
    param_init['C']['resolution'] = 2
    init_DEM_and_network(BV, param_init)
    param_init['C']['computeGlobal'] = False



if detect_dike:
    '''trouve les digues par sur élévation à la moyenne glissante de l altitude
    '''
    area_min = 10000 # aire maximale des zones en m2
    output_name = os.path.join(work_directory, 'constraint_lines_' + str(area_min) + 'm2')
    width_range =  [10,20]#np.arange(5,15,10)
    height_range = [0.4]#np.arange(0.4,1.4,0.2)
    _, nstructure, multi_polygones = find_dikes([],param_init['I']['DEM_path'],areaMin = area_min , output_name = output_name,method ='avg',
        width_range =width_range ,height_range = height_range,res =None )
    gdf_polygone = smooth_polygone(multi_polygones, 1000, dilate_distance=dilate_distance,
                                   save_path=work_directory)


    if False:#detect_subbassin or os.path.exists(path.join(param_init.work_path, 'boundary_area.shp')):

        gdf_crest = gpd.read_file(path.join(param_init.work_path, 'boundary_area.shp'))
        gdf_polygone.drop_duplicates(subset='geometry')
        gdf_crest.drop_duplicates(subset='geometry')
        gdf_crest = detect_crest_from_area_acc(gdf_polygone, gdf_crest, dilate_distance=dilate_distance,
                                   save_path=work_directory)
        gdf_intersection = detect_crest_from_area_acc(gdf_polygone, gdf_crest)
        gdf_intersection.to_file(os.path.join(param_init.work_path, 'crest_lines_from_area.shp'))



if is_mesh_dike:
    maxWidth = 50
    step_long = 5*size_mesh
    step_lat = 2

    multi_polygones = os.path.join(param_init.work_path, 'polygone_smooth.shp')
    gdf_polygones = gpd.read_file(multi_polygones )
    #maillage des digues####################
    gdf_dike = mesh_dike(param_init['I']['DEM_path'], gdf_polygones,param_init['C']['DEM_CRS_ID'],maxWidth,step_long = step_long, step_lat = step_lat ,save_path = work_directory )
    gdf_dike.to_file(os.path.join(param_init.work_path, 'crest_lines.shp'))
    write_i2s(gdf_dike, param_init.work_path, size_mesh, 'crest_lines')



if detect_and_correct_bank:
    '''Détecte les berges sur le DEM par une sur cote, ne garde que celle pour une largeur plus grande qie mesh_size
    '''

    param_init['XS']['creation_step'] = size_mesh*2
    param_init['XS']['width'] = 50
    #parametres à garder fixes
    param_init['XS']['number_of_points'] = 50
    param_init['H']['dxlat'] =-1
    param_init['H']['himposed'] = 0.5
    param_init['H']['createBanks']= True
    param_init['H']['createBanksMethods']= 'Himposed'


    bank_file_init = only_banks(param_init, size_mesh)
    gdf_bank_init = gpd.read_file(bank_file_init)
    #gdf_bank = enforce_min_distance(gdf_bank_init, size_mesh, iterations=10, step=0.1)
    multi_polygones = os.path.join(param_init.work_path, 'polygone_smooth.shp')
    gdf_polygones = gpd.read_file(multi_polygones )
    crest_file = os.path.join(param_init.work_path, 'crest_lines.shp')
    gdf_crest = gpd.read_file(crest_file)

    gdf_bank = enforce_min_distance(
        gdf_bank_init,
        size_mesh,
        iterations=1,
        step=1,
        #gdf_polygone=gdf_polygones,
        gdf_lignes_buffer=gdf_crest,
        use_buffer_distance=True
    )

    gdf_bank.to_file(os.path.join(param_init.work_path, 'banks_lines_OK.shp'))
    write_i2s(gdf_bank,param_init.work_path, size_mesh,'bank_lines')

    print( 'Bank file has been created and checked')

if find_centerline:
    bank_file = os.path.join(param_init.work_path, 'banks_lines_all_banks.shp')
    gdf_bank = gpd.read_file(bank_file)

    centerlines, outline = sort_lines_and_centerline(gdf_bank,size_mesh*2)
    gdf_centerline = gpd.GeoDataFrame(geometry=centerlines)
    gdf_centerline.to_file(os.path.join(param_init.work_path, 'centerline_from_bank.shp'))
    write_i2s(gdf_centerline, param_init.work_path, size_mesh, 'centerline')



if adjust_polygone:
    '''Bouge les berges pour avoir la distance d'une maille netre les lignes de contraintes
    '''

    bank_file = os.path.join(param_init.work_path, 'banks_lines_OK.shp')
    gdf_bank= gpd.read_file(bank_file)
    #gdf_bank = enforce_min_distance(gdf_bank_init, size_mesh, iterations=10, step=0.1)
    multi_polygones = os.path.join(param_init.work_path, 'polygone_smooth.shp')
    gdf_polygones = gpd.read_file(multi_polygones )
    crest_file = os.path.join(param_init.work_path, 'crest_lines.shp')
    gdf_crest = gpd.read_file(crest_file)

    gdf_open_polygone = extract_lines_outside_buffers(gdf_polygones, gdf_bank ,gdf_crest,size_mesh)
    gdf_open_polygone.to_file(os.path.join(param_init.work_path, 'open_polygone.shp'))
    write_i2s(gdf_open_polygone, param_init.work_path, size_mesh, 'open_polygone')



if export_line_for_bluekenue:
    '''exporte les berges au format i2s pour faire des maillages structurés dans les lits
    '''
    separate_bank_to_i2s(shift_bank_name,work_directory ,shp =False)



if is_density:
    '''
    créer une carte de densité au format bluekenue poiur indiquer la taille de maille minimale dans les zones de digues
    '''
    ###création des lignes pour l'aide à la création du maillage sous bluekenue######
    multi_polygones = os.path.join(param_init.work_path, 'open_polygone.shp')
    gdf_polygone = gpd.read_file(multi_polygones)

    bank_file = os.path.join(param_init.work_path, 'banks_lines_OK.shp')
    gdf_bank= gpd.read_file(bank_file)
    density_map_from_poly(gdf_polygone,size_mesh, dilate_distance,boundary_line = boundary_line,gdf_bank =gdf_bank,save_path =work_directory)
    print('density map created')

if is_mesh_channel:
    maxWidth = 50
    step_long = 20
    nb_trans_point = 5
    #maillage des lit mineurs####################
    bank_file = os.path.join(param_init.work_path, 'banks_lines_OK.shp')
    gdf_bank= gpd.read_file(bank_file)
    gdf_channel = mesh_channel(bank_file, 'EPSG:2154',maxWidth,step_long, nb_trans_point,method = 'cephee', save_path = work_directory )
    channel_file = os.path.join(param_init.work_path, 'channel.shp')
    gdf_channel.to_file(channel_file)


if create_mesh:
    import shapefile
    import triangle
    import matplotlib.pyplot as plt

    import gmsh
    from shapely.geometry import Polygon, LineString, Point


    # Charger plusieurs shapefiles de contraintes
    files = [os.path.join(param_init.work_path,"banks_lines_OK.shp"),
             os.path.join(param_init.work_path,"centerline_from_bank.shp"),
             os.path.join(param_init.work_path,"crest_lines.shp"),
             os.path.join(param_init.work_path,"open_polygone.shp")]
    gdf_contraintes = gpd.GeoDataFrame(pd.concat([gpd.read_file(f) for f in files], ignore_index=True))

    # ------------------------------------------------------
    # Paramètres
    # ------------------------------------------------------
    lc = 200  # Taille globale de maille
    lc_ref = 10  # Taille locale (raffinement)
    mesh_file = "maillage.msh"

    zones_shp = os.path.join(param_init.work_path,"bank_density.shp")  # Polygones ou points

    # ------------------------------------------------------
    # Initialiser Gmsh
    # ------------------------------------------------------
    gmsh.initialize()
    gmsh.model.add("maillage_shp")

    # ------------------------------------------------------
    # 1. Charger contour extérieur
    # ------------------------------------------------------
    gdf_contour = gpd.read_file('/Users/cassan/PycharmProjects/cephee/tmp/stOmer/data/test_gmesh.shp')
    poly: Polygon = gdf_contour.geometry.iloc[0]

    # Créer points et lignes du contour extérieur
    pts = []
    for (x, y) in poly.exterior.coords[:-1]:  # on ignore le dernier car = premier
        pts.append(gmsh.model.geo.addPoint(x, y, 0, lc))

    lines = []
    for i in range(len(pts)):
        lines.append(gmsh.model.geo.addLine(pts[i], pts[(i + 1) % len(pts)]))

    # Boucle fermée pour la surface extérieure
    ext_loop = gmsh.model.geo.addCurveLoop(lines)

    # Si trous intérieurs (holes)
    int_loops = []
    for interior in poly.interiors:
        pts_hole = [gmsh.model.geo.addPoint(x, y, 0, lc) for (x, y) in interior.coords[:-1]]
        lines_hole = [gmsh.model.geo.addLine(pts_hole[i], pts_hole[(i + 1) % len(pts_hole)]) for i in
                      range(len(pts_hole))]
        int_loops.append(gmsh.model.geo.addCurveLoop(lines_hole))

    # Créer la surface avec trous éventuels
    surf = gmsh.model.geo.addPlaneSurface([ext_loop] + int_loops)

    # ------------------------------------------------------
    # 2. Contraintes internes (lignes)
    # ------------------------------------------------------

    for line in gdf_contraintes.geometry:
        if isinstance(line, LineString):
            for (x, y, *_) in line.coords:
                pts_c = [gmsh.model.geo.addPoint(x, y, 0, lc) for (x, y, *_) in line.coords]
            #pts_c = [gmsh.model.geo.addPoint(x, y, 0, lc) for (x, y) in line.coords]
            lns_c = [gmsh.model.geo.addLine(pts_c[i], pts_c[i + 1]) for i in range(len(pts_c) - 1)]
            for l in lns_c:
                gmsh.model.geo.addCurveLoop([l])
                # forcer la ligne à apparaître dans le maillage
                gmsh.model.geo.synchronize()
                gmsh.model.mesh.embed(1, [l], 2, surf)

    # ------------------------------------------------------
    # 3. Zones de raffinement
    # ------------------------------------------------------
    gdf_zones = gpd.read_file(zones_shp)
    zone_entities = []
    for geom in gdf_zones.geometry:
        if isinstance(geom, Polygon):
            pts_z = [gmsh.model.geo.addPoint(x, y, 0, lc) for (x, y) in geom.exterior.coords[:-1]]
            lines_z = [gmsh.model.geo.addLine(pts_z[i], pts_z[(i + 1) % len(pts_z)]) for i in range(len(pts_z))]
            loop_z = gmsh.model.geo.addCurveLoop(lines_z)
            surf_z = gmsh.model.geo.addPlaneSurface([loop_z])
            zone_entities.append(surf_z)
        elif isinstance(geom, Point):
            p = gmsh.model.geo.addPoint(geom.x, geom.y, 0, lc)
            zone_entities.append(p)

    gmsh.model.geo.synchronize()

    # Définir un champ Distance basé sur ces zones
    if zone_entities:
        f1 = gmsh.model.mesh.field.add("Distance")
        gmsh.model.mesh.field.setNumbers(f1, "FacesList", zone_entities)

        f2 = gmsh.model.mesh.field.add("Threshold")
        gmsh.model.mesh.field.setNumber(f2, "IField", f1)
        gmsh.model.mesh.field.setNumber(f2, "LcMin", lc_ref)
        gmsh.model.mesh.field.setNumber(f2, "LcMax", lc)
        gmsh.model.mesh.field.setNumber(f2, "DistMin", 0.1)
        gmsh.model.mesh.field.setNumber(f2, "DistMax", 0.3)

        gmsh.model.mesh.field.setAsBackgroundMesh(f2)

    # ------------------------------------------------------
    # 4. Génération du maillage
    # ------------------------------------------------------
    gmsh.model.geo.synchronize()
    gmsh.model.mesh.generate(2)
    gmsh.write(mesh_file)

    print(f"✅ Maillage généré : {mesh_file}")

    gmsh.finalize()