nb1-calculations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate additional variables and store them in the COAMPS-TC, drifter, and IBTrACS datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import Tuple\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy\n",
    "import xarray as xr\n",
    "\n",
    "from configure import read_stored_variable\n",
    "from src import buoy_accessor, coamps, kinematics, geodesy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read stored variables from `io.pynb`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "%run 'nb0-datasets.ipynb'\n",
    "%run -i configure.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df = read_stored_variable('drifter_df')\n",
    "ian_spotter_coamps_df = read_stored_variable('ian_spotter_coamps_df')\n",
    "fiona_spotter_coamps_df = read_stored_variable('fiona_spotter_coamps_df')\n",
    "coamps_ds = read_stored_variable('coamps_ds')\n",
    "ibtracs_df = read_stored_variable('ibtracs_df')\n",
    "bathymetry_ds = read_stored_variable('bathymetry_ds')\n",
    "time_slice_full = read_stored_variable('time_slice_full')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## COAMPS-TC calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Local coordinate system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each time step, convert the longitude and latitude grid to a local Cartesian coordinate system centered at the minimum pressure.  The local coordinate system will be used to colocate the buoy with respect to the hurricane center."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def localize_coamps_grid(\n",
    "    longitude: np.ndarray,\n",
    "    latitude: np.ndarray,\n",
    "    pressure: np.ndarray\n",
    ") -> Tuple[np.ndarray, np.ndarray, float, float]:\n",
    "    \"\"\"Localize a COAMPS-TC grid defined in `longitude` and `latitude` to a\n",
    "    cartesian coordinate system centered at the minimum value of `pressure`.\n",
    "\n",
    "    Args:\n",
    "        longitude (np.ndarray): longitude array with shape (x,)\n",
    "        latitude (np.ndarray):  latitude array with shape (y,)\n",
    "        pressure (np.ndarray): pressure array with shape (x,y)\n",
    "\n",
    "    Returns:\n",
    "        Tuple: meshed x- and y-grids with shape (x,y) and the coordinates\n",
    "        at the xy-origin (coordinates of the the minimum pressure).\n",
    "    \"\"\"\n",
    "    longitude_grid, latitude_grid = np.meshgrid(longitude, latitude)\n",
    "    origin_indices = np.unravel_index(pressure.argmin(), pressure.shape)\n",
    "    origin = (longitude_grid[origin_indices], latitude_grid[origin_indices])\n",
    "    x_grid, y_grid = geodesy.lonlat_to_xy(longitude_grid, latitude_grid, origin)\n",
    "    return x_grid, y_grid, origin[0], origin[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "coamps_ds['x'], coamps_ds['y'], coamps_ds['x-center'], coamps_ds['y-center']\\\n",
    "    = xr.apply_ufunc(\n",
    "    localize_coamps_grid,\n",
    "    coamps_ds.longitude,\n",
    "    coamps_ds.latitude,\n",
    "    coamps_ds.mslp,\n",
    "    input_core_dims=[[\"longitude\"], [\"latitude\"], [\"latitude\", \"longitude\"]],\n",
    "    output_core_dims=[[\"latitude\", \"longitude\"], [\"latitude\", \"longitude\"], [], []],\n",
    "    vectorize=True\n",
    ")\n",
    "\n",
    "coamps_ds['x'].attrs['long_name'] = 'x-distance from minimum pressure'\n",
    "coamps_ds['x'].attrs['units'] = 'm'\n",
    "coamps_ds['y'].attrs['long_name'] = 'x-distance from minimum pressure'\n",
    "coamps_ds['y'].attrs['units'] = 'm'\n",
    "coamps_ds['x-center'].attrs['long_name'] = 'longitude of minimum pressure'\n",
    "coamps_ds['y-center'].attrs['long_name'] = 'latitude of minimum pressure'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpolate local coordinates onto a regular grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Localizing the longitude and latitude grid to Cartesian coordinates (the previous step) produces irregular xy-grids. To optimize colocation, the wind field is interpolated onto a regular (rectangular) grid and stored in a new Dataset, `coamps_xy_ds`.  **Note:** this may take a minute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def regrid_irregular_field(field_irr, x_irr, y_irr, x_reg, y_reg):\n",
    "    \"\"\"Interpolate a field defined at irregular points onto a rectangular grid.\n",
    "\n",
    "    Args:\n",
    "        field_irr (np.ndarray): field at irregular points with shape (m, n)\n",
    "        x_irr (np.ndarray): irregular x-coordinates with shape (m, n)\n",
    "        y_irr (np.ndarray): irregular y-coordinates with shape (m, n)\n",
    "        x_reg (np.ndarray): regular x-coordinates with shape (x)\n",
    "        y_reg (np.ndarray): regular y-coordinates with shape (y)\n",
    "\n",
    "    Returns:\n",
    "        np.ndarray: field interpolated onto the regular grid with shape (x, y).\n",
    "    \"\"\"\n",
    "    x_grid, y_grid = np.meshgrid(x_reg, y_reg)\n",
    "    points_irr = np.c_[x_irr.ravel(), y_irr.ravel()]\n",
    "    field_reg = scipy.interpolate.griddata(points_irr,\n",
    "                                           field_irr.ravel(),\n",
    "                                           (x_grid, y_grid),\n",
    "                                           method='linear')\n",
    "    return field_reg\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create x and y arrays which will define the rectangular grid (in units of m).\n",
    "# The grid is centered on the storm and aligned with north.  The goal is to\n",
    "# limit distortion and achieve similar resolution to the original grid, 0.025\n",
    "# degrees of latitude (about 2775 m).\n",
    "x_arr = np.arange(-750*10**3, 750*10**3, 2775)  # 0.025 deg\n",
    "y_arr = np.arange(-750*10**3, 750*10**3, 2775)\n",
    "\n",
    "# Initialize the Dataset.  Variables will be added after.\n",
    "coamps_xy_ds = xr.Dataset(\n",
    "    coords=dict(\n",
    "        time=coamps_ds['time'],\n",
    "        x=([\"x\"], x_arr),\n",
    "        y=([\"y\"], y_arr),\n",
    "    ),\n",
    "    attrs=dict(description=\"Storm-relative COAMPS-TC fields.\"),\n",
    ")\n",
    "\n",
    "# Interpolate the wind field magnitude and components and store in the Dataset.\n",
    "#TODO: can this be sped up?\n",
    "for field in ['ws', 'u', 'v']:\n",
    "    coamps_xy_ds[field] = xr.apply_ufunc(\n",
    "        regrid_irregular_field,\n",
    "        coamps_ds[field],\n",
    "        coamps_ds['x'],\n",
    "        coamps_ds['y'],\n",
    "        x_arr,\n",
    "        y_arr,\n",
    "        input_core_dims=[[\"latitude\", \"longitude\"],\n",
    "                         [\"latitude\", \"longitude\"],\n",
    "                         [\"latitude\", \"longitude\"],\n",
    "                         [\"x\"],\n",
    "                         [\"y\"]],\n",
    "        output_core_dims=[[\"x\", \"y\"]],\n",
    "        exclude_dims=set((\"x\", \"y\")),  # dimensions allowed to change size.\n",
    "        vectorize=True,\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IBTrACS calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cast the 34 kt, 50 kt, and 64 kt wind speed radii, max winds, and min pressure to floats.  These will be used for future plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ibtracs_column_to_float(ibtracs_da):\n",
    "    ibtracs_da_as_float = (ibtracs_da\n",
    "                           .replace(' ', np.nan)\n",
    "                           .apply(float))\n",
    "    return ibtracs_da_as_float"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "for speed in ['34', '50', '64']:\n",
    "    for dir in ['NE', 'NW', 'SW', 'SE']:\n",
    "        col_name = f'USA_R{speed}_{dir}'\n",
    "        ibtracs_df[col_name] = ibtracs_column_to_float(ibtracs_df[col_name])\n",
    "\n",
    "ibtracs_df['USA_WIND'] = ibtracs_column_to_float(ibtracs_df['USA_WIND'])\n",
    "ibtracs_df['USA_PRES'] = ibtracs_column_to_float(ibtracs_df['USA_PRES'])\n",
    "ibtracs_df['USA_RMW'] = ibtracs_column_to_float(ibtracs_df['USA_RMW'])\n",
    "ibtracs_df['STORM_SPEED'] = ibtracs_column_to_float(ibtracs_df['STORM_SPEED'])\n",
    "ibtracs_df['STORM_DIR'] = ibtracs_column_to_float(ibtracs_df['STORM_DIR'])\n",
    "ibtracs_df['LON'] = ibtracs_column_to_float(ibtracs_df['LON'])\n",
    "ibtracs_df['LAT'] = ibtracs_column_to_float(ibtracs_df['LAT'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Drifter calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate drifter (or buoy) variables.  Most variables derive from the scalar spectra or GPS data reported by the buoy.  Drifter datasets are stored as DataFrames and rely on methods in the `drifter_df.buoy` defined in `src/` (for reuse and to improve notebook readability).  Accessors extend Pandas DataFrame objects with custom methods, properties, etc. (https://pandas.pydata.org/docs/development/extending.html).  Less routine calculations use functions defined directly in this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Depth colocation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine the depth (positive down) at each drifter position by colocating with the GEBCO elevation grid.  Any negative values (on land) should be screened."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def colocate_bathrymetry(drifter_df, bathy_ds):\n",
    "    latitude_sort_indices = np.searchsorted(bathy_ds['lat'].values, drifter_df['latitude'])\n",
    "    longitude_sort_indices = np.searchsorted(bathy_ds['lon'].values, drifter_df['longitude'])\n",
    "    depth = -1.0 * bathy_ds['elevation'].values[latitude_sort_indices, longitude_sort_indices]\n",
    "    depth[depth < 0] = np.NaN\n",
    "    return depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df['depth'] = colocate_bathrymetry(drifter_df, bathymetry_ds)\n",
    "drifter_df = drifter_df.query('depth > 0')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Drift speed and direction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate drift speed and direction from GPS positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df.loc[:, ['drift_speed', 'drift_direction']] = drifter_df.groupby(\n",
    "    level='id',\n",
    "    as_index=True,\n",
    "    group_keys=False\n",
    ").apply(\n",
    "    lambda df: df.buoy.drift_speed_and_direction(),\n",
    ")\n",
    "\n",
    "ian_spotter_coamps_df.loc[:, ['drift_speed', 'drift_direction']] = ian_spotter_coamps_df.groupby(\n",
    "    level='id',\n",
    "    as_index=True,\n",
    "    group_keys=False\n",
    ").apply(\n",
    "    lambda df: df.buoy.drift_speed_and_direction(),\n",
    ")\n",
    "\n",
    "fiona_spotter_coamps_df.loc[:, ['drift_speed', 'drift_direction']] = fiona_spotter_coamps_df.groupby(\n",
    "    level='id',\n",
    "    as_index=True,\n",
    "    group_keys=False\n",
    ").apply(\n",
    "    lambda df: df.buoy.drift_speed_and_direction(),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trim spectral variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trim the frequency arrays of all buoys to the upper cutoff of 0.5 Hz.\n",
    "(The last frequency bin of the Spotter buoy data contains all the energy from 0.49805 Hz to 0.6543 Hz and needs to be trimmed to produce an accurate mean square slope estimate.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trim_spectral_vars(drifter_df):\n",
    "    valid_frequencies = drifter_df['frequency'] <= 0.5\n",
    "    spectral_vars = ['energy_density', 'a1', 'a2', 'b1', 'b2', 'frequency']\n",
    "    for var in spectral_vars:\n",
    "        drifter_df[var] = drifter_df[var][valid_frequencies]\n",
    "    return drifter_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df = drifter_df.apply(trim_spectral_vars, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unify microSWIFT frequency arrays (the arrays have slight differences due to floating point errors)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def unify_microswift_frequencies(drifter_df):\n",
    "    frequencies = drifter_df['frequency']\n",
    "    drifter_df.loc[:, 'frequency'] = (drifter_df['frequency']\n",
    "                                      .apply(lambda f: frequencies.iloc[0]))\n",
    "    return drifter_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "is_microswift = drifter_df['drifter_type'] == 'microswift'\n",
    "drifter_df.loc[is_microswift] = unify_microswift_frequencies(drifter_df.loc[is_microswift])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wave direction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate frequency-dependent wave direction from the directional moments as:\n",
    "\\begin{equation*}\n",
    "    \\mod\\left(270\\degree - \\frac{180\\degree}{\\pi} \\arctan_2(b1, a1), 360\\degree \\right)\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df['direction'] = drifter_df.buoy.wave_direction()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean square slope"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the spectral mean square slope. In the frequency domain, this is proportional to the fourth moment of the spectrum:\n",
    "\\begin{equation*}\n",
    "\\mathrm{mss} = \\int \\frac{E(f) (2 \\pi f)^4}{g^2} df\n",
    "\\end{equation*}\n",
    "Where $E(f)$ is the energy density spectrum as a function of frequency, $f$, and $g$ is the acceleration of gravity.  Using deep water dispersion, $\\omega^2 = gk$ (with $\\omega=2\\pi f$), and with $E(f) \\propto a^2$, this is proportional to the product of amplitude and wavenumber squared, $(ak)^2$, or slope squared.\n",
    "\n",
    "We can also calculate mean square slope using wavenumber calculated from the dispersion relationship (the Florida shelf is in relatively shallow waters, thus the deep water assumption may not hold for the lower frequencies in the spectrum).  Frequency-based mean square slope (which assumes deep water) and wavenumber-based mean square slope (corrected for intermediate to shallow water) will be compared in a subsequent notebook.\n",
    "\n",
    "Mean square slope is calculated from the wavenumber spectrum as:\n",
    "\\begin{equation*}\n",
    "\\mathrm{mss} = \\int E(k) k^2 dk\n",
    "\\end{equation*}\n",
    "Where $k$ are the wavenumbers calculated from the (intrinsic) inverse dispersion relationship:\n",
    "\\begin{equation*}\n",
    "\\omega^2 = gk \\tanh(kh)\n",
    "\\end{equation*}\n",
    "which is valid for all water depths, and the wavenumber spectrum $E(k)$ can be calculated from the frequency spectrum as:\n",
    "\\begin{equation*}\n",
    "E(k) = E(\\omega)\\frac{d\\omega}{dk}\n",
    "\\end{equation*}\n",
    "with\n",
    "\\begin{equation*}\n",
    "E(\\omega) = E(f)\\frac{df}{d\\omega}\n",
    "\\end{equation*}\n",
    "where $d\\omega / dk = c_g$ (where $c_g$ is the group velocity) is the Jacobian which maps changes in angular frequency to changes in wavenumber and $df / d\\omega = 1/(2\\pi)$ is the Jacobian which maps changes in frequency to changes in angular frequency."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transform the frequency spectra to wavenumber spectra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df[['energy_density_wn', 'wavenumber']] = drifter_df.buoy.wavenumber_energy_density(var_rtol=0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the wavenumber-based mean square slope."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df['mean_square_slope_wn'] = drifter_df.buoy.wavenumber_mean_square_slope()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the frequency-based mean square slope (which assumes deep water) for comparison.  Since Spotters have a small amount of windage, this is the mean square slope in the \"observed\" reference frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df['mean_square_slope_observed'] = drifter_df.buoy.mean_square_slope()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Moment-weighted means"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate mean square slope weighted wave direction and spread from the mss-weighted directional moments:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{gather*}\n",
    "    \\overline{a}_1 = \\frac{\\int a_1 (f) f^4 E(f)  df}{\\int f^4 E(f) df}\n",
    "    \\qquad\n",
    "    \\overline{b}_1 = \\frac{\\int b_1 (f) f^4 E(f)  df}{\\int f^4 E(f) df}\n",
    "\\end{gather*}\n",
    "\n",
    "and $a_1(f)$ and $b_1(f)$ are the first two directional moments which represent wave motion along the N-S and E-W principle axes at each frequency, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Any weighted direction needs to be done from moments, since the direction wraps\n",
    "drifter_df['mss_weighted_a1'] = drifter_df.buoy.moment_weighted_mean(column='a1', n=4)\n",
    "drifter_df['mss_weighted_b1'] = drifter_df.buoy.moment_weighted_mean(column='b1', n=4)\n",
    "drifter_df['mss_weighted_direction'] = drifter_df.buoy.wave_direction(a1_col='mss_weighted_a1', b1_col='mss_weighted_b1')\n",
    "drifter_df['mss_weighted_directional_spread'] = drifter_df.buoy.wave_directional_spread(\n",
    "    a1_col='mss_weighted_a1',\n",
    "    b1_col='mss_weighted_b1'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distance from eye"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the distance and heading from the eye of the storm to each buoy.  This requires interpolation of the best track data onto the buoy times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_time_interpolation(\n",
    "    interpolator: scipy.interpolate.interp1d,\n",
    "    query_times: pd.DatetimeIndex,\n",
    ") -> Tuple[np.ndarray, ...]:\n",
    "    \"\"\" Interpolate a time series of data at query times (datetimes).\"\"\"\n",
    "    numeric_query_times = pd.to_numeric(query_times)\n",
    "    return interpolator(numeric_query_times).T\n",
    "\n",
    "\n",
    "def storm_distances(\n",
    "    distance: np.ndarray,\n",
    "    bearing_deg: np.ndarray,\n",
    ") -> Tuple[np.ndarray, np.ndarray]:\n",
    "    \"\"\" Return storm east and north components from distance and bearing.\"\"\"\n",
    "    bearing_rad = np.deg2rad(bearing_deg)\n",
    "    distance_east = distance * np.sin(bearing_rad)\n",
    "    distance_north = distance * np.cos(bearing_rad)\n",
    "    return distance_east, distance_north\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct a 1-dimensonal linear interpolation object to interpolate the best track latitude and longitude (`coords`) at new query times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = ibtracs_df[['LON', 'LAT']]\n",
    "direction = ibtracs_df['STORM_DIR']\n",
    "direction_unwrap = np.unwrap(direction, period=360)\n",
    "storm_speed = ibtracs_df['STORM_SPEED'] # kt\n",
    "radius_max_winds = ibtracs_df['USA_RMW'] # nm\n",
    "max_wind_speed = ibtracs_df['USA_WIND']# kt\n",
    "\n",
    "times = pd.to_numeric(ibtracs_df.index.tz_localize(None).to_numpy())\n",
    "coords_interpolator = scipy.interpolate.interp1d(times, coords, fill_value=np.nan, axis=0)\n",
    "direction_interpolator = scipy.interpolate.interp1d(times, direction_unwrap, fill_value=np.nan, axis=0)\n",
    "storm_speed_interpolator = scipy.interpolate.interp1d(times, storm_speed, fill_value=np.nan, axis=0)\n",
    "radius_max_winds_interpolator = scipy.interpolate.interp1d(times, radius_max_winds, fill_value=np.nan, axis=0)\n",
    "max_wind_speed_interpolator = scipy.interpolate.interp1d(times, max_wind_speed, fill_value=np.nan, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add the `ibtracs_df` columns to the buoy DataFrame.  These will be used to compute the distance and heading."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "storm_positions = get_time_interpolation(coords_interpolator, drifter_df.index.get_level_values(level='time'))\n",
    "drifter_df['storm_longitude'], drifter_df['storm_latitude'] = storm_positions\n",
    "\n",
    "storm_headings = get_time_interpolation(direction_interpolator, drifter_df.index.get_level_values(level='time'))\n",
    "drifter_df['storm_heading'] = storm_headings % 360\n",
    "\n",
    "storm_speeds = get_time_interpolation(storm_speed_interpolator, drifter_df.index.get_level_values(level='time'))\n",
    "drifter_df['storm_speed_kt'] = storm_speeds\n",
    "\n",
    "storm_radius_max_winds = get_time_interpolation(radius_max_winds_interpolator, drifter_df.index.get_level_values(level='time'))\n",
    "drifter_df['storm_radius_max_wind_nmi'] = storm_radius_max_winds\n",
    "\n",
    "storm_max_wind_speeds = get_time_interpolation(max_wind_speed_interpolator, drifter_df.index.get_level_values(level='time'))\n",
    "drifter_df['storm_max_wind_speed_kt'] = storm_max_wind_speeds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the distance and bearing from the storm to the buoy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df['storm_distance'], drifter_df['storm_back_bearing'] \\\n",
    "    = geodesy.great_circle_pairwise(longitude_a=drifter_df['storm_longitude'],\n",
    "                                    latitude_a=drifter_df['storm_latitude'],\n",
    "                                    longitude_b=drifter_df['longitude'],\n",
    "                                    latitude_b=drifter_df['latitude'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the east and north distance from the storm in kilometers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df['storm_distance_east'], drifter_df['storm_distance_north'] \\\n",
    "            = storm_distances(distance=drifter_df['storm_distance'],\n",
    "                              bearing_deg=drifter_df['storm_back_bearing'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trapping criterion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kudryavtsev et al. (2021) derive a criterion for wave trapping (``resonance'' between the group velocity and storm translation speed) based on scaling for self-similar wave growth.  Resonance occurs when the nondimensonal radius of max wind\n",
    "\n",
    "\\begin{equation*}\n",
    "    \\tilde{R}_m = \\frac{g R_m}{u^2_m}\n",
    "\\end{equation*}\n",
    "\n",
    "and nondimensional critical fetch length\n",
    "\n",
    "\\begin{equation*}\n",
    "    \\tilde{L}_{cr}^{m} = -c_\\alpha^{-1/q} \\frac{q}{1+q} \\left(\\frac{u_m}{2 V}\\right)^{1/q}\n",
    "\\end{equation*}\n",
    "\n",
    "match ($\\tilde{R}_m \\approx \\tilde{L}_{cr}^{m}$).  When $\\tilde{R}_m > \\tilde{L}_{cr}^{m}$ the tropical cyclone is slower than the optimal speed for resonance, and when $\\tilde{R}_m < \\tilde{L}_{cr}^{m}$ it is faster.\n",
    "Here $R_m$ is the storm's dimensional radius of max wind, $u_m$ is the storm's max wind speed, and $g$ is the acceleration of gravity.  In the nondimensional critical fetch length, $u_m$ is the maximum wind speed, the $V$ is the storm's translation speed and $q$ and $c_\\alpha$ are constants related to wave growth.  This equation is evaluated using $u_r = u_m$ with $q=-1/4$ and $c_\\alpha = 11.8$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def kudryavtsev_trapping_criterion(\n",
    "    radius_max_winds_nmi,\n",
    "    max_wind_speed_kt,\n",
    "    storm_speed_kt,\n",
    "    c_alpha = 11.8,\n",
    "    q = -0.25,\n",
    "):\n",
    "    \"\"\"\n",
    "    Return the ratio of dimensionless radius of max wind to critical fetch.\n",
    "    \"\"\"\n",
    "    g = 9.81\n",
    "\n",
    "    # Unit conversions to SI\n",
    "    radius_max_winds_m = nmi_to_m(radius_max_winds_nmi)\n",
    "    max_wind_speed_mps = knots_to_mps(max_wind_speed_kt)\n",
    "    storm_speed_mps = knots_to_mps(storm_speed_kt)\n",
    "\n",
    "    dimless_radius_max_wind = g * radius_max_winds_m / max_wind_speed_mps**2\n",
    "    c_cr = -c_alpha**(-1/q) * q/(1 + q)\n",
    "    dimless_critical_fetch = c_cr * (max_wind_speed_mps / (2 * storm_speed_mps))**(1/q)\n",
    "    return dimless_radius_max_wind / dimless_critical_fetch\n",
    "\n",
    "def knots_to_mps(knots):\n",
    "    \"\"\" Convert knots to meters per second. \"\"\"\n",
    "    return knots / 1.944\n",
    "\n",
    "def nmi_to_m(nmi):\n",
    "    \"\"\" Convert nautical miles to meters. \"\"\"\n",
    "    return nmi * 1852"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate kudryavtsev trapping criterion\n",
    "drifter_df['storm_trapping_criterion'] = kudryavtsev_trapping_criterion(\n",
    "    radius_max_winds_nmi=drifter_df['storm_radius_max_wind_nmi'],\n",
    "    max_wind_speed_kt=drifter_df['storm_max_wind_speed_kt'],\n",
    "    storm_speed_kt=drifter_df['storm_speed_kt'],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Storm-centric coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Locate the buoys in a coordinate system that is rotated into the storm (storm heading aligned with y-axis)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rotate_xy(\n",
    "    x: np.ndarray,\n",
    "    y: np.ndarray,\n",
    "    theta: np.ndarray\n",
    ") -> Tuple[np.ndarray, np.ndarray]:\n",
    "    \"\"\"\n",
    "    Perform 2D coordinate rotation of points (x, y).\n",
    "\n",
    "    Args:\n",
    "        x (np.ndarray): x coordinate to rotate\n",
    "        y (np.ndarray): y coordinate to rotate\n",
    "        theta (np.ndarray): rotation angles in [rad]\n",
    "\n",
    "    Returns:\n",
    "        Tuple[np.ndarray, np.ndarray]: rotated x and y coordinates\n",
    "    \"\"\"\n",
    "    x_rot = x*np.cos(theta) - y*np.sin(theta)\n",
    "    y_rot = x*np.sin(theta) + y*np.cos(theta)\n",
    "    return x_rot, y_rot\n",
    "\n",
    "def normalize_storm_distance(distance_km, rmw_nmi):\n",
    "    rmw_km = rmw_nmi * 1.852\n",
    "    return distance_km / rmw_km\n",
    "\n",
    "def storm_distance_magnitude(distance_x, distance_y):\n",
    "    return np.sqrt(distance_x**2 + distance_y**2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df['storm_distance_x'], drifter_df['storm_distance_y'] = \\\n",
    "        rotate_xy(x=drifter_df['storm_distance_east'],\n",
    "                  y=drifter_df['storm_distance_north'],\n",
    "                  theta=np.deg2rad(drifter_df['storm_heading']))\n",
    "\n",
    "drifter_df['storm_distance_x_normalized'] = normalize_storm_distance(\n",
    "    drifter_df['storm_distance_x'],\n",
    "    drifter_df['storm_radius_max_wind_nmi'],\n",
    ")\n",
    "\n",
    "drifter_df['storm_distance_y_normalized'] = normalize_storm_distance(\n",
    "    drifter_df['storm_distance_y'],\n",
    "    drifter_df['storm_radius_max_wind_nmi'],\n",
    ")\n",
    "\n",
    "drifter_df['storm_distance_xy_normalized_magnitude'] = storm_distance_magnitude(\n",
    "    drifter_df['storm_distance_x_normalized'],\n",
    "    drifter_df['storm_distance_y_normalized'],\n",
    ")\n",
    "\n",
    "drifter_df['storm_frame_azimuth'] = (drifter_df['storm_back_bearing'] - drifter_df['storm_heading']) % 360\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Storm position"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Classify the quadrant of each buoy based on its bearing relative to the center and heading of the storm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "quadrant_definitions = {\n",
    "    'front_right': (0, 90),\n",
    "    'back_right': (90, 180),\n",
    "    'back_left': (180, 270),\n",
    "    'front_left': (270, 360),\n",
    "}\n",
    "\n",
    "def classify_storm_position(\n",
    "    relative_bearing: np.ndarray,\n",
    "    position_definitions: dict[str, tuple] = quadrant_definitions\n",
    ") -> np.ndarray:\n",
    "    \"\"\" Classify the relative bearing of a point with respect to the storm. \"\"\"\n",
    "    positions = np.empty(relative_bearing.shape, dtype=object)\n",
    "    for label, limits in position_definitions.items():\n",
    "        in_position = np.logical_and(relative_bearing >= limits[0],\n",
    "                                     relative_bearing < limits[1])\n",
    "        positions[in_position] = label\n",
    "    return positions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df['storm_quadrant'] = classify_storm_position(drifter_df['storm_back_bearing'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### COAMPS-TC wind colocation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine the surface wind speed at each drifter position by colocating with COAMPS-TC model output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "TEMPORAL_TOLERANCE = np.timedelta64(60, 'm')\n",
    "\n",
    "def colocate_drifter_and_coamps(\n",
    "    drifter_df: pd.DataFrame,\n",
    "    coamps_xy_ds: xr.Dataset\n",
    ") -> dict:\n",
    "    \"\"\" Colocate COAMPS-TC wind speed at each position in drifter_df.\"\"\"\n",
    "    wind_speed = {\n",
    "        'ws': [],\n",
    "        'u': [],\n",
    "        'v': [],\n",
    "    }\n",
    "    for field in ['ws', 'u', 'v']:\n",
    "        buoy = {\n",
    "            'time': drifter_df.index.to_numpy(dtype='datetime64[m]'),\n",
    "            'latitude': drifter_df['storm_distance_north'].values,\n",
    "            'longitude': drifter_df['storm_distance_east'].values,\n",
    "        }\n",
    "        wind_speed[field] = coamps.match_model_and_buoy_by_interpolation(\n",
    "            buoy,\n",
    "            model={\n",
    "                'time': np.array(coamps_xy_ds['time'].values, dtype='datetime64[m]'),\n",
    "                'latitude': coamps_xy_ds['y'].values * 10**(-3),\n",
    "                'longitude': coamps_xy_ds['x'].values * 10**(-3),\n",
    "                'field': coamps_xy_ds[field].values,\n",
    "            },\n",
    "            temporal_tolerance=TEMPORAL_TOLERANCE,\n",
    "            bounds_error=False,\n",
    "        )\n",
    "    return wind_speed\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_ids = drifter_df.index.get_level_values(level='id').unique()\n",
    "for drifter_id in drifter_ids:\n",
    "    wind_speed = colocate_drifter_and_coamps(drifter_df.loc[drifter_id],\n",
    "                                             coamps_xy_ds)\n",
    "    drifter_df.loc[drifter_id, 'wind_speed'] = wind_speed['ws']\n",
    "    drifter_df.loc[drifter_id, 'wind_speed_u'] = wind_speed['u']\n",
    "    drifter_df.loc[drifter_id, 'wind_speed_v'] = wind_speed['v']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wind-wave alignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Directional alignment between the wind and the waves (wind-wave alignment) is calculated using COAMPS-TC 10-m wind speeds and bulk wave direction measured by the buoys. This can also be thought of as the wave heading relative to wind direction.\n",
    "Here we define wind-wave alignment as the smallest angle resulting from a difference between the wind and wave directions\n",
    "\n",
    "\\begin{equation*}\n",
    "\\Delta\\theta\n",
    "=\n",
    "\\begin{cases}\n",
    "    \\theta_{\\textrm{wind}} - \\theta_{\\textrm{wave}}, & \\text{if } -180\\degree \\leq \\theta_{\\textrm{wind}} - \\theta_{\\textrm{wave}} \\leq 180\\degree\\\\\n",
    "    \\theta_{\\textrm{wind}} - \\theta_{\\textrm{wave}} -  360\\degree, & \\text{if } \\theta_{\\textrm{wind}} - \\theta_{\\textrm{wave}} > 180\\degree\\\\\n",
    "    \\theta_{\\textrm{wind}} - \\theta_{\\textrm{wave}} +  360\\degree, & \\text{if } \\theta_{\\textrm{wind}} - \\theta_{\\textrm{wave}} < -180\\degree\\\\\n",
    "\\end{cases}\n",
    "\\end{equation*}\n",
    "\n",
    "or more concisely as\n",
    "\n",
    "\\begin{equation*}\n",
    "\\Delta\\theta\n",
    "= (\\theta_{\\textrm{wind}} - \\theta_{\\textrm{wave}} + 180\\degree \\mod \\;\\;\\: 360\\degree)  - 180\\degree\n",
    "\\end{equation*}\n",
    "\n",
    "where $\\theta_{\\textrm{wind}}$ and $\\theta_{\\textrm{wave}}$ are in the ``going to'' convention and $\\mathrm{mod}$ is the modulo operation.\n",
    "This definition ensures, for example, that wind heading towards 0$\\degree$ (north) and waves heading to 270$\\degree$ (west) results in an alignment of 90$\\degree$ and not $-$270$\\degree$.\n",
    "The signed wind-wave alignment is $\\in[-180\\degree, 180\\degree]$.\n",
    "Wind-wave alignment is calculated using an energy-weighted wave direction from the buoy's energy-weighted directional moments \\cite{kuik_method_1988}\n",
    "\n",
    "\\begin{equation}\n",
    "    \\theta_\\mathrm{wave} = \\left(270\\degree - \\frac{180\\degree}{\\pi} \\mathrm{arctan2}(\\overline{b}_1, \\overline{a}_1) \\right) \\mod \\;\\;\\: 360\\degree\n",
    "\\end{equation}\n",
    "\n",
    "where\n",
    "\n",
    "\\begin{gather*}\n",
    "    \\overline{a}_1 = \\frac{\\int a_1 (f) E(f) df}{\\int E(f) df}\n",
    "    \\qquad\n",
    "    \\overline{b}_1 = \\frac{\\int b_1 (f) E(f) df}{\\int E(f) df}\n",
    "\\end{gather*}\n",
    "\n",
    "and $a_1(f)$ and $b_1(f)$ are the first two directional moments which represent wave motion along the N-S and E-W principle axes at each frequency, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wind_wave_alignment(wind_dir_from: np.ndarray, wave_dir_from: np.ndarray):\n",
    "    \"\"\" Compute alignment as the angle between wave and wind vectors.\n",
    "\n",
    "    Compute the alignment of the waves relative to the wind. The\n",
    "    difference will yield a combination of angles that form either the\n",
    "    major or minor arc between the wind and wave unit vectors, however this\n",
    "    returns the signed minor arc angle on [-180 180] only.\n",
    "\n",
    "    Args:\n",
    "        wind_dir_from (np.ndarray): wind direction in \"coming from\"\n",
    "            (meteorological) convention\n",
    "        wave_dir_from (np.ndarray): wave direction in \"coming from\"\n",
    "            (meteorological) convention\n",
    "\n",
    "    Returns:\n",
    "        np.ndarray: alignment on [-180 180] deg\n",
    "    \"\"\"\n",
    "    wave_dir_to = kinematics.coming_to_going(wave_dir_from)\n",
    "    wind_dir_to = kinematics.coming_to_going(wind_dir_from)\n",
    "    alignment_full = wind_dir_to - wave_dir_to\n",
    "    return (alignment_full + 180) % 360 - 180\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate wind direction from the components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df['wind_direction'] = coamps.wind_direction(drifter_df['wind_speed_u'], drifter_df['wind_speed_v'])\n",
    "\n",
    "ian_spotter_coamps_df['wind_direction'] = coamps.wind_direction(ian_spotter_coamps_df['wind_speed_u'], ian_spotter_coamps_df['wind_speed_v'])\n",
    "\n",
    "fiona_spotter_coamps_df['wind_direction'] = coamps.wind_direction(fiona_spotter_coamps_df['wind_speed_u'], fiona_spotter_coamps_df['wind_speed_v'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the wind-wave alignment as the angle between the wind and wave vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df['mean_wind_wave_alignment'] = wind_wave_alignment(drifter_df['wind_direction'],\n",
    "                                                             drifter_df['mean_direction'])\n",
    "drifter_df['mean_wind_wave_alignment_abs'] = drifter_df['mean_wind_wave_alignment'].abs()\n",
    "\n",
    "drifter_df['peak_wind_wave_alignment'] = wind_wave_alignment(drifter_df['wind_direction'],\n",
    "                                                             drifter_df['peak_direction'])\n",
    "drifter_df['peak_wind_wave_alignment_abs'] = drifter_df['peak_wind_wave_alignment'].abs()\n",
    "\n",
    "drifter_df['mss_weighted_wind_wave_alignment'] = wind_wave_alignment(drifter_df['wind_direction'],\n",
    "                                                                     drifter_df['mss_weighted_direction'])\n",
    "drifter_df['mss_weighted_wind_wave_alignment_abs'] = drifter_df['mss_weighted_wind_wave_alignment'].abs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "ian_spotter_coamps_df['mean_wind_wave_alignment'] = wind_wave_alignment(ian_spotter_coamps_df['wind_direction'],\n",
    "                                                                        ian_spotter_coamps_df['mean_direction'])\n",
    "ian_spotter_coamps_df['mean_wind_wave_alignment_abs'] = ian_spotter_coamps_df['mean_wind_wave_alignment'].abs()\n",
    "\n",
    "fiona_spotter_coamps_df['mean_wind_wave_alignment'] = wind_wave_alignment(fiona_spotter_coamps_df['wind_direction'],\n",
    "                                                                        fiona_spotter_coamps_df['mean_direction'])\n",
    "fiona_spotter_coamps_df['mean_wind_wave_alignment_abs'] = fiona_spotter_coamps_df['mean_wind_wave_alignment'].abs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The structure of a storm will contain regions of cross, aligned, and opposing swell (see Figure 3 Holthuijsen et al. 2012). These are regions in which the swell and locally-generated wind sea (and thus the wind direction) are crossing over one another (cross), nearly in-line (aligned), or are in entirely opposite directions (opposing). Holthuijsen et al. define these quantitatively as:\n",
    "- aligned swell: swell traveling within 45 degrees from the wind direction\n",
    "- cross swell: swell traveling within 45 degrees from the normal to the wind direction\n",
    "- opposing swell: swell traveling within 45 degrees from the opposing direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "def classify_wind_wave_alignment(alignment):\n",
    "    \"\"\" Classify wind-wave alignment based on Holthuijsen et al. 2012. \"\"\"\n",
    "    classifications = np.empty(alignment.shape, dtype=object)\n",
    "\n",
    "    alignment_abs = np.abs(alignment)\n",
    "    aligned = alignment_abs <= 45\n",
    "    crossing = np.logical_and(alignment_abs > 45, alignment_abs <= 135)\n",
    "    opposing = np.logical_and(alignment_abs > 135, alignment_abs <= 180)\n",
    "\n",
    "    classifications[aligned] = 'aligned'\n",
    "    classifications[crossing] = 'crossing'\n",
    "    classifications[opposing] = 'opposing'\n",
    "\n",
    "    return classifications\n",
    "\n",
    "cat_map = {\n",
    "    'aligned': 0,\n",
    "    'crossing': 1,\n",
    "    'opposing': 2,\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df['mean_wind_wave_alignment_cat'] = classify_wind_wave_alignment(drifter_df['mean_wind_wave_alignment'])\n",
    "drifter_df['mean_wind_wave_alignment_int'] = drifter_df['mean_wind_wave_alignment_cat'].map(cat_map)\n",
    "drifter_df['peak_wind_wave_alignment_cat'] = classify_wind_wave_alignment(drifter_df['peak_wind_wave_alignment'])\n",
    "drifter_df['peak_wind_wave_alignment_int'] = drifter_df['peak_wind_wave_alignment_cat'].map(cat_map)\n",
    "drifter_df['mss_weighted_wind_wave_alignment_cat'] = classify_wind_wave_alignment(drifter_df['mss_weighted_wind_wave_alignment'])\n",
    "drifter_df['mss_weighted_wind_wave_alignment_int'] = drifter_df['mss_weighted_wind_wave_alignment_cat'].map(cat_map)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wave age\n",
    "\n",
    "Wave age is defined as the ratio of the wave phase speed, typically defined using the peak period, to the wind speed. It is a proxy for the magnitude of wind forcing and wave growth. \n",
    "\\begin{equation*}\n",
    "\\textrm{wave age} \\equiv \\frac{c_p}{U_{10}}\n",
    "\\end{equation*}\n",
    "where the phase speed, $c_p=\\omega/k$, is the wave's angular frequency divided by its wavenumber. Using the bulk peak period to compute the phase speed in the deep water limit, it is \n",
    "\\begin{equation*}\n",
    "c_p = \\frac{g}{2\\pi} T_p\n",
    "\\end{equation*}\n",
    "For low wave age values, the waves are strongly forced and growth rates, normalized by the wave frequency used in the phase speed calculation,\n",
    "are high. For values of wave age around unity the growth rates are very small and may change sign as unity is exceeded (Melville, 2004)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def phase_speed(period):\n",
    "    \"\"\" Calculate phase speed from wave period. \"\"\"\n",
    "    return 9.81 * period /(2*np.pi)\n",
    "\n",
    "def wave_age(phase_speed, wind_speed):\n",
    "    \"\"\" Calculate wave age from phase speed and wind speed. \"\"\"\n",
    "    return phase_speed / wind_speed\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifter_df['mean_phase_speed'] = phase_speed(drifter_df['mean_period'])\n",
    "drifter_df['peak_phase_speed'] = phase_speed(drifter_df['peak_period'])\n",
    "drifter_df['peak_wave_age'] = wave_age(drifter_df['peak_phase_speed'], drifter_df['wind_speed'])\n",
    "drifter_df['mean_wave_age'] = wave_age(drifter_df['mean_phase_speed'], drifter_df['wind_speed'])\n",
    "drifter_df['peak_inverse_wave_age'] = drifter_df['peak_wave_age']**(-1)\n",
    "drifter_df['mean_inverse_wave_age'] = drifter_df['mean_wave_age']**(-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DataFrame by drifter type"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Split the drifter DataFrames by type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_drifter_type_df(drifter_df, drifter_type):\n",
    "    \"\"\" Return a DataFrame containing only the specified drifter type. \"\"\"\n",
    "    is_drifter_type = drifter_df['drifter_type'] == drifter_type\n",
    "    return drifter_df.loc[is_drifter_type, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "microswift_df = get_drifter_type_df(drifter_df, 'microswift')\n",
    "spotter_df = get_drifter_type_df(drifter_df, 'spotter')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Store"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Store the updated variables for use in subsequent notebooks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "%store microswift_df\n",
    "%store spotter_df\n",
    "%store ian_spotter_coamps_df\n",
    "%store fiona_spotter_coamps_df\n",
    "%store coamps_ds\n",
    "%store coamps_xy_ds\n",
    "%store ibtracs_df\n",
    "%store swan_currents_ds\n",
    "%store swan_no_currents_ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "Holthuijsen, L. H., Powell, M. D., & Pietrzak, J. D. (2012). Wind and waves in extreme hurricanes. J. Geophys. Res., 117(C9), C09003. https://doi.org/10.1029/2012JC007983\n",
    "\n",
    "Kudryavtsev, V., Yurovskaya, M., & Chapron, B. (2021). Self‐Similarity of Surface Wave Developments Under Tropical Cyclones. Journal of Geophysical Research: Oceans, 126(4), e2020JC016916. https://doi.org/10.1029/2020JC016916\n",
    "\n",
    "\n",
    "Melville, W. K. (2004). Wave slope and wave age effects in measurements of electromagnetic bias. Journal of Geophysical Research, 109(C7), C07018. https://doi.org/10.1029/2002JC001708\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "NHCI_geo_python3-9",
   "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.3"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}