{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "978b8621-e050-47c3-8e66-807cbac823e6",
   "metadata": {},
   "source": [
    "# Basic use for assessing storm climatology in a dataset\n",
    "Here, we examplify usage of `huracanpy` with the dataset of TC in ERA-20C detected by the TRACK algorithm.\n",
    "This is meant to show an example of workflow. Please refer to specific parts of the documentation to learn about each part (e.g. loading, plotting, etc.) in more detail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0757f0c8-9a56-4f49-b08b-63768585fb98",
   "metadata": {},
   "outputs": [],
   "source": [
    "import huracanpy\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import cartopy.crs as ccrs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9ea5a70-a092-4d42-9f6b-378c6fca0728",
   "metadata": {},
   "source": [
    "## Read the file\n",
    "`huracanpy`'s `load` function can handle different track file types. Here, the data is available as a netcdf file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b4ceb72-dbdf-441a-93e2-98807f7fda81",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = huracanpy.load(\n",
    "    huracanpy.example_ERA20C_file\n",
    ")\n",
    "data.psl.attrs[\"units\"] = \"hPa\"  # Fixing misspelled pressure unit\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d79c309-719b-4f28-96b6-6d90831d8032",
   "metadata": {},
   "source": [
    "## Add useful information\n",
    "After loading, you can add various useful information for the analysis (basin, season, category...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e81a21a-3f89-488c-a621-6f966328e4ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the accessor's add_ function to add the info you want\n",
    "## SSHS category\n",
    "data = data.hrcn.add_saffir_simpson_category(wind_name = \"wind_speed_10m\")\n",
    "## Presure category\n",
    "data = data.hrcn.add_pressure_category(slp_name = \"psl\")\n",
    "## Season \n",
    "data = data.hrcn.add_season()\n",
    "# More info are available, in particular geographical ones, but we do not need them for this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fad57ab-5bd3-4cb3-8592-b76e4f5fa0cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a14b9f64-3d86-4f45-bee0-28d1c6ce70fc",
   "metadata": {},
   "source": [
    "# Check the content of the file\n",
    "`huracanpy` provide a coarse plotting function that you can use for checking what is in your data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b894434f-7f28-4993-bbab-91ef1630c988",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Basic plot of the data points\n",
    "data.hrcn.plot_tracks(intensity_var_name = \"wind_speed_10m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4db9aa6d-a806-4956-81d4-7bd653eeffb3",
   "metadata": {},
   "source": [
    "## Climatological metrics\n",
    "You can compute basic statistics: frequency, TC days, ACE. Here shown as yearly averages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83c93284-c843-4017-b627-a2009978424e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Frequency (Number of track per year)\n",
    "data.track_id.hrcn.nunique() / data.season.hrcn.nunique() \n",
    "# number of unique tracks / number of unique season"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6923e97b-496d-4267-a37f-f5c8c61a2ec6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TCD (Accumulated duration of storms per year)\n",
    "data.hrcn.get_track_duration().sum() / 24 / data.season.hrcn.nunique()\n",
    "# Compute duration per track, convert to days / number of unique season"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29c21041-3beb-425e-8e8e-c1624616d862",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ACE per year\n",
    "data = data.hrcn.add_ace(wind_name = \"wind_speed_10m\")\n",
    "data.ace.groupby(data.season).sum().mean()\n",
    "# Get ace for each point, sum by season and average over the seasons\n",
    "# NB: By default, huracanpy computes ACE only for points with wind above 34 knots"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "392197d9-0ad1-480f-bb11-5e65d78e42ba",
   "metadata": {},
   "source": [
    "## Variability\n",
    "With xarray's grouping functionnalities, you can show variations of these statistics. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "785c5b89-de05-47ff-8e06-be88e39eb8dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Interannual\n",
    "\n",
    "fig, axs = plt.subplots(3, sharex=True)\n",
    "# Frequency\n",
    "# In this case, it is easier to go through a dataframe\n",
    "data.to_dataframe().groupby(\"season\").track_id.nunique().plot(ax=axs[0]) \n",
    "axs[0].set_ylabel(\"Number of tracks\")\n",
    "# TCD\n",
    "data.groupby(\"season\").apply(lambda s: s.hrcn.get_track_duration().sum()).plot(ax=axs[1])\n",
    "axs[1].set_ylabel(\"TC days\")\n",
    "# ACE\n",
    "data.groupby(\"season\").sum().ace.plot(ax=axs[2])\n",
    "axs[2].set_ylabel(\"ACE\")\n",
    "\n",
    "for ax in axs:\n",
    "    ax.set_ylim(0)\n",
    "    ax.set_xlabel(\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c20dfe44-ebc7-482b-9c48-3007acfc51ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Seasonal\n",
    "\n",
    "gen = data.hrcn.get_gen_vals()  # Extract the point of genesis for each track\n",
    "(\n",
    "    gen.groupby(\"time.month\").count().lon  # compute number of genesis points per month\n",
    "    / gen.season.hrcn.nunique()  # Normalize by number of season\n",
    ").plot(marker=\"o\")  # plot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "266059b6-22e5-442b-984c-f5e8974aafa8",
   "metadata": {},
   "source": [
    "## Track statistics\n",
    "You can also compute track-level statistics such as duration and lifetime maximum intensity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1843fe4b-299e-43d7-975b-e83e7019c2de",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, sharey=True, figsize=[15, 5])\n",
    "\n",
    "# Duration\n",
    "data.hrcn.get_track_duration().plot.hist(ax=axs[0])\n",
    "\n",
    "# Maximum wind speed\n",
    "data.wind_speed_10m.groupby(data.track_id).max().plot.hist(ax=axs[1])\n",
    "\n",
    "# Minimum SLP\n",
    "data.psl.groupby(data.track_id).min().plot.hist(ax=axs[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c0e3e2e-6a4f-4f6b-a431-0cc794a44be2",
   "metadata": {},
   "source": [
    "## Lifecycles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdc25eb7-4654-4ea2-8ab9-030967d96a2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute times from extremum\n",
    "data[\"time_from_slp_min\"] = data.hrcn.get_time_from_apex(intensity_var_name = \"psl\", stat = \"min\")\n",
    "data[\"time_from_wind_max\"] = data.hrcn.get_time_from_apex(intensity_var_name = \"wind_speed_10m\",)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "214a3d42-ff7e-4074-9569-205e50a0c610",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot with seaborn\n",
    "fig, axs = plt.subplots(2, sharex=True)\n",
    "# SLP lifecycle\n",
    "sns.lineplot(data=data, x=\"time_from_slp_min\", y=\"psl\", ax=axs[0])\n",
    "# Wind lifecycle\n",
    "sns.lineplot(data=data, x=\"time_from_wind_max\", y=\"wind_speed_10m\", ax=axs[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cdd26fe-c54a-438e-8ff6-feb8287fdaed",
   "metadata": {},
   "source": [
    "## Track density"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e6bf09a-8e4d-4c8e-973e-f580e6f13aa4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Density of all points\n",
    "d = data.hrcn.get_density(n_seasons = data.season.hrcn.nunique(), bin_size = 5)\n",
    "huracanpy.plot.density(d, \n",
    "                      cbar_kwargs=dict(label=\"Number of points per 5x5° box per year\"),\n",
    "                      )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bac93d64-4d0d-40db-a4be-f307f05a6369",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Density of all points\n",
    "d = gen.hrcn.get_density(n_seasons = data.season.hrcn.nunique(), bin_size = 10)\n",
    "huracanpy.plot.density(d, \n",
    "                      cbar_kwargs=dict(label=\"Number of genesis per 10x10° box per year\"),\n",
    "                      )"
   ]
  }
 ],
 "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.11.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}