{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "name": "#%% Script to explore Mustuque data\n"
    }
   },
   "outputs": [],
   "source": [
    "# Alexander Rey, December 2021\n",
    "# 1 RBR, 2 Nortek Eco, GFS winds\n",
    "\n",
    "import os\n",
    "\n",
    "import pandas as pd\n",
    "import geopandas as gp\n",
    "import numpy as np\n",
    "import math\n",
    "import contextily as ctx\n",
    "import cmocean.cm as cmo\n",
    "\n",
    "from xarray.backends import NetCDF4DataStore\n",
    "import xarray as xr\n",
    "\n",
    "from datetime import datetime, timedelta\n",
    "from netCDF4 import num2date\n",
    "from metpy.units import units\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from siphon.catalog import TDSCatalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "outputs": [],
   "source": [
    "def find_time_var(var, time_basename='time'):\n",
    "    for coord_name in var.coords:\n",
    "        if coord_name.startswith(time_basename):\n",
    "            return var.coords[coord_name]\n",
    "    raise ValueError('No time variable found for ' + var.name)\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%  Helper function for finding proper time variable\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "outputs": [],
   "source": [
    "eco_path_1  = \"//srv-ott3/Projects/13539.101 L'Ansecoy Bay, Mustique/03_Data/03_Field/01_September 2021 trip/02_Nortek ECO/Eco61_20210910184606.csv\"\n",
    "eco_path_2  = \"//srv-ott3/Projects/13539.101 L'Ansecoy Bay, Mustique/03_Data/03_Field/01_September 2021 trip/02_Nortek ECO/Eco214_20210902120703.csv\"\n",
    "\n",
    "eco1_pd     = pd.read_csv(eco_path_1, sep=',', header=[28], index_col=0,\n",
    "                          parse_dates=False)\n",
    "\n",
    "eco2_pd     = pd.read_csv(eco_path_2, sep=',', header=[28], index_col=0,\n",
    "                          parse_dates=False)\n",
    "\n",
    "# Drop first row\n",
    "eco1_pd.drop(eco1_pd.index[0], inplace=True)\n",
    "eco2_pd.drop(eco2_pd.index[0], inplace=True)\n",
    "\n",
    "eco1_pd.index = pd.to_datetime(eco1_pd.index, format='%m/%d/%Y %I:%M:%S %p')\n",
    "eco2_pd.index = pd.to_datetime(eco2_pd.index, format='%m/%d/%Y %I:%M:%S %p')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Import Nortek Eco data\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "outputs": [],
   "source": [
    "rbr_path  = \"C:/Users/arey/files/Projects/Mustique/041279_20211203_1541Ruskin.xlsx\"\n",
    "rbr_pd    = pd.read_excel(rbr_path, sheet_name='Wave', parse_dates=True, header=1, index_col=0)\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Import RBR data\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "outputs": [],
   "source": [
    "minTime  = eco1_pd.index[0]\n",
    "metDate  = minTime - timedelta(\n",
    "                            hours = minTime.hour % 6,\n",
    "                            minutes=minTime.minute,\n",
    "                            seconds=minTime.second,\n",
    "                            microseconds=minTime.microsecond)\n",
    "\n",
    "GFS_Lat = 12.890861\n",
    "GFS_Lon = -61.1791"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Setup GFS Import\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "'gfs_4_20211205_0600_000.grb2'",
     "output_type": "error",
     "traceback": [
      "\u001B[1;31m---------------------------------------------------------------------------\u001B[0m",
      "\u001B[1;31mKeyError\u001B[0m                                  Traceback (most recent call last)",
      "\u001B[1;32m~\\AppData\\Local\\Temp/ipykernel_18196/2112783647.py\u001B[0m in \u001B[0;36m<module>\u001B[1;34m\u001B[0m\n\u001B[0;32m     18\u001B[0m \u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m     19\u001B[0m     \u001B[1;31m# Generate URLs for specific grib file\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[1;32m---> 20\u001B[1;33m     \u001B[0mbest_ds\u001B[0m         \u001B[1;33m=\u001B[0m \u001B[0mbest_gfs\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mdatasets\u001B[0m\u001B[1;33m[\u001B[0m\u001B[1;34m'gfs_4_'\u001B[0m\u001B[1;33m+\u001B[0m\u001B[0mdate\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mstrftime\u001B[0m\u001B[1;33m(\u001B[0m\u001B[1;34m'%Y%m%d_%H%M'\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m+\u001B[0m\u001B[1;34m'_000.grb2'\u001B[0m\u001B[1;33m]\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0m\u001B[0;32m     21\u001B[0m \u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m     22\u001B[0m     \u001B[1;31m# Format the query parameters\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n",
      "\u001B[1;32m~\\Anaconda3\\envs\\BairdBase\\lib\\site-packages\\siphon\\catalog.py\u001B[0m in \u001B[0;36m__getitem__\u001B[1;34m(self, item)\u001B[0m\n\u001B[0;32m     34\u001B[0m         \u001B[1;32mtry\u001B[0m\u001B[1;33m:\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m     35\u001B[0m             \u001B[0mitem\u001B[0m \u001B[1;33m+\u001B[0m \u001B[1;34m''\u001B[0m  \u001B[1;31m# Raises if item not a string\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[1;32m---> 36\u001B[1;33m             \u001B[1;32mreturn\u001B[0m \u001B[0msuper\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mIndexableMapping\u001B[0m\u001B[1;33m,\u001B[0m \u001B[0mself\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0m__getitem__\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mitem\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0m\u001B[0;32m     37\u001B[0m         \u001B[1;32mexcept\u001B[0m \u001B[0mTypeError\u001B[0m\u001B[1;33m:\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m     38\u001B[0m             \u001B[1;32mreturn\u001B[0m \u001B[0mlist\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mself\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mvalues\u001B[0m\u001B[1;33m(\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m[\u001B[0m\u001B[0mitem\u001B[0m\u001B[1;33m]\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n",
      "\u001B[1;31mKeyError\u001B[0m: 'gfs_4_20211205_0600_000.grb2'"
     ]
    }
   ],
   "source": [
    "var = ['u-component_of_wind_height_above_ground', 'v-component_of_wind_height_above_ground']\n",
    "\n",
    "wndu_1d = []\n",
    "wndv_1d = []\n",
    "time_1d = []\n",
    "\n",
    "# Set times to download\n",
    "startdate  = metDate - timedelta(hours=18)\n",
    "enddate    = metDate + timedelta(days=90)\n",
    "\n",
    "date_list  = pd.date_range(startdate, enddate, freq='6H')\n",
    "\n",
    "# Loop through dates\n",
    "for date in date_list:\n",
    "    # Base URL for 0.5 degree GFS data\n",
    "    best_gfs = TDSCatalog('https://www.ncei.noaa.gov/thredds/catalog/model-gfs-g4-anl-files/' +\n",
    "                        date.strftime('%Y%m') + '/' + date.strftime('%Y%m%d') + '/' + 'catalog.xml')\n",
    "\n",
    "    # Generate URLs for specific grib file\n",
    "    best_ds         = best_gfs.datasets['gfs_4_'+date.strftime('%Y%m%d_%H%M')+'_000.grb2']\n",
    "\n",
    "    # Format the query parameters\n",
    "    ncss  = best_ds.subset()\n",
    "    query = ncss.query()\n",
    "\n",
    "    # Extract data from specific point\n",
    "    query.lonlat_point(GFS_Lon, GFS_Lat).time(date)\n",
    "    query.accept('netcdf')\n",
    "    query.variables(var[0], var[1])\n",
    "    query.vertical_level(10)\n",
    "\n",
    "    data = ncss.get_data(query)\n",
    "    data = xr.open_dataset(NetCDF4DataStore(data), drop_variables='height_above_ground4')\n",
    "\n",
    "    wndu_3d = data[var[0]]\n",
    "    wndv_3d = data[var[1]]\n",
    "\n",
    "    # Read the individual point (with units) and append to the list\n",
    "    wndu_1d.append(wndu_3d.metpy.unit_array.squeeze())\n",
    "    wndv_1d.append(wndv_3d.metpy.unit_array.squeeze())\n",
    "    time_1d.append(find_time_var(wndu_3d))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% GFS Met Import\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'list' object has no attribute 'squeeze'",
     "output_type": "error",
     "traceback": [
      "\u001B[1;31m---------------------------------------------------------------------------\u001B[0m",
      "\u001B[1;31mAttributeError\u001B[0m                            Traceback (most recent call last)",
      "\u001B[1;32m~\\AppData\\Local\\Temp/ipykernel_18196/2073007942.py\u001B[0m in \u001B[0;36m<module>\u001B[1;34m\u001B[0m\n\u001B[0;32m     19\u001B[0m \u001B[0mmet_wdir\u001B[0m \u001B[1;33m=\u001B[0m  \u001B[0mnp\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mdegrees\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mnp\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0marctan2\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mwndu_np\u001B[0m\u001B[1;33m,\u001B[0m \u001B[0mwndv_np\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m)\u001B[0m \u001B[1;33m%\u001B[0m \u001B[1;36m360\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m     20\u001B[0m \u001B[1;33m\u001B[0m\u001B[0m\n\u001B[1;32m---> 21\u001B[1;33m \u001B[0mtime_1d\u001B[0m  \u001B[1;33m=\u001B[0m \u001B[0mtime_1d\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0msqueeze\u001B[0m\u001B[1;33m(\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0m",
      "\u001B[1;31mAttributeError\u001B[0m: 'list' object has no attribute 'squeeze'"
     ]
    }
   ],
   "source": [
    "time_np = np.asanyarray(time_1d).squeeze()\n",
    "wndu_np = np.empty(shape=(len(wndu_1d), 1), dtype=float)\n",
    "for x in range(0, len(wndu_1d)):\n",
    "    wndu_np[x] = wndu_1d[x].m.item(0)\n",
    "\n",
    "wndv_np = np.empty(shape=(len(wndv_1d), 1), dtype=float)\n",
    "for x in range(0, len(wndv_1d)):\n",
    "    wndv_np[x] = wndv_1d[x].m.item(0)\n",
    "\n",
    "time_np = np.empty(shape=(len(time_1d), 1), dtype=float)\n",
    "for x in range(0, len(time_1d)):\n",
    "    time_np[x] = time_1d[x].values[0]\n",
    "\n",
    "wndu_np = wndu_np.squeeze()\n",
    "wndv_np = wndv_np.squeeze()\n",
    "time_np = time_np.squeeze()\n",
    "\n",
    "met_wind =  np.sqrt(wndu_np ** 2 + wndv_np **2)\n",
    "met_wdir =  np.degrees(np.arctan2(wndu_np, wndv_np)) % 360"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Format GFS\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 648x936 with 5 Axes>",
      "image/png": "[encoded data, 175292 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 648x936 with 5 Axes>",
      "image/png": "[encoded data, 135908 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 648x936 with 5 Axes>",
      "image/png": "[encoded data, 129972 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 648x936 with 5 Axes>",
      "image/png": "[encoded data, 139612 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 648x936 with 5 Axes>",
      "image/png": "[encoded data, 140772 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 648x936 with 5 Axes>",
      "image/png": "[encoded data, 134972 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 648x936 with 5 Axes>",
      "image/png": "[encoded data, 132520 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 648x936 with 5 Axes>",
      "image/png": "[encoded data, 134416 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 648x936 with 5 Axes>",
      "image/png": "[encoded data, 138784 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 648x936 with 5 Axes>",
      "image/png": "[encoded data, 134028 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 648x936 with 5 Axes>",
      "image/png": "[encoded data, 109836 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 648x936 with 5 Axes>",
      "image/png": "[encoded data, 95628 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Loop through different dates\n",
    "dateSrt = datetime(2021, 9, 10)\n",
    "dateEnd = datetime(2021, 12, 6)\n",
    "datePlt = datetime(2021, 12, 5)\n",
    "pltCnt  = 0\n",
    "\n",
    "while  datePlt<(dateEnd+timedelta(weeks=0)):\n",
    "    fig, axes = plt.subplots(nrows=5, ncols=1, figsize=(9, 13), sharex=True)\n",
    "\n",
    "    axes[0].plot(pd.to_datetime(time_np), met_wind, color='tab:blue')\n",
    "    axes[0].set_title('Wind Speed [m/s]')\n",
    "    axes[0].set_ylabel('Wind Speed [m/s]')\n",
    "\n",
    "    rbr_pd['Significant wave height'].plot(ax=axes[1], color='tab:orange')\n",
    "    axes[1].set_title('RBR: Significant wave height [m]')\n",
    "    axes[1].set_ylabel('Significant wave height [m]')\n",
    "\n",
    "    rbr_pd['Depth'].plot(ax=axes[2], color='tab:green')\n",
    "    axes[2].set_title('RBR: Depth [m]')\n",
    "    axes[2].set_ylim(3.8, 4.6)\n",
    "    axes[2].set_ylabel('Depth [m]')\n",
    "\n",
    "    eco1_pd['Upper speed'].astype(float).plot(ax=axes[3], color='tab:red')\n",
    "    # eco1_pd['Upper speed'].astype(float).rolling('21600s').mean().plot(ax=axes[2])\n",
    "    axes[3].set_title('Nearshore Eco: Upper speed [m/s]')\n",
    "    axes[3].set_ylabel('Upper speed [m/s]')\n",
    "\n",
    "    eco2_pd['Upper speed'].astype(float).plot(ax=axes[4], color='tab:purple')\n",
    "    # eco2_pd['Upper speed'].astype(float).rolling('21600s').mean().plot(ax=axes[3])\n",
    "    axes[4].set_title('Offshore Eco: Upper speed [m/s]')\n",
    "    axes[4].set_ylabel('Upper speed [m/s]')\n",
    "\n",
    "    fig.patch.set_facecolor('white')\n",
    "    fig.autofmt_xdate()\n",
    "    axes[4].set_xlabel('')\n",
    "\n",
    "    axes[3].sharey(axes[4])\n",
    "\n",
    "    plt.setp(axes, xlim=(dateSrt, datePlt))\n",
    "\n",
    "    if pltCnt == 0:\n",
    "        # Set axis ticks\n",
    "        minor_ticks = pd.date_range(datetime(2021, 9, 10), datetime(2021, 12, 3), freq='7D')\n",
    "        for axIDX in range(0, 5):\n",
    "            axes[axIDX].set_xticks(minor_ticks, minor=True)\n",
    "            axes[axIDX].grid(which='minor', alpha=0.2)\n",
    "\n",
    "        saveName = 'C:/Users/arey/files/Projects/Mustique/Figures/Combined/Combined_Surface_Full.pdf'\n",
    "    else:\n",
    "        minor_ticks = pd.date_range(dateSrt, datePlt, freq='1D')\n",
    "        for axIDX in range(0, 5):\n",
    "            axes[axIDX].set_xticks(minor_ticks, minor=True)\n",
    "            axes[axIDX].grid(which='minor', alpha=0.2)\n",
    "            axes[axIDX].grid(which='major', alpha=0.2)\n",
    "\n",
    "        saveName = 'C:/Users/arey/files/Projects/Mustique/Figures/Combined/Combined_Surface_Week_' + str(pltCnt) + '.pdf'\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "    fig.savefig(saveName, bbox_inches='tight')\n",
    "\n",
    "    pltCnt  = pltCnt + 1\n",
    "    dateSrt = dateSrt + timedelta(weeks=1)\n",
    "    datePlt = dateSrt + timedelta(weeks=1)\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Plot Met Data\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "pycharm-78fc648f",
   "language": "python",
   "display_name": "PyCharm (Mustique)"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}