{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Retieving wind drift factor from trajectory\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from datetime import datetime, timedelta\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport cmocean\nfrom opendrift.models.oceandrift import OceanDrift\nfrom opendrift.models.physics_methods import wind_drift_factor_from_trajectory, distance_between_trajectories, skillscore_liu_weissberg"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "A very simple drift model is: current + wind_drift_factor*wind\nwhere wind_drift_factor for surface drift is typically\n0.033 if Stokes drift included, and 0.02 in addition to Stokes drift.\nThis example illustrates how a best-fit wind_drift_factor\ncan be calculated from an observed trajectory, using two different methods.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First we simulate a synthetic drifter trajectory\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ot = OceanDrift(loglevel=50)\not.add_readers_from_list([ot.test_data_folder() +\n    '16Nov2015_NorKyst_z_surface/norkyst800_subset_16Nov2015.nc',\n    ot.test_data_folder() + '16Nov2015_NorKyst_z_surface/arome_subset_16Nov2015.nc'], lazy=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Adding some horizontal diffusivity as \"noise\"\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ot.set_config('drift:horizontal_diffusivity', 10)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Using a wind_drift_factor of 0.33 i.e. drift is current + 3.3% of wind speed\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ot.seed_elements(lon=4, lat=60, number=1, time=ot.env.readers[list(ot.env.readers)[0]].start_time,\n        wind_drift_factor=0.033)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ot.run(duration=timedelta(hours=12), time_step=600)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Secondly, calculating the wind_drift_factor which reproduces the \"observed\" trajectory with minimal difference\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "drifter_lons = ot.history['lon'][0]\ndrifter_lats = ot.history['lat'][0]\ndrifter_times = ot.get_time_array()[0]\ndrifter={'lon': drifter_lons, 'lat': drifter_lats,\n        'time': drifter_times, 'linewidth': 2, 'color': 'b', 'label': 'Synthetic drifter'}\n\no = OceanDrift(loglevel=50)\no.add_readers_from_list([o.test_data_folder() +\n    '16Nov2015_NorKyst_z_surface/norkyst800_subset_16Nov2015.nc',\n    o.test_data_folder() + '16Nov2015_NorKyst_z_surface/arome_subset_16Nov2015.nc'], lazy=False)\nt = o.env.get_variables_along_trajectory(variables=['x_sea_water_velocity', 'y_sea_water_velocity', 'x_wind', 'y_wind'],\n        lons=drifter_lons, lats=drifter_lats, times=drifter_times)\n\nwind_drift_factor, azimuth = wind_drift_factor_from_trajectory(t)\n\no.seed_elements(lon=4, lat=60, number=1, time=ot.env.readers[list(ot.env.readers)[0]].start_time,\n                wind_drift_factor=0.033)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "New simulation, this time without diffusivity/noise\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "o.run(duration=timedelta(hours=12), time_step=600)\n\no.plot(fast=True, legend=True, drifter=drifter)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The mean retrieved wind_drift_factor is 0.036, slichtly higher than the value 0.033 used for the simulation.\nThe difference is due to the \"noise\" added by horizontal diffusion.\nNote that the retieved wind_drift_factor is linked to the wind and current used for the retrieval,\nother forcing datasets will yeld different value of the wind_drift_factor.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(wind_drift_factor.mean())\n\nplt.hist(wind_drift_factor, label='Retrieved wind_drift_factor')\nplt.axvline(x=0.033, label='Actual wind_drift_factor of 0.033', color='k')\nplt.axvline(x=wind_drift_factor.mean(), label='Mean retieved wind_drift_factor of %.3f' % wind_drift_factor.mean(), color='r')\nplt.ylabel('Number')\nplt.xlabel('Wind drift factor  [fraction]')\nplt.legend(loc='lower left')\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "A polar 2D histogram showing also the azimuth offset direction of the retrieved wind drift factor.\nSee Sutherland et al. (2020), https://doi.org/10.1175/JTECH-D-20-0013.1\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "wmax = wind_drift_factor.max()\nwbins = np.arange(0, wmax+.005, .005)\nabins = np.linspace(-180, 180, 30)\nhist, _, _ = np.histogram2d(azimuth, wind_drift_factor, bins=(abins, wbins))\nA, W = np.meshgrid(abins, wbins)\nfig, ax = plt.subplots(subplot_kw=dict(projection='polar'))\nax.set_theta_zero_location('N', offset=0)\nax.set_theta_direction(-1)\npc = ax.pcolormesh(np.radians(A), W, hist.T, cmap=cmocean.cm.dense)\nplt.arrow(np.pi, wmax, 0, -wmax, width=.015, facecolor='k', zorder=100,\n          head_width=.8, lw=2, head_length=.005, length_includes_head=True)\nplt.text(np.pi, wmax*.5, ' Wind direction', fontsize=18)\nplt.grid()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Alternative method, using skillscore\n\nHere we calculate trajectories for a range of wind_drift_factors,\nand calculate the Liu-Weissberg skillscore for each.\nThe optimized wind_drift_factor corresponds to the maximum skillscore.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "o = OceanDrift(loglevel=50)\no.add_readers_from_list([o.test_data_folder() +\n    '16Nov2015_NorKyst_z_surface/norkyst800_subset_16Nov2015.nc',\n    o.test_data_folder() + '16Nov2015_NorKyst_z_surface/arome_subset_16Nov2015.nc'], lazy=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Calulating trajectories for 100 different wind_drift_factors between 0 and 0.05\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "wdf = np.linspace(0.0, 0.05, 100)\no.seed_elements(lon=4, lat=60, time=ot.env.readers[list(ot.env.readers)[0]].start_time,\n                wind_drift_factor=wdf, number=len(wdf))\no.run(duration=timedelta(hours=12), time_step=600)\no.plot(linecolor='wind_drift_factor', drifter=drifter)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plotting and finding the wind_drift_factor which maximises the skillscore\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "skillscore = o.skillscore_trajectory(drifter_lons, drifter_lats, drifter_times, tolerance_threshold=1)\nind = np.argmax(skillscore)\nplt.plot(wdf, skillscore)\nplt.xlabel('Wind drift factor  [fraction]')\nplt.ylabel('Liu-Weissberg skillscore')\nplt.title('Maximum skillscore %.3f for wdf=%.3f' % (skillscore[ind], wdf[ind]))\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We see that using skillscore from the full trajectories gives a better estimate\nthan calculating the average of the wind_drift_factor from \none position to the next (i.e. polar histogram above).\nThis is even more clear if increasing the diffusivity (i.e. noise) above from 10 m2/s to 200 m2/s:\nThe histogram method then gives 0.071, which is much to high (true is 0.033), and the histogram is noisy.\nThe skillscore method is still robust, and gives a `wind_drift_factor` of 0.036, only slightly too high.\n\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}