{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6b06ec72-4a81-421e-8628-b3684f6198b7",
   "metadata": {},
   "source": [
    "# Studying a set of tracks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98f00647-d338-4bea-b551-511758fa39c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import huracanpy\n",
    "import numpy as np\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c39bbf8-757d-40ae-b80a-18d3af3d7fd2",
   "metadata": {},
   "source": [
    "## Load tracks\n",
    "Here we show an example with a CSV file that is embedded within HuracanPy.\n",
    "HuracanPy supports many track file formats, see [huracanpy.load guide](../load.ipynb) for more details.\n",
    "\n",
    "Load the ERA5 1996 TC tracks with `huracanpy.load`. \n",
    "These are tracks detected by TempestExtremes in ERA5 for the year 1996 and are embedded within HuracanPy as an example.\n",
    "Here the file extension is '.csv', the function will automatically recognise how to open it.\n",
    "\n",
    "The tracks are loaded as an xarray.Dataset, with one dimension \"record\" corresponding to each point.\n",
    "Variables indicate position in space and time, as well as additional attributes such as maximum wind speed and minimum slp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c9386e3-4cb6-4f51-88c6-aadaec541c31",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "file = huracanpy.example_year_file\n",
    "print(file.split(\"/\")[-1])\n",
    "\n",
    "tracks = huracanpy.load(file)\n",
    "tracks"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13fcb4af-c3f6-4353-b772-630f0fae6852",
   "metadata": {},
   "source": [
    "## Adding info to the tracks\n",
    "HuracanPy has several function to add useful information to the tracks (for full list, see [huracanpy.info](../../api/info.rst)). Here for example we add basin and Saffir-Simpson hurrican scale category information.]\n",
    "\n",
    "### Add basin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c999ed6f-d9af-40e9-be12-f29bc027dd49",
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks = tracks.hrcn.add_basin()\n",
    "tracks.basin"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68a94a8b-8b1a-4436-975b-7d52b35ef6d2",
   "metadata": {},
   "source": [
    "#### Show distribution of TC points among basins\n",
    "(calling seaborn function, works better with categorical labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0b81f23-0c0a-4636-a2b8-3548892f76b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.countplot(tracks.basin)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcca0209-e0ca-4380-be3c-e024bbce8de6",
   "metadata": {},
   "source": [
    "### Add Saffir-Simpson and pressure-based categories\n",
    "Note: in ERA5 data, wind is stored in wind10 in m/s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c087c73e-6caf-4d1c-97af-89a840063330",
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks = tracks.hrcn.add_saffir_simpson_category(wind_name=\"wind10\", wind_units=\"m s-1\")\n",
    "tracks = tracks.hrcn.add_pressure_category(slp_name=\"slp\", slp_units=\"Pa\")\n",
    "tracks[[\"saffir_simpson_category\", \"pressure_category\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01a995ac-1f56-4a19-9f75-50aec2e66b38",
   "metadata": {},
   "source": [
    "#### Show distribution of TC points among categories\n",
    "(using xarray's built-in function)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9f93b37-b280-458b-9330-68fbca67d3f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.arange(-1.5, 5.5 + 1)\n",
    "tracks.saffir_simpson_category.plot.hist(bins=bins, alpha=0.5)\n",
    "tracks.pressure_category.plot.hist(bins=bins, alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6b5f0ce-3124-43b3-8f3d-38647cba0217",
   "metadata": {},
   "source": [
    "## Plotting\n",
    "HuracanPy embeds basic plotting functions, which are mainly meant for having a preliminary look at your data. In particular here we show how to plot the track points themselves, and track density. The [example gallery](../../examples/index.rst) displays nice plots made from HuracanPy and the associated scripts.\n",
    "### Plotting the tracks\n",
    "Plot ERA5 tracks colored by wind intensity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6efa21e4-f03e-4124-8228-6c66cc2a3847",
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks.hrcn.plot_tracks(intensity_var_name=\"wind10\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae30c01a-2721-4d4c-a9a6-6c65d0177a2b",
   "metadata": {},
   "source": [
    "### Plotting track density\n",
    "You can plot the track density directly with `plot_density`, which is based on a simple 2D histogram of TC points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f41cec0d-ae16-43b8-b0e7-037ae6376af4",
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks.hrcn.plot_density()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2953228-eeec-4926-a9cb-19e8e814dd19",
   "metadata": {},
   "source": [
    "You can also get the underlying density matrix with `get_density` and then use it to make you own plots in your favourite way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c7368a1-e5fa-41ab-8ad2-d068a9c7cc60",
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks.hrcn.get_density()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11280ed3-0270-4c71-a20a-099a773eec54",
   "metadata": {},
   "source": [
    "### Plotting genesis points\n",
    "`get_gen_vals` allows you to subset only the genesis points in an efficient way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9239e0ed-5c07-4531-8585-c3f1cc5618f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "gen_points = tracks.hrcn.get_gen_vals()\n",
    "gen_points"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11215983-84ee-48c6-941b-20679cc06bae",
   "metadata": {},
   "source": [
    "If you use `plot_tracks` on these, you can display only the genesis points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbfbe423-976c-464e-b255-bca2aadf2eb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "gen_points.hrcn.plot_tracks()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef3b9157-c798-41da-958e-f14021c9afdc",
   "metadata": {},
   "source": [
    "## Compute statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43219cdb-5e4f-4ce1-9f37-7cd05355d280",
   "metadata": {},
   "source": [
    "### Number of cyclones\n",
    "Count number of unique track ids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0891efaf-8ff1-4751-93a8-190c485b29e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks.track_id.hrcn.nunique()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbde1b28-0840-4d87-8103-0a6128bdd1ae",
   "metadata": {},
   "source": [
    "### Cyclone duration & TC days\n",
    "Get the duration for each track.\n",
    "The result is an `xarray.Dataset` with \"track_id\" as a coordinate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e7c6c67-cc01-404d-843c-c3b63dc62a27",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "TC_duration = tracks.hrcn.get_track_duration()\n",
    "TC_duration"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ed0039a-8a7e-447b-8cc1-c18f402129d3",
   "metadata": {},
   "source": [
    "Compute the total number of TC days.\n",
    "Sum all the durations (and divide by 24 because durations are in hours)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2f170ad-fc36-4b74-a15d-edf1a6f9911e",
   "metadata": {},
   "outputs": [],
   "source": [
    "TC_duration.sum() / 24"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7216f5f7-c943-4cdc-951e-48027ec95e16",
   "metadata": {},
   "source": [
    "### Cyclone Intensity"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afb889aa-3b27-43de-8a3f-6d1958cf3b84",
   "metadata": {},
   "source": [
    "There are two ways to obtain the lifetime maximum intensity (LMI) of each tracks\n",
    "\n",
    "1. Use `get_apex_vals`, which returns the subset of points only at specified LMI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9db8f886-6ce1-4221-bfa8-9873a82d4464",
   "metadata": {},
   "outputs": [],
   "source": [
    "lmi_points = tracks.hrcn.get_apex_vals(var_name=\"wind10\")\n",
    "lmi_points"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c444bdb9-9d11-4b72-9a10-0ccdf0559b31",
   "metadata": {},
   "source": [
    "2. Compute lifetime maximum intensity per track with xarray's `groupby`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fde2d7e-e691-4882-a9bd-d72db83ba1d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "lmi_wind = tracks.wind10.groupby(tracks.track_id).max()\n",
    "lmi_wind"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "665bcca8-2bd8-4c45-8a2b-0648117ddba9",
   "metadata": {},
   "source": [
    "You can then plot the LMI distribution using xarray's built-in plot function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "662bd037-0c1c-4256-a765-65f17d6119b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "lmi_wind.plot.hist()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95fa6ade-b25d-4c1a-b655-a996fe0fb108",
   "metadata": {},
   "source": [
    "### ACE\n",
    "Accumulated cyclone energy (ACE) is a commonly used measure of cyclone activity that combines the energy and duration of cyclones.\n",
    "\n",
    "#### Compute ACE for each individual track point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef7b7f31-4799-459a-a65d-8d0828da9074",
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks = tracks.hrcn.add_ace(wind_name=\"wind10\", wind_units=\"m s-1\")\n",
    "tracks.ace"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70c17332-47b6-4c03-926b-06eaefb67098",
   "metadata": {},
   "source": [
    "#### Compute total ACE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e88dc197-8167-437a-9d91-d9a5cd219b68",
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks.ace.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b772c181-dd5f-4ad4-b8bb-bb0a23120611",
   "metadata": {},
   "source": [
    "## Compositing lifecycle\n",
    "Add time from apex variable to be able to superimpose all the tracks centered on apex.\n",
    "Here we use minimum pressure as the apex point\n",
    "\n",
    "### Add time from minimum pressure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "616d40ca-9e14-4f44-a21b-3903dd4c5026",
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks = tracks.hrcn.add_time_from_apex(intensity_var_name=\"slp\", stat=\"min\")\n",
    "tracks.time_from_apex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ad10480-b4b8-4749-a98d-811dbb3d93aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks.time_from_apex / np.timedelta64(1, \"h\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0241bc41-79ad-4db0-a887-02fb91546a11",
   "metadata": {},
   "source": [
    "### Plot composite SLP lifecycle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "268f293b-601d-408b-a987-ffb639784415",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert time_from_apex to hours\n",
    "tracks[\"time_from_apex\"] = tracks.time_from_apex / np.timedelta64(1, \"h\")\n",
    "\n",
    "# Use xarray's where to mask points too far away from apex (48 hours away)\n",
    "tracks_close_to_apex = tracks.where(np.abs(tracks.time_from_apex) <= 48, drop=True)\n",
    "\n",
    "# Seaborn lineplot allows for drawing composites with uncertainty range\n",
    "# x-axis is time from apex\n",
    "# y-axis is slp, converted to hPa\n",
    "sns.lineplot(x=tracks_close_to_apex.time_from_apex, y=tracks_close_to_apex.slp / 100)"
   ]
  }
 ],
 "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}