main.py

mail@pastecode.io avatar
unknown
plain_text
2 years ago
28 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):
    nc = NetCDFFile(path)
    lat = nc.variables['latitude'][:]
    lon = nc.variables['longitude'][:]
    time = nc.variables['time'][:]
    var = nc.variables[var][:]
    nc.close()
    print(var)
    map = Basemap(width=5000000, height=3500000,
                  resolution='l', projection='cyl',
                  llcrnrlon=lon.min(), llcrnrlat=lat.min(), urcrnrlon=lon.max(), urcrnrlat=lat.max(), lat_0=lat.mean(),
                  lon_0=lon.mean())
    lons, lats = np.meshgrid(lon, lat)
    xi, yi = map(lons, lats)

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

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

    jour = 26
    heure = 10
    date_exact = jour * 24 + heure
    cs = map.pcolor(xi, yi, np.squeeze(var[date_exact, :, :]))
    cbar = map.colorbar(cs, location='bottom', pad="10%")
    cbar.set_label(unite)
    plt.title(titre)


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


def plot_array(array, x, y):
    lon = met_to_deg(x)
    lat = met_to_deg(y)

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

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


def input_fuzz(variable, input_to_fuzz):
    fuzz_input = []
    for i in variable.terms.keys():
        fuzz_input.append(fuzz.interp_membership(variable.universe, variable.terms[i].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)')
    time = 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])

    time['Young'] = fuzz.zmf(time.universe, 3, 6)
    time['Middle'] = fuzz.gaussmf(time.universe, 6, 1)
    time['Old'] = fuzz.smf(time.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()
    time.view()

    SF1_rules = [ctrl.Rule(temperature['Low'] & d_speed['--'] & time['Young'], lifespan['Dying+']),
                 ctrl.Rule(temperature['Low'] & d_speed['--'] & time['Middle'], lifespan['Dying+']),
                 ctrl.Rule(temperature['Low'] & d_speed['--'] & time['Old'], lifespan['Dying+']),
                 ctrl.Rule(temperature['Low'] & d_speed['-'] & time['Young'], lifespan['Dying']),
                 ctrl.Rule(temperature['Low'] & d_speed['-'] & time['Middle'], lifespan['Dying+']),
                 ctrl.Rule(temperature['Low'] & d_speed['-'] & time['Old'], lifespan['Dying+']),
                 ctrl.Rule(temperature['Low'] & d_speed['St'] & time['Young'], lifespan['Dying']),
                 ctrl.Rule(temperature['Low'] & d_speed['St'] & time['Middle'], lifespan['Dying']),
                 ctrl.Rule(temperature['Low'] & d_speed['St'] & time['Old'], lifespan['Dying']),
                 ctrl.Rule(temperature['Low'] & d_speed['+'] & time['Young'], lifespan['St']),
                 ctrl.Rule(temperature['Low'] & d_speed['+'] & time['Middle'], lifespan['St']),
                 ctrl.Rule(temperature['Low'] & d_speed['+'] & time['Old'], lifespan['St']),
                 ctrl.Rule(temperature['Low'] & d_speed['++'] & time['Young'], lifespan['Living']),
                 ctrl.Rule(temperature['Low'] & d_speed['++'] & time['Middle'], lifespan['Living']),
                 ctrl.Rule(temperature['Low'] & d_speed['++'] & time['Old'], lifespan['Living']),
                 ctrl.Rule(temperature['Ok'] & d_speed['--'] & time['Young'], lifespan['Dying']),
                 ctrl.Rule(temperature['Ok'] & d_speed['--'] & time['Middle'], lifespan['Dying']),
                 ctrl.Rule(temperature['Ok'] & d_speed['--'] & time['Old'], lifespan['Dying']),
                 ctrl.Rule(temperature['Ok'] & d_speed['-'] & time['Young'], lifespan['Dying']),
                 ctrl.Rule(temperature['Ok'] & d_speed['-'] & time['Middle'], lifespan['Dying']),
                 ctrl.Rule(temperature['Ok'] & d_speed['-'] & time['Old'], lifespan['Dying']),
                 ctrl.Rule(temperature['Ok'] & d_speed['St'] & time['Young'], lifespan['Living']),
                 ctrl.Rule(temperature['Ok'] & d_speed['St'] & time['Middle'], lifespan['St']),
                 ctrl.Rule(temperature['Ok'] & d_speed['St'] & time['Old'], lifespan['Dying']),
                 ctrl.Rule(temperature['Ok'] & d_speed['+'] & time['Young'], lifespan['Living']),
                 ctrl.Rule(temperature['Ok'] & d_speed['+'] & time['Middle'], lifespan['Living']),
                 ctrl.Rule(temperature['Ok'] & d_speed['+'] & time['Old'], lifespan['Living']),
                 ctrl.Rule(temperature['Ok'] & d_speed['++'] & time['Young'], lifespan['Living+']),
                 ctrl.Rule(temperature['Ok'] & d_speed['++'] & time['Middle'], lifespan['Living+']),
                 ctrl.Rule(temperature['Ok'] & d_speed['++'] & time['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(time, 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 i in range(0, len(SF1_rules)):  # on initialise à 0 pour chaque règle
        csq_SF1[str(SF1_rules[i].consequent)] = 0
    for i in range(0, len(SF1_rules)):
        csq_i = str(SF1_rules[i].consequent)
        csq_SF1[csq_i] = max(csq_SF1[csq_i], declenchement[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 i in range(0, len(SF2_rules)):  # on initialise à 0 pour chaque règle
        csq_SF2[str(SF2_rules[i].consequent)] = 0
    for i in range(0, len(SF2_rules)):
        csq_i = str(SF2_rules[i].consequent)
        csq_SF2[csq_i] = max(csq_SF2[csq_i], declenchement[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)


# https://stackoverflow.com/questions/11352047/finding-moving-average-from-data-points-in-python

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):
    if axis == 0:  # latitude
        for i in range(0, len(axis_array)):
            array[i, :] = movingaverage(array[i, :], smooth_parameter)
    elif axis == 1:  # longitude
        for i in range(0, len(axis_array)):
            array[:, i] = movingaverage(array[:, i], smooth_parameter)
    return array


def psi_compute(uwind, vwind, t):
    uwind = np.squeeze(uwind[t, :, :])
    vwind = np.squeeze(vwind[t, :, :])  # var[temps,latitude,longitude]

    smooth = 5
    intx = scipy.integrate.cumtrapz(vwind, lon, axis=1, initial=0)[0]
    inty = scipy.integrate.cumtrapz(uwind, lat, axis=0, initial=0)

    intx_av = np.zeros(np.shape(intx))
    inty_av = np.zeros(np.shape(inty))
    for i in range(0, len(lat)):
        inty_av[i, :] = movingaverage(inty[i, :], smooth)

    psi1 = intx - inty_av

    intx = scipy.integrate.cumtrapz(vwind, lon, axis=1, initial=0)
    inty = scipy.integrate.cumtrapz(uwind, lat, axis=0, initial=0)[:, 0][:, None]

    intx_av = np.zeros(np.shape(intx))

    for i in range(0, len(lon)):
        intx_av[:, i] = movingaverage(intx[:, i], smooth)

    psi2 = intx_av - inty
    psi = 0.5 * (psi1 + psi2)
    plot_array(psi,lon,lat)
    return psi


def cart_to_polar(x0, y0, x1, y1):  # fonctionne 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 psi_r1_compute(r, r_map, angle_map, psi):
    study = np.isclose(r_map, r, atol=20000)  # On prend tous les points à une distance r du TC ; tol en m
    # Revoir la tolérance
    angle_psi = np.where(study, angle_map,
                         0)  # On prend les valeurs d'angle à une distance r du TC, on met 0 aux autres
    psi_r = np.where(study, psi,
                     np.nan)  # On va stocker suivant theta = [0,2pi] les différentes valeurs de psi sur le rayon r
    # plot_array(psi_r,lon,lat)

    return psi_r


def psi_TC_compute(x0, y0, study_univers):
    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 = 0
    i = 0  # Compteur
    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 = psi_r1_compute(r1, r_map, angle_map, psi)
        to_add = np.nanmean(psi_r1) + psi_inf
        if not stable:
            to_add = np.nanmean(psi_r1) + psi_inf
        else:
            to_add = psi_TC_r[k]

        psi_TC_r.append(to_add)

        # 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 i > 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 = i

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

    return psi_TC_r, psi_TC_map


def xy_to_indx_lonlat(x, y, lon, lat):
    idx_lon = np.where(np.isclose(lon, x, atol=25000))[-1]  # sachant que la résolution est de env 25000 km
    idx_lat = np.where(np.isclose(lat, y, atol=25000))[-1]  # au pire on a deux éléments
    if len(idx_lon) > 1:
        idx_lon = np.squeeze(idx_lon)[-1]
    if len(idx_lat) > 1:
        idx_lat = np.squeeze(idx_lat)[-1]
    return idx_lon, idx_lat


def met_to_deg(x):  # à préciser
    return x / 111000


def deg_to_met(x):  # à préciser
    return x * 111000


def position_TC(jour, heure, path):
    nc = NetCDFFile(path)
    date_exact = jour * 24 + heure
    x = scipy.integrate.cumtrapz(vwind[date_exact, :, :], lon, axis=1, initial=0)
    y = scipy.integrate.cumtrapz(uwind[date_exact, :, :], lat, axis=0, initial=0)
    psi_faux = y - x

    latitude = np.array(nc.variables['latitude'][:])
    longitude = np.array(nc.variables['longitude'][:])

    # On trouve la vorticité max à cette date
    a = np.unravel_index(np.argmax(psi_faux), psi_faux.shape)

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


def taille_TC(jour, heure, seuil, path):
    nc = NetCDFFile(path)
    date_exact = jour * 24 + heure

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

    vorticity = np.array(nc.variables['vo'][:])
    latitude = np.array(nc.variables['latitude'][:])
    longitude = np.array(nc.variables['longitude'][:])

    # 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(vorticity[date_exact, :, 0:74]), vorticity[date_exact, :, 0:74].shape)
    lon_TC = vo_max[1]
    lat_TC = vo_max[0]

    s = seuil
    vorticity[vorticity < 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 i in range(93):
        if vorticity[
            date_exact, i, lon_TC] > 0:  # Taking the first non zero value index as the start point of the TC latitude
            lat_vo_min = i
            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 i in range(lat_vo_min, 93):
            if vorticity[
                date_exact, i, lon_TC] == 0:  # Taking the last non zero value index as the end point of the TC latitude
                lat_vo_max = i - 1
                break

    # Hovering the longitude, fixing latitude of current position

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

    if lon_vo_min == 0:
        lon_vo_max = 0
    else:
        for i in range(lon_vo_min, 93):
            if vorticity[
                date_exact, lat_TC, i] == 0:  # Taking the last non zero value index as the end point of the TC longitude
                lon_vo_max = i - 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')

    return longueur_latitudinale * longueur_longitudinale


def parametre_coriolis(latitude):
    latitude = np.radians(met_to_deg(latitude))

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


def get_fut_pos(x_TC, y_TC, psi):
    # On dépasse de +20, nécessaire car on balaie selon le rayon, faut sortir du carré
    univers = np.arange(x_TC, max(lon) + 2000000, 5000)

    psi_TC_r, psi_TC_xy = psi_TC_compute(x_TC, y_TC, univers)  # psi selon r, et psi selon x;y

    psi_TC_on_psi0 = [-x / psi_TC_r[0] for x in psi_TC_r]
    # TROUVER K
    univers = (univers - min(univers))
    study_r = int(np.round(0.05*int(np.where(np.isclose(max(psi_TC_on_psi0),psi_TC_on_psi0))[-1][-1])))

    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)) #faire avec exp


    myline = np.linspace(0, 700000, 500)
    L = myline[np.where(np.isclose(mymodel(myline), psi_TC_on_psi0[-1], atol=0.2))[-1][-1]]

    K = (L ** 2) / 4
    myline = np.linspace(0, L, 500)
    psi_large = psi - psi_TC_xy

    f = np.expand_dims(lat, axis=1) * np.ones(np.shape(psi_large))  # paramètre de Coriolis
    f = parametre_coriolis(f)

    delta = lat[0] - lat[1]  # le même mais négatif en longitude

    n = scipy.ndimage.laplace(psi_large) / (delta ** 2) + f

    # axis 0 = dérivée selon les lignes ; axis 1 = dérivée selon les colonnes

    # d'après leur définitino le - n'est pas là, mais jsp pk il le faut pour qu'on est le bon sesn du vent
    v_large = np.array(np.gradient(psi_large, lon, axis=1))
    u_large = np.array(-1) * np.array(np.gradient(psi_large, lat, axis=0))

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

    smooth = 5
    u_large = smooth_array(u_large, lon, smooth, axis=1)
    v_large = smooth_array(v_large, lat, smooth, axis=0)
    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)

    a = 1
    # Pour dans x heures :
    t_1 = 12 * 60 * 60
    dt = 1 * 60 * 60

    C = 0.7156
    nm = 1.2830

    x_TC_fut = x_TC
    y_TC_fut = y_TC

    # A Faire : un Mask avec le TC, et faire une moyenne de u/v

    for i in range(0, t_1, int(np.round(dt))):

        TC_lat = np.radians(met_to_deg(y_TC_fut))
        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_fut, y_TC_fut, lon, lat)
        Cx0 = a * u_large[idx_y, idx_x] - K * n_grad_y[idx_y, idx_x]  # en m/s, normalement
        Cy0 = a * v_large[idx_y, idx_x] + K * n_grad_x[idx_y, idx_x]  # en m/s, normalement

        x_TC_fut = x_TC_fut + Cx0 * dt#* m
        y_TC_fut = y_TC_fut + Cy0 * dt#* m


        print("dx", Cx0 * dt, a * u_large[idx_y, idx_x], - K * n_grad_y[idx_y, idx_x])
        print("dy", Cy0 * dt, a * v_large[idx_y, idx_x], K * n_grad_x[idx_y, idx_x],m)



    #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(myline, np.exp(mymodel(myline)))
    #plt.plot(univers,psi_TC_on_psi0-min(psi_TC_on_psi0)+0.01)
    #plt.plot(univers,np.grandient(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)
    #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(x_TC), met_to_deg(y_TC), ms=10, marker='o', markeredgecolor="green")
    #plt.plot(met_to_deg(x_TC_fut), met_to_deg(y_TC_fut), 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_fut, y_TC_fut


jour = 28  # Commence au jour zéro
heure = 2
t = jour * 24 + heure

data_path = 'D:/TT/UV/SY10/KATRINA.nc'

nc = NetCDFFile(data_path)
lat = nc.variables['latitude'][:]
lon = nc.variables['longitude'][:]
time = nc.variables['time'][:]
uwind = nc.variables['u'][:]
vwind = nc.variables['v'][:]
vorticity = nc.variables['vo'][:]

lon = deg_to_met(lon)
lat = deg_to_met(lat)

lons, lats = np.meshgrid(lon, lat)

nc.close()

psi = psi_compute(uwind, vwind, t)  # psi c'est la stream function

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

# plot_array(psi, lon, lat)
# plot_array(psi, lon, lat)
# 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))

x_TC_fut, y_TC_fut = get_fut_pos(x_TC, y_TC, psi)  # prévision de la position future du cyclone

# x_TC_fut2, y_TC_fut2 = get_fut_pos(x_TC_fut, y_TC_fut, K, psi)

# plt.plot(x_TC_fut2, y_TC_fut2, ms=10, marker='o', markeredgecolor="red")


# plt.plot(univers+np.array(-x_TC), fit_k)


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

plt.show()