{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot Antarctic Tidal Currents\n",
    "=============================\n",
    "\n",
    "Demonstrates plotting hourly tidal currents around Antarctica\n",
    "\n",
    "OTIS format tidal solutions provided by Ohio State University and ESR  \n",
    "- http://volkov.oce.orst.edu/tides/region.html  \n",
    "- https://www.esr.org/research/polar-tide-models/list-of-polar-tide-models/\n",
    "- ftp://ftp.esr.org/pub/datasets/tmd/  \n",
    "\n",
    "Finite Element Solution (FES) provided by AVISO  \n",
    "- https://www.aviso.altimetry.fr/en/data/products/auxiliary-products/global-tide-fes.html\n",
    "\n",
    "#### Python Dependencies\n",
    " - [numpy: Scientific Computing Tools For Python](https://www.numpy.org)  \n",
    " - [scipy: Scientific Tools for Python](https://www.scipy.org/)  \n",
    " - [pyproj: Python interface to PROJ library](https://pypi.org/project/pyproj/)  \n",
    " - [netCDF4: Python interface to the netCDF C library](https://unidata.github.io/netcdf4-python/)  \n",
    " - [matplotlib: Python 2D plotting library](http://matplotlib.org/)  \n",
    " - [cartopy: Python package designed for geospatial data processing](https://scitools.org.uk/cartopy/docs/latest/)  \n",
    "\n",
    "#### Program Dependencies\n",
    "\n",
    "- `arguments.py`: load the nodal corrections for tidal constituents  \n",
    "- `astro.py`: computes the basic astronomical mean longitudes  \n",
    "- `crs.py`: Coordinate Reference System (CRS) routines  \n",
    "- `io.model.py`: retrieves tide model parameters for named tide models\n",
    "- `io.OTIS.py`: extract tidal harmonic constants from OTIS tide models  \n",
    "- `io.ATLAS.py`: extract tidal harmonic constants from ATLAS netcdf models  \n",
    "- `io.FES.py`: extract tidal harmonic constants from FES tide models  \n",
    "- `predict.py`: predict tidal values using harmonic constants  \n",
    "- `time.py`: utilities for calculating time operations\n",
    "- `utilities.py`: download and management utilities for files\n",
    "\n",
    "This notebook uses Jupyter widgets to set parameters for calculating the tidal maps.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pyproj\n",
    "import datetime\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "matplotlib.rcParams['axes.linewidth'] = 2.0\n",
    "matplotlib.rcParams[\"animation.html\"] = \"jshtml\"\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "import cartopy.crs as ccrs\n",
    "from IPython.display import HTML\n",
    "\n",
    "# import tide programs\n",
    "import pyTMD.io\n",
    "import pyTMD.predict\n",
    "import pyTMD.tools\n",
    "import pyTMD.utilities\n",
    "import timescale.time\n",
    "\n",
    "# autoreload\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set parameters  for program\n",
    "\n",
    "- Model directory  \n",
    "- Tide model  \n",
    "- Date to run  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# available model list\n",
    "model_list = sorted(pyTMD.io.model.global_current() + pyTMD.io.model.antarctic_current())\n",
    "# display widgets for setting directory and model\n",
    "TMDwidgets = pyTMD.tools.widgets()\n",
    "TMDwidgets.model.options = model_list\n",
    "TMDwidgets.model.value = 'CATS2008'\n",
    "TMDwidgets.VBox([\n",
    "    TMDwidgets.directory,\n",
    "    TMDwidgets.model,\n",
    "    TMDwidgets.atlas,\n",
    "    TMDwidgets.compress,\n",
    "    TMDwidgets.datepick\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Setup tide model parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get model parameters\n",
    "model = pyTMD.io.model(TMDwidgets.directory.value,\n",
    "    format=TMDwidgets.atlas.value,\n",
    "    compressed=TMDwidgets.compress.value\n",
    "   ).current(TMDwidgets.model.value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Setup coordinates for calculating tidal currents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create an image around Antarctica\n",
    "xlimits = [-560.*5e3,560.*5e3]\n",
    "ylimits = [-560.*5e3,560.*5e3]\n",
    "spacing = [20e3,-20e3]\n",
    "# x and y coordinates\n",
    "x = np.arange(xlimits[0],xlimits[1]+spacing[0],spacing[0])\n",
    "y = np.arange(ylimits[1],ylimits[0]+spacing[1],spacing[1])\n",
    "xgrid,ygrid = np.meshgrid(x,y)\n",
    "# x and y dimensions\n",
    "nx = int((xlimits[1]-xlimits[0])/spacing[0])+1\n",
    "ny = int((ylimits[0]-ylimits[1])/spacing[1])+1\n",
    "# convert image coordinates from polar stereographic to latitude/longitude\n",
    "crs1 = pyproj.CRS.from_epsg(3031)\n",
    "crs2 = pyproj.CRS.from_epsg(4326)\n",
    "transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True)\n",
    "lon,lat = transformer.transform(xgrid.flatten(), ygrid.flatten())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate tide map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert from calendar date to days relative to Jan 1, 1992 (48622 MJD)\n",
    "YMD = TMDwidgets.datepick.value\n",
    "ts = timescale.time.Timescale().from_calendar(YMD.year, YMD.month,\n",
    "    YMD.day, hour=np.arange(24))\n",
    "\n",
    "# save tide currents\n",
    "tide = {}\n",
    "# iterate over u and v currents\n",
    "for TYPE in model.type:\n",
    "    # read tidal constants and interpolate to grid points\n",
    "    if model.format in ('OTIS','ATLAS','TMD3'):\n",
    "        amp,ph,D,c = pyTMD.io.OTIS.extract_constants(lon, lat, model.grid_file,\n",
    "            model.model_file['u'], model.projection, type=TYPE, crop=True,\n",
    "            method='spline', grid=model.format)\n",
    "        DELTAT = np.zeros_like(ts.tide)\n",
    "    elif (model.format == 'netcdf'):\n",
    "        amp,ph,D,c = pyTMD.io.ATLAS.extract_constants(lon, lat, model.grid_file,\n",
    "            model.model_file[TYPE], type=TYPE, crop=True, method='spline',\n",
    "            scale=model.scale, compressed=model.compressed)\n",
    "        DELTAT = np.zeros_like(ts.tide)\n",
    "    elif (model.format == 'GOT'):\n",
    "        amp,ph,c = pyTMD.io.GOT.extract_constants(lon, lat, model.model_file[TYPE],\n",
    "            crop=True, method='spline', scale=model.scale,\n",
    "            compressed=model.compressed)\n",
    "        # delta time (TT - UT1)\n",
    "        DELTAT = ts.tt_ut1\n",
    "    elif (model.format == 'FES'):\n",
    "        amp,ph = pyTMD.io.FES.extract_constants(lon, lat, model.model_file[TYPE],\n",
    "            type=TYPE, version=model.version, crop=True, method='spline',\n",
    "            scale=model.scale, compressed=model.compressed)\n",
    "        c = model.constituents\n",
    "        # delta time (TT - UT1)\n",
    "        DELTAT = ts.tt_ut1\n",
    "\n",
    "    # calculate complex phase in radians for Euler's\n",
    "    cph = -1j*ph*np.pi/180.0\n",
    "    # calculate constituent oscillation\n",
    "    hc = amp*np.exp(cph)\n",
    "\n",
    "    # allocate for tide current map calculated every hour\n",
    "    tide[TYPE] = np.ma.zeros((ny,nx,24))\n",
    "    for hour in range(24):\n",
    "        # predict tidal elevations at time and infer minor corrections\n",
    "        TIDE = pyTMD.predict.map(ts.tide[hour], hc, c, deltat=DELTAT[hour],\n",
    "            corrections=model.format)\n",
    "        MINOR = pyTMD.predict.infer_minor(ts.tide[hour], hc, c,\n",
    "            deltat=DELTAT[hour], corrections=model.format)\n",
    "        # add major and minor components and reform grid\n",
    "        tide[TYPE][:,:,hour] = np.reshape((TIDE+MINOR),(ny,nx))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create animation of hourly tidal oscillation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# output Antarctic Tidal Current Animation\n",
    "projection = ccrs.Stereographic(central_longitude=0.0,\n",
    "    central_latitude=-90.0,true_scale_latitude=-71.0)\n",
    "# figure axis and image objects\n",
    "ax1,im = ({},{})\n",
    "fig, (ax1['u'],ax1['v']) = plt.subplots(num=1, ncols=2,\n",
    "    figsize=(11.5,7), subplot_kw=dict(projection=projection))\n",
    "vmin = np.min([tide['u'].min(),tide['v'].min()])\n",
    "vmax = np.max([tide['u'].max(),tide['v'].max()])\n",
    "extent = (xlimits[0],xlimits[1],ylimits[0],ylimits[1])\n",
    "for TYPE,ax in ax1.items():\n",
    "    # plot tidal currents\n",
    "    im[TYPE] = ax.imshow(np.zeros((ny,nx)),\n",
    "        interpolation='nearest', vmin=vmin, vmax=vmax,\n",
    "        transform=projection, extent=extent, origin='upper',\n",
    "        animated=True)\n",
    "    # add high resolution cartopy coastlines\n",
    "    ax.coastlines('10m')\n",
    "    # set x and y limits\n",
    "    ax.set_xlim(xlimits)\n",
    "    ax.set_ylim(ylimits)\n",
    "    # stronger linewidth on frame\n",
    "    ax.spines['geo'].set_linewidth(2.0)\n",
    "    ax.spines['geo'].set_capstyle('projecting')\n",
    "\n",
    "# Add colorbar with a colorbar axis\n",
    "# Add an axes at position rect [left, bottom, width, height]\n",
    "cbar_ax = fig.add_axes([0.085, 0.075, 0.83, 0.035])\n",
    "# extend = add extension triangles to upper and lower bounds\n",
    "# options: neither, both, min, max\n",
    "cbar = fig.colorbar(im['u'], cax=cbar_ax, extend='both',\n",
    "    extendfrac=0.0375, drawedges=False, orientation='horizontal')\n",
    "# rasterized colorbar to remove lines\n",
    "cbar.solids.set_rasterized(True)\n",
    "# Add label to the colorbar\n",
    "cbar.ax.set_title(f'{model.name} Tidal Velocity', fontsize=13,\n",
    "    rotation=0, y=-1.65, va='top')\n",
    "cbar.ax.set_xlabel('cm/s', fontsize=13, rotation=0, va='center')\n",
    "cbar.ax.xaxis.set_label_coords(1.075, 0.5)\n",
    "# ticks lines all the way across\n",
    "cbar.ax.tick_params(which='both', width=1, length=18,\n",
    "    labelsize=13, direction='in')\n",
    "\n",
    "# add title (date and time)\n",
    "ttl = fig.suptitle(None, y=0.97, fontsize=13)\n",
    "# adjust subplot within figure\n",
    "fig.subplots_adjust(left=0.02,right=0.98,bottom=0.1,top=0.98,wspace=0.04)\n",
    "           \n",
    "# animate each map\n",
    "def animate_maps(hour):\n",
    "    # set map data iterating over u and v currents\n",
    "    for TYPE in model.type:\n",
    "        im[TYPE].set_data(tide[TYPE][:,:,hour])\n",
    "    # set title\n",
    "    args = (YMD.year,YMD.month,YMD.day,hour)\n",
    "    ttl.set_text('{0:4d}-{1:02d}-{2:02d}T{3:02d}:00:00'.format(*args))\n",
    "\n",
    "# set animation\n",
    "anim = animation.FuncAnimation(fig, animate_maps, frames=24)\n",
    "%matplotlib inline\n",
    "HTML(anim.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.4 64-bit",
   "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.10.12"
  },
  "vscode": {
   "interpreter": {
    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}