{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Using wind measurements\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from datetime import datetime, timedelta, date\nfrom opendrift.readers import reader_netCDF_CF_generic\nfrom opendrift.models.openoil import OceanDrift\nfrom opendrift.readers import reader_timeseries\nfrom opendrift.readers.operators import ops\nimport opendrift.models.physics_methods as phy\nimport xarray as xr\nimport numpy as np\nimport pyproj\nimport matplotlib.pyplot as plt\nimport matplotlib as mpl\n\n# Preparation of figsize, dpi, fontsize and font style\n# Could also be done using mpl style sheet\n# https://matplotlib.org/stable/users/explain/customizing.html\nmpl.rcParams['figure.figsize'] = (5, 5)\nmpl.rcParams['savefig.dpi'] = 300\nmpl.rcParams['font.family'] = \"sans-serif\"\nmpl.rcParams['font.sans-serif'] = \"DejaVu Sans\"\nmpl.rcParams['font.size'] = 22\nmpl.rcParams['mathtext.fontset'] = \"dejavusans\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This example shows how to use wind or currents measurments \ncombined with an ocean model to improove the drift\ncalculation arround the measurment station.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Fisrt, we import two times the same model, because one\nwill be used for the usual simulation, and one will\nbe used for the measurement combined simulation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reader_norkyst = reader_netCDF_CF_generic.Reader('https://thredds.met.no/thredds/dodsC/sea/norkyst800m/1h/aggregate_be')\nreader_norkyst2 = reader_netCDF_CF_generic.Reader('https://thredds.met.no/thredds/dodsC/sea/norkyst800m/1h/aggregate_be')\n\n#Verification\nprint(reader_norkyst)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We then import wind measurement data,\nthat will be stored in a timeserie reader.\nDon't forget to give the lon and lat of the \nmeasurement station.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "time_list = [datetime(2024, 5, 27, 6, 00, 00) + i*timedelta(minutes=10) for i in range(7)] # Manual import for this example\nwind_speed_array = np.array([1.9, 3.0, 3.5, 3.4, 3.9, 2.8, 2.5])\nwind_speed_angle = np.array([22, 25, 24, 19, 17, 30, 16])\n\n\n\n# Alternative if data are in a csv file\n\n# data_wind_bridge = np.genfromtxt(\"Simulation_data/data_wind_21_mai_2024.txt\", dtype = '<U10', skip_header=1)\n# temporal_time = np.array([data_wind_bridge[i, 0] + \"T\" + data_wind_bridge[i, 1] for i in range(len(data_wind_bridge))])\n# time_list = [datetime.fromisoformat(temporal_time[i]) for i in range(len(temporal_time))]\n# wind_speed_array = data_wind_bridge[:, 3]\n# wind_speed_angle = data_wind_bridge[:, 5]\n# wind_speed_array = wind_speed_array.astype(float)\n# wind_speed_angle = wind_speed_angle.astype(float)\n\n\n#Calculation of wind_x and wind_y\nwind_speed_angle = np.deg2rad(wind_speed_angle)\nwind_x = wind_speed_array * np.cos(wind_speed_angle)\nwind_y = wind_speed_array * np.sin(wind_speed_angle)\n\n#Creation of the wind measurement reader\nlon, lat = 5.17, 60.37\nparameter_value_map = {'time': time_list, 'wind_speed': wind_speed_array, 'x_wind': wind_x, 'y_wind':  wind_y, 'lon':lon, 'lat':lat}\nwind_obs_reader = reader_timeseries.Reader(parameter_value_map)\n\n#Preparing future plots\ntext = [{'s': '* Station', 'x': lon, 'y': lat, 'fontsize': 20, 'color': 'k', 'backgroundcolor': 'white', 'bbox': dict(facecolor='none', edgecolor='none', alpha=0.4), 'zorder': 1000}]\n\n#Verification\nprint(wind_obs_reader)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we combine the second norkyst reader \nwith the measurement one, using a gaussian \ncombined reader.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "c = reader_norkyst2.combine_gaussian(wind_obs_reader, std = 1000)\nprint(c)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Configuration of the simulations\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "time_start = datetime(2024, 5, 27, 6, 00, 00)\ntime_end = time_start + timedelta(hours = 2) \ntime = [time_start, time_end]\ntime_step = 150\ntime_step_output = 150\n\n# Seeding configuration\nlon = 5.187\nlat = 60.385\nradius = 0\nnumber = 1\nwind_drift_factor = 0.03"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Usual simulation\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "o = OceanDrift(loglevel=20)\no.add_reader([reader_norkyst])\n\n# o.set_config('drift:current_uncertainty', .1)\n# o.set_config('drift:wind_uncertainty', 1)\n\no.seed_elements(lon=lon, lat=lat, radius=radius, number=number, time=time_start, wind_drift_factor=wind_drift_factor)\no.run(end_time=time_end, time_step=time_step, time_step_output=time_step_output)\no.plot(fast = False, background=['x_sea_water_velocity', 'y_sea_water_velocity'], legend=['Norkyst only', 'Gaussian measurement'], buffer = .023, markersize = 150, linewidth = 3, title = \"\", xlocs = mpl.ticker.MaxNLocator(5), ylocs = mpl.ticker.MaxNLocator(5), clabel = r\"Wind speed $\\mathrm{(m.s^{-1})}$\", cpad = 0.08, text=text)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Measurement modified model simulation\nand comparison\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "o2 = OceanDrift(loglevel=20)\no2.set_config('drift:max_speed', 10) #Necessary for the moment because of a bug to be solved\no2.add_reader([c, reader_norkyst])\n\n#o2.set_config('drift:current_uncertainty', .1)\n#o2.set_config('drift:wind_uncertainty', 1)\n\no2.seed_elements(lon=lon, lat=lat, radius=radius, number=number, time=time_start, wind_drift_factor=wind_drift_factor)\no2.run(end_time=time_end, time_step=time_step, time_step_output=time_step_output)\n\no.plot(fast = False, compare=o2, background=['x_sea_water_velocity', 'y_sea_water_velocity'], legend=['Norkyst only', 'Gaussian measurement'], buffer = .023, markersize = 150, linewidth = 3, title = \"\", xlocs = mpl.ticker.MaxNLocator(5), ylocs = mpl.ticker.MaxNLocator(5), clabel = r\"Wind speed $\\mathrm{(m.s^{-1})}$\", cpad = 0.08, text=text)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here, we generate more particles and look\nat the skillscore impact of the measurment \ncombined simulation\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "radius = 100\nnumber = 100\n  \no = OceanDrift(loglevel=20)\no.add_reader([reader_norkyst])\no.seed_elements(lon=lon, lat=lat, radius=radius, number=number, time=time_start, wind_drift_factor=wind_drift_factor)\no.run(end_time=time_end, time_step=time_step, time_step_output=time_step_output)\n\no2 = OceanDrift(loglevel=20)\no2.set_config('drift:max_speed', 10) #Necessary for the moment because of a bug to be solved\no2.add_reader([c, reader_norkyst])\no2.seed_elements(lon=lon, lat=lat, radius=radius, number=number, time=time_start, wind_drift_factor=wind_drift_factor)\no2.run(end_time=time_end, time_step=time_step, time_step_output=time_step_output)\n\no.plot(fast = False, compare=o2, background=['x_sea_water_velocity', 'y_sea_water_velocity'], legend=['Norkyst only', 'Gaussian measurement'], buffer = .023, markersize = 70, linewidth = 1, title = \"\", xlocs = mpl.ticker.MaxNLocator(5), ylocs = mpl.ticker.MaxNLocator(5), clabel = r\"Wind speed $\\mathrm{(m.s^{-1})}$\", cpad = 0.08, text=text)\n\nlon_drifters_1, lat_drifters_1 = o.get_lonlats()\nlon_drifters_2, lat_drifters_2 = o2.get_lonlats()\n\nskillscores = phy.skillscore_liu_weissberg(lon_drifters_1, lat_drifters_1, lon_drifters_2, lat_drifters_2)\n\nplt.figure()\nplt.hist(skillscores, bins = 100, range = (0, 1), facecolor = 'none', edgecolor = 'C0')\nplt.xlabel(\"Skillscore\")\nplt.ylabel(\"Histogram (percentages)\")\nplt.show()\n\nplt.figure()\nplt.hist(skillscores, bins = 100, range = (0.9, 0.95), facecolor = 'none', edgecolor = 'C0')\nplt.xlabel(\"Skillscore\")\nplt.ylabel(\"Histogram (percentages)\")\nplt.show()"
      ]
    }
  ],
  "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
}