Untitled

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 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

x_TC, y_TC = posTC_compute(t0, tmax, want_to_plot=False)  # En mètres

jour = 24
heure = 0
t0 = jour * 24 + heure
tmax = t0 + 7 * 24