Main

mail@pastecode.io avatar
unknown
plain_text
2 years ago
48 kB
3
Indexable
Never
import numpy as np

import matplotlib.pyplot as plt
import scipy

from netCDF4 import Dataset as NetCDFFile
from mpl_toolkits.basemap import Basemap
import skfuzzy as fuzz
from skfuzzy import control as ctrl

from geopy import distance


# FAIRE GAFFE, LES DONNES NETCDF SONT : [ORDONNEE;ABSCISSE]

# mettre les notations de Kauffman sur le rapport / dire que la T norm c'est min / Bibliographie,... TOUT METTRE
# dire ce qu'on aurait pu faire si on avait plus de temps

# longitude / latitude / time / d = divergence / r = humidity / t = temperature / u=uwind / v=vwind / vo=vortcity

def plotNC(path, var, unite, titre, date_plot):
    nc_plot = NetCDFFile(path)
    lat_plot = nc_plot.variables['latitude'][:]
    lon_plot = nc_plot.variables['longitude'][:]
    var = nc_plot.variables[var][:]
    nc_plot.close()
    print(var)
    map_plot = Basemap(width=5000000, height=3500000,
                       resolution='l', projection='cyl',
                       llcrnrlon=lon_plot.min(), llcrnrlat=lat_plot.min(), urcrnrlon=lon_plot.max(),
                       urcrnrlat=lat_plot.max(), lat_0=lat_plot.mean(),
                       lon_0=lon_plot.mean())
    lons_plot, lats_plot = np.meshgrid(lon_plot, lat_plot)
    xi, yi = map_plot(lons_plot, lats_plot)

    map_plot.drawmapboundary(fill_color='aqua')
    map_plot.fillcontinents(color='coral', lake_color='aqua')
    map_plot.drawcoastlines()

    parallels = np.arange(lat_plot.min(), lat_plot.max(), 5.)  # make latitude lines every 5 degrees from xx to xx
    meridians = np.arange(lon_plot.min(), lat_plot.min(), 5.)  # make longitude lines every 5 degrees from xx to xx
    map_plot.drawparallels(parallels, labels=[1, 0, 0, 0], fontsize=10)
    map_plot.drawmeridians(meridians, labels=[0, 0, 0, 1], fontsize=10)

    cs = map_plot.pcolor(xi, yi, np.squeeze(var[date_plot, :, :]))
    cbar = map_plot.colorbar(cs, location='bottom', pad="10%")
    cbar.set_label(unite)
    plt.title(titre)


def plotTC(path, date_plot):
    plt.subplot(3, 3, 1)
    plotNC(path, 'r', '%', 'Humidité relative', date_plot)
    plt.plot(-82, 24, ms=10, marker="o", markeredgecolor="red")
    plt.subplot(3, 3, 2)
    plotNC(path, 'd', 'jsp', 'Divergence', date_plot)
    plt.subplot(3, 3, 3)
    plotNC(path, 'vo', 'jsp', 'Vorticité', date_plot)
    # plt.subplot(3, 3, 4)
    # plotNC(path, 't', 'Kelvin', 'Température')
    plt.subplot(3, 3, 5)
    plotNC(path, 'u', 'm/s', 'Uwind', date_plot)  # U wind = composante horizontontale
    plt.subplot(3, 3, 6)
    plotNC(path, 'v', 'm/s', 'Vwind', date_plot)  # V wind = composante verticale


def plot_array(array, x, y):
    lon_plot = met_to_deg(x)
    lat_plot = met_to_deg(y)

    map_plot = Basemap(width=5000000, height=3500000,
                       resolution='l', projection='cyl',
                       llcrnrlon=lon_plot.min(), llcrnrlat=lat_plot.min(), urcrnrlon=lon_plot.max(),
                       urcrnrlat=lat_plot.max(), lat_0=lat_plot.mean(),
                       lon_0=lon_plot.mean())
    lons, lats = np.meshgrid(lon_plot, lat_plot)
    xi, yi = map_plot(lons, lats)

    map_plot.drawmapboundary(fill_color='aqua')
    map_plot.fillcontinents(color='coral', lake_color='aqua')
    map_plot.drawcoastlines()
    cs = map_plot.pcolor(xi, yi, array)
    cbar = map_plot.colorbar(cs, location='bottom', pad="10%")


def input_fuzz(variable, input_to_fuzz):
    fuzz_input = []
    for vari in variable.terms.keys():
        fuzz_input.append(fuzz.interp_membership(variable.universe, variable.terms[vari].mf, input_to_fuzz))
    return fuzz_input


def IRR_2var(SF_rules, input_fuzz1, input_fuzz2):
    IRR = np.empty([len(SF_rules), 2])
    k = -1
    for j in range(0, len(SF_rules)):
        if (j % len(input_fuzz2)) == 0:
            k += 1
        IRR[j, 0] = input_fuzz1[k]
        IRR[j, 1] = input_fuzz2[j % len(input_fuzz2)]

    return IRR


def IRR_3var(SF_rules, input_fuzz1, input_fuzz2, input_fuzz3):
    IRR = np.empty([len(SF_rules), 3])

    k = -1
    m = -1
    for j in range(0, len(SF_rules)):
        if (j % len(input_fuzz3)) == 0:
            k += 1
            if (k % len(input_fuzz2)) == 0:
                m += 1
        IRR[j, 0] = input_fuzz1[m]
        IRR[j, 1] = input_fuzz2[j % len(input_fuzz2)]
        IRR[j, 2] = input_fuzz3[j % len(input_fuzz3)]
    return IRR


def SF1_compute(input_temperature, input_d_speed, input_time):
    pas = 0.01
    temperature = ctrl.Antecedent(np.arange(263, 323, pas), 'Temperature (K)')
    d_speed = ctrl.Antecedent(np.arange(-100, 100, pas), 'Variation of wind speed (m/s)')
    timeTC = ctrl.Antecedent(np.arange(0, 15, pas), 'Time since TC formation (days)')
    lifespan = ctrl.Consequent(np.arange(0, 15, pas), 'Lifespan')

    temperature['Low'] = fuzz.trapmf(temperature.universe, [263, 263, 293, 303])
    temperature['Ok'] = fuzz.trapmf(temperature.universe, [293, 303, 323, 323])

    d_speed['--'] = fuzz.trapmf(d_speed.universe, [-100, -100, -60, -20])
    d_speed['-'] = fuzz.trimf(d_speed.universe, [-60, -30, 0])
    d_speed['St'] = fuzz.trimf(d_speed.universe, [-30, 0, 20])
    d_speed['+'] = fuzz.trimf(d_speed.universe, [0, 30, 60])
    d_speed['++'] = fuzz.trapmf(d_speed.universe, [20, 60, 100, 100])

    timeTC['Young'] = fuzz.zmf(timeTC.universe, 3, 6)
    timeTC['Middle'] = fuzz.gaussmf(timeTC.universe, 6, 1)
    timeTC['Old'] = fuzz.smf(timeTC.universe, 6, 9)

    lifespan['Dying+'] = fuzz.trapmf(lifespan.universe, [0, 0, 3, 5])
    lifespan['Dying'] = fuzz.trimf(lifespan.universe, [3, 5, 7])
    lifespan['St'] = fuzz.trimf(lifespan.universe, [5, 7, 10])
    lifespan['Living'] = fuzz.trimf(lifespan.universe, [7, 10, 13])
    lifespan['Living+'] = fuzz.trapmf(lifespan.universe, [10, 13, 15, 15])

    lifespan.view()
    temperature.view()
    d_speed.view()
    timeTC.view()

    SF1_rules = [ctrl.Rule(temperature['Low'] & d_speed['--'] & timeTC['Young'], lifespan['Dying+']),
                 ctrl.Rule(temperature['Low'] & d_speed['--'] & timeTC['Middle'], lifespan['Dying+']),
                 ctrl.Rule(temperature['Low'] & d_speed['--'] & timeTC['Old'], lifespan['Dying+']),
                 ctrl.Rule(temperature['Low'] & d_speed['-'] & timeTC['Young'], lifespan['Dying']),
                 ctrl.Rule(temperature['Low'] & d_speed['-'] & timeTC['Middle'], lifespan['Dying+']),
                 ctrl.Rule(temperature['Low'] & d_speed['-'] & timeTC['Old'], lifespan['Dying+']),
                 ctrl.Rule(temperature['Low'] & d_speed['St'] & timeTC['Young'], lifespan['Dying']),
                 ctrl.Rule(temperature['Low'] & d_speed['St'] & timeTC['Middle'], lifespan['Dying']),
                 ctrl.Rule(temperature['Low'] & d_speed['St'] & timeTC['Old'], lifespan['Dying']),
                 ctrl.Rule(temperature['Low'] & d_speed['+'] & timeTC['Young'], lifespan['St']),
                 ctrl.Rule(temperature['Low'] & d_speed['+'] & timeTC['Middle'], lifespan['St']),
                 ctrl.Rule(temperature['Low'] & d_speed['+'] & timeTC['Old'], lifespan['St']),
                 ctrl.Rule(temperature['Low'] & d_speed['++'] & timeTC['Young'], lifespan['Living']),
                 ctrl.Rule(temperature['Low'] & d_speed['++'] & timeTC['Middle'], lifespan['Living']),
                 ctrl.Rule(temperature['Low'] & d_speed['++'] & timeTC['Old'], lifespan['Living']),
                 ctrl.Rule(temperature['Ok'] & d_speed['--'] & timeTC['Young'], lifespan['Dying']),
                 ctrl.Rule(temperature['Ok'] & d_speed['--'] & timeTC['Middle'], lifespan['Dying']),
                 ctrl.Rule(temperature['Ok'] & d_speed['--'] & timeTC['Old'], lifespan['Dying']),
                 ctrl.Rule(temperature['Ok'] & d_speed['-'] & timeTC['Young'], lifespan['Dying']),
                 ctrl.Rule(temperature['Ok'] & d_speed['-'] & timeTC['Middle'], lifespan['Dying']),
                 ctrl.Rule(temperature['Ok'] & d_speed['-'] & timeTC['Old'], lifespan['Dying']),
                 ctrl.Rule(temperature['Ok'] & d_speed['St'] & timeTC['Young'], lifespan['Living']),
                 ctrl.Rule(temperature['Ok'] & d_speed['St'] & timeTC['Middle'], lifespan['St']),
                 ctrl.Rule(temperature['Ok'] & d_speed['St'] & timeTC['Old'], lifespan['Dying']),
                 ctrl.Rule(temperature['Ok'] & d_speed['+'] & timeTC['Young'], lifespan['Living']),
                 ctrl.Rule(temperature['Ok'] & d_speed['+'] & timeTC['Middle'], lifespan['Living']),
                 ctrl.Rule(temperature['Ok'] & d_speed['+'] & timeTC['Old'], lifespan['Living']),
                 ctrl.Rule(temperature['Ok'] & d_speed['++'] & timeTC['Young'], lifespan['Living+']),
                 ctrl.Rule(temperature['Ok'] & d_speed['++'] & timeTC['Middle'], lifespan['Living+']),
                 ctrl.Rule(temperature['Ok'] & d_speed['++'] & timeTC['Old'], lifespan['Living+']),
                 ]
    SF1_ctrl = ctrl.ControlSystem(SF1_rules)
    SF1 = ctrl.ControlSystemSimulation(SF1_ctrl)

    input_temperature_fuzz = input_fuzz(temperature, input_temperature)
    input_d_speed_fuzz = input_fuzz(d_speed, input_d_speed)
    input_time_fuzz = input_fuzz(timeTC, input_time)

    print(input_temperature_fuzz)
    print(input_d_speed_fuzz)
    print(input_time_fuzz)
    IRR = IRR_3var(SF1_rules, input_temperature_fuzz, input_d_speed_fuzz, input_time_fuzz)

    declenchement = IRR.min(axis=1)
    print(IRR)
    print(declenchement)

    csq_SF1 = {}  # variable de type dictionnaire
    for rule_i in range(0, len(SF1_rules)):  # on initialise à 0 pour chaque règle
        csq_SF1[str(SF1_rules[rule_i].consequent)] = 0
    for rule_i in range(0, len(SF1_rules)):
        csq_i = str(SF1_rules[rule_i].consequent)
        csq_SF1[csq_i] = max(csq_SF1[csq_i], declenchement[rule_i])  # Comme dans le TP6

    print("Conséquences SF1 :", csq_SF1)

    # return SF1


def SF2_compute(input_humidity, input_TC_size):
    pas = 0.01
    humidity = ctrl.Antecedent(np.arange(0, 100, pas), 'humidity')  # Humidity (%)'
    TC_size = ctrl.Antecedent(np.arange(0, 350, pas), 'TC_size')  # TC size (km)
    TC_size_variation = ctrl.Consequent(np.arange(0, 60, pas), 'TC_var')  # Variation of TC size for next day (km)

    humidity['--'] = fuzz.trapmf(humidity.universe, [0, 0, 20, 40])  # 20%
    humidity['-'] = fuzz.trimf(humidity.universe, [20, 40, 60])  # 40%
    humidity['+'] = fuzz.trimf(humidity.universe, [40, 60, 80])  # 60%
    humidity['++'] = fuzz.trapmf(humidity.universe, [60, 80, 100, 100])  # 80%

    TC_size['Little+'] = fuzz.trapmf(TC_size.universe, [0, 0, 75, 120])
    TC_size['Little'] = fuzz.trimf(TC_size.universe, [75, 150, 200])
    # TC_size['Normal'] = fuzz.trimf(TC_size.universe, [120, 150, 200])
    TC_size['Huge'] = fuzz.trimf(TC_size.universe, [150, 200, 250])
    TC_size['Huge+'] = fuzz.trapmf(TC_size.universe, [200, 250, 350, 350])

    # 10 13 15 20 25 50

    TC_size_variation['Very Low'] = fuzz.trapmf(TC_size_variation.universe, [0, 0, 10, 13])  # 10
    TC_size_variation['Low'] = fuzz.trimf(TC_size_variation.universe, [10, 13, 15])  # 13
    TC_size_variation['Quite Low'] = fuzz.trimf(TC_size_variation.universe, [13, 15, 20])  # 15
    TC_size_variation['Quite High'] = fuzz.trimf(TC_size_variation.universe, [15, 20, 25])  # 20
    TC_size_variation['High'] = fuzz.trimf(TC_size_variation.universe, [20, 25, 50])  # 25
    TC_size_variation['Very High'] = fuzz.trapmf(TC_size_variation.universe, [25, 50, 50, 50])  # 50

    humidity.view()
    TC_size.view()
    TC_size_variation.view()
    SF2_rules = [ctrl.Rule(humidity['--'] & TC_size['Little+'], TC_size_variation['Very Low']),
                 ctrl.Rule(humidity['--'] & TC_size['Little'], TC_size_variation['Quite Low']),
                 ctrl.Rule(humidity['--'] & TC_size['Huge'], TC_size_variation['Quite Low']),
                 ctrl.Rule(humidity['--'] & TC_size['Huge+'], TC_size_variation['Quite Low']),
                 ctrl.Rule(humidity['-'] & TC_size['Little+'], TC_size_variation['Low']),
                 ctrl.Rule(humidity['-'] & TC_size['Little'], TC_size_variation['Quite Low']),
                 ctrl.Rule(humidity['-'] & TC_size['Huge'], TC_size_variation['Quite Low']),
                 ctrl.Rule(humidity['-'] & TC_size['Huge+'], TC_size_variation['Quite Low']),
                 ctrl.Rule(humidity['+'] & TC_size['Little+'], TC_size_variation['Quite High']),
                 ctrl.Rule(humidity['+'] & TC_size['Little'], TC_size_variation['Quite High']),
                 ctrl.Rule(humidity['+'] & TC_size['Huge'], TC_size_variation['Quite High']),
                 ctrl.Rule(humidity['+'] & TC_size['Huge+'], TC_size_variation['Quite High']),
                 ctrl.Rule(humidity['++'] & TC_size['Little+'], TC_size_variation['High']),
                 ctrl.Rule(humidity['++'] & TC_size['Little'], TC_size_variation['Very High']),
                 ctrl.Rule(humidity['++'] & TC_size['Huge'], TC_size_variation['Very High']),
                 ctrl.Rule(humidity['++'] & TC_size['Huge+'], TC_size_variation['Very Low'])]

    SF2_ctrl = ctrl.ControlSystem(SF2_rules)
    SF2 = ctrl.ControlSystemSimulation(SF2_ctrl)

    input_TC_size_fuzz = input_fuzz(TC_size, input_TC_size)
    input_humidity_fuzz = input_fuzz(humidity, input_humidity)

    IRR = IRR_2var(SF2_rules, input_humidity_fuzz, input_TC_size_fuzz)
    print(IRR)

    declenchement = IRR.min(axis=1)
    print(declenchement)
    csq_SF2 = {}  # variable de type dictionnaire
    for rule_i in range(0, len(SF2_rules)):  # on initialise à 0 pour chaque règle
        csq_SF2[str(SF2_rules[rule_i].consequent)] = 0
    for rule_i in range(0, len(SF2_rules)):
        csq_i = str(SF2_rules[rule_i].consequent)
        csq_SF2[csq_i] = max(csq_SF2[csq_i], declenchement[rule_i])  # Comme dans le TP6

    print("Conséquences SF2 :", csq_SF2)

    # SF2.input['humidity'] = input_humidity
    # SF2.input['TC_size'] = input_TC_size

    # print(SF2.output['TC_var'])
    # TC_size_variation.view(sim=SF2)


def movingaverage(interval, window_size):
    window = np.ones(int(window_size)) / float(window_size)
    return np.convolve(interval, window, 'same')


def smooth_array(array, axis_array, smooth_parameter, axis):
    # Fortement inspiré de https://stackoverflow.com/questions/11352047/finding-moving-average-from-data-points-in-python
    # Objectif : Appliquer une convolution de manière à lisser les données
    # L'élément numéro x prendra la valeur de la moyenne
    if axis == 0:  # latitude
        for row in range(0, len(axis_array)):
            array[row, :] = movingaverage(array[row, :], smooth_parameter)
    elif axis == 1:  # longitude
        for col in range(0, len(axis_array)):
            array[:, col] = movingaverage(array[:, col], smooth_parameter)
    return array


def psi_compute(u, v, smooth_parameter):
    # Intégrales de U et V pour déterminer psi, à l'instant t
    # Cf.Formule
    intx = scipy.integrate.cumtrapz(v, lon, axis=1, initial=0)
    inty = scipy.integrate.cumtrapz(u, lat, axis=0, initial=0)

    return (inty - intx) / 2

    # return 0.5 * (psi1 + psi2)


def cart_to_polar(x0, y0, x1, y1):
    # Convertit des coordonnées cartésiennes (x,y) en polaires (r,theta), dans le repère (o,x0,y0)
    # Cette fonction fonctionne aussi avec des matrices
    x0 = met_to_deg(x0)
    y0 = met_to_deg(y0)
    x1 = met_to_deg(x1)
    y1 = met_to_deg(y1)
    x_ref = x1 - x0
    y_ref = y1 - y0
    x_ref, y_ref = np.meshgrid(x_ref, y_ref)
    r = np.sqrt(deg_to_met(x_ref) ** 2 + deg_to_met(y_ref) ** 2)
    angle = np.arctan2(x_ref, y_ref) - np.radians(90)

    return r, angle


def array_r1_compute(r, r_map, angle_map, array):
    study = np.isclose(r_map, r, atol=50000)  # On prend tous les points à une distance r du centre de r1 ; tol en m
    # angle_psi = np.where(study, angle_map,0) # On prend les valeurs d'angle à une distance r du TC, on met 0 aux autres

    array_r = np.where(study, array,
                       np.nan)  # On va stocker suivant theta = [0,2pi] les différentes valeurs du champ sur le rayon r

    return array_r


def psi_TC_compute(x0, y0, study_univers, psi_t):
    r_map, angle_map = cart_to_polar(x0, y0, lon, lat)
    psi_TC_r = []  # Liste : psi(r)
    psi_TC_map = np.zeros(np.shape(angle_map))  # Map : psi(x;y)
    k = 0
    psi_inf = - np.mean(psi_t)
    compteur = 0
    stable = False
    # On se place ds le centre du TC, et on balaie selon x pour construire psi(r)
    for x1 in study_univers:
        y1 = y0
        r1, angle1 = cart_to_polar(x0, y0, x1, y1)

        psi_r1 = array_r1_compute(r1, r_map, angle_map, psi_t)

        if not stable:
            to_add = np.nanmean(psi_r1) + psi_inf
        else:
            to_add = psi_TC_r[k]
        if compteur < 4:
            psi_TC_r.append(to_add)
        else:
            psi_TC_r.append((to_add + sum(psi_TC_r[-3:])) / 4)
        # là où il n'y a pas de NaN dans psi_r1 on met mean(psi_r1)+psi_inf, sinon on met psi_TC_map (= on change rien)
        mask = np.invert(np.isnan(psi_r1))
        psi_TC_map = np.where(mask, to_add, psi_TC_map)

        # Test de Stabilité
        if compteur > 30:
            val_past = psi_TC_r[-30] / psi_TC_r[0]
            val_now = to_add / psi_TC_r[0]
            if abs(val_past - val_now) < 0.015:
                stable = True
                k = compteur

        # Si Stable, alors on fixe la valeur à 0
        compteur += 1

    return psi_TC_r, psi_TC_map


def xy_to_indx_lonlat(x, y, lon_map, lat_map):
    # On prend l'index ou (x,y)=(lon,lat), avec une tolérance de 25 km,
    # sachant que l'écart entre 2 valeurs de (lon,lat) est d'env 25 km aussi
    idx_lon = np.where(np.isclose(lon_map, x, atol=25000))[-1][-1]  # sachant que la résolution est de env 25000 km
    idx_lat = np.where(np.isclose(lat_map, y, atol=25000))[-1][-1]
    return idx_lon, idx_lat


def met_to_deg(x):
    return x / 111000


def deg_to_met(x):
    return x * 111000


def position_TC(jour_TC, heure_TC, path, range_grid):
    nc_TC = NetCDFFile(path)
    date = jour_TC * 24 + heure_TC

    u = uwind[date, :, :]
    v = vwind[date, :, :]
    vor = vorticity[date, :, :]

    x = scipy.integrate.cumtrapz(v, lon, axis=1, initial=0)
    y = scipy.integrate.cumtrapz(u, lat, axis=0, initial=0)
    psi_vor = (y - x)
    psi_vor = psi_vor / np.min(psi_vor) + vor / np.min(vor)

    latitude = np.array(nc_TC.variables['latitude'][:])
    longitude = np.array(nc_TC.variables['longitude'][:])
    # On trouve la vorticité max à cette date, dans une certaine zone
    psi_vor = np.where(range_grid, psi_vor, float('nan'))
    a = np.unravel_index(np.nanargmin(psi_vor), psi_vor.shape)

    # Puis on envoie les coordonnées
    # print(f"(Lon,Lat) = {latitude[a[0]], longitude[a[1]]}")

    nc_TC.close()

    return latitude[a[0]], longitude[a[1]]


def taille_TC(jour_TC, heure_TC, seuil, path):
    nc_TC = NetCDFFile(path)
    date_exact = jour_TC * 24 + heure_TC

    # np.shape(vorticity) = (744,93,117)
    # np.shape(latitude) = (93,)
    # np.shape(longitude) = (117,)

    vo = np.array(nc_TC.variables['vo'][:])

    # print(longitude[80]) gives -74.0

    # Limit longitude study up to -74.0 to the right to ignore high vorticity where there is no cyclone; can be
    # observed on the .gif

    vo_max = np.unravel_index(np.argmax(vo[date_exact, :, 0:74]), vo[date_exact, :, 0:74].shape)
    lon_TC = vo_max[1]
    lat_TC = vo_max[0]

    s = seuil
    vo[vo < s] = 0

    # Calculating min and max lon and lat

    # Hovering the latitude, fixing longitude of current position

    lat_vo_min = 0  # Initialising values
    lat_vo_max = 0  # Initialising values
    for j in range(93):
        if vo[
            date_exact, j, lon_TC] > 0:  # Taking the first non zero value index as the start point of the TC latitude
            lat_vo_min = j
            break

    if lat_vo_min == 0:  # If lat_vo_min==0 it means that all values were below s
        lat_vo_max = 0
    else:
        for j in range(lat_vo_min, 93):
            if vo[
                date_exact, j, lon_TC] == 0:  # Taking the last non zero value index as the end point of the TC latitude
                lat_vo_max = j - 1
                break

    # Hovering the longitude, fixing latitude of current position

    lon_vo_min = 0  # Initialising values
    lon_vo_max = 0  # Initialising values
    for j in range(117):
        if vo[
            date_exact, lat_TC, j] > 0:  # Taking the first non zero valueindex as the start point of the TC longitude
            lon_vo_min = j
            break

    if lon_vo_min == 0:
        lon_vo_max = 0
    else:
        for j in range(lon_vo_min, 93):
            if vo[
                date_exact, lat_TC, j] == 0:  # Taking the last non zero value index as the end point of the TC longitude
                lon_vo_max = j - 1
                break

    # print(f"LAT : Min:{lat_vo_min},Max:{lat_vo_max}")
    # print(f"LON : Min{lon_vo_min},Max:{lon_vo_max}")

    # ---------------------------------------------------------------------------------

    lat1, lon1, lat2, lon2, R = lat_vo_min, lon_vo_min, lat_vo_max, lon_vo_max, 6373.0

    # Distance between lon_vo_min and lon_vo_max, lat_TC constant
    coordinates_from = [lat_TC, lon1]
    coordinates_to = [lat_TC, lon2]

    distance_geopy = distance.distance(coordinates_from, coordinates_to).km
    distance_geopy_great_circle = distance.great_circle(coordinates_from, coordinates_to).km

    longueur_longitudinale = (distance_geopy + distance_geopy_great_circle) / 2
    # print('Longueur longitudinale', longueur_longitudinale)

    # Distance between lat_vo_min and lat_vo_max, lon_TC constant

    coordinates_from = [lat1, lon_TC]
    coordinates_to = [lat2, lon_TC]

    distance_geopy = distance.distance(coordinates_from, coordinates_to).km
    distance_geopy_great_circle = distance.great_circle(coordinates_from, coordinates_to).km

    longueur_latitudinale = (distance_geopy + distance_geopy_great_circle) / 2

    # print('Longueur latitudinale', longueur_latitudinale)

    # print('AIRE', longueur_latitudinale*longueur_longitudinale,'M2')

    nc_TC.close()
    return longueur_latitudinale * longueur_longitudinale


def parametre_coriolis(latitude):
    # Cf.Formule
    latitude = np.radians(met_to_deg(latitude))

    return 2 * 0.72921 * (10 ** (-4)) * np.sin(latitude)


def L_compute(psi_TC_on_psi0, univers):
    # Nous allons faire une régression sur les premiers éléments de la courbe, selon la fonction exponentielle
    # Ensuite nous allons déterminer l'abscisse (donc L) comme dit dans le rapport
    # La régression s'arrêtera à l'abscisse ou la dérivée seconde de la courbe devient négative, c.-à-d. quand
    # la dérivée première atteint son maximale.
    psi_grad = np.gradient(psi_TC_on_psi0)
    study_r = np.argmax(psi_grad)
    # Pour faire une régression exponentielle, on passe par la fonction log, cela donnera une fonction affine
    # On additionne donc PSI à son min ne pas avoir de valeurs négatives.

    mymodel = np.poly1d(
        np.polyfit(univers[0:study_r], np.log(psi_TC_on_psi0[0:study_r] - min(psi_TC_on_psi0) + 0.01), 1))

    myline = np.arange(0, 5000000, univers[1] - univers[0])

    L = univers[
        np.where(np.isclose(np.exp(mymodel(myline)), psi_TC_on_psi0[-50] - min(psi_TC_on_psi0) + 0.01, atol=1))[-1][
            -1]]

    # myline = np.arange(0, L, -pas) ; Pour l'affichage

    return L


def K_compute(psi_TC_r, univers):
    ############## 1) TROUVER K #################
    # La fonction de lissage requiert un passage en 2D ; ce passage n'est que temporaire

    psi_TC_r = smooth_array(np.expand_dims(psi_TC_r, axis=1), [1], smooth_parameter=10, axis=1)
    psi_TC_r = np.squeeze(psi_TC_r)
    # Le lissage entraîne du "bruit" aux valeurs extrêmes, que l'on retire
    remove_ind = np.arange(-5, 5, 1)
    psi_TC_r = np.delete(psi_TC_r, remove_ind)
    univers = np.delete(univers, remove_ind)

    psi_TC_on_psi0 = [-x / psi_TC_r[0] for x in psi_TC_r]
    if psi_TC_on_psi0[5] > psi_TC_on_psi0[50]:  # Si la courbe descend, on la met à l'opposée
        psi_TC_on_psi0 = [-x for x in psi_TC_on_psi0]
    # L'univers, selon x, est négatif. On le passe positif pour faciliter la régression, nécessaire pour trouver K
    univers = univers - min(univers)
    L = L_compute(psi_TC_on_psi0, univers)

    return (L ** 2) / 4


def get_fut_pos(x0_TC, y0_TC, psi_t, t_1h, var):
    # Plusieurs manipulations vont nous permettre de déterminer les différentes variables nécessaires
    # au calcul de la vélocité du TC
    # Il sera aussi nécessaire de lisser les données à plusieurs reprises
    # ---------------------------------------------------------------------
    # Détermination de PSI_TC, à l'aide de la formule
    # PSI_TC est utile pour : 1) Trouver K ; 2) Déterminer PSI_LARGE?
    # Nous allons le déterminer avec les coordonnées polaires
    # L'idée est de longer selon l'axe x, et de calculer PSI_TC sur le cercle de rayon x?

    # L'Univers est pour définir l'espace de travail de x
    pas = 3000
    univers = np.arange(x0_TC, max(lon) + 2000000,
                        pas)  # On va au delà de lon, nécessaire car on travaille avec des cercles, faut sortir du carré

    # On veut PSI_TC selon R (une liste), et PSI_TC selon (x,y) (un champ 2D)
    psi_TC_r, psi_TC_xy = psi_TC_compute(x0_TC, y0_TC, univers, psi_t)

    K = K_compute(psi_TC_r, univers)

    ####### 2) PSI_LARGE ########
    psi_large = (psi_t - psi_TC_xy)

    # PSI_lARGE est nécessaire pour :
    # 1) Déterminer U et V étant le gradient de PSI_LARGE
    # 2) Déterminer la vorticité absolue n (il faut aussi le paramètre de Coriolis f)
    f = np.expand_dims(lat, axis=1) * np.ones(np.shape(psi_large))
    f = parametre_coriolis(f)

    delta = lat[0] - lat[1]  # Résolution de nos données
    n = scipy.ndimage.laplace(psi_large) / (delta ** 2) + 1 * f

    u_TC = np.gradient(psi_TC_xy, lat, axis=0)
    v_TC = -np.gradient(psi_TC_xy, lon, axis=1)

    v_large = vwind - v_TC  # axis 1 = dérivée selon les colonnes
    u_large = uwind - u_TC

    n_grad_y = np.array(np.gradient(n, np.flip(lat), axis=0))
    n_grad_x = np.array(np.gradient(n, lon, axis=1))

    # On lisse les données
    smooth = 5

    n_grad_y = smooth_array(n_grad_y, lon, smooth, axis=1)
    n_grad_x = smooth_array(n_grad_x, lat, smooth, axis=0)

    # On a tous les paramètres nécessaires au calcul de la vélocité, et donc de la position future du cyclone
    # La vélocité dépend des champs n et PSI. Ces derniers n'étant pas uniformes, la vélocité n'est pas la même partout
    # Il est donc nécessaire d'étudier la vélocité à petit pas, avec un pas de temps dt
    # Soit t_1 (heures) l'instant futur du cyclone
    # Nous travaillons en unités SI

    t_1s = t_1h * 60 * 60
    dts = 0.5 * 60 * 60  # Soit 1 heure

    C = 0.7156
    nm = 1.2830

    x_TC_out = x0_TC  # On prépare la sortie
    y_TC_out = y0_TC

    a = 1.5

    # Idée pour le flou : Ici r_TC c'est la taille max ; Mais ça sera un nb flou.
    # On trace la moyenne de U,V,N en fonction du rayon
    # Pour chaque valeur du rayon on multiplie par le degré de vérité de r_TC
    r_TC = 400000
    r_min = 0  # Pour éviter de prendre les valeurs centrales, qui sont pas 100% fiables
    # Pour chaque pas de temps dt, allant de l'instant t0 (maintenant) à l'instant t_1

    for k in range(0, t_1s, int(np.round(dts))):

        r_map, angle_map = cart_to_polar(x_TC_out, y_TC_out, lon, lat)
        u_large_r1 = []
        v_large_r1 = []
        n_grad_y_r1 = []
        n_grad_x_r1 = []

        r_univers = np.arange(r_min, r_TC, 10000)

        for r in r_univers:
            u_large_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, u_large)))
            v_large_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, v_large)))
            n_grad_y_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, n_grad_y)))
            n_grad_x_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, n_grad_x)))

        # FLOU
        # Ces lignes là seront à changer pour le flou
        u = np.mean(u_large_r1)
        v = np.mean(v_large_r1)
        ny = np.mean(n_grad_y)
        nx = np.mean(n_grad_x)
        #

        TC_lat = np.radians(met_to_deg(y_TC_out))  # Détermination de la latitude du cyclone, en rad
        m = (C * (np.cos(TC_lat)) ** (-1)) * (np.tan(0.25 * (np.pi - 2 * TC_lat))) ** nm

        # idx_x, idx_y = xy_to_indx_lonlat(x_TC_out, y_TC_out, lon, lat)  # On récupère les index selon [lon,lat] du TC

        Cx0 = a * u - K * ny  # en m/s
        Cy0 = a * v + K * nx  # en m/s

        x_TC_out = x_TC_out + Cx0 * dts * m
        y_TC_out = y_TC_out + Cy0 * dts * m

    want_to_plot = False
    if want_to_plot:
        plt.subplot(2, 2, 1)
        # plot_array(psi, lon, lat)
        # plt.plot(met_to_deg(x_TC), met_to_deg(y_TC), ms=4, marker='o', markeredgecolor="green")
        plt.subplot(2, 2, 2)
        # plt.plot(univers,psi_TC_on_psi0-min(psi_TC_on_psi0)+0.01)
        # plt.plot(univers,np.gradient(psi_TC_on_psi0,univers))
        # plt.plot(univers, smooth_array(np.expand_dims(np.gradient(psi_TC_on_psi0),1),[1],2,axis=1))

        plot_array(psi_TC_xy, lon, lat)

        plt.subplot(2, 2, 3)
        # plt.plot(myline, np.exp(mymodel(myline)))
        # plt.plot(univers, psi_TC_on_psi0 - min(psi_TC_on_psi0) + 0.01)
        # plot_array(vorticity[t, :, :], lon, lat)
        # plt.plot(met_to_deg(x_TC), met_to_deg(y_TC), ms=10, marker='o', markeredgecolor="green")
        # plt.streamplot(met_to_deg(lon),np.flip(met_to_deg(lat)),u_large,v_large )
        plt.subplot(2, 2, 4)
        plot_array(n_grad_y, lon, lat)
        plt.plot(met_to_deg(x0_TC), met_to_deg(y0_TC), ms=10, marker='o', markeredgecolor="green")
        plt.plot(met_to_deg(x_TC_out), met_to_deg(y_TC_out), ms=10, marker='o', markeredgecolor="red")
        y_v, x_v = position_TC(jour, heure + 12, data_path)
        plt.plot(x_v, y_v, ms=10, marker='x', markeredgecolor="yellow")

        # plt.figure() #Le Cyclone circule le long des lignes
        # plt.streamplot(met_to_deg(lon), met_to_deg(np.flip(lat)), np.flip(u_large, 0), np.flip(v_large, 0),
        #              density=5,linewidth=0.3,color=np.flip(psi_large,axis=0))
        # plt.streamplot(met_to_deg(lon), met_to_deg(np.flip(lat)), np.flip(n_grad_x, 0), np.flip(n_grad_y, 0),
        #               density=4, linewidth=0.3, color=np.flip(n, axis=0))

    return x_TC_out, y_TC_out


def optimisation():
    # Fonction vite fait pour le fun, où on teste des paramètres et ça nous affiche l'erreur
    # Y a qqs paramètres dans les CAF que j'ai pu affiner comme ça
    jour_1 = 3
    h_1 = 0
    t1 = t0 + 24 * jour_1 + h_1  # en heures
    dt = 6  # ttes les dt heures

    for s in range(80, 150, 10):
        for s2 in [2, 10]:
            instant = 0
            instant_s = 0
            x_TC_fut = 0 * np.arange(t0, t1 + dt, dt)
            y_TC_fut = 0 * np.arange(t0, t1 + dt, dt)

            y_TC, x_TC = position_TC(jour, heure, data_path)
            x_TC = deg_to_met(x_TC)
            y_TC = deg_to_met(y_TC)

            x_TC_fut[0] = x_TC
            y_TC_fut[0] = y_TC

            err = np.zeros(np.shape(x_TC_fut))
            err[0] = 0
            for ti in range(t0, t1, dt):
                psi = psi_compute(uwind, vwind, ti, 8)

                y_TC, x_TC = position_TC(0, ti, data_path)
                x_TC, y_TC = deg_to_met(x_TC), deg_to_met(y_TC)

                err[instant] = np.sqrt(((x_TC_fut[instant] - x_TC) ** 2) + ((y_TC_fut[instant] - y_TC) ** 2))

                x_TC_fut[instant + 1], y_TC_fut[instant + 1] = get_fut_pos(x_TC, y_TC,
                                                                           psi, 24,
                                                                           [s,
                                                                            s2])  # prévision de la position future du cyclone
                instant += 1

            err = met_to_deg(err)
            print("s =", s, "s2= ", s2, "MEAN", np.mean(err), "; MAX", np.max(err), "; MIN", np.min(err[1:]))
            instant_s += 1

    plt.subplot(2, 1, 1)
    plt.plot(met_to_deg(x_TC_fut), met_to_deg(y_TC_fut), ms=3, marker='x', markeredgecolor="red")
    plot_array(psi, lon, lat)

    plt.subplot(2, 1, 2)
    plt.plot(err)


def traj_pred(x0, y0, psi, t1, dt):
    # Prédiction d'une trajectoire, à partir de l'emplacement actuel du TC, et des champs
    # Ressort la prédiction au temps t0+t1
    # Les points sont espacés d'un temps dt

    time_step_traj = np.arange(0, t1, dt)
    x_TC_fut = np.zeros(np.shape(time_step_traj))
    y_TC_fut = np.zeros(np.shape(time_step_traj))
    x_TC_fut[0] = x0
    y_TC_fut[0] = y0

    instant = 0
    for temps in time_step_traj:
        x_TC_fut[instant], y_TC_fut[instant] = get_fut_pos(x0, y0, psi, temps, 0)
        instant += 1

    return x_TC_fut, y_TC_fut


def posTC_compute(t0, t_end, want_to_plot):
    # Ressort la position du TC en fonction du temps sur notre temps d'étude [t0,t_end]
    y_TC = [0 for t in range(0, t0, 1)]  # Avant t0, la position du cyclone n'est pas définie
    x_TC = [0 for t in range(0, t0, 1)]
    # Nous allons utiliser les champs PSI et de Vorticité pour trouver l'emplacement du cyclone
    # A leur valeur maximale se trouve le cyclone. Cependant cela n'est pas toujours vrai, car la zone d'étude est grande,
    # d'autres perturbations peuvent s'apparenter un à un cyclone.
    # L'idée est d'utiliser une meshgrid pour afficher la recherche du point, à partir de l'emplacement du point
    # à l'instant précédent
    range_grid = np.ones(np.shape(np.meshgrid(lon, lat)))[0, :, :]  # Grille pour filtrer la recherche du xy_TC

    if want_to_plot:
        plot_array(vorticity[t0 + 30, :, :], lon, lat)

    for t in range(t0, t_end, 1):
        #
        y_TC.append(deg_to_met(position_TC(0, t, data_path, range_grid)[0]))
        x_TC.append(deg_to_met(position_TC(0, t, data_path, range_grid)[1]))

        # Mise à jour de la grille de recherche
        range_lon = np.isclose(met_to_deg(lon), met_to_deg(x_TC[-1]), atol=2)
        range_lat = np.isclose(met_to_deg(lat), met_to_deg(y_TC[-1]), atol=2)

        range_grid = np.meshgrid(range_lon, range_lat)
        # Nous souhaitons les emplacements où les valeurs en lon ET en lat sont True
        range_grid = range_grid[0] * range_grid[1]

        if want_to_plot:
            plt.plot(met_to_deg(x_TC[-1]), met_to_deg(y_TC[-1]), ms=5, marker='x', markeredgecolor='red')
            plt.show(block=False)
            plt.pause(0.001)

    return x_TC, y_TC  # En degrés


def xy_pred(x0, y0, u_large, v_large, n_grad_x, n_grad_y, r_TC, t_1s, dts):

    x_TC_out = x0  # On prépare la sortie
    y_TC_out = y0
    # A partir des champs et de la position du TC, nous pouvons prédire sa trajectoire pour dans un temps t_1s
    # La vitesse varie dans le champ, nous ferons donc des "petits pas" de temps dts, pour arriver à l'instant final
    # Retourne la position finale (x;y)
    # Cependant il n'est pas intéressant de prendre les valeurs de vitesse en 1 point : on s'intéresse à une zone
    # de r mètres autour du cyclone. On prendra la vitesse moyenne dans cette zone
    # _____________
    # Divers paramètres pour le calcul :
    r_min = 20000
    C = 0.7156
    nm = 1.2830
    a = 1.5

    for k in range(0, t_1s, int(np.round(dts))):

        r_map, angle_map = cart_to_polar(x_TC_out, y_TC_out, lon, lat)

        u_large_r1 = []
        v_large_r1 = []
        n_grad_y_r1 = []
        n_grad_x_r1 = []

        r_univers = np.arange(r_min, r_TC, 10000)

        for r in r_univers:
            u_large_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, u_large)))
            v_large_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, v_large)))
            n_grad_y_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, n_grad_y)))
            n_grad_x_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, n_grad_x)))

        # FLOU
        # Ces lignes là seront à changer pour le flou
        u = np.mean(u_large_r1)
        v = np.mean(v_large_r1)
        ny = np.mean(n_grad_y)
        nx = np.mean(n_grad_x)
        #

        TC_lat = np.radians(met_to_deg(y_TC_out))  # Détermination de la latitude du cyclone, en rad
        # Calcul du paramètre m. Celui-ci corrige 'l'erreur' lié au fait que nous travaillons sur une sphère, et donc
        # que les distances sont courbées comparées à une étude sur un plan 2D
        m = (C * (np.cos(TC_lat)) ** (-1)) * (np.tan(0.25 * (np.pi - 2 * TC_lat))) ** nm

        # idx_x, idx_y = xy_to_indx_lonlat(x_TC_out, y_TC_out, lon, lat)  # On récupère les index selon [lon,lat] du TC

        Cx0 = a * u  # - K * ny  # Vélocité selon x au centre du cyclone, en m/s
        Cy0 = a * v  # + K * nx  # Vélocité selon y au centre du cyclone, en m/s

        # On actualise la position
        x_TC_out = x_TC_out + Cx0 * dts * m
        y_TC_out = y_TC_out + Cy0 * dts * m

    return x_TC_out, y_TC_out


def fuzzy_gradient(fuzz_array, delta, axis):  # TRAPEZES, donc matrice 3D de dim(4,lat,lon)
    # Nos nombres flous sont de la forme [a,b,c,d]
    # Ils sont semblables à des matrices classiques en 2D, sauf qu'il y a une 3è dimensions pour [a,b,c,d]
    # Méthode des différences finies, mais en flou
    # L'intérêt de créer une fonction à part est de gérer l'inversion lorsque que l'on soustrait un IFT à un autre
    # Rappel : A-B = (a1-b4,a2-b3,a3-b2,a4-b1)
    # Pour cela on crée directement une matrice B = np.flip(A)=(a4,a3,a2,a1)
    # Ainsi on peut directement exécuter A-B pour calculer la dérivée

    grad = np.zeros(np.shape(fuzz_array))
    fuzz_array_B = np.flip(fuzz_array, axis=2)

    if axis == 1:
        for col in range(0, np.shape(fuzz_array)[0] - 1):
            for row in range(0, np.shape(fuzz_array)[1] - 1):
                if col == 0:
                    grad[0, row, :] = (fuzz_array[0, row + 1, :] - fuzz_array_B[0, row, :]) / delta
                elif col == np.shape(fuzz_array)[0] - 1:
                    grad[-1, row, :] = (fuzz_array[-1, row, :] - fuzz_array_B[-1, row - 1, :]) / delta
                else:
                    grad[col, row, :] = (fuzz_array[col, row + 1, :] - fuzz_array_B[col, row, :]) / (delta)

    if axis == 0:
        for col in range(0, np.shape(fuzz_array)[0] - 1):
            for row in range(0, np.shape(fuzz_array)[1] - 1):
                if row == 0:
                    grad[col, 0, :] = (fuzz_array[col + 1, 0, :] - fuzz_array_B[col, 0, :]) / delta
                elif row == np.shape(fuzz_array)[1] - 1:
                    grad[col, -1, :] = (fuzz_array[col, -1, :] - fuzz_array_B[col - 1, -1, :]) / delta
                else:
                    grad[col, row, :] = (fuzz_array[col + 1, row, :] - fuzz_array_B[col, row, :]) / (delta)
    return grad


def fuzz_CAF1(uwind, vwind, x_TC, y_TC):
    # Fuzzification des matrices, avec des NFT
    # Soit x le nombre net, alors cela revient à l'ajout d'une dimension à la matrice, de valeur : [x+a,x+b,x+c,x+d]
    a, b, c, d = [-0.1, 0, 0, 0.1]

    uwind = uwind[:, :, np.newaxis] # Extension de la matrice en 3D
    uwind = np.dstack((uwind + a, uwind + b, uwind + c, uwind + d))  # La matrice est maintenant fuzzifiée

    a, b, c, d = [-0.1, 0, 0, 0.1]
    vwind = vwind[:, :, np.newaxis]
    vwind = np.dstack((vwind + a, vwind + b, vwind + c, vwind + d))  # La matrice est maintenant fuzzifiée

    dNFT = [0, 0, 0, 0]
    x_TC = np.array(x_TC + dNFT)
    y_TC = np.array(y_TC + dNFT)
    return uwind, vwind, x_TC, y_TC


def CAF1_compute(uwind, vwind, x_TC, y_TC):
    # Calcul PSI, PSI_TC
    # ---------
    # PSI : Ce dernier est calculé à partir d'une intégrale, donc d'une somme et d'une division par un nombre net
    # U et V sont des IFT de la forme (a,b,c,d)
    # Nous pouvons calculer l'intégrale pour a,b,c et d, 1 à 1

    psi = np.zeros(np.shape(uwind))
    for i in range(0, 4):
        psi[:, :, i] = psi_compute(uwind[:, :, i], vwind[:, :, i], 6)
    # ----------
    # PSI_TC : Comme au dessus, PSI_TC se base sur une moyenne, donc par le calcul d'additions et d'une division par un nombre net
    # Nous pouvons calculer directement PSI_TC pour les 4 termes qui composent notre IFT
    # Nous allons le déterminer avec les coordonnées polaires
    # L'idée est de longer selon l'axe x, et de calculer PSI_TC sur le cercle de rayon x

    res = 3000
    psi_TC_xy = np.zeros(np.shape(psi))
    psi_TC_r = np.zeros((res, 4))

    K = np.zeros(4)
    for i in range(0, 4):
        univers = np.linspace(x_TC[i], max(lon) + 2500000, res)
        # On va au delà de lon, nécessaire car on travaille avec des cercles, faut sortir du carré

        # L'Univers est pour définir l'espace de travail de x
        # On veut PSI_TC selon R (une liste), et PSI_TC selon (x,y) (un champ 2D)
        psi_TC_r[:, i], psi_TC_xy[:, :, i] = psi_TC_compute(x_TC[i], y_TC[i], univers, psi[:, :, i])
        K[i] = K_compute(psi_TC_r[:, i], univers)

    # Le paramètre K nous sera aussi nécessaire

    psi_large = (psi - np.flip(psi_TC_xy))

    return psi, psi_large, psi_TC_xy, K


def CAF2_compute(psi, psi_TC, x_TC, y_TC, r_TC, temps_etude):
    f = np.expand_dims(lat, axis=1) * np.ones(np.shape(psi[:, :, 0]))
    f = parametre_coriolis(f)
    f = np.dstack(([f, f, f, f]))

    delta = lat[1] - lat[0]  # Résolution de nos données
    n = scipy.ndimage.laplace(psi_large) / (delta ** 2) + 1 * f

    u_TC = fuzzy_gradient(psi_TC, delta=lat[1] - lat[0], axis=0)
    v_TC = -np.flip(fuzzy_gradient(psi_TC, delta=lon[1] - lon[0], axis=1))

    v_large = vwind[:, :, :] - np.flip(v_TC, axis=2)  # La négation du nombre flou demande le np.flip()
    u_large = uwind[:, :, :] - np.flip(u_TC, axis=2)

    n_grad_y = np.array(np.gradient(n, np.flip(lat), axis=0))
    n_grad_x = np.array(np.gradient(n, lon, axis=1))

    plot_array(u_TC[:, :, 0], lon, lat)
    plt.show()
    plot_array(u_TC[:, :, 1], lon, lat)
    plt.show()
    plot_array(u_TC[:, :, 2], lon, lat)
    plt.show()

    # On lisse les données
    # smooth = 5

    # n_grad_y = smooth_array(n_grad_y, lon, smooth, axis=1)
    # n_grad_x = smooth_array(n_grad_x, lat, smooth, axis=0)

    # On a tous les paramètres nécessaires au calcul de la vélocité, et donc de la position future du cyclone
    # La vélocité dépend des champs n et PSI. Ces derniers n'étant pas uniformes, la vélocité n'est pas la même partout
    # Il est donc nécessaire d'étudier la vélocité à petit pas, avec un pas de temps dt
    # Soit t_1 (heures) l'instant futur du cyclone
    # Nous travaillons en unités SI

    t_1s = temps_etude * 60 * 60
    dts = 0.5 * 60 * 60  # Soit 1 heure

    # Idée pour le flou : Ici r_TC c'est la taille max ; Mais ça sera un nb flou.
    # On trace la moyenne de U,V,N en fonction du rayon
    # Pour chaque valeur du rayon on multiplie par le degré de vérité de r_TC
    # Pour chaque pas de temps dt, allant de l'instant t0 (maintenant) à l'instant t_1
    x = np.zeros(4)
    y = np.zeros(4)
    for i in range(0, 4):
        x[i], y[i] = xy_pred(x_TC[i], y_TC[i], u_large[:, :, i], v_large[:, :, i], n_grad_x[:, :, i], n_grad_y[:, :, i],
                             r_TC, t_1s, dts)
    print(x, y)
    return x, y

# Récupération des données
data_path = 'D:/KATRINA.nc' # Données selon la pression, pour la vitesse
nc = NetCDFFile(data_path)
uwind = np.mean(nc.variables['u'][:, (1, 3, 6, 9), :, :], axis=1)  # Moyenne de la vitesse sur les différents
vwind = np.mean(nc.variables['v'][:, (1, 3, 6, 9), :, :], axis=1)  # niveaux de pression, pour une meilleure précision
nc.close()

data_path = 'D:/TT/UV/SY10/KATRINA.nc' # Données à niveau de pression précis
nc = NetCDFFile(data_path)
lat = deg_to_met(nc.variables['latitude'][:])  # En degré, qu'on convertit en mètre
lon = deg_to_met(nc.variables['longitude'][:])
time = nc.variables['time'][:]
# uwind = nc.variables['u'][:]
# vwind = nc.variables['v'][:]
vorticity = nc.variables['vo'][:]
nc.close()

# Les calculs se font en unités SI (mètres, secondes,...). Pour l'affichage, nous convertissons cependant l'abscisse
# et l'ordonnée en mètre
# plt.streamplot(met_to_deg(lon), met_to_deg(np.flip(lat)), np.flip(uwind[t,:,:], 0), np.flip(vwind[t,:,:], 0),
#               density=5,linewidth=0.3,color=np.flip(psi, axis=0))

# Définition de l'instant 0
if data_path == 'D:/TT/UV/SY10/WILMA.nc':
    jour = 16
    heure = 18
    t0 = jour * 24 + heure
    tmax = t0 + 9 * 24
if data_path == 'D:/TT/UV/SY10/FERN.nc':
    jour = 0
    heure = 0
    t0 = jour * 24 + heure
    tmax = t0 + 4 * 24
if data_path == 'D:/TT/UV/SY10/KATRINA.nc':
    jour = 24
    heure = 0
    t0 = jour * 24 + heure
    tmax = t0 + 7 * 24

x_TC, y_TC = posTC_compute(t0, tmax, want_to_plot=False)  # En mètres
uwind = uwind[t0, :, :]
vwind = vwind[t0, :, :]

#_________________________
# Prédit toutes les dt_traj heures, la traj de (t0) à (t0+t1). Tous les points d'une trajectoire sont espacés d'un
# temps dt
study = True
dt_traj = 12 # Temps entre chaque trajectoire
t1 = 72  # Durée d'une trajectoire
dt = 12  # Durée entre les points d'une trajectoire
if study:
    ci = 0
    for t0 in range(t0, tmax - 24, dt_traj):  # Le domaine temporel d'étude
        print(tmax - t0)
        psi = psi_compute(uwind, vwind, 6)

        # Prediction de (xy) dans un temps t1 à partir de l'instant t0
        x, y = traj_pred(x_TC[t0], y_TC[t0], psi, t1, dt)

        # Affichage
        c = ['yellow', 'red', 'blue', 'green', 'purple', 'brown', 'red']
        plt.plot(met_to_deg(x), met_to_deg(y), ms=5, marker='x', markeredgecolor=c[ci % 7])
        plt.plot(met_to_deg(x_TC[t0]), met_to_deg(y_TC[t0]), ms=7, marker='o', markeredgecolor='black')
        ci += 1
        # plt.show(block=False) #Permet l'affichage pendant que le programme fonctionne
        # plt.pause(0.001)


#________________________
fuzz = False
if fuzz:
    # Fuzzification uwind vwind
    # Créer une dimension pour les coordonnées des nombres trapézoidaux, de forme [a,b,c,d]
    uwind, vwind, x_TC, y_TC = fuzz_CAF1(uwind, vwind, x_TC[t0], y_TC[t0])

    psi, psi_large, psi_TC, K = CAF1_compute(uwind, vwind, x_TC, y_TC)

    x, y = CAF2_compute(psi, psi_TC, x_TC, y_TC, r_TC=300000, temps_etude=t1)

    plot_array(psi[:, :, 0], lon, lat)
    plt.plot(met_to_deg(x), met_to_deg(y), ms=5, marker='x', markeredgecolor='red')
    plt.show()

# SF1_compute(298, -36.5, 8.1) # temperature ; d_speed ; âge
# SF2_compute(54, 224)  # humidité ; taille du cyclone

plt.show()