{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d373c1cc-f372-4076-a9a0-df6654063333",
   "metadata": {},
   "source": [
    "# Example workflows\n",
    "In this notebook, we demonstrate the type of analysis workflow you can build using HuracanPy. To know more about the available options and functions, please go through the [huracanpy.load guide](index.rst), and/or browse the [API documentation](../api/index.rst). The [example gallery](../examples/index.rst) gallery also provides useful demonstrations.\n",
    "\n",
    "We show three workflows: \n",
    "\n",
    "1. Studying a specific cyclone\n",
    "2. Studying a set of tracks\n",
    "3. Comparing a set of detected/modelled tracks to an observationnal reference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98f00647-d338-4bea-b551-511758fa39c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import huracanpy\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8ec02c1-fd58-46b6-9211-e3ac1b36a2fd",
   "metadata": {},
   "source": [
    "## 1. Studying a specific cyclone\n",
    "In this example, we want to study hurricane Wilma (the deepest Atlantic hurricane on record)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b458a5ab-e8dc-4a94-9a9d-51078981e617",
   "metadata": {},
   "source": [
    "### 1a. Load IBTrACS and subset the specific hurricane\n",
    "Two subsets of IBTrACS are embedded within HuracanPy: WMO and JTWC. \n",
    "You can also retrieve the full and last IBTrACS file from the online website. \n",
    "Default behavior is loading the embedded WMO subset. For more information, see [huracanpy.load guide](load.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "999c3ab4-a444-4b4a-bbe6-c36c9b19fade",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we load the WMO subset. This raises a warning that reminds you of the main caveats.\n",
    "ib = huracanpy.load(source=\"ibtracs\")\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.\n",
    "ib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "faf70c53-ce52-4bf1-8b02-7e263d05751d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Wilma corresponds to index 2005289N18282, so we subset this storm. There are two ways of doing this:\n",
    "# 1. Use warray's where\n",
    "Wilma = ib.where(ib.track_id == \"2005289N18282\", drop=True)\n",
    "# 2. Use huracanpy's sel_id method (more efficient and shorter, but does the same thing)\n",
    "# Note: the `.hrcn` is called an accessor, and allows you to call HuracanPy functions as methods on the xarray objects.\n",
    "Wilma = ib.hrcn.sel_id(\"2005289N18282\")\n",
    "# The Wilma object contains only the data for Wilma:\n",
    "Wilma"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee08d4f1-e18f-4820-810d-0005b92f9105",
   "metadata": {},
   "source": [
    "### 1b. Add category info\n",
    "You can add the Saffir-Simpson and/or the pressure category of Wilma to the tracks (for full list of available info, see [huracanpy.info](../api/info.rst))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99f5e05f-77f1-48af-85e4-e4dbeea3ea7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add Saffir-Simpson Category\n",
    "Wilma = Wilma.hrcn.add_saffir_simpson_category(wind_name=\"wind\", wind_units=\"knots\")\n",
    "Wilma.saffir_simpson_category  # This is stored in the `saffir_simpson_category` variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59c13d3d-bdaf-4f32-a2b8-e9fb9912c7a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add pressure category\n",
    "Wilma = Wilma.hrcn.add_pressure_category(slp_name=\"slp\")\n",
    "Wilma.pressure_category"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8e42735-daac-4a59-9d3e-92db5a6b46e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note: Most of the accessor methods have a get_* and an add_* version.\n",
    "# get_ return the values of what you ask for as a DataArray, while add_ adds it directly to the dataset with a default name.\n",
    "# In the previous case, we could have called get_pressure_category\n",
    "Wilma.hrcn.get_pressure_category(slp_name=\"slp\")\n",
    "# we could then save it as a variable, and potentially add it to the dataset separately"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58a151fe-4f12-48f9-8f7b-92e8e3eda692",
   "metadata": {},
   "source": [
    "### 1c. Plot the track and its evolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3a5cca2-4899-487e-abea-1244ea1aa4b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the track on a map, colored by Saffir-Simpson category\n",
    "Wilma.hrcn.plot_tracks(\n",
    "    intensity_var_name=\"saffir_simpson_category\", scatter_kws={\"palette\": \"turbo\"}\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "318f54d6-f974-4ca1-8876-1b46978a0043",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot intensity time series using matplotlib\n",
    "fig, axs = plt.subplots(2, sharex=True)\n",
    "axs[0].plot(Wilma.time, Wilma.wind)\n",
    "axs[1].plot(Wilma.time, Wilma.slp)\n",
    "axs[0].set_ylabel(\"Wind / kn\")\n",
    "axs[1].set_ylabel(\"SLP / hPa\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "997d6782-d0ea-4274-98d9-bd1930a259c4",
   "metadata": {},
   "source": [
    "### 1d. Calculate properties "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcdec2df-6d47-4b3c-9c94-ee6d7237694c",
   "metadata": {},
   "source": [
    "#### Duration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe2d0d35-2782-4080-b0e9-acc67a3c5d70",
   "metadata": {},
   "outputs": [],
   "source": [
    "Wilma.hrcn.get_track_duration()  # Note duration is in h"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98b4cd18-d260-491b-96ed-061991d1c141",
   "metadata": {},
   "source": [
    "#### ACE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60aafd0b-1f7f-46e2-bdec-933c44d60863",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Compute ACE for each point\n",
    "Wilma = Wilma.hrcn.add_ace(wind_units=\"knots\")\n",
    "Wilma.ace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b06d515b-c77c-4296-81d3-28e86ec38816",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot cumulated ACE\n",
    "plt.plot(Wilma.time, Wilma.ace.cumsum())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0362378e-149e-4e02-975a-7e123399ff51",
   "metadata": {},
   "source": [
    "#### Translation speed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcf8bcc5-2059-4aa7-8efc-94f52a8a1144",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute translation speed\n",
    "Wilma = Wilma.hrcn.add_translation_speed()\n",
    "# Plot translation speed against latitude\n",
    "plt.plot(Wilma.lat, Wilma.translation_speed)\n",
    "plt.xlabel(\"Latitude / °\")\n",
    "plt.ylabel(\"Translation speed / m/s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4eee8b7e-7b71-4528-b411-4b91db62ead5",
   "metadata": {},
   "source": [
    "#### Intensification rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ea0f5ef-d7f4-4a69-b37e-b7a437713c79",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add intensification rate in wind and pressure\n",
    "Wilma = Wilma.hrcn.add_rate(var_name=\"wind\")\n",
    "Wilma = Wilma.hrcn.add_rate(var_name=\"slp\")\n",
    "# NB: The rates will be in unit/s, where unit is the unit of the variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2baf0771-a278-4b93-ad1c-81ec83eabcbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot intensity time series\n",
    "fig, axs = plt.subplots(2, sharex=True)\n",
    "axs[0].plot(Wilma.time, Wilma.rate_wind * 3600)  # Convert to kn/h\n",
    "axs[1].plot(Wilma.time, Wilma.rate_slp * 3600)  # Convert to hPa/h\n",
    "axs[0].set_ylabel(\"kn/h\")\n",
    "axs[1].set_ylabel(\"hPa/h\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b06ec72-4a81-421e-8628-b3684f6198b7",
   "metadata": {},
   "source": [
    "## 2. Studying a set of tracks"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c39bbf8-757d-40ae-b80a-18d3af3d7fd2",
   "metadata": {},
   "source": [
    "### 2a. Loading data\n",
    "Here we show an example with a csv file that is embedded within HuracanPy for example. HuracanPy supports many track files format, see [huracanpy.load guide](load.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c9386e3-4cb6-4f51-88c6-aadaec541c31",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading the ERA5 1996 TC tracks\n",
    "## The tracks detected by TempestExtremes in ERA5 for the year 1996 are embedded within the package as an example.\n",
    "file = huracanpy.example_year_file\n",
    "print(file)\n",
    "## Load the tracks with huracanpy.load.\n",
    "## Here the file extension is '.csv', the function will automatically recognise how to open it.\n",
    "tracks = huracanpy.load(file)\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.\n",
    "tracks"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13fcb4af-c3f6-4353-b772-630f0fae6852",
   "metadata": {},
   "source": [
    "### 2b. 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 SSHS category information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c999ed6f-d9af-40e9-be12-f29bc027dd49",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add basin\n",
    "tracks = tracks.hrcn.add_basin()  # Add basin attribute\n",
    "tracks.basin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0b81f23-0c0a-4636-a2b8-3548892f76b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show distribution of TC points among basins (calling seaborn function, works better with categorical labels)\n",
    "sns.countplot(tracks.basin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c087c73e-6caf-4d1c-97af-89a840063330",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add SSHS and pressure categories\n",
    "tracks = tracks.hrcn.add_saffir_simpson_category(wind_name=\"wind10\", wind_units=\"m s-1\")\n",
    "tracks = tracks.hrcn.add_pressure_category(\n",
    "    slp_name=\"slp\",\n",
    ")\n",
    "## (In ERA5 data, wind is stored in wind10 in m/s)\n",
    "tracks[[\"saffir_simpson_category\", \"pressure_category\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9f93b37-b280-458b-9330-68fbca67d3f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show distribution of TC points among categories (using xarray's built-in function)\n",
    "tracks.saffir_simpson_category.plot.hist(\n",
    "    bins=[-1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5], alpha=0.5\n",
    ")\n",
    "tracks.pressure_category.plot.hist(\n",
    "    bins=[-1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5], alpha=0.5\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6b5f0ce-3124-43b3-8f3d-38647cba0217",
   "metadata": {},
   "source": [
    "### 2c. 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. You can learn more in the [huracanpy.plot guide](plot.ipynb). The [example gallery](../examples/index.rst)  also displays nice plots made from HuracanPy and the associated scripts.\n",
    "#### Plotting the tracks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6efa21e4-f03e-4124-8228-6c66cc2a3847",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot ERA5 tracks colored by wind intensity\n",
    "tracks.hrcn.plot_tracks(\n",
    "    intensity_var_name=\"wind10\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae30c01a-2721-4d4c-a9a6-6c65d0177a2b",
   "metadata": {},
   "source": [
    "#### Plotting track density"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f41cec0d-ae16-43b8-b0e7-037ae6376af4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can plot the track density directly with `plot_density`, which is based on a simple 2D histogram of TC points\n",
    "tracks.hrcn.plot_density()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c7368a1-e5fa-41ab-8ad2-d068a9c7cc60",
   "metadata": {},
   "outputs": [],
   "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\n",
    "tracks.hrcn.get_density()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11280ed3-0270-4c71-a20a-099a773eec54",
   "metadata": {},
   "source": [
    "#### Plotting genesis points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9239e0ed-5c07-4531-8585-c3f1cc5618f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# `get_gen_vals` allows you to subset only the genesis points in an efficient way\n",
    "gen_points = tracks.hrcn.get_gen_vals()\n",
    "gen_points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbfbe423-976c-464e-b255-bca2aadf2eb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you use `plot_tracks` on these, you can display only the genesis points.\n",
    "gen_points.hrcn.plot_tracks()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef3b9157-c798-41da-958e-f14021c9afdc",
   "metadata": {},
   "source": [
    "### 2d. Compute statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43219cdb-5e4f-4ce1-9f37-7cd05355d280",
   "metadata": {},
   "source": [
    "#### Number of cyclones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0891efaf-8ff1-4751-93a8-190c485b29e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks.track_id.hrcn.nunique()  # Count number of unique track ids"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbde1b28-0840-4d87-8103-0a6128bdd1ae",
   "metadata": {},
   "source": [
    "#### Cyclones duration & TC days"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e7c6c67-cc01-404d-843c-c3b63dc62a27",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Get the duration for each track\n",
    "TC_duration = tracks.hrcn.get_track_duration()\n",
    "TC_duration  # xarray.Dataset with track_id as dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2f170ad-fc36-4b74-a15d-edf1a6f9911e",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Compute the total number of TC days\n",
    "## Sum all the durations (and divide by 24 because durations are in hours)\n",
    "TC_duration.sum() / 24"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7216f5f7-c943-4cdc-951e-48027ec95e16",
   "metadata": {},
   "source": [
    "#### Cyclone Intensity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9db8f886-6ce1-4221-bfa8-9873a82d4464",
   "metadata": {},
   "outputs": [],
   "source": [
    "# There are two ways to obtain the lifetime maximum intensity (LMI) of each tracks\n",
    "## 1. Use `get_apex_vals`, which return the subset of points only as specified LMI\n",
    "tracks.hrcn.get_apex_vals(varname=\"wind10\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fde2d7e-e691-4882-a9bd-d72db83ba1d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "## 2. Compute lifetime maximum intensity per track with xarray's groupby\n",
    "LMI_wind = tracks.wind10.groupby(tracks.track_id).max()\n",
    "LMI_wind  # xarray.Dataset with track_id as dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "662bd037-0c1c-4256-a765-65f17d6119b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can then plot the LMI distribution using xarray's built-in plot function.\n",
    "LMI_wind.plot.hist()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82f009c6-8d92-4cdd-8af4-c72752848d07",
   "metadata": {},
   "source": [
    "#### ACE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef7b7f31-4799-459a-a65d-8d0828da9074",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute ACE for each point\n",
    "tracks = tracks.hrcn.add_ace(wind_name=\"wind10\", wind_units=\"m s**-1\")\n",
    "tracks.ace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e88dc197-8167-437a-9d91-d9a5cd219b68",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Compute total ACE\n",
    "tracks.ace.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b772c181-dd5f-4ad4-b8bb-bb0a23120611",
   "metadata": {},
   "source": [
    "### 2e. Compositing lifecycle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "616d40ca-9e14-4f44-a21b-3903dd4c5026",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add time from apex variable to be able to superimpose all the tracks centered on apex\n",
    "tracks = tracks.hrcn.add_time_from_apex(\n",
    "    intensity_var_name=\"slp\", stat=\"min\"\n",
    ")  # Add time from minimum pressure\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": "code",
   "execution_count": null,
   "id": "268f293b-601d-408b-a987-ffb639784415",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot composite SLP lifecycle\n",
    "## Convert time_from_apex to hours\n",
    "tracks[\"time_from_apex\"] = tracks.time_from_apex / np.timedelta64(1, \"h\")\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",
    "## Seaborn lineplot allows for drawing composites with uncertainty range\n",
    "sns.lineplot(\n",
    "    x=tracks_close_to_apex.time_from_apex,  # x-axis is time from apex\n",
    "    y=tracks_close_to_apex.slp / 100,\n",
    ")  # y-axis is slp, converted to hPa"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aef8b2fb-d9b0-47e8-a70e-1024c7e6e2b1",
   "metadata": {},
   "source": [
    "## 3. Comparing two datasets\n",
    "In this part, we compare the set of 1996 tracks above to IBTrACS which we use as reference. \n",
    "To start with, note that for all that was shown above, you can superimpose several sets and therefore compare several sources/models/trackers/etc. Below we show specific functions for matching tracks and computing detection scores."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53a94bdb-6d94-47c9-9296-5826986ea0cb",
   "metadata": {},
   "source": [
    "### 3a. Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efcb9ebc-bf0e-4432-b803-e53925b4fcdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remember IBTrACS is stored in the `ib` object from the first part above.\n",
    "# Here we subset the 1996 tracks with xarray's where method:\n",
    "ib_1996 = ib.where(ib.time.dt.year == 1996, drop=True)\n",
    "ib_1996"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd786725-709d-4d66-9fc2-3f964117fb3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The tracks from ERA5 are stored in `tracks`. For clarity, we name it `ERA5` from now:\n",
    "ERA5 = tracks.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9c764fc-9f16-49b1-abcb-08faed69e8e3",
   "metadata": {},
   "source": [
    "### 3b. Superimposing several sets on one plot\n",
    "To start with, note that for all that was shown above, you can superimpose several sets and therefore compare several sources/models/trackers/etc. Here we only show one example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73f68f44-7745-4852-afbc-b0e212903a07",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute LMI for both sets\n",
    "LMI_wind_ib = ib_1996.wind.groupby(ib_1996.track_id).max()\n",
    "LMI_wind_ib = LMI_wind_ib / 1.94  # Convert kn to m/s\n",
    "LMI_wind_ERA5 = ERA5.wind10.groupby(ERA5.track_id).max()\n",
    "# Plot both histograms\n",
    "LMI_wind_ib.plot.hist(\n",
    "    bins=[10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65], color=\"k\", label=\"IBTrACS\"\n",
    ")\n",
    "LMI_wind_ERA5.plot.hist(\n",
    "    bins=[10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65], label=\"ERA5\", alpha=0.8\n",
    ")\n",
    "# Labels\n",
    "plt.legend()\n",
    "plt.xlabel(\"Lifetime maximum wind speed / m/s\")\n",
    "plt.ylabel(\"Number of tracks\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da63e862-5675-4277-bc20-e7e9d9dd5168",
   "metadata": {},
   "source": [
    "### 3c. Matching tracks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c66913f0-cbe9-41b2-a740-645e74942788",
   "metadata": {},
   "outputs": [],
   "source": [
    "matches = huracanpy.assess.match([ERA5, ib_1996], names=[\"ERA5\", \"IBTrACS\"])\n",
    "matches  # each row is a pair of tracks that matched, with both ids, the number of time steps and the mean distance between the tracks over their matching period."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80c90c22-a6ad-4219-9efa-9ba062d7f57d",
   "metadata": {},
   "source": [
    "### 3d. Computing scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73d70abf-7ba0-472d-bf25-f5c47d1775fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Probability of detection (POD) : Proportion of observed tracks that are found in ERA5.\n",
    "huracanpy.assess.pod(matches, ref=ib_1996, ref_name=\"IBTrACS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31c17f9d-e830-4dd3-8084-f09430f58b2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# False alarm rate (FAR) : Proportion of detected tracks that were not observed\n",
    "huracanpy.assess.far(matches, detected=ERA5, detected_name=\"ERA5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01fc6748-fef2-4cb4-9b77-b9d91b122cbe",
   "metadata": {},
   "source": [
    "### 3e. Venn diagrams\n",
    "Venn diagrams are a convenient way to show the overlap between two datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d719ad71-830c-472d-906e-bfe044b390c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "huracanpy.plot.venn([ERA5, ib_1996], matches, labels=[\"ERA5\", \"IBTrACS\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3a2dcca-7026-4b1e-99f6-e867e2e4150c",
   "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.11.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}