ERA5_storms_CPS_diagram.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Suivi de tempêtes dans les réanalyses ERA5 - Diagramme de phase\n",
    "\n",
    "**Auteur : Frédéric FERRY - Ecole Nationale de la Météorologie (Météo France)**\n",
    "\n",
    "Le diagramme de phase des cyclones (Cyclone Phase Space Diagram : moe.met.fsu.edu/~rhart/papers-hart/2003Hart.pdf) est un outil diagnostique utile permettant de caractériser la structure thermodynamique des cyclones tropicaux et extra-tropicaux. Il se base sur l'analyse de la symétrie thermique du système dans la basse troposphère et sur sa structure à cœur chaud ou à cœur froid dans la basse et la haute troposphère.\n",
    "\n",
    "Les paramètres utilisés pour construire ce diagnostic permettent de décrire la structure à cœur chaud d'un cyclone tropical, la structure à cœur froid d'un cyclone extra-tropical ainsi que ses stades de développement (formation, intensification, occlusion, affaiblissement, développement de séclusions chaudes). Cet outil permet aussi d'analyser la structure hybride des cyclones subtropicaux et de suivre les transtions tropicales de cyclones subtropicaux ou transitions extra-tropicales de cyclones tropicaux.\n",
    "\n",
    "Le présent calepin permet :\n",
    "- Le suivi d'une tempête à partir de réanalyses ERA5 de Pmer et de géopotentiel et d'un fichier .txt de tracking obtenu grâce à un programme de tracking automatique (non fourni ici)\n",
    "- Le calcul pas à pas des paramètres de symétrie thermique et de cœur chaud/cœur froid du diagrammme de phase\n",
    "- Le tracé du diagramme de phase tel que disponible sur http://moe.met.fsu.edu/cyclonephase/\n",
    "\n",
    "Le diagramme de phase est ici illustré avec le cas de la transition extra-tropicale du cyclone tropical Ophelia (octobre 2017)\n",
    "\n",
    "Pour faire fonctionner ce calepin il est nécessaire de télécharger les données ERA5 au format netcdf, au pas de temps horaire et sur domaine limité (-100W-50E, 0-90N) :\n",
    "- https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=form\n",
    "- https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-pressure-levels?tab=form"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import xarray as xr\n",
    "import netCDF4\n",
    "\n",
    "import math\n",
    "import numpy as np\n",
    "from scipy.ndimage import gaussian_filter\n",
    "\n",
    "from cartopy import config\n",
    "from cartopy.util import add_cyclic_point\n",
    "import cartopy.feature as cfeature\n",
    "import cartopy.crs as ccrs\n",
    "from cartopy.mpl.geoaxes import GeoAxes\n",
    "from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter\n",
    "import cartopy.io.img_tiles as cimgt\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import AxesGrid\n",
    "import matplotlib.collections as collections\n",
    "\n",
    "import IPython.display as IPdisplay, matplotlib.font_manager as fm\n",
    "from PIL import Image\n",
    "import glob\n",
    "\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "from tqdm import tqdm\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Fonctions pour les graphiques."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_maxmin_points(lon, lat, data, extrema, nsize, symbol, color='k',\n",
    "                       plotValue=True, transform=None):\n",
    "\n",
    "    from scipy.ndimage import maximum_filter, minimum_filter\n",
    "\n",
    "    if (extrema == 'max'):\n",
    "        data_ext = maximum_filter(data, nsize, mode='nearest')\n",
    "    elif (extrema == 'min'):\n",
    "        data_ext = minimum_filter(data, nsize, mode='nearest')\n",
    "    else:\n",
    "        raise ValueError('Value for hilo must be either max or min')\n",
    "\n",
    "    mxy, mxx = np.where(data_ext == data)\n",
    "\n",
    "    for i in range(len(mxy)):\n",
    "        #print('Low '+str(i+1))\n",
    "        #print(lon[mxy[i], mxx[i]], lat[mxy[i], mxx[i]], int(data[mxy[i], mxx[i]]))\n",
    "        ax.text(lon[mxy[i], mxx[i]], lat[mxy[i], mxx[i]], symbol, color=color, size=12,\n",
    "                clip_on=True, horizontalalignment='center', verticalalignment='center',\n",
    "                transform=transform)\n",
    "        ax.text(lon[mxy[i], mxx[i]], lat[mxy[i], mxx[i]],\n",
    "                '\\n' + str(int(data[mxy[i], mxx[i]])),\n",
    "                color=color, size=10, clip_on=True, fontweight='bold',\n",
    "                horizontalalignment='center', verticalalignment='top', transform=transform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_animation(gif_filepath):\n",
    "    from PIL import Image\n",
    "    import os\n",
    "    from IPython.display import Image as IPImage\n",
    "    from IPython.display import display\n",
    "    import time\n",
    "    \n",
    "    image_folder = './anim/'+storm+'/' # répertoire contenant les fichiers PNG\n",
    "    output_file = gif_filepath # nom du fichier de sortie\n",
    "    animation_speed = 0.9 # vitesse de l'animation en secondes\n",
    "    \n",
    "    # Liste tous les fichiers PNG dans le répertoire image_folder\n",
    "    files = sorted(os.listdir(image_folder))\n",
    "    image_files = [f for f in files if f.endswith('.png')]\n",
    "    \n",
    "    # Ouvre chaque fichier PNG et ajoute l'image à une liste\n",
    "    images = []\n",
    "    for filename in image_files:\n",
    "        img = Image.open(os.path.join(image_folder, filename))\n",
    "        images.append(img)\n",
    "    \n",
    "    # Crée l'animation GIF\n",
    "    images[0].save(output_file, save_all=True, append_images=images[1:], duration=int(animation_speed*1000), loop=0)\n",
    "    # Affiche l'animation GIF dans Jupyter\n",
    "    with open(output_file,'rb') as f:\n",
    "        display(IPImage(data=f.read(), format='png'))\n",
    "    # Efface les fichiers PNG\n",
    "    for filename in image_files:\n",
    "        os.remove(image_folder+filename)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Fonction pour le calcul de distances sur la sphère."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def haversine(lon1, lat1, lon2, lat2):\n",
    "    # convert decimal degrees to radians \n",
    "    lon1 = np.deg2rad(lon1)\n",
    "    lon2 = np.deg2rad(lon2)\n",
    "    lat1 = np.deg2rad(lat1)\n",
    "    lat2 = np.deg2rad(lat2)\n",
    "\n",
    "    # haversine formula \n",
    "    dlon = lon2 - lon1 \n",
    "    dlat = lat2 - lat1 \n",
    "    a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2\n",
    "    c = 2 * np.arcsin(np.sqrt(a)) \n",
    "    r = 6371\n",
    "    return c * r"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Fonction pour le calcul d'angles sur la sphère."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def get_bearing(lat1, lon1, lat2, lon2):\n",
    "    dLon = (lon2 - lon1)\n",
    "    x = math.cos(math.radians(lat2)) * math.sin(math.radians(dLon))\n",
    "    y = math.cos(math.radians(lat1)) * math.sin(math.radians(lat2)) - math.sin(math.radians(lat1)) * math.cos(math.radians(lat2)) * math.cos(math.radians(dLon))\n",
    "    brng = np.arctan2(x,y)\n",
    "    brng = np.degrees(brng)\n",
    "    brng = (brng +360) % 360\n",
    "    return brng"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Trajectoires de toutes les tempêtes"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Liste des tempêtes dont le fichier de tracking automatique est disponible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "storms=['Nov1982', 'Oct1987', 'Herta', 'Viviane', 'Braer', 'Lothar', 'Martin', 'Klaus', 'Xynthia',\n",
    "        'Joachim', 'Zeus', 'Ophelia', 'Eleanor', 'Zorbas', 'Lorenzo', 'Alex', 'Ciaran']"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Tracé des trajectoires de toutes les tempêtes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latS=25\n",
    "latN=70\n",
    "lonW=-75\n",
    "lonE=30\n",
    "bounds = [(lonW, lonE, latS, latN)]\n",
    "projection=ccrs.NearsidePerspective(central_longitude=(lonW+lonE)/2, central_latitude=(latS+latN)/2)\n",
    "                             \n",
    "fig = plt.figure(figsize=(15., 10.))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=projection)\n",
    "ax.set_title('A selection of storms tracked from ERA5 hourly MSLP data at ~30km resolution',loc='center',fontsize=14)\n",
    "ax.gridlines(draw_labels=True, color='gray', alpha=0.8, linestyle='-')\n",
    "#ax.stock_img()\n",
    "LAND = cfeature.NaturalEarthFeature('physical', 'land', '10m',edgecolor='face',facecolor=cfeature.COLORS['land'],linewidth=.1)\n",
    "ax.add_feature(LAND)\n",
    "OCEAN = cfeature.NaturalEarthFeature('physical', 'ocean', '10m',edgecolor='face',facecolor=cfeature.COLORS['water'],linewidth=.1)\n",
    "ax.add_feature(OCEAN)\n",
    "ax.set_extent(*bounds, crs=ccrs.PlateCarree())\n",
    "\n",
    "for storm in storms:\n",
    "    liste_time,liste_lon,liste_lat,liste_pres=np.loadtxt('./txt/'+storm+'.txt',skiprows=0, dtype='U13,float,float,int',unpack=True)\n",
    "    date1=str(liste_time[0])\n",
    "    date2=str(liste_time[-1])\n",
    "    ax.plot(liste_lon,liste_lat, transform=ccrs.PlateCarree(), label=storm+' '+date1+'-'+date2)\n",
    "    ax.text(liste_lon[0]+1, liste_lat[0], storm, verticalalignment='center', horizontalalignment='right',transform=ccrs.PlateCarree())\n",
    "\n",
    "ax.legend(loc=\"upper center\", ncol=3)\n",
    "plt.show()\n",
    "\n",
    "figname='./figs/tracks2'\n",
    "fig.savefig(figname+'.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Suivi d'une tempête"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Choix de la tempête."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_menu():\n",
    "    print ('1 -- Tempête de Novembre 1982' )\n",
    "    print (\"2 -- 'Ourangan' d'octobre 1987\")\n",
    "    print ('3 -- Herta - Février 1990' )\n",
    "    print ('4 -- Viviane - Février 1990' )\n",
    "    print ('5 -- Braer - Janvier 1993' )\n",
    "    print ('6 -- Lothar - Décembre 1999' )\n",
    "    print ('7 -- Martin - Décembre 1999' )\n",
    "    print ('8 -- Klaus - Janvier 2009' )\n",
    "    print ('9 -- Xynthia - Février 2010' )\n",
    "    print ('10 -- Joachim - Décembre 2011' )\n",
    "    print ('11 -- Zeus - Mars 2017' )\n",
    "    print ('12 -- Ophelia - Octobre 2017 (transition extra tropicale)' )\n",
    "    print ('13 -- Eleanor - Janvier 2018' )\n",
    "    print ('14 -- Zorbas - Septembre 2018 (medicane)')\n",
    "    print ('15 -- Lorenzo - Octobre 2019 (transition extra tropicale)' )\n",
    "    print ('16 -- Alex - Octobre 2020' )\n",
    "    print ('17 -- Ciaran - Octobre-novembre 2023' )\n",
    "\n",
    "print_menu()\n",
    "\n",
    "option = int(input('Enter number of the desired storm : ')) \n",
    "if option == 1:\n",
    "    storm='Nov1982'\n",
    "elif option == 2:\n",
    "    storm='Oct1987'\n",
    "elif option == 3:\n",
    "    storm='Herta'\n",
    "elif option == 4:\n",
    "    storm='Viviane'\n",
    "elif option == 5:\n",
    "    storm='Braer'\n",
    "elif option == 6:\n",
    "    storm='Lothar'\n",
    "elif option == 7:\n",
    "    storm='Martin'\n",
    "elif option == 8:\n",
    "    storm='Klaus'\n",
    "elif option == 9:\n",
    "    storm='Xynthia'\n",
    "elif option == 10:\n",
    "    storm='Joachim'\n",
    "elif option == 11:\n",
    "    storm='Zeus'\n",
    "elif option == 12:\n",
    "    storm='Ophelia'\n",
    "elif option == 13:\n",
    "    storm='Eleanor'\n",
    "elif option == 14:\n",
    "    storm='Zorbas'\n",
    "elif option == 15:\n",
    "    storm='Lorenzo'\n",
    "elif option == 16:\n",
    "    storm='Alex'\n",
    "elif option == 17:\n",
    "    storm='Ciaran'\n",
    "\n",
    "else:\n",
    "    print('Invalid option. Please enter a number between 1 and 17.')\n",
    "\n",
    "dir_figs='./figs/'+storm+'/CPS/'\n",
    "dir_anim='./anim/'+storm+'/'\n",
    "\n",
    "if not os.path.exists(dir_figs):\n",
    "    os.mkdir(dir_figs)\n",
    "    \n",
    "if not os.path.exists(dir_anim):\n",
    "    os.mkdir(dir_anim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latS=30\n",
    "latN=70\n",
    "lonW=-60\n",
    "lonE=30\n",
    "\n",
    "if storm=='Braer':\n",
    "    latS=40\n",
    "    latN=80\n",
    "    lonW=-80\n",
    "    lonE=20\n",
    "if storm=='Ciaran':\n",
    "    latS=30\n",
    "    latN=70\n",
    "    lonW=-70\n",
    "    lonE=30\n",
    "if storm=='Lorenzo':\n",
    "    latS=25\n",
    "    latN=60\n",
    "    lonW=-60\n",
    "    lonE=20\n",
    "if storm=='Zorbas':\n",
    "    latS=30\n",
    "    latN=45\n",
    "    lonW=5\n",
    "    lonE=35\n",
    "    \n",
    "bounds = [(lonW, lonE, latS, latN)]"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Ouverture du fichier texte de tracking de la tempête choisie."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "liste_time,liste_lon,liste_lat,liste_pres=np.loadtxt('./txt/'+storm+'.txt',\n",
    "        skiprows=0, dtype='U13,float,float,int',unpack=True)\n",
    "print(liste_time)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Changement du pas de temps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interval=int(input(\"Enter the desired time step (3h or 6h advised): \"))\n",
    "\n",
    "liste_time=liste_time[::interval]\n",
    "liste_lon=liste_lon[::interval]\n",
    "liste_lat=liste_lat[::interval]\n",
    "liste_pres=liste_pres[::interval]\n",
    "\n",
    "print(liste_time)\n",
    "print(liste_time.shape)\n",
    "\n",
    "date1=str(liste_time[0])\n",
    "date2=str(liste_time[-1])"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Tracé"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def plot_background(ax):\n",
    "    ax.coastlines(\"10m\", zorder=3, color='grey')\n",
    "    ax.gridlines()\n",
    "    ax.set_xticks(np.linspace(-180, 180, 37), crs=ccrs.PlateCarree())\n",
    "    ax.set_yticks(np.linspace(-90, 90, 37), crs=ccrs.PlateCarree())\n",
    "    lon_formatter = LongitudeFormatter(zero_direction_label=True)\n",
    "    lat_formatter = LatitudeFormatter()\n",
    "    ax.xaxis.set_major_formatter(lon_formatter)\n",
    "    ax.yaxis.set_major_formatter(lat_formatter)\n",
    "    return(ax)\n",
    "\n",
    "projection=ccrs.PlateCarree()\n",
    "\n",
    "vmin = np.min(liste_pres)\n",
    "vmax = np.max(liste_pres)\n",
    "\n",
    "fig = plt.figure(figsize=(15., 15.))\n",
    "ax = fig.add_subplot(111, projection=projection)\n",
    "ax.set_title(storm+' storm track from '+liste_time[0]+' to '+liste_time[-1],loc='center',fontsize=14)\n",
    "\n",
    "#ax.stock_img()\n",
    "LAND = cfeature.NaturalEarthFeature('physical', 'land', '10m',edgecolor='face',facecolor=cfeature.COLORS['land'],linewidth=.1)\n",
    "ax.add_feature(LAND)\n",
    "OCEAN = cfeature.NaturalEarthFeature('physical', 'ocean', '10m',edgecolor='face',facecolor=cfeature.COLORS['water'],linewidth=.1)\n",
    "ax.add_feature(OCEAN)\n",
    "# Create a Stamen terrain background instance.\n",
    "#stamen_terrain = cimgt.Stamen('terrain-background')\n",
    "# Add the Stamen data at zoom level 8.\n",
    "#ax.add_image(stamen_terrain, 8)\n",
    "plot_background(ax)\n",
    "ax.set_extent(*bounds, crs=ccrs.PlateCarree())\n",
    "\n",
    "ax.plot(liste_lon,liste_lat, color='k', transform=ccrs.PlateCarree())\n",
    "for i in range(len(liste_time)):\n",
    "    ax.scatter(liste_lon[i],liste_lat[i], c=liste_pres[i],\n",
    "               cmap='plasma', vmin=vmin, vmax=vmax, s=50,\n",
    "               label=liste_pres[i],transform=ccrs.PlateCarree(), zorder=10)\n",
    "\n",
    "plt.legend(ncol=2)\n",
    "plt.show()\n",
    "figname=dir_figs+storm+'_era_track'\n",
    "fig.savefig(figname+'.png')\n",
    "\n",
    "fig = plt.figure(figsize=(15., 8.))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_title(storm+' MSLP evolution from '+liste_time[0]+' to '+liste_time[-1],loc='center',fontsize=14)\n",
    "ax.scatter(liste_time,liste_pres, color='red')\n",
    "ax.plot(liste_time,liste_pres, color='red')\n",
    "\n",
    "plt.xlabel('Date', fontsize=12)\n",
    "plt.ylabel('Mean Sea Level Pressure (hPa)', fontsize=12)\n",
    "locs, labels = plt.xticks()\n",
    "plt.setp(labels, rotation=90)\n",
    "plt.grid()\n",
    "ax.tick_params(axis='x', labelsize=10)\n",
    "\n",
    "plt.show()\n",
    "figname=dir_figs+storm+'_era_slp_evolution'\n",
    "fig.savefig(figname+'.png')"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Ouverture des données de Pmer et Z."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_data=\"E:/BIG_DATA/ERA5_storms/data/\"\n",
    "\n",
    "f1    = xr.open_mfdataset(dir_data+storm+\"/msl*.nc\")\n",
    "f2    = xr.open_dataset(dir_data+storm+\"/z.nc\")\n",
    "print(f1)\n",
    "print(f2)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Choix du pas de temps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interval=int(input(\"Enter the desired time interval to match with track data : \"))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Ouverture des données de Pmer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f1    = xr.open_mfdataset(dir_data+storm+\"/msl*.nc\").sel(\n",
    "    time=slice(date1,date2)).sel(\n",
    "    latitude=slice(latN,latS)).sel(\n",
    "    longitude=slice(lonW,lonE))\n",
    "\n",
    "mslp = f1['msl'][::interval]/100\n",
    "lat  = mslp.latitude.values\n",
    "time  = mslp.time.values\n",
    "lon  = mslp.longitude.values\n",
    "\n",
    "print(mslp)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Ouverture des données de Z."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f2    = xr.open_dataset(dir_data+storm+\"/z.nc\").sel(\n",
    "    time=slice(date1,date2)).sel(\n",
    "    latitude=slice(latN,latS)).sel(\n",
    "    longitude=slice(lonW,lonE))\n",
    "\n",
    "lev = f2.level.values\n",
    "\n",
    "z = f2['z'][::interval]/9.81\n",
    "print(z)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Grille 2D lat/lon (pour la fonction min/max)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lat  = mslp.latitude.values\n",
    "lon  = mslp.longitude.values\n",
    "lon_grid, lat_grid = np.meshgrid(lon, lat)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Définition des couches verticales."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(lev)\n",
    "lower_layer = int(input('Enter 1 for 900-600 hPa lower layer and 2 for 925-700 hPa lower layer : ')) \n",
    "if (lower_layer==1):\n",
    "    listlev1=[900, 800, 700, 600]\n",
    "if (lower_layer==2):\n",
    "    listlev1=[925, 875, 825, 775, 700]\n",
    "higher_layer = int(input('Enter 1 for 600-300 hPa higher layer and 2 for 700-400 hPa higher layer : ')) \n",
    "if (higher_layer==1):\n",
    "    listlev2=[600, 500, 400, 300]\n",
    "if (higher_layer==2):\n",
    "    listlev2=[700, 600, 500, 400]"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Calcul des épaisseurs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(listlev1)\n",
    "print(z)\n",
    "thickness=z.sel(level=listlev1[-1])-z.sel(level=listlev1[0])\n",
    "print(thickness)\n",
    "print(np.min(thickness))\n",
    "print(np.max(thickness))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Cartes Pmer + tracking + et épaisseurs en basses couches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "mslp_levels = np.arange(900,1072,2)\n",
    "#dz_levels = np.arange(3000,3525,25)\n",
    "dz_levels=np.linspace(np.min(thickness),np.max(thickness),20)\n",
    "\n",
    "for i in tqdm(range(len(time))):\n",
    "    #print(str(time[i])[0:13])\n",
    "    fig = plt.figure(figsize=(17., 12.))\n",
    "    ax = fig.add_subplot(1, 1, 1, projection=projection)\n",
    "    ax.set_title('Storm '+storm+' - MSLP tracking and low level thickness',loc='left',fontsize=14)\n",
    "    ax.set_title(str(time[i])[0:13],loc='right',fontsize=14)\n",
    "    plot_background(ax)\n",
    "    ax.set_extent(*bounds, crs=ccrs.PlateCarree())\n",
    "    \n",
    "    c1 = ax.contour(lon, lat, mslp[i,:,:], levels=mslp_levels, \n",
    "                    colors=\"black\", linewidths=1, transform=ccrs.PlateCarree())\n",
    "    ax.clabel(c1,fmt='%4.1i',fontsize=10)\n",
    "    ax.scatter(liste_lon[i],liste_lat[i], c='green', transform=ccrs.PlateCarree())\n",
    "    ax.plot(liste_lon[0:i+1],liste_lat[0:i+1], c='red', marker='+', transform=ccrs.PlateCarree())\n",
    "    plot_maxmin_points(lon_grid, lat_grid, mslp[i,:,:], 'min', 25,\n",
    "                       symbol='L', color='b', transform=ccrs.PlateCarree())\n",
    "    cf = ax.contourf(lon, lat, thickness[i,:,:], levels=dz_levels, \n",
    "                     cmap='jet', transform=ccrs.PlateCarree())\n",
    "    cb = fig.colorbar(cf, orientation='horizontal', aspect=65, shrink=0.5, pad=0.1)\n",
    "    cb.set_label('m')\n",
    "    #plt.show()\n",
    "    figname=dir_anim+'/MSL_thickness_tracking_'+str(time[i])[0:13]\n",
    "    fig.savefig(figname+'.png',bbox_inches='tight')\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "gif_filepath = dir_anim+'/MSL_thickness_tracking.gif'\n",
    "make_animation(gif_filepath)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Cartes Pmer + tracking + et épaisseurs 900-600hPa avec découpage dans un rayon de 500km autour du centre de la tempête."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "max_dist = 500      # max distance in km\n",
    "\n",
    "for i in tqdm(range(len(time))):\n",
    "    #print(str(time[i])[0:13])\n",
    "    data_slp = mslp[i,:,:]\n",
    "    data_thickness = thickness[i,:,:]\n",
    "        \n",
    "    # Center coordinates\n",
    "    clat = liste_lat[i]\n",
    "    clon = liste_lon[i] \n",
    "        \n",
    "    # Calculate distance between center and all other lat/lon pairs\n",
    "    distance = haversine(lon[:,np.newaxis], lat, clon, clat)\n",
    "    distance2=np.transpose(distance)\n",
    "    # Mask distance array where distance > max_dist\n",
    "    distance_m = np.ma.masked_greater(distance2, max_dist)\n",
    "    # Mask the data array based on the distance mask\n",
    "    data_slp_m = np.ma.masked_where(distance2 > max_dist, data_slp)\n",
    "    data_thickness_m = np.ma.masked_where(distance2 > max_dist, data_thickness)\n",
    "\n",
    "    fig = plt.figure(figsize=(17., 12.))\n",
    "    ax = fig.add_subplot(1, 1, 1, projection=projection)\n",
    "    ax.set_title('Storm '+storm+' - MSLP tracking and low level thickness',loc='left',fontsize=14)\n",
    "    ax.set_title(str(time[i])[0:13],loc='right',fontsize=14)\n",
    "    plot_background(ax)\n",
    "    ax.set_extent(*bounds, crs=ccrs.PlateCarree())\n",
    "    # MSLP in contours and min max\n",
    "    c1 = ax.contour(lon, lat, data_slp_m, levels=mslp_levels,\n",
    "                    colors=\"black\", linewidths=1, transform=ccrs.PlateCarree())\n",
    "    ax.clabel(c1,fmt='%4.1i',fontsize=10)\n",
    "    ax.scatter(liste_lon[i],liste_lat[i], c='green', transform=ccrs.PlateCarree())\n",
    "    ax.plot(liste_lon[0:i+1],liste_lat[0:i+1], c='red', marker='+', transform=ccrs.PlateCarree())\n",
    "    cf = ax.contourf(lon, lat, data_thickness_m, levels=dz_levels, \n",
    "                     cmap='jet', transform=ccrs.PlateCarree())\n",
    "    cb = fig.colorbar(cf, orientation='horizontal', aspect=65, shrink=0.5, pad=0.1)\n",
    "    cb.set_label('m')\n",
    "    #plt.show()\n",
    "    figname=dir_anim+'/MSL_tracking_zoom_'+str(time[i])[0:13]\n",
    "    fig.savefig(figname+'.png',bbox_inches='tight')\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "gif_filepath = dir_anim+'/MSL_tracking_zoom.gif'\n",
    "make_animation(gif_filepath)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Choix définitif du rayon max pour le calcul des paramètres du CPS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_dist = int(input('Enter the value of the maximum radius (in km) to compute the CPS diagrams : ')) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calcul du paramètre de symétrie thermique du CPS\n",
    "\n",
    "$$B=h*(\\overline{Z_{600}-Z_{900}})_R-\\overline{Z_{600}-Z_{900}})_L$$\n",
    "\n",
    "h = hemisphere (1 = NH, -1 = SH).\n",
    "\n",
    "$$(\\overline{Z_{600}-Z_{900}})_R$$ : mean 900-600hPa thickness in semicircle right of motion\n",
    "\n",
    "$$(\\overline{Z_{600}-Z_{900}})_L$$ : mean 900-600hPa thickness in semicircle left of motion"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Sélection des données."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print('Available dates : ', liste_time[0:-1])\n",
    "\n",
    "date = input('Enter a date (YYYY-MM-DDTXX) : ')\n",
    "data = thickness.sel(time=date)\n",
    "data_slp = mslp.sel(time=date)\n",
    "print(data)\n",
    "\n",
    "idate = list(liste_time).index(date)\n",
    "print('index of date : ', idate)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Angle de la trajectoire du système (par rapport à l'instant T+1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "clat = liste_lat[idate]\n",
    "clon = liste_lon[idate]\n",
    "ang=get_bearing(clat, clon, liste_lat[idate+1], liste_lon[idate+1])"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Calcul de l'angle entre le centre et toutes les paires de lat/lon du domaine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "qq_ang_all=np.zeros([len(lon),len(lat)],dtype='f')\n",
    "for x in range(0,len(lon)):\n",
    "    for y in range(0,len(lat)):\n",
    "        qq_ang_all[x,y]=get_bearing(clat, clon, lat[y], lon[x])     \n",
    "qq_ang_all=np.transpose(qq_ang_all)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Création des cadrans d'épaisseurs à droite et à gauche de la trajectoire. On distingue selon que le centre se déplace vers le NE et SE ou le NW et le SW."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "Zl=np.zeros([len(lat),len(lon)],dtype='f')\n",
    "Zr=np.zeros([len(lat),len(lon)],dtype='f')\n",
    "\n",
    "for jlat in tqdm(range(0,len(lat)-1)):\n",
    "    for jlon in range(0,len(lon)-1):\n",
    "        #set values along great circle line to missing (neither left nor right)\n",
    "        if qq_ang_all[jlat,jlon] == ang:\n",
    "            Zl[jlat,jlon] = 0\n",
    "            Zr[jlat,jlon] = 0\n",
    "        #for storm travel angles in quadrants 1 and 2 (NE and SE)\n",
    "        elif (ang >= 0 and ang < 180):\n",
    "            if (qq_ang_all[jlat,jlon] > ang and qq_ang_all[jlat,jlon]  < ang+180):\n",
    "                Zl[jlat,jlon] = 0\n",
    "                Zr[jlat,jlon] = data[jlat,jlon]\n",
    "            else:\n",
    "                Zr[jlat,jlon] = 0\n",
    "                Zl[jlat,jlon] = data[jlat,jlon]\n",
    "        #for storm travel angles in quadrants 3 and 4 (NW and SW)\n",
    "        elif (ang >= 180 and ang < 360):\n",
    "            if (qq_ang_all[jlat,jlon] > ang-180 and qq_ang_all[jlat,jlon] < ang):\n",
    "                Zr[jlat,jlon] = 0\n",
    "                Zl[jlat,jlon] = data[jlat,jlon]\n",
    "            else:\n",
    "                Zl[jlat,jlon] = 0\n",
    "                Zr[jlat,jlon] = data[jlat,jlon]"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "On ne conserve les données que dans un rayon défini autour du centre. On remplace les zeros par des NaN sinon problème sur le calcul de la moyenne."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Calculate distance between center and all other lat/lon pairs\n",
    "distance = haversine(lon[:,np.newaxis], lat, clon, clat)\n",
    "distance2=np.transpose(distance)\n",
    "# Mask distance array where distance > max_dist\n",
    "distance_m = np.ma.masked_greater(distance2, max_dist)\n",
    "# Mask the data array based on the distance mask\n",
    "data_slp_m = np.ma.masked_where(distance2 > max_dist, data_slp)\n",
    "Zr_m = np.ma.masked_where(distance2 > max_dist, Zr)\n",
    "Zl_m = np.ma.masked_where(distance2 > max_dist, Zl)\n",
    "\n",
    "Zr_m[Zr_m == 0] = np.nan\n",
    "Zl_m[Zl_m == 0] = np.nan"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Vérification sur une carte du cadran droite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "levels=np.arange(0,370,10)\n",
    "\n",
    "fig = plt.figure(figsize=(17., 12.))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=projection)\n",
    "fig.suptitle('Storm '+storm+' - MSLP, low level thickness and right quadrant',fontsize=16)\n",
    "ax.set_title(liste_time[idate],fontsize=14)\n",
    "plot_background(ax)\n",
    "ax.set_extent(*bounds, crs=ccrs.PlateCarree())\n",
    " \n",
    "c0 = ax.contour(lon, lat, data_slp_m, mslp_levels,\n",
    "                colors=\"black\", linewidths=1, transform=ccrs.PlateCarree())\n",
    "ax.clabel(c0,fmt='%4.1i',fontsize=10)\n",
    "cf1 = ax.contourf(lon, lat, Zr_m, dz_levels,\n",
    "                 cmap='jet', transform=ccrs.PlateCarree())\n",
    "cb = fig.colorbar(cf1, orientation='horizontal', aspect=65, shrink=0.5, pad=0.1)\n",
    "cb.set_label('m')\n",
    "c3 = ax.contour(lon, lat, qq_ang_all, levels,\n",
    "                colors=\"grey\", linewidths=0.5, alpha=0.8, transform=ccrs.PlateCarree())\n",
    "ax.clabel(c3,fmt='%4.1i',fontsize=10)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "figname=dir_figs+'/Thermal_symmetry0'\n",
    "fig.savefig(figname+'.png',bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Vérification sur une carte du cadran gauche."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "levels=np.arange(0,370,10)\n",
    "\n",
    "fig = plt.figure(figsize=(17., 12.))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=projection)\n",
    "fig.suptitle('Storm '+storm+' - MSLP, low level thickness and left quadrant',fontsize=16)\n",
    "ax.set_title(liste_time[idate],fontsize=14)\n",
    "plot_background(ax)\n",
    "ax.set_extent(*bounds, crs=ccrs.PlateCarree())\n",
    " \n",
    "c0 = ax.contour(lon, lat, data_slp_m, mslp_levels,\n",
    "                colors=\"black\", linewidths=1, transform=ccrs.PlateCarree())\n",
    "ax.clabel(c0,fmt='%4.1i',fontsize=10)\n",
    "cf2 = ax.contourf(lon, lat, Zl_m, dz_levels,\n",
    "                 cmap='jet', transform=ccrs.PlateCarree())\n",
    "cb = fig.colorbar(cf2, orientation='horizontal', aspect=65, shrink=0.5, pad=0.1)\n",
    "cb.set_label('m')\n",
    "c3 = ax.contour(lon, lat, qq_ang_all, levels,\n",
    "                colors=\"grey\", linewidths=0.5, alpha=0.8, transform=ccrs.PlateCarree())\n",
    "ax.clabel(c3,fmt='%4.1i',fontsize=10)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "figname=dir_figs+'/Thermal_symmetry1'\n",
    "fig.savefig(figname+'.png',bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Calcul du paramètre de symétrie thermique B."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# cyclone thermal symmetry parameter using classic average\n",
    "Zr_mean=np.nanmean(Zr_m)\n",
    "Zl_mean=np.nanmean(Zl_m)\n",
    "B=Zr_mean-Zl_mean\n",
    "print('Right quadrant thickness : ', Zr_mean)\n",
    "print('Left quadrant thickness : ', Zl_mean)\n",
    "\n",
    "if clat<0:\n",
    "    B=-B"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Vérification sur une carte du calcul du paramètre de symétrie thermique B."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "levels=np.arange(0,370,10)\n",
    "\n",
    "fig = plt.figure(figsize=(17., 12.))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=projection)\n",
    "fig.suptitle('Storm '+storm+' - MSLP, low level thickness and right/left quadrants',fontsize=16)\n",
    "ax.set_title(liste_time[idate]+' : B = '+str(int(B)),fontsize=14)\n",
    "plot_background(ax)\n",
    "ax.set_extent(*bounds, crs=ccrs.PlateCarree())\n",
    " \n",
    "c0 = ax.contour(lon, lat, data_slp_m, mslp_levels,\n",
    "                colors=\"black\", linewidths=1, transform=ccrs.PlateCarree())\n",
    "ax.clabel(c0,fmt='%4.1i',fontsize=10)\n",
    "cf1 = ax.contourf(lon, lat, Zr_m, dz_levels,\n",
    "                 cmap='jet', transform=ccrs.PlateCarree())\n",
    "cf2 = ax.contourf(lon, lat, Zl_m, dz_levels,\n",
    "                 cmap='jet', transform=ccrs.PlateCarree())\n",
    "cb = fig.colorbar(cf1, orientation='horizontal', aspect=65, shrink=0.5, pad=0.1)\n",
    "cb.set_label('m')\n",
    "c3 = ax.contour(lon, lat, qq_ang_all, levels,\n",
    "                colors=\"grey\", linewidths=0.5, alpha=0.8, transform=ccrs.PlateCarree())\n",
    "ax.clabel(c3,fmt='%4.1i',fontsize=10)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "figname=dir_figs+'/Thermal_symmetry'\n",
    "fig.savefig(figname+'.png',bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Calcul du paramètre de symétrie thermique B pour tous les temps et réalisation des cartes (attention, la dernière date est absente car on a pas l'instant d'après pour calculer l'angle de la trajectoire)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "B_series=[]\n",
    "\n",
    "for i in tqdm(range(len(time)-1)):\n",
    "    #print(liste_time[i])\n",
    "    data = thickness[i,:,:]\n",
    "    data_slp = mslp[i,:,:]\n",
    "\n",
    "    # Center coordinates\n",
    "    clat = liste_lat[i]\n",
    "    clon = liste_lon[i]\n",
    "\n",
    "    # Calculate bearing between center and all other lat/lon pairs\n",
    "    qq_ang_all=np.zeros([len(lon),len(lat)],dtype='f')\n",
    "    for x in range(0,len(lon)):\n",
    "        for y in range(0,len(lat)):\n",
    "            qq_ang_all[x,y]=get_bearing(clat, clon, lat[y], lon[x])     \n",
    "    qq_ang_all=np.transpose(qq_ang_all)\n",
    "\n",
    "    # Angle of travel of the center\n",
    "    ang=get_bearing(clat, clon, liste_lat[i+1], liste_lon[i+1])\n",
    "\n",
    "    # Define right/left cadrans\n",
    "    Zl=np.zeros([len(lat),len(lon)],dtype='f')\n",
    "    Zr=np.zeros([len(lat),len(lon)],dtype='f')\n",
    "\n",
    "    for jlat in range(0,len(lat)-1):\n",
    "        for jlon in range(0,len(lon)-1):\n",
    "            #set values along great circle line to missing (neither left nor right)\n",
    "            if qq_ang_all[jlat,jlon] == ang:\n",
    "                Zl[jlat,jlon] = 0\n",
    "                Zr[jlat,jlon] = 0\n",
    "            #for storm travel angles in quadrants 1 and 2 (NE and SE)\n",
    "            elif (ang >= 0 and ang < 180):\n",
    "                if (qq_ang_all[jlat,jlon] > ang and qq_ang_all[jlat,jlon]  < ang+180):\n",
    "                    Zl[jlat,jlon] = 0\n",
    "                    Zr[jlat,jlon] = data[jlat,jlon]\n",
    "                else:\n",
    "                    Zr[jlat,jlon] = 0\n",
    "                    Zl[jlat,jlon] = data[jlat,jlon]\n",
    "            #for storm travel angles in quadrants 3 and 4 (NW and SW)\n",
    "            elif (ang >= 180 and ang < 360):\n",
    "                if (qq_ang_all[jlat,jlon] > ang-180 and qq_ang_all[jlat,jlon] < ang):\n",
    "                    Zr[jlat,jlon] = 0\n",
    "                    Zl[jlat,jlon] = data[jlat,jlon]\n",
    "                else:\n",
    "                    Zl[jlat,jlon] = 0\n",
    "                    Zr[jlat,jlon] = data[jlat,jlon]\n",
    "\n",
    "    # Calculate distance between center and all other lat/lon pairs\n",
    "    distance = haversine(lon[:,np.newaxis], lat, clon, clat)\n",
    "    distance2=np.transpose(distance)\n",
    "    # Mask distance array where distance > max_dist\n",
    "    distance_m = np.ma.masked_greater(distance2, max_dist)\n",
    "    # Mask the data array based on the distance mask\n",
    "    data_slp_m = np.ma.masked_where(distance2 > max_dist, data_slp)\n",
    "    Zr_m = np.ma.masked_where(distance2 > max_dist, Zr)\n",
    "    Zl_m = np.ma.masked_where(distance2 > max_dist, Zl)\n",
    "    \n",
    "    Zr_m[Zr_m == 0] = np.nan\n",
    "    Zl_m[Zl_m == 0] = np.nan\n",
    "\n",
    "    # cyclone thermal symmetry parameter using weighted average\n",
    "    Zr_mean=np.nanmean(Zr_m)\n",
    "    Zl_mean=np.nanmean(Zl_m)\n",
    "    B=Zr_mean-Zl_mean\n",
    "    \n",
    "    if clat<0:\n",
    "        B=-B\n",
    "    #print('B = ',B)\n",
    "    \n",
    "    B_series.append(B)\n",
    "\n",
    "    fig = plt.figure(figsize=(17., 12.))\n",
    "    ax = fig.add_subplot(1, 1, 1, projection=projection)\n",
    "    fig.suptitle('Storm '+storm+' - MSLP, low level thickness and right/left quadrants',fontsize=16)\n",
    "    ax.set_title(liste_time[i]+' : B = '+str(int(B)),fontsize=14)\n",
    "    plot_background(ax)\n",
    "    ax.set_extent(*bounds, crs=ccrs.PlateCarree())\n",
    "\n",
    "    c0 = ax.contour(lon, lat, data_slp_m, mslp_levels,\n",
    "                    colors=\"black\", linewidths=1, transform=ccrs.PlateCarree())\n",
    "    ax.clabel(c0,fmt='%4.1i',fontsize=10)\n",
    "    cf1 = ax.contourf(lon, lat, Zr_m, dz_levels,\n",
    "                     cmap='jet', transform=ccrs.PlateCarree())\n",
    "    cf2 = ax.contourf(lon, lat, Zl_m, dz_levels,\n",
    "                     cmap='jet', transform=ccrs.PlateCarree())\n",
    "    cb = fig.colorbar(cf1, orientation='horizontal', aspect=65, shrink=0.5, pad=0.1)\n",
    "    cb.set_label('m')\n",
    "    c3 = ax.contour(lon, lat, qq_ang_all, levels,\n",
    "                    colors=\"grey\", linewidths=0.5, alpha=0.8, transform=ccrs.PlateCarree())\n",
    "    ax.clabel(c3,fmt='%4.1i',fontsize=10)\n",
    "\n",
    "    #plt.show()\n",
    "    figname=dir_anim+'/Thermal_symmetry_'+liste_time[i]\n",
    "    fig.savefig(figname+'.png',bbox_inches='tight')\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "gif_filepath = dir_anim+'/Thermal_symmetry.gif'\n",
    "make_animation(gif_filepath)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Tracé de la série temporelle du paramètre B."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(B_series)\n",
    "\n",
    "fig = plt.figure(figsize=(15., 8.))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_title('Storm '+storm+' - low level thermal asymmetry parameter (B)'\n",
    "             ,loc='left',fontsize=14)\n",
    "ax.set_title(liste_time[0]+'-'+liste_time[-2],loc='right',fontsize=14)\n",
    "ax.plot(time[0:-1], B_series, label='B')\n",
    "plt.axhline(y=10, color='r')\n",
    "\n",
    "figname=dir_figs+'/Thermal_symmetry_timeseries'\n",
    "fig.savefig(figname+'.png',bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calcul du paramètre de coeur chaud/coeur froid du CPS\n",
    "\n",
    "Approximation of the height gradient on an isobaric surface (i.e. magnitude of the geostrophic wind) :\n",
    "\n",
    "$$\\Delta{Z}=Z_{max}-Z_{min}$$\n",
    "\n",
    "Magnitude of the cyclone thermal wind measured as the vertical profile of cyclone geopotential height gradient (cyclone strength) :\n",
    "\n",
    "$$-V_T^L=\\frac{\\partial{\\Delta Z}}{\\partial{ln(P)}}_{900}^{600}$$\n",
    "\n",
    "$$-V_T^U=\\frac{\\partial{\\Delta Z}}{\\partial{ln(P)}}_{600}^{300}$$"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Selection de la date."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print('Available dates : ', liste_time[0:-1])\n",
    "\n",
    "date = input('Enter a date (YYYY-MM-DDTXX) : ')\n",
    "data = thickness.sel(time=date)\n",
    "data_slp = mslp.sel(time=date)\n",
    "print(data)\n",
    "\n",
    "idate = list(liste_time).index(date)\n",
    "print('index of date : ', idate)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Calcul de deltaZ = Zmax - Zmin pour tous les niveaux et dans chaque couche (pour un seul temps)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "clat = liste_lat[idate]\n",
    "clon = liste_lon[idate]\n",
    "\n",
    "deltaZ1=[]\n",
    "deltaZ2=[]\n",
    "\n",
    "#low layer loop\n",
    "for plev in listlev1:\n",
    "\n",
    "    ilev = list(lev).index(plev)\n",
    "    data_z1 = z[idate,ilev,:,:]\n",
    "\n",
    "    # Calculate distance between center and all other lat/lon pairs\n",
    "    distance = haversine(lon[:,np.newaxis], lat, clon, clat)\n",
    "    distance2=np.transpose(distance)\n",
    "    # Mask distance array where distance > max_dist\n",
    "    distance_m = np.ma.masked_greater(distance2, max_dist)\n",
    "    # Mask the data array based on the distance mask\n",
    "    data_z1_m = np.ma.masked_where(distance2 > max_dist, data_z1)\n",
    "    zmin=np.min(data_z1_m)\n",
    "    zmax=np.max(data_z1_m)\n",
    "    delta_z1=zmax-zmin\n",
    "    deltaZ1.append(int(delta_z1))\n",
    "\n",
    "#high layer loop\n",
    "for plev in listlev2:\n",
    "\n",
    "    ilev = list(lev).index(plev)\n",
    "    data_z2 = z[idate,ilev,:,:]\n",
    "\n",
    "    # Calculate distance between center and all other lat/lon pairs\n",
    "    distance = haversine(lon[:,np.newaxis], lat, clon, clat)\n",
    "    distance2=np.transpose(distance)\n",
    "    # Mask distance array where distance > max_dist\n",
    "    distance_m = np.ma.masked_greater(distance2, max_dist)\n",
    "    # Mask the data array based on the distance mask\n",
    "    data_z2_m = np.ma.masked_where(distance2 > max_dist, data_z2)\n",
    "    zmin=np.min(data_z2_m)\n",
    "    zmax=np.max(data_z2_m)\n",
    "    delta_z2=zmax-zmin\n",
    "    deltaZ2.append(int(delta_z2))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Vérification du calcul de la pente de la surface deltaZ (couche supérieure, dernier niveau)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(zmax)\n",
    "print(zmin)\n",
    "\n",
    "fig = plt.figure(figsize=(17., 12.))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=projection)\n",
    "fig.suptitle('Storm '+storm+' - Geopotential at '+str(plev)+ 'hPa',fontsize=16)\n",
    "ax.set_title(liste_time[idate]+ ' : zmax-zmin = '+str(int(delta_z2)),fontsize=14)\n",
    "plot_background(ax)\n",
    "ax.set_extent(*bounds, crs=ccrs.PlateCarree())\n",
    " \n",
    "cf = ax.contourf(lon, lat, data_z2_m,\n",
    "                     cmap='jet', transform=ccrs.PlateCarree())\n",
    "cb = fig.colorbar(cf, orientation='horizontal', aspect=65, shrink=0.5, pad=0.1)\n",
    "cb.set_label('m')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "figname=dir_figs+'/deltaz_map'\n",
    "fig.savefig(figname+'.png',bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Calcul de ln(P) pour chaque couche."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(\"Lower layer vertical levels : \", listlev1)\n",
    "print(\"Lower layer deltaZ : \", deltaZ1)\n",
    "print(\"Higher layer vertical levels : \", listlev2)\n",
    "print(\"Higher layer deltaZ : \", deltaZ2)\n",
    "\n",
    "listlev11=[n * 100 for n in listlev1]\n",
    "lnP1=np.log(listlev11)\n",
    "\n",
    "listlev22=[n * 100 for n in listlev2]\n",
    "lnP2=np.log(listlev22)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Calcul des paramètres de coeur chaud/coeur froid (-VTU, -VTL) pour tous les temps. (attention, on ne prend pas la dernière date)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "VTL_series=[]\n",
    "VTU_series=[]\n",
    "\n",
    "for i in tqdm(range(len(time)-1)):\n",
    "    #print(liste_time[i])\n",
    "\n",
    "    # Center coordinates\n",
    "    clat = liste_lat[i]\n",
    "    clon = liste_lon[i]\n",
    "    deltaZ1=[]\n",
    "    deltaZ2=[]\n",
    "\n",
    "    for plev in listlev1:\n",
    "\n",
    "        ilev = list(lev).index(plev)\n",
    "        data_z1 = z[i,ilev,:,:]\n",
    "        # Calculate distance between center and all other lat/lon pairs\n",
    "        distance = haversine(lon[:,np.newaxis], lat, clon, clat)\n",
    "        distance2=np.transpose(distance)\n",
    "        # Mask distance array where distance > max_dist\n",
    "        distance_m = np.ma.masked_greater(distance2, max_dist)\n",
    "        # Mask the data array based on the distance mask\n",
    "        data_z1_m = np.ma.masked_where(distance2 > max_dist, data_z1)\n",
    "        zmin=np.min(data_z1_m)\n",
    "        zmax=np.max(data_z1_m)\n",
    "        delta_z1=zmax-zmin\n",
    "        deltaZ1.append(int(delta_z1))\n",
    "\n",
    "    for plev in listlev2:\n",
    "\n",
    "        ilev = list(lev).index(plev)\n",
    "        data_z2 = z[i,ilev,:,:]\n",
    "\n",
    "        # Calculate distance between center and all other lat/lon pairs\n",
    "        distance = haversine(lon[:,np.newaxis], lat, clon, clat)\n",
    "        distance2=np.transpose(distance)\n",
    "        # Mask distance array where distance > max_dist\n",
    "        distance_m = np.ma.masked_greater(distance2, max_dist)\n",
    "        # Mask the data array based on the distance mask\n",
    "        data_z2_m = np.ma.masked_where(distance2 > max_dist, data_z2)\n",
    "        zmin=np.min(data_z2_m)\n",
    "        zmax=np.max(data_z2_m)\n",
    "        delta_z2=zmax-zmin\n",
    "        deltaZ2.append(int(delta_z2))\n",
    "\n",
    "    coef1=(deltaZ1[-1]-deltaZ1[0])/(lnP1[-1]-lnP1[0])\n",
    "    VTL_series.append(int(coef1))\n",
    "\n",
    "    coef2=(deltaZ2[-1]-deltaZ2[0])/(lnP2[-1]-lnP2[0])\n",
    "    VTU_series.append(int(coef2))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Tracé des séries temporelles des paramètres -VTU et -VTL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(VTU_series)\n",
    "print(VTL_series)\n",
    "\n",
    "fig = plt.figure(figsize=(15., 8.))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_title('Storm '+storm+' - Low level thermal wind parameter'\n",
    "             ,loc='left',fontsize=14)\n",
    "ax.set_title(liste_time[0]+'-'+liste_time[-2],loc='right',fontsize=14)\n",
    "ax.plot(time[0:-1], VTL_series, label='-$VT_T^L$')\n",
    "\n",
    "plt.axhline(y=0, color='r')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "figname=dir_figs+'/VTL_timeseries'\n",
    "fig.savefig(figname+'.png',bbox_inches='tight')\n",
    "\n",
    "fig = plt.figure(figsize=(15., 8.))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_title('Storm '+storm+' - High level thermal wind parameter'\n",
    "             ,loc='left',fontsize=14)\n",
    "ax.set_title(liste_time[0]+'-'+liste_time[-2],loc='right',fontsize=14)\n",
    "ax.plot(time[0:-1], VTU_series, label='-$VT_T^U$')\n",
    "\n",
    "plt.axhline(y=0, color='r')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "figname=dir_figs+'/VTU_timeseries'\n",
    "fig.savefig(figname+'.png',bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Tracé des diagrammes des phases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "time_str=[x for x in range(len(time)-1)]\n",
    "date_str=[x for x in range(len(time)-1)]\n",
    "for i in range(len(time)-1):\n",
    "    time_str[i] = str(time[i])\n",
    "    date_str[i] = time_str[i][5:13]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "x1=-600\n",
    "x2=300\n",
    "y1=-50\n",
    "y2=140\n",
    "\n",
    "if storm=='Zorbas':\n",
    "    x1=-200\n",
    "    x2=300\n",
    "    y1=-50\n",
    "    y2=60\n",
    "    \n",
    "fig = plt.figure(figsize=(15., 8.))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_title('Storm '+storm,loc='left',fontsize=14)\n",
    "ax.set_title(liste_time[0]+'-'+liste_time[-2],loc='right',fontsize=14)\n",
    "plt.axhline(y=10, color='grey', linewidth=10,zorder=1)\n",
    "plt.axvline(x=0, color='grey', linewidth=10,zorder=1)\n",
    "plt.plot(VTL_series, B_series)\n",
    "plt.scatter(VTL_series,  B_series, c='k',zorder=2)\n",
    "plt.scatter(VTL_series[0], B_series[0], c='green',zorder=2)\n",
    "plt.scatter(VTL_series[-1], B_series[-1], c='red',zorder=2)\n",
    "plt.xlabel(\"$-V_T^L$ [low level thermal wind]\", labelpad=20, fontsize = 16)\n",
    "plt.ylabel(\"B [low level thickness asymmetry]\", fontsize = 16)\n",
    "\n",
    "plt.text(x1-35,-15,'Symmetric', rotation=90., color='red', fontsize = 14)\n",
    "plt.text(x1-35, 20,'Asymmetric', rotation=90., color='red', fontsize = 14)\n",
    "plt.text(-100, -30,'Cold core', rotation=0., color='blue', fontsize = 14)\n",
    "plt.text(100, -30,'Warm core', rotation=0., color='red', fontsize = 14)\n",
    "plt.xlim(x1, x2)\n",
    "plt.ylim(y1, y2)\n",
    "    \n",
    "for i,type in enumerate(date_str):\n",
    "    x = VTL_series[i]\n",
    "    y = B_series[i]\n",
    "    plt.text(x, y, type, fontsize=9)\n",
    "    \n",
    "xrange = [(x1, -x1)]\n",
    "yrange1 = (y1, y1+110)\n",
    "\n",
    "c1 = collections.BrokenBarHCollection(xrange, yrange1, facecolor='blue', alpha=0.2)\n",
    "ax.add_collection(c1)\n",
    "xrange = [(x1, -x1)]\n",
    "yrange1 = (10, 10+120)\n",
    "c1 = collections.BrokenBarHCollection(xrange, yrange1, facecolor='blue', alpha=0.2)\n",
    "ax.add_collection(c1)\n",
    "\n",
    "xrange = [(0, 0+x2)]\n",
    "yrange1 = (y1, y1+110)\n",
    "\n",
    "c1 = collections.BrokenBarHCollection(xrange, yrange1, facecolor='red', alpha=0.2)\n",
    "ax.add_collection(c1)\n",
    "xrange = [(0, 0+x2)]\n",
    "yrange1 = (10, 10+120)\n",
    "c1 = collections.BrokenBarHCollection(xrange, yrange1, facecolor='red', alpha=0.2)\n",
    "ax.add_collection(c1)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "figname=dir_figs+'/CPS1'\n",
    "fig.savefig(figname+'.png',bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if storm=='Zorbas':\n",
    "    x1=-300\n",
    "    x2=300\n",
    "    y1=-300\n",
    "    y2=300\n",
    "\n",
    "else:\n",
    "    x1=-600\n",
    "    x2=300\n",
    "    y1=-600\n",
    "    y2=300\n",
    "\n",
    "fig = plt.figure(figsize=(15., 8.))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_title('Storm '+storm,loc='left',fontsize=14)\n",
    "ax.set_title(liste_time[0]+'-'+liste_time[-2],loc='right',fontsize=14)\n",
    "plt.axhline(y=0, color='grey', linewidth=10,zorder=1)\n",
    "plt.axvline(x=0, color='grey', linewidth=10,zorder=1)\n",
    "plt.plot(VTL_series, VTU_series)\n",
    "plt.scatter(VTL_series, VTU_series, c='k',zorder=2)\n",
    "plt.scatter(VTL_series[0], VTU_series[0], c='green',zorder=2)\n",
    "plt.scatter(VTL_series[-1], VTU_series[-1], c='red',zorder=2)\n",
    "\n",
    "plt.xlabel(\"$-V_T^L$ [Low level thermal wind]\", fontsize = 16)\n",
    "plt.ylabel(\"$-V_T^U$ [High level thermal wind]\", fontsize = 16)\n",
    "plt.text(-200, y1+20,'Deep cold core', rotation=0., color='blue', fontsize = 14)\n",
    "plt.text(90, 270,'Deep warm core', rotation=0., color='red', fontsize = 14)\n",
    "plt.text(90, y1+20,'Shallow warm core', rotation=0., color='orange', fontsize = 14)\n",
    "\n",
    "plt.xlim(x1, x2)\n",
    "plt.ylim(y1, y2)\n",
    "\n",
    "for i,type in enumerate(date_str):\n",
    "    x = VTL_series[i]\n",
    "    y = VTU_series[i]\n",
    "    plt.text(x, y, type, fontsize=9)\n",
    "\n",
    "xrange = [(0, 0+x2)]\n",
    "yrange1 = (0, 0+y2)\n",
    "c1 = collections.BrokenBarHCollection(xrange, yrange1, facecolor='r', alpha=0.2)\n",
    "ax.add_collection(c1)\n",
    "xrange = [(-600, 600)]\n",
    "yrange1 = (-600, 600)\n",
    "c1 = collections.BrokenBarHCollection(xrange, yrange1, facecolor='b', alpha=0.2)\n",
    "ax.add_collection(c1)\n",
    "xrange = [(x1, -x1)]\n",
    "yrange1 = (0, 300)\n",
    "c1 = collections.BrokenBarHCollection(xrange, yrange1, facecolor='orange', alpha=0.2)\n",
    "ax.add_collection(c1)\n",
    "xrange = [(0, x2)]\n",
    "yrange1 = (y1, -y1)\n",
    "c1 = collections.BrokenBarHCollection(xrange, yrange1, facecolor='orange', alpha=0.2)\n",
    "ax.add_collection(c1)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "figname=dir_figs+'/CPS2'\n",
    "fig.savefig(figname+'.png',bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}