{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot Antarctic Tide Range\n",
    "=========================\n",
    "\n",
    "Demonstrates plotting the tidal amplitudes 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",
    "Global Tide Model (GOT) solutions provided by Richard Ray at GSFC  \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",
    "- `convert_crs.py`: convert points to and from Coordinates Reference Systems  \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",
    "\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 numpy as np\n",
    "import matplotlib\n",
    "matplotlib.rcParams['axes.linewidth'] = 2.0\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import ipywidgets as widgets\n",
    "\n",
    "# import tide programs\n",
    "import pyTMD.io\n",
    "import pyTMD.tools\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  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# available model list\n",
    "model_list = sorted(pyTMD.io.model.global_ocean() + pyTMD.io.model.antarctic_ocean())\n",
    "# display widgets for setting directory and model\n",
    "TMDwidgets = pyTMD.tools.widgets()\n",
    "TMDwidgets.model.options = model_list\n",
    "TMDwidgets.model.value = 'GOT4.10'\n",
    "widgets.VBox([\n",
    "    TMDwidgets.directory,\n",
    "    TMDwidgets.model,\n",
    "    TMDwidgets.atlas,\n",
    "    TMDwidgets.compress\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",
    "   ).elevation(TMDwidgets.model.value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Setup coordinates for calculating tides"
   ]
  },
  {
   "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": [
    "#### Infer minor amplitudes from the major constituents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def infer_minor_amplitudes(zmajor,constituents):\n",
    "    # number of constituents\n",
    "    npts,nc = np.shape(zmajor)    \n",
    "    cindex = ['q1','o1','p1','k1','n2','m2','s2','k2']\n",
    "    # re-order zmajor to correspond to cindex\n",
    "    z8 = np.ma.zeros((npts,8))\n",
    "    ni = 0\n",
    "    for i,c in enumerate(cindex):\n",
    "        j = [j for j,val in enumerate(constituents) if val == c]\n",
    "        if j:\n",
    "            j1, = j\n",
    "            z8[:,i] = zmajor[:,j1]\n",
    "            ni += 1\n",
    "    # list of minor constituents\n",
    "    minor = ['2q1','sigma1','rho1','m1','m1','chi1','pi1','phi1','theta1',\n",
    "        'j1','oo1','2n2','mu2','nu2','lambda2','l2','l2','t2']\n",
    "    # only add minor constituents that are not on the list of major values\n",
    "    minor_flag = [m not in constituents for m in minor]\n",
    "    # estimate minor constituents\n",
    "    zmin = np.zeros((npts,18))\n",
    "    zmin[:,0] = 0.263*z8[:,0] - 0.0252*z8[:,1]# 2Q1\n",
    "    zmin[:,1] = 0.297*z8[:,0] - 0.0264*z8[:,1]# sigma1\n",
    "    zmin[:,2] = 0.164*z8[:,0] + 0.0048*z8[:,1]# rho1\n",
    "    zmin[:,3] = 0.0140*z8[:,1] + 0.0101*z8[:,3]# M1\n",
    "    zmin[:,4] = 0.0389*z8[:,1] + 0.0282*z8[:,3]# M1\n",
    "    zmin[:,5] = 0.0064*z8[:,1] + 0.0060*z8[:,3]# chi1\n",
    "    zmin[:,6] = 0.0030*z8[:,1] + 0.0171*z8[:,3]# pi1\n",
    "    zmin[:,7] = -0.0015*z8[:,1] + 0.0152*z8[:,3]# phi1\n",
    "    zmin[:,8] = -0.0065*z8[:,1] + 0.0155*z8[:,3]# theta1\n",
    "    zmin[:,9] = -0.0389*z8[:,1] + 0.0836*z8[:,3]# J1\n",
    "    zmin[:,10] = -0.0431*z8[:,1] + 0.0613*z8[:,3]# OO1\n",
    "    zmin[:,11] = 0.264*z8[:,4] - 0.0253*z8[:,5]# 2N2\n",
    "    zmin[:,12] = 0.298*z8[:,4] - 0.0264*z8[:,5]# mu2\n",
    "    zmin[:,13] = 0.165*z8[:,4] + 0.00487*z8[:,5]# nu2\n",
    "    zmin[:,14] = 0.0040*z8[:,5] + 0.0074*z8[:,6]# lambda2\n",
    "    zmin[:,15] = 0.0131*z8[:,5] + 0.0326*z8[:,6]# L2\n",
    "    zmin[:,16] = 0.0033*z8[:,5] + 0.0082*z8[:,6]# L2\n",
    "    zmin[:,17] = 0.0585*z8[:,6]# t2\n",
    "    # only add minor constituents that are not major\n",
    "    return np.where(minor_flag, np.abs(zmin), 0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate tide map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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, model.projection, type=model.type,\n",
    "        method='spline', grid=model.format)\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=model.type, method='spline',\n",
    "        scale=model.scale, compressed=model.compressed)\n",
    "elif (model.format == 'GOT'):\n",
    "    amp,ph,c = pyTMD.io.GOT.extract_constants(lon, lat, model.model_file,\n",
    "        method='spline', scale=model.scale,\n",
    "        compressed=model.compressed)\n",
    "elif (model.format == 'FES'):\n",
    "    amp,ph = pyTMD.io.FES.extract_constants(lon, lat, model.model_file,\n",
    "        type=model.type, version=model.version, method='spline',\n",
    "        scale=model.scale, compressed=model.compressed)\n",
    "    c = model.constituents\n",
    "\n",
    "# calculate minor constituent amplitudes\n",
    "minor_amp = infer_minor_amplitudes(amp,c)\n",
    "# calculate sum of major and minor amplitudes\n",
    "tide_range = np.sum(amp,axis=1) + np.sum(minor_amp,axis=1)\n",
    "# convert from meters to centimeters\n",
    "tide_cm = 100.0*np.reshape(tide_range,(ny,nx))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create plot of tidal range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot Antarctic tide range\n",
    "projection = ccrs.Stereographic(central_longitude=0.0,\n",
    "    central_latitude=-90,true_scale_latitude=-71.0)\n",
    "fig, ax = plt.subplots(num=1, figsize=(9,8),\n",
    "    subplot_kw=dict(projection=projection))\n",
    "# plot tide height\n",
    "vmin,vmax = (0, np.max(tide_cm))\n",
    "extent = (xlimits[0],xlimits[1],ylimits[0],ylimits[1])\n",
    "im = ax.imshow(tide_cm, interpolation='nearest',\n",
    "    vmin=vmin, vmax=vmax, transform=projection,\n",
    "    extent=extent, origin='upper')\n",
    "# add high resolution cartopy coastlines\n",
    "ax.coastlines('10m')\n",
    "\n",
    "# Add colorbar and adjust size\n",
    "# pad = distance from main plot axis\n",
    "# extend = add extension triangles to upper and lower bounds\n",
    "# options: neither, both, min, max\n",
    "# shrink = percent size of colorbar\n",
    "# aspect = lengthXwidth aspect of colorbar\n",
    "cbar = plt.colorbar(im, ax=ax, pad=0.025, extend='max',\n",
    "    extendfrac=0.0375, shrink=0.85, aspect=22.5, drawedges=False)\n",
    "# rasterized colorbar to remove lines\n",
    "cbar.solids.set_rasterized(True)\n",
    "# Add label to the colorbar\n",
    "cbar.ax.set_ylabel(f'{model.name} Tide Range', labelpad=10, fontsize=13)\n",
    "cbar.ax.set_title('cm', fontsize=13, va='bottom')\n",
    "# ticks lines all the way across\n",
    "cbar.ax.tick_params(which='both', width=1, length=20,\n",
    "    labelsize=13, direction='in')\n",
    "\n",
    "# set x and y limits\n",
    "ax.set_xlim(xlimits)\n",
    "ax.set_ylim(ylimits)\n",
    "\n",
    "# stronger linewidth on frame\n",
    "ax.spines['geo'].set_linewidth(2.0)\n",
    "ax.spines['geo'].set_capstyle('projecting')\n",
    "# adjust subplot within figure\n",
    "fig.subplots_adjust(left=0.02,right=0.98,bottom=0.05,top=0.98)\n",
    "# show the plot\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.6 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
}