{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check Tide Map\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",
    "- 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](https://matplotlib.org/)  \n",
    " - [ipyleaflet: Jupyter / Leaflet bridge enabling interactive maps](https://github.com/jupyter-widgets/ipyleaflet)  \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",
    "- `io.ATLAS.py`: extract tidal harmonic constants from netcdf models  \n",
    "- `io.GOT.py`: extract tidal harmonic constants from GSFC GOT models  \n",
    "- `io.FES.py`: extract tidal harmonic constants from FES tide models  \n",
    "- `convert_ll_xy.py`: converts lat/lon points to and from projected coordinates \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",
    "import IPython.display\n",
    "\n",
    "# import tide programs\n",
    "import pyTMD.io\n",
    "import pyTMD.time\n",
    "import pyTMD.tools\n",
    "import pyTMD.convert_ll_xy\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.ocean_elevation())\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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# default coordinates to use\n",
    "LAT,LON = (32.86710263,-117.25750387)\n",
    "m = pyTMD.tools.leaflet(center=(LAT,LON), zoom=12,\n",
    "    zoom_control=True, marker_control=True)\n",
    "# show map\n",
    "m.map"
   ]
  },
  {
   "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)\n",
    "    \n",
    "# read tidal constants and interpolate to grid points\n",
    "if model.format in ('OTIS','ATLAS','ESR'):\n",
    "    # if reading a single OTIS solution\n",
    "    xi,yi,hz,mz,iob,dt = pyTMD.io.OTIS.read_otis_grid(model.grid_file)\n",
    "elif (model.format == 'netcdf'):\n",
    "    # if reading a netCDF OTIS atlas solution\n",
    "    xi,yi,hz = pyTMD.io.ATLAS.read_netcdf_grid(model.grid_file,\n",
    "        compressed=model.compressed, type='z')\n",
    "    # invert bathymetry mask\n",
    "    mz = np.invert(hz.mask)\n",
    "elif (model.format == 'GOT'):\n",
    "    # if reading a NASA GOT solution\n",
    "    hc,xi,yi,c = pyTMD.io.GOT.read_ascii_file(model.model_file[0],\n",
    "        compressed=model.compressed)\n",
    "    # invert tidal constituent mask\n",
    "    mz = np.invert(hc.mask)\n",
    "elif (model.format == 'FES'):\n",
    "    # if reading a FES netCDF solution\n",
    "    hc,xi,yi = pyTMD.io.FES.read_netcdf_file(model.model_file[0],\n",
    "        compressed=model.compressed, type='z', version=model.version)\n",
    "    # invert tidal constituent mask\n",
    "    mz = np.invert(hc.mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_coordinates(sender):\n",
    "    # leaflet location\n",
    "    LAT,LON = np.copy(m.marker.location)\n",
    "    # verify longitudes\n",
    "    LON = m.wrap_longitudes(LON)\n",
    "    # adjust dimensions of input coordinates to be iterable\n",
    "    LON = np.atleast_1d(LON)\n",
    "    LAT = np.atleast_1d(LAT)\n",
    "    # read tidal constants and interpolate to grid points\n",
    "    if model.format in ('OTIS','ATLAS','ESR'):\n",
    "        # if reading a single OTIS solution\n",
    "        xi,yi,hz,mz,iob,dt = pyTMD.io.OTIS.read_otis_grid(model.grid_file)\n",
    "        # convert coordinate systems of input latitude and longitude\n",
    "        x,y = pyTMD.convert_ll_xy(np.atleast_1d(LON), np.atleast_1d(LAT),\n",
    "            model.projection, 'F')\n",
    "        # adjust longitudinal convention of input latitude and longitude\n",
    "        # to fit tide model convention (if global)\n",
    "        if (np.min(x) < np.min(xi)) & (model.projection == '4326'):\n",
    "            lt0, = np.nonzero(x < 0)\n",
    "            x[lt0] += 360.0\n",
    "        if (np.max(x) > np.max(xi)) & (model.projection == '4326'):\n",
    "            gt180, = np.nonzero(x > 180)\n",
    "            x[gt180] -= 360.0\n",
    "    elif (model.format == 'netcdf'):\n",
    "        # if reading a netCDF OTIS atlas solution\n",
    "        # adjust longitudinal convention of input latitude and longitude\n",
    "        # to fit tide model convention\n",
    "        x,y = np.copy([LON,LAT]).astype(np.float64)\n",
    "        lt0, = np.nonzero(x < 0)\n",
    "        x[lt0] += 360.0\n",
    "    elif (model.format == 'GOT'):\n",
    "        # if reading a NASA GOT solution\n",
    "        # adjust longitudinal convention of input latitude and longitude\n",
    "        # to fit tide model convention\n",
    "        x,y = np.copy([LON,LAT]).astype(np.float64)\n",
    "        lt0, = np.nonzero(x < 0)\n",
    "        x[lt0] += 360.0\n",
    "    elif (model.format == 'FES'):\n",
    "        # if reading a FES netCDF solution\n",
    "        # adjust longitudinal convention of input latitude and longitude\n",
    "        # to fit tide model convention\n",
    "        x,y = np.copy([LON,LAT]).astype(np.float64)\n",
    "        lt0, = np.nonzero(x < 0)\n",
    "        x[lt0] += 360.0\n",
    "    # update plot coordinates\n",
    "    m.point.set_xdata(x)\n",
    "    m.point.set_ydata(y)\n",
    "    # refresh plot\n",
    "    IPython.display.display(m.figure)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "# check coordinates on tide grid\n",
    "m.figure,ax = plt.subplots(num=1, figsize=(8.25,5.25))\n",
    "ax.imshow(mz, interpolation='nearest',\n",
    "    extent=(xi.min(),xi.max(),yi.min(),yi.max()),\n",
    "    origin='lower', cmap='gray')\n",
    "m.point, = ax.plot([],[],'r*')\n",
    "update_coordinates(None)\n",
    "# no ticks on the x and y axes\n",
    "ax.get_xaxis().set_ticks([])\n",
    "ax.get_yaxis().set_ticks([])\n",
    "# stronger linewidth on frame\n",
    "[i.set_linewidth(2.0) for i in ax.spines.values()]\n",
    "# adjust subplot within figure\n",
    "m.figure.tight_layout()\n",
    "IPython.display.clear_output(wait=True)\n",
    "m.marker.observe(update_coordinates)"
   ]
  },
  {
   "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
}