{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot ATLAS Compact\n",
    "==================\n",
    "\n",
    "Check if a given point is within a tide model domain\n",
    "\n",
    "OTIS format tidal solutions provided by Ohio State University and ESR  \n",
    "- http://volkov.oce.orst.edu/tides/region.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](https://matplotlib.org/)  \n",
    "\n",
    "#### Program Dependencies\n",
    "\n",
    "- `convert_ll_xy.py`: convert lat/lon points to and from projected coordinates  \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",
    "- `predict.py`: predict tidal values using harmonic constants  \n",
    "- `time.py`: utilities for calculating time operations\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": [
    "from __future__ import print_function\n",
    "\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import ipywidgets as widgets\n",
    "\n",
    "import pyTMD.time\n",
    "import pyTMD.tools\n",
    "import pyTMD.io\n",
    "\n",
    "# autoreload\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# available model list\n",
    "model_list = sorted(pyTMD.io.model.ATLAS_compact())\n",
    "# display widgets for setting directory and model\n",
    "TMDwidgets = pyTMD.tools.widgets()\n",
    "TMDwidgets.model.options = model_list\n",
    "TMDwidgets.model.value = 'TPXO8-atlas'\n",
    "TMDwidgets.compress.value = False\n",
    "widgets.VBox([\n",
    "    TMDwidgets.directory,\n",
    "    TMDwidgets.model,\n",
    "    TMDwidgets.compress,\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get model parameters\n",
    "model = pyTMD.io.model(TMDwidgets.directory.value,\n",
    "    compressed=TMDwidgets.compress.value\n",
    "   ).elevation(TMDwidgets.model.value)\n",
    "# read each constituent\n",
    "constituents,nc = pyTMD.read_tide_model.read_constituents(model.model_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if reading a global solution with localized solutions\n",
    "x0,y0,hz0,mz0,iob,dt,pmask,local = pyTMD.io.OTIS.read_atlas_grid(model.grid_file)\n",
    "xi,yi,hz = pyTMD.io.OTIS.combine_atlas_model(x0,y0,hz0,pmask,local,VARIABLE='depth')\n",
    "mz = pyTMD.io.OTIS.create_atlas_mask(x0,y0,mz0,local,VARIABLE='depth')\n",
    "# resample global solution to 2 arc-minute solution\n",
    "x30,y30,hz30 = pyTMD.io.OTIS.interpolate_atlas_model(x0,y0,hz0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "# plot the ATLAS mask\n",
    "fig,ax1 = plt.subplots(num=1, figsize=(8.25,5.25), dpi=120)\n",
    "ax1.imshow(mz, interpolation='nearest',\n",
    "    extent=(xi.min(),xi.max(),yi.min(),yi.max()),\n",
    "    vmin=0, vmax=1, origin='lower', cmap='gray_r', alpha=0.5)\n",
    "ax1.imshow(pmask, interpolation='nearest',\n",
    "    extent=(x0.min(),x0.max(),y0.min(),y0.max()),\n",
    "    vmin=0, vmax=1, origin='lower', cmap='Purples', alpha=0.5)\n",
    "# no ticks on the x and y axes\n",
    "ax1.get_xaxis().set_ticks([])\n",
    "ax1.get_yaxis().set_ticks([])\n",
    "# stronger linewidth on frame\n",
    "[i.set_linewidth(2.0) for i in ax1.spines.values()]\n",
    "# adjust subplot within figure\n",
    "fig.subplots_adjust(left=0.02,right=0.98,bottom=0.05,top=0.98)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# percent difference between grids\n",
    "percent = 100.0*(hz30 - hz)/hz30\n",
    "# plot the percent difference between ATLAS depth\n",
    "fig,ax2 = plt.subplots(num=2, figsize=(8.25,5.25), dpi=120)\n",
    "im = ax2.imshow(percent, interpolation='nearest',\n",
    "    extent=(xi.min(),xi.max(),yi.min(),yi.max()),\n",
    "    vmin=-40, vmax=40, origin='lower', cmap='PRGn')\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=ax2, pad=0.025, extend='both',\n",
    "    extendfrac=0.0375, orientation='horizontal', shrink=0.925,\n",
    "    aspect=22, drawedges=False)\n",
    "# rasterized colorbar to remove lines\n",
    "cbar.solids.set_rasterized(True)\n",
    "# Add label to the colorbar\n",
    "cbar.ax.set_xlabel(f'{model.name} Bathymetry Differences', fontsize=13)\n",
    "cbar.ax.set_ylabel('%', fontsize=13, rotation=0)\n",
    "cbar.ax.yaxis.set_label_coords(1.0325, 0.15)\n",
    "# ticks lines all the way across\n",
    "cbar.ax.tick_params(which='both', width=1, length=23,\n",
    "    labelsize=13, direction='in')\n",
    "# no ticks on the x and y axes\n",
    "ax2.get_xaxis().set_ticks([])\n",
    "ax2.get_yaxis().set_ticks([])\n",
    "# stronger linewidth on frame\n",
    "[i.set_linewidth(2.0) for i in ax2.spines.values()]\n",
    "# adjust subplot within figure\n",
    "fig.subplots_adjust(left=0.02,right=0.98,bottom=0.05,top=0.98)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# total amplitude difference for all constituents\n",
    "diff = np.zeros_like(hz30)\n",
    "power = np.zeros_like(hz30)\n",
    "# for each constituent\n",
    "for i,c in enumerate(constituents):\n",
    "    # if reading a global solution with localized solutions\n",
    "    z0,zlocal = pyTMD.io.OTIS.read_atlas_elevation(model.model_file,i,c)\n",
    "    xi,yi,z = pyTMD.io.OTIS.combine_atlas_model(x0,y0,z0,pmask,zlocal,VARIABLE='z')\n",
    "    # resample global solution to 2 arc-minute solution\n",
    "    x30,y30,z30 = pyTMD.io.OTIS.interpolate_atlas_model(x0,y0,z0)\n",
    "    # add to total amplitude difference\n",
    "    diff += (z30.real - z.real)**2 + (z30.imag - z.imag)**2\n",
    "    power += z30.real**2 + z30.imag**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the percent difference\n",
    "percent = 100.0*np.sqrt(diff/power)\n",
    "# plot the percent between ATLAS tidal amplitudes\n",
    "fig,ax3 = plt.subplots(num=3, figsize=(8.25,5.25), dpi=120)\n",
    "im = ax3.imshow(percent, interpolation='nearest',\n",
    "    extent=(xi.min(),xi.max(),yi.min(),yi.max()),\n",
    "    vmin=0, vmax=100, origin='lower', cmap='BuPu')\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=ax3, pad=0.025, extend='max',\n",
    "    extendfrac=0.0375, orientation='horizontal', shrink=0.925,\n",
    "    aspect=22, drawedges=False)\n",
    "# rasterized colorbar to remove lines\n",
    "cbar.solids.set_rasterized(True)\n",
    "# Add label to the colorbar\n",
    "cbar.ax.set_xlabel(f'{model.name} Tide Height Differences', fontsize=13)\n",
    "cbar.ax.set_ylabel('%', fontsize=13, rotation=0)\n",
    "cbar.ax.yaxis.set_label_coords(1.0325, 0.15)\n",
    "# ticks lines all the way across\n",
    "cbar.ax.tick_params(which='both', width=1, length=23,\n",
    "    labelsize=13, direction='in')\n",
    "# no ticks on the x and y axes\n",
    "ax3.get_xaxis().set_ticks([])\n",
    "ax3.get_yaxis().set_ticks([])\n",
    "# stronger linewidth on frame\n",
    "[i.set_linewidth(2.0) for i in ax3.spines.values()]\n",
    "# adjust subplot within figure\n",
    "fig.subplots_adjust(left=0.02,right=0.98,bottom=0.05,top=0.98)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
  },
  "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.10.6 (main, Nov 14 2022, 16:10:14) [GCC 11.3.0]"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}