Fonctions autres
unknown
plain_text
2 years ago
35 kB
4
Indexable
import numpy as np
import matplotlib.pyplot as plt
import scipy
from netCDF4 import Dataset as NetCDFFile
from mpl_toolkits.basemap import Basemap
from scipy.optimize import curve_fit
from geopy import distance
TC_name = "KATRINA"
path = 'D:TT/UV/SY10/'
# Nos données sont divisées en 3 fichiers
# -- 1 FICHIER DEEP MEAN LAYER (DML) : La vitesse à chaque niveau de pression.
# Nous ferons une moyenne de cette vitesse sur tous les niveaux de pression.
# Travailler sur plusieurs niveaux de pression est plus efficace.
# -- 1 FICHIER PRESSION (P) : Il s'agit de la pression à une altitude égale à la surface de la mer
# -- 1 FICHIER STANDARD : Contenant le reste des informations nécessaires (température et vorticité),
# ici à un niveau de pression précis.
# Nous aurions pu aussi choisir un fichier DML pour la pression et vorticité de manière à travailler sur
# tous les niveaux de pression, mais cela demande une trop grande capacité de stockage.
# Chaque variable est en 3D, 1 pour le temps, 1 pour la longitude (axe x) et 1 pour la latitude (axe y)
# L'indice du temps est en Heures, et commence à t=0. Une matrice comporte des données allant du 1er au Dernier jour
# du mois
# Pour accéder à une matrice au 2è jours du mois à 5h du matin, on pose donc t = (2-1) * 24 + 5 = 29
# On pose t0 l'instant t où le cyclone apparaît, et tmax l'instant t où le cyclone n'est plus.
if TC_name == "WILMA":
data_path_DML = path + 'WILMA_DML.nc' # Données selon la pression, pour la vitesse
data_path = path + "WILMA.nc" # Données à niveau de pression précis
pressure_path = path + "WILA_P.n" # Pression au niveau de la mer
jour = 17 # Jour d'apparition du cyclone
heure = 18 # Heure d'apparition
t_TC0 = (jour - 1) * 24 + heure # Instant t d'apparition
tmax = t_TC0 + 9 * 24 # Instant t de disparition
if TC_name == 'FERN':
data_path_DML = path + "FERN_DML.nc"
data_path = path + "FERN.nc"
pressure_path = path + "FERN_P.nc"
jour = 1
heure = 0
t_TC0 = (jour - 1) * 24 + heure
tmax = t_TC0 + 4 * 24
if TC_name == 'KATRINA':
data_path_DML = path + "KATRINA_DML.nc"
data_path = path + "KATRINA.nc"
pressure_path = path + "KATRINA_P.nc"
jour = 25
heure = 14
t_TC0 = (jour - 1) * 24 + heure
tmax = t_TC0 + 6 * 24 - 13
if TC_name == 'CHARLEY':
data_path_DML = path + "CHARLEY_DML.nc"
data_path = path + "CHARLEY.nc"
pressure_path = path + "CHARLEY_P.nc"
jour = 10
heure = 0
t_TC0 = (jour - 1) * 24 + heure
tmax = t_TC0 + 4 * 24
def met_to_deg(x): # Conversion d'une valeur métrique en degrés
return x / 111000
def deg_to_met(x): # Conversion d'une valeur degrés en mètres
return x * 111000
def parametre_coriolis(latitude): # Calcul du paramètre de Coriolis pour latitude donnée
# Cf.Formule
latitude = np.radians(met_to_deg(latitude))
return 2 * 0.72921 * (10 ** (-4)) * np.sin(latitude)
def plotNC(path, var, unite, titre, date_plot):
# Affiche une variable quelconque d'un fichier .nc
nc_plot = NetCDFFile(path)
lat_plot = nc_plot.variables['latitude'][:]
lon_plot = nc_plot.variables['longitude'][:]
var = nc_plot.variables[var][:]
nc_plot.close()
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):
# Affiche toutes les informations nécessaires en une fois
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):
# Affiche une matrice 2D selon la longitude x et la latitude y
# Les coordonnées sont en degrés
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 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):
# Fonction calculant 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 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):
# 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
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=30000) # 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, is_fuzzy):
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_TC sera soustrait par psi_inf. En flou, la soustraction nécessite un np.flip, que l'on traitera après.
# Dans un premier temps nous ne calculerons donc que la moyenne
# En net, on peut effectuer la soustraction directement.
if is_fuzzy:
psi_inf = 0
else:
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 < 100:
psi_TC_r.append(to_add)
else:
psi_TC_r.append((to_add + sum(psi_TC_r[-1:])) / 2)
# 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 > 100:
val_past = psi_TC_r[-10] / psi_TC_r[0]
val_now = to_add / psi_TC_r[0]
if abs(val_past - val_now) < 0.005:
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):
# Permet à partir d'une position (x,y) d'obtenir les indices des matrices (lon_map),(lat_map) correspondantes
# 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 position_TC(t, range_grid):
u = uwind[t, :, :]
v = vwind[t, :, :]
vor = vorticity[t, :, :]
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)
# 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]]}")
return lat[a[0]], lon[a[1]]
def taille_TC(t, seuil):
lon_TC = idx_x_TC[t]
lat_TC = idx_y_TC[t]
vo = np.where(seuil,vorticity[t, :, :],0)
vo[vo < 0.3 * np.max(vo)] = 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(len(lat)):
if vo[
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, len(lat)):
if vo[
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(len(lon)):
if vo[
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, len(lon)):
if vo[
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[lat_vo_min], lon[lon_vo_min], lat[lat_vo_max], lon[lon_vo_max], 6373.0
# Maintenant on travaille avec des degrés, et non plus des index
lat1 = met_to_deg(lat1)
lon1 = met_to_deg(lon1)
lat2 = met_to_deg(lat2)
lon2 = met_to_deg(lon2)
lat_TC = met_to_deg(y_TC[t])
lon_TC = met_to_deg(x_TC[t])
# 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)/2
def get_fut_pos(x0_TC, y0_TC, psi_t, t_1h, var, uwind, vwind):
# 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, is_fuzzy=False)
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) + 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 = 50000
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, 15000)
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(t0):
# 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(30000, 1000000, 10000):
for s2 in [0.2, 0.5]:
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)
x_TC_fut[0] = x_TC[t0]
y_TC_fut[0] = y_TC[t0]
err = np.zeros(np.shape(x_TC_fut))
err[0] = 0
for ti in range(t0, t1, dt):
psi = psi_compute(uwind[t0, :, :], vwind[t0, :, :])
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, u, v):
# 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:
print(temps)
x_TC_fut[instant], y_TC_fut[instant] = get_fut_pos(x0, y0, psi, temps, 0, u, v)
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)]
idx_x_TC = [0 for t in range(0, t0, 1)]
idx_y_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):
#
xy_pred = position_TC(t, range_grid)
x_pred = xy_pred[1]
y_pred = xy_pred[0]
idx_pred = xy_to_indx_lonlat(x_pred, y_pred, lon, lat)
y_TC.append(xy_pred[0])
x_TC.append(xy_pred[1])
idx_x_TC.append(idx_pred[0])
idx_y_TC.append(idx_pred[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, idx_x_TC, idx_y_TC # En Mètres
# Calcule le module de la vitesse en tout point à tout moment
def calcul_pth(t,study_grid):
pth = np.sqrt(np.power(uwind[t,:,:], 2) + np.power(vwind[t,:,:], 2))
pth = np.where(study_grid,pth,0)
return (pth)
# Donne en sortie la vitesse du TC (utilise la matrice créée juste avant avec calcul_pth)
def vitesse_TC(t, study_grid):
res = np.max(calcul_pth(t,study_grid))
return (res) # Pythagore
# Donne la sortie la variation de vitesse du TC, à distance d'une heure)
def variation_vitesse_TC(t,study_grid):
rep = vitesse_TC(t, calcul_pth(t,study_grid)) - vitesse_TC(t-12,calcul_pth(t,study_grid))
return (rep)
def relation_pression_taille_TC(t, study_grid,idx_x, idx_y):
t_list = []
taille_list = []
pressure_list = []
for i in range(t_TC0,tmax):
t_list.append(i)
pressure_list.append(pressure[i, idx_y, idx_x]) # METTRE AU NIVEAU DU TC
taille_list.append(taille_TC(i, study_grid))
fig, ax1 = plt.subplots()
color = 'tab:red'
ax1.set_xlabel('Temps')
ax1.set_ylabel('Pression (Pa)', color=color)
ax1.plot(t_list, pressure_list, color=color)
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis
color = 'tab:blue'
ax2.set_ylabel('Taille (Km)', color=color) # we already handled the x-label with ax1
ax2.plot(t_list, taille_list, color=color)
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout() # otherwise the right y-label is slightly clipped
plt.show()
def relation_pression_vitesse_TC(t, study_grid,idx_x, idx_y):
t_list = []
taille_list = []
pressure_list = []
for i in range(t_TC0,tmax):
t_list.append(i)
pressure_list.append(pressure[i, idx_y, idx_x])
taille_list.append(vitesse_TC(i, study_grid))
fig, ax1 = plt.subplots()
color = 'tab:red'
ax1.set_xlabel('Temps')
ax1.set_ylabel('Pression (Pa)', color=color)
ax1.plot(t_list, pressure_list, color=color)
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis
color = 'tab:blue'
ax2.set_ylabel('Vitesse (Km/h)', color=color) # we already handled the x-label with ax1
ax2.plot(t_list, taille_list, color=color)
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout() # otherwise the right y-label is slightly clipped
plt.show()
def get_angle(x1, y1, x0, y0):
# Entrée en mètre
r, angle = cart_to_polar(x0, y0, x1, y1)
angle = np.degrees(angle[0][0])
# plot_array(u[date_exact,:,:],lon,lat)
# plt.plot(met_to_deg(x0), met_to_deg(y0), ms=5, marker='o', markeredgecolor='red')
# plt.plot(met_to_deg(x1), met_to_deg(y1), ms=8, marker='x', markeredgecolor='green')
# plt.show()
return np.round(angle, decimals=1)
def regression_trajectoire_TC(t,want_to_plot):
# Données entrainement
x = [met_to_deg(x_TC[i]) for i in range(t_TC0,t)]
y = [met_to_deg(y_TC[i]) for i in range(t_TC0,t)]
if x[-1] > x[-5]: #Si le cyclone va vers la droite
sens = 1
else:
sens = -1
x_pred = [met_to_deg(lon[i]) for i in range(0,len(lon))]
plot_array(vorticity[t, :, :], lon, lat)
def reg_poly(x, y):
def objective(x, a, b, c):
return a * (x ** 2) + b * x + c
params = np.polyfit(x,y,2)
a = params[0]
b = params[1]
c = params[2]
y_test = []
for pt in x_pred:
y_test.append(objective(pt, a, b, c))
x1 = x_pred[idx_x_TC[t]+(sens*15)]
y1 = objective(x_pred[idx_x_TC[t]+(sens*15)],a,b,c)
x0 = met_to_deg(x_TC[t])
y0 = met_to_deg(y_TC[t])
angle_poly = get_angle(x1,y1,x0,y0) + 180
if want_to_plot:
y_test = np.array(y_test)
plt.scatter(x_pred, y_test)
plt.scatter(x1,y1)
plt.show()
erreur = []
for i in range(0,len(x)):
e = (y[i] - objective(x[i],a,b,c)) ** 2
erreur.append(np.sqrt(e))
print(f"Erreur poly en moyenne:{np.mean(erreur)}")
return np.mean(erreur),angle_poly
def reg_lin(x, y):
def objective(x, a, b):
return a * x + b
params = curve_fit(objective, x, y)
a, b = params[0]
y_test = []
for pt in x_pred:
y_test.append(objective(pt, a, b))
x1 = x_pred[idx_x_TC[t]+(sens*15)]
y1 = objective(x_pred[idx_x_TC[t]+(sens*15)],a,b)
x0 = met_to_deg(x_TC[t])
y0 = met_to_deg(y_TC[t])
angle_lin = get_angle(x1,y1,x0,y0) + 180
if want_to_plot:
plt.scatter(x_pred, y_test)
plt.scatter(x1, y1)
plt.show()
erreur = []
for i in range(0,len(x)):
e = (y[i] - objective(x[i],a,b)) ** 2
erreur.append(np.sqrt(e))
print(f"Erreur lin en moyenne:{np.mean(erreur)}")
return np.mean(erreur),angle_lin
err_poly,angle_poly = reg_poly(x,y)
err_lin, angle_lin = reg_lin(x,y)
if err_poly < err_lin:
return angle_poly
else:
return angle_lin
# Convertir une position polaire en position cartésienne
def polar_to_cart(r, angle):
angle = np.radians(angle)
x = r * np.cos(angle)
y = r * np.sin(angle)
return x, y
# 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
# Estimation en Net
def study():
t0 = t_TC0
dt_traj = 24 # Temps entre chaque trajectoire
t1 = 48 # Durée d'une trajectoire
dt = 40 # Durée entre les points d'une trajectoire
ci = 0
plot_array(vorticity[t0, :, :], lon, lat)
for t0 in range(t0, tmax - 72, dt_traj): # Le domaine temporel d'étude
print("Heures restantes:", tmax - t0)
psi = psi_compute(uwind[t0, :, :], vwind[t0, :, :])
# 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, uwind[t0, :, :], vwind[t0, :, :])
# 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)
nc_DML = NetCDFFile(data_path_DML)
# uwind = np.mean(nc_DML.variables['u'][:, (1, 3, 6, 9), :, :], axis=1) # Moyenne de la vitesse sur les différents
# vwind = np.mean(nc_DML.variables['v'][:, (1, 3, 6, 9), :, :], axis=1) # niveaux de pression, pour une meilleure précision
nc_DML.close()
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'][:]
temperature = nc.variables['t'][:] # Prendre la température à la surface de la mer??
humidity = nc.variables['r'][:]
uwind = nc.variables['u'][:]
vwind = nc.variables['v'][:]
vorticity = nc.variables['vo'][:]
nc.close()
# pressure = np.array(NetCDFFile(pressure_path).variables['sp'][:])
x_TC, y_TC, idx_x_TC, idx_y_TC = posTC_compute(t_TC0, tmax, want_to_plot=False) # En mètres
Editor is loading...