{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d4b7b6b-1126-49ec-b8b4-337321000fe8",
   "metadata": {
    "editable": true,
    "nbsphinx": "hidden",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "np.set_printoptions(threshold=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fec6b138-477f-4332-8582-682bb3c849c4",
   "metadata": {},
   "source": [
    "# Unit-Aware Calculations\n",
    "Where possible, HuracanPy functions are doing unit-aware calculations.\n",
    "This is achieved using functionality from [metpy and the metpy accessor](https://unidata.github.io/MetPy/latest/tutorials/xarray_tutorial.html#units).\n",
    "\n",
    "Typically, for metpy unit-aware calculations to work, the units must be included in the attributes of the variable. e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a540f017-6d7d-4a64-95b2-9289c434ce43",
   "metadata": {},
   "outputs": [],
   "source": [
    "import huracanpy\n",
    "\n",
    "tracks = huracanpy.load(huracanpy.example_csv_file)\n",
    "\n",
    "# Add units to some variables\n",
    "tracks.slp.attrs[\"units\"] = \"Pa\"\n",
    "tracks.wind10.attrs[\"units\"] = \"m s-1\"\n",
    "tracks.lon.attrs[\"units\"] = \"degrees\"\n",
    "tracks.lat.attrs[\"units\"] = \"degrees\"\n",
    "\n",
    "# Use a unit-aware huracanpy calculation\n",
    "dv_dt = tracks.hrcn.get_rate(var_name=\"wind10\")\n",
    "print(dv_dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bac8262-6993-4f0f-b24a-533dc47ef98c",
   "metadata": {},
   "source": [
    "## MetPy functionality\n",
    "\n",
    "When the track data has units in the attributes you can use functionalities from the `metpy` accessor\n",
    "(see [the metpy docs](https://unidata.github.io/MetPy/latest/tutorials/xarray_tutorial.html#units) for more details).\n",
    "\n",
    "### Unit Conversion\n",
    "For example, if we want to convert the intensification to more readable units we can\n",
    "use the metpy unit conversion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a078a1a9-e5aa-4f9c-82af-ad8036edb1f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(dv_dt.metpy.convert_units(\"m s-1 day-1\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdb152f6-9c78-47f9-a8fa-087d2f1ae50c",
   "metadata": {},
   "source": [
    "### dequantify and quantify\n",
    "\n",
    "Note that, the returned `Dataarray` contains a `pint.Quantity` instead of a numpy array.\n",
    "This is a numpy-like array that enables the unit-aware calculations. However it can\n",
    "sometimes result in strange behaviour. The `metpy` accessor provides a method for\n",
    "undoing this (and HuracanPy functions will do this automatically)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84bdfcc6-a8b3-4f3a-983d-670aa19da2e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(dv_dt.metpy.convert_units(\"m s-1 day-1\").metpy.dequantify())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a085500-aaaf-46d1-a8ff-ed5d732e91aa",
   "metadata": {},
   "source": [
    "Using `dequantify` has reverted to a numpy array with the units in the attributes\n",
    "\n",
    "Be aware that fairly trivial operations will lose these units if they are not done\n",
    "in a unit-aware way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddebc148-c7be-466c-8038-777c0c524c2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(dv_dt + 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47b834f8-0f9f-4a60-833f-e3c696698903",
   "metadata": {},
   "source": [
    "If you want to do your own unit-aware calculations, you will need to use the `metpy` methods, e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2cb6b68-6951-46bd-9b6b-e88eb40c237d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from metpy.units import units\n",
    "\n",
    "print((dv_dt.metpy.quantify() + 0 * units(\"m s-2\")).metpy.dequantify())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6a7f241-5210-4e07-866a-bdeab59e8cd0",
   "metadata": {},
   "source": [
    "## HuracanPy functionality\n",
    "The following functions support unit-aware calculations\n",
    "\n",
    "- [huracanpy.calc.delta](../api/_autosummary/huracanpy.calc.delta.rst)\n",
    "- [huracanpy.calc.rate](../api/_autosummary/huracanpy.calc.rate.rst)\n",
    "- [huracanpy.calc.distance](../api/_autosummary/huracanpy.calc.distance.rst)\n",
    "- [huracanpy.calc.translation_speed](../api/_autosummary/huracanpy.calc.translation_speed.rst)\n",
    "- [huracanpy.info.category](../api/_autosummary/huracanpy.info.category.rst)\n",
    "- [huracanpy.tc.ace](../api/_autosummary/huracanpy.tc.ace.rst)\n",
    "- [huracanpy.tc.pace](../api/_autosummary/huracanpy.tc.pace.rst)\n",
    "- [huracanpy.tc.beta_drift](../api/_autosummary/huracanpy.tc.beta_drift.rst)\n",
    "\n",
    "### Assumed units\n",
    "Unless you have a netCDF file which includes units, or manually add the units, your track variables will not have units.\n",
    "To account for this, HuracanPy functions that use units will assume that your variables have default units with the following convention\n",
    "\n",
    "- Distance in metres\n",
    "- Speed in metres per second\n",
    "- Pressure in hectoPascals (with checks on reasonable values)\n",
    "- Latitude/longitude in degrees\n",
    "\n",
    "e.g. calculating accumulated cyclone energy (ACE), the input is assumed to be in metres per second and converted to knots\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d67ed5f1-e6fd-415e-9afc-d37b0d69df4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reload the tracks without units\n",
    "tracks = huracanpy.load(huracanpy.example_csv_file)\n",
    "\n",
    "# Calculate ACE without specifying units\n",
    "ace = tracks.hrcn.get_ace(wind_name=\"wind10\", sum_by=\"track_id\")\n",
    "\n",
    "# Result has units of knots**2\n",
    "print(ace)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8f3dbfe-b71b-403a-b6e9-e837a74fab42",
   "metadata": {},
   "source": [
    "The units can also be specified directly as a function argument"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e2057c4-d08f-474b-9b00-28601c561c46",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tell the function that the input wind is in metres per second (same result)\n",
    "ace = tracks.hrcn.get_ace(wind_name=\"wind10\", sum_by=\"track_id\", wind_units=\"m s-1\")\n",
    "print(ace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc79a277-6b8b-4185-a3f9-d6b2b0cf4138",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tell the function that the input wind is in kilometres per second for a wildly\n",
    "# inaccurate result\n",
    "ace = tracks.hrcn.get_ace(wind_name=\"wind10\", sum_by=\"track_id\", wind_units=\"km s-1\")\n",
    "print(ace)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3451463f-c5ee-4f52-9b29-83f94ad6cb6e",
   "metadata": {},
   "source": [
    "In some cases, the default will check that the values have a reasonable magnitude for\n",
    "the assumed units. e.g. calculated tropical-cyclone category from pressure the units\n",
    "of the track are in Pascals but assumed to be hectoPascals if unspecified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fb159f4-6016-475e-8ea4-6043b77aa79a",
   "metadata": {},
   "outputs": [],
   "source": [
    "category = tracks.hrcn.get_pressure_category()\n",
    "print(category)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "170b45f2-1a7e-4739-8c2d-c1fb587dca21",
   "metadata": {},
   "source": [
    "If you want to override this check, you can specify the units, either through the attributes or the function arguments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3aa2f03e-b883-4248-b827-4f5d758f25aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "category = tracks.hrcn.get_pressure_category(slp_units=\"hPa\")\n",
    "print(category)"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}