{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7c2a896d-f8ff-43f9-92b5-3e3c66b01315",
   "metadata": {},
   "source": [
    "# Pressure-Wind Relation\n",
    "\n",
    "In this notebook we show three examples of plotting using `huracanpy.plot.pressure_wind_relation`.\n",
    "The first two plots use the JTWC version of IBTrACS in huracanpy because it has more\n",
    "consistent definitions of winds\n",
    "- The first figure shows a plot of maximum intensity for all the TCs in the JTWC dataset\n",
    "- The second figure plots separate colors for each basin. This shows the consistency of North Atlantic and East Pacific basins, in which JTWC define wind speeds a 1-minute sustained winds, whereas other basins use longer time periods\n",
    "- The third figure is a repeat of the second figure, but using the WMO version of IBTrACS. In this case the wind definitions are different for each basin and we see less consistent pressure-wind relations across basins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2dab350-3564-48ec-84fd-0fa6a9182d0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import huracanpy\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3caa3690-7b25-4de2-8c7c-8573d4ae571b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use JTWC for more consistent definitions of wind speed\n",
    "tracks = huracanpy.load(source=\"ibtracs\", ibtracs_subset=\"jtwc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6759741e-33cf-421b-88b1-f1c50b44e1d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find each track at maximum intensity as a tropical cyclone\n",
    "# Subset to tropical part of tracks without nans in pressure and wind\n",
    "tracks = tracks.isel(\n",
    "    record=np.where(\n",
    "        np.isin(tracks.status, [\"TD\", \"TS\", \"TY\", \"ST\", \"HU\", \"HR\"])\n",
    "        & ~np.isnan(tracks.slp)\n",
    "        & ~np.isnan(tracks.wind)\n",
    "    )[0]\n",
    ")\n",
    "\n",
    "# Get maximum intensity (by wind)\n",
    "tracks_mi = tracks.hrcn.get_apex_vals(\"wind\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c446f6e9-f3db-4c16-b522-8e7b5927c9e0",
   "metadata": {},
   "source": [
    "## Single plot\n",
    "Note that we are setting `wind_units=\"knots\"` and `wind_convention=\"1min\"` to match the\n",
    "data in the IBTrACS JTWC subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed5fdaf6-d423-44f8-9cf7-0fa2d570036c",
   "metadata": {},
   "outputs": [],
   "source": [
    "huracanpy.plot.pressure_wind_relation(\n",
    "    tracks_mi.slp, tracks_mi.wind, wind_units=\"knots\", wind_convention=\"1min\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab93c96a-5b1d-4b2f-bf34-b32cfe5f57a8",
   "metadata": {},
   "source": [
    "## Overlaying Multiple Plots\n",
    "The `pressure_wind_relation` function returns the grid\n",
    "(generated by :py:class:`seaborn.JointGrid`) and the bins used for pressure and wind.\n",
    "To overlay more data onto the same figure pass the `grid` back to\n",
    "`pressure_wind_relation`. The example below also passes the bins back so that the\n",
    "histograms use the same binning.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "251013a9-2e46-4fd1-ab5c-3464cfd9c1de",
   "metadata": {},
   "outputs": [],
   "source": [
    "grid, bins_pressure, bins_wind = None, None, None\n",
    "for basin, _tracks in tracks_mi.groupby(\"basin\"):\n",
    "    try:\n",
    "        grid, bins_pressure, bins_wind = huracanpy.plot.pressure_wind_relation(\n",
    "            _tracks.slp,\n",
    "            _tracks.wind,\n",
    "            wind_units=\"knots\",\n",
    "            wind_convention=\"1min\",\n",
    "            bins_pressure=bins_pressure,\n",
    "            bins_wind=bins_wind,\n",
    "            grid=grid,\n",
    "            label=basin,\n",
    "        )\n",
    "    except np.linalg.LinAlgError:\n",
    "        print(f\"Basin {basin} quadratic fit not possible\")\n",
    "\n",
    "grid.ax_joint.legend(ncol=2)\n",
    "grid.fig.suptitle(\"Pressure-wind relation by basin. IBTrACS-JTWC\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74416d46-fd3f-41d7-9ffb-c273191a1d75",
   "metadata": {},
   "source": [
    "## Using the IBTrACS WMO subset\n",
    "Here we still set `wind_units=\"knots\"`, but set `wind_convention=\"10min\"`. This is the\n",
    "default for the wind convention. The WMO subset has a mix of different conventions for\n",
    "the wind speeds, and you can see from the figure the pressure-wind fits are more\n",
    "separated than the JTWC versions above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bf868e3-a389-4c41-b2d3-83d9736f4d5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Re-do by basin plot for IBTrACS-WMO where each basin uses the local convention\n",
    "# for wind speeds\n",
    "tracks = huracanpy.load(source=\"ibtracs\", ibtracs_subset=\"wmo\")\n",
    "\n",
    "# Find each track at maximum intensity\n",
    "# No \"status\" variable included with WMO subset, so include all points without nans in\n",
    "# pressure and wind\n",
    "tracks = tracks.isel(record=np.where(~np.isnan(tracks.slp) & ~np.isnan(tracks.wind))[0])\n",
    "\n",
    "# Get maximum intensity (by wind)\n",
    "tracks_mi = tracks.hrcn.get_apex_vals(\"wind\")\n",
    "\n",
    "grid, bins_pressure, bins_wind = None, None, None\n",
    "for basin, _tracks in tracks_mi.groupby(\"basin\"):\n",
    "    # Unlike JTWC, this does not fail for South Atlantic\n",
    "    # This is because there are some extra non-tropical cyclones not filtered out\n",
    "    grid, bins_pressure, bins_wind = huracanpy.plot.pressure_wind_relation(\n",
    "        _tracks.slp,\n",
    "        _tracks.wind,\n",
    "        wind_units=\"knots\",\n",
    "        wind_convention=\"10min\",\n",
    "        bins_pressure=bins_pressure,\n",
    "        bins_wind=bins_wind,\n",
    "        grid=grid,\n",
    "        label=basin,\n",
    "    )\n",
    "\n",
    "\n",
    "# Adjust the xlimits\n",
    "# This needs to be done because the fit for the South Atlantic is based on three tracks\n",
    "# and pushes the xlimit to unreasonable values by default\n",
    "grid.ax_joint.set_xlim(0, tracks_mi.wind.values.max())\n",
    "grid.ax_joint.legend(ncol=2)\n",
    "grid.fig.suptitle(\"Pressure-wind relation by basin. IBTrACS-WMO\")\n",
    "plt.tight_layout()"
   ]
  }
 ],
 "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
}