fix: allow variable case for Doodson number formalisms

doc/source/notebooks/Plot-Tide-Form-Factor.ipynb

@@ -0,0 +1,340 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plot Tide Form Factor\n",
+    "======================\n",
+    "\n",
+    "This ({nb-download}`notebook <Plot-Tide-Form-Factor.ipynb>`) demonstrates plotting tidal form factors for classifying tides\n",
+    "\n",
+    "- Daily tidal form factors for determining the dominant species of a region using the classifications from [Courtier (1938)](https://journals.lib.unb.ca/index.php/ihr/article/download/27428/1882520184). The dominant species classifications do have limitations as pointed out by [Amin (1986)](https://journals.lib.unb.ca/index.php/ihr/article/download/23443/27218/0)\n",
+    "- Monthly tidal form factors for semi-diurnal species from [Byun and Hart](https://doi.org/10.5194/os-16-965-2020)\n",
+    "\n",
+    "OTIS format tidal solutions provided by Oregon State University and ESR  \n",
+    "- [http://volkov.oce.orst.edu/tides/region.html](http://volkov.oce.orst.edu/tides/region.html) \n",
+    "- [https://www.esr.org/research/polar-tide-models/list-of-polar-tide-models/](https://www.esr.org/research/polar-tide-models/list-of-polar-tide-models/)\n",
+    "- [ftp://ftp.esr.org/pub/datasets/tmd/](ftp://ftp.esr.org/pub/datasets/tmd/)  \n",
+    "\n",
+    "Global Tide Model (GOT) solutions provided by Richard Ray at GSFC  \n",
+    "- [https://earth.gsfc.nasa.gov/geo/data/ocean-tide-models](https://earth.gsfc.nasa.gov/geo/data/ocean-tide-models)\n",
+    "\n",
+    "Finite Element Solution (FES) provided by AVISO  \n",
+    "- [https://www.aviso.altimetry.fr/en/data/products/auxiliary-products/global-tide-fes.html](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",
+    "- `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.GOT.py`: extract tidal harmonic constants from GOT tide 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 numpy as np\n",
+    "import matplotlib\n",
+    "matplotlib.rcParams['axes.linewidth'] = 2.0\n",
+    "import matplotlib.pyplot as plt\n",
+    "import matplotlib.colors as colors\n",
+    "import cartopy.crs as ccrs\n",
+    "import ipywidgets\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.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",
+    "TMDwidgets.VBox([\n",
+    "    TMDwidgets.directory,\n",
+    "    TMDwidgets.model,\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",
+    "    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 a global image\n",
+    "xlimits = [-180,180]\n",
+    "ylimits = [-90, 90]\n",
+    "spacing = [0.25, 0.25]\n",
+    "# x and y coordinates\n",
+    "x = np.arange(xlimits[0],xlimits[1]+spacing[0],spacing[0])\n",
+    "y = np.arange(ylimits[0],ylimits[1]+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[1]-ylimits[0])/spacing[1])+1\n",
+    "# flatten latitude and longitude to arrays\n",
+    "lon,lat =  xgrid.flatten(), ygrid.flatten()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Calculate tidal amplitudes and phases"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# read tidal constants and interpolate to grid points\n",
+    "if model.format in ('OTIS','ATLAS-compact','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, crop=True,\n",
+    "        method='spline', grid=model.file_format)\n",
+    "elif (model.format == 'ATLAS-netcdf'):\n",
+    "    amp,ph,D,c = pyTMD.io.ATLAS.extract_constants(lon, lat, model.grid_file,\n",
+    "        model.model_file, type=model.type, crop=True, method='spline',\n",
+    "        scale=model.scale, compressed=model.compressed)\n",
+    "elif model.format in ('GOT-ascii', 'GOT-netcdf'):\n",
+    "    amp,ph,c = pyTMD.io.GOT.extract_constants(lon, lat, model.model_file,\n",
+    "        grid=model.file_format, crop=True, method='spline',\n",
+    "        scale=model.scale, compressed=model.compressed)\n",
+    "elif (model.format == 'FES-netcdf'):\n",
+    "    amp,ph = pyTMD.io.FES.extract_constants(lon, lat, model.model_file,\n",
+    "        type=model.type, version=model.version, crop=True,\n",
+    "        method='spline', scale=model.scale, compressed=model.compressed)\n",
+    "    c = model.constituents"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Calculate tidal form factors\n",
+    "\n",
+    "Courtier form factor:\n",
+    "Ratios between major diurnal tides and major semi-diurnal tides\n",
+    "\n",
+    "- F: < 0.25: Semi-diurnal\n",
+    "- F: 0.25 - 1.5: Mixed predominantly semi-diurnal\n",
+    "- F: 1.5 - 3.0: Mixed predominantly diurnal\n",
+    "- F: > 3.0: Diurnal\n",
+    "\n",
+    "Byut-Hart form factor:\n",
+    "Ratios between semi-diurnal tides for monthly tidal envelopes\n",
+    "\n",
+    "- E: < 0.8: Spring-Neap\n",
+    "- E: 0.8 - 1.0: Mixed predominantly Spring-Neap\n",
+    "- E: 1.0 - 1.15: Mixed predominantly Perigean-Apogean\n",
+    "- E: > 2.0: Perigean-Apogean\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "TMDwidgets.form_factor = ipywidgets.Dropdown(\n",
+    "    options=['Courtier','Byun-Hart'],\n",
+    "    value='Courtier',\n",
+    "    description='Factor:',\n",
+    "    disabled=False,\n",
+    "    style=TMDwidgets.style,\n",
+    ")\n",
+    "display(TMDwidgets.form_factor)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# find constituents for tidal form factors\n",
+    "k1 = c.index('k1')\n",
+    "o1 = c.index('o1')\n",
+    "m2 = c.index('m2')\n",
+    "s2 = c.index('s2')\n",
+    "n2 = c.index('n2')\n",
+    "# select form factor\n",
+    "if TMDwidgets.form_factor.value == 'Courtier':\n",
+    "    # tidal form factor from Courtier\n",
+    "    factor = np.reshape((amp[:,k1] + amp[:,o1])/(amp[:,m2] + amp[:,s2]), (ny,nx))\n",
+    "    boundary = np.array([0.0, 0.25, 1.5, 3.0, 5.0])\n",
+    "    ticklabels = ['Semi-Diurnal', 'Mixed SD', 'Mixed D', 'Diurnal']\n",
+    "    longname = 'Tide Species Classification'\n",
+    "elif TMDwidgets.form_factor.value == 'Byun-Hart':\n",
+    "    # semi-diurnal form factor from Byun and Hart\n",
+    "    factor = np.reshape((amp[:,m2] + amp[:,n2])/(amp[:,m2] + amp[:,s2]), (ny,nx))\n",
+    "    boundary = np.array([0.0, 0.8, 1.0, 1.15, 2.0])\n",
+    "    ticklabels = ['Spring-Neap', 'Mixed S-N', 'Mixed P-A', 'Perigean-Apogean']\n",
+    "    longname = 'Semi-Diurnal Classification'\n",
+    "# calculate ticks for labels\n",
+    "ticks = 0.5*(boundary[1:] + boundary[:-1])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Create plot of tidal form factors"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# cartopy transform for Equirectangular Projection\n",
+    "projection = ccrs.PlateCarree()\n",
+    "# create figure axis\n",
+    "fig, ax = plt.subplots(num=1, figsize=(5.5,3.5),\n",
+    "    subplot_kw=dict(projection=projection))\n",
+    "# create boundary norm\n",
+    "norm = colors.BoundaryNorm(boundary, ncolors=256)\n",
+    "# plot tidal form factor\n",
+    "extent = (xlimits[0],xlimits[1],ylimits[0],ylimits[1])\n",
+    "im = ax.imshow(factor, interpolation='nearest',\n",
+    "    norm=norm, cmap='plasma', transform=projection,\n",
+    "    extent=extent, origin='lower')\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, extend='neither',\n",
+    "    extendfrac=0.0375, orientation='horizontal', pad=0.025,\n",
+    "    shrink=0.90, 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_title(longname, fontsize=13,\n",
+    "    rotation=0, y=-2.0, va='top')\n",
+    "# Set the tick levels for the colorbar\n",
+    "cbar.set_ticks(ticks=ticks, labels=ticklabels)\n",
+    "\n",
+    "# axis = equal\n",
+    "ax.set_aspect('equal', adjustable='box')\n",
+    "# set x and y limits\n",
+    "ax.set_xlim(xlimits)\n",
+    "ax.set_ylim(ylimits)\n",
+    "\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",
+    "ax.spines['geo'].set_linewidth(2.0)\n",
+    "ax.spines['geo'].set_capstyle('projecting')\n",
+    "\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": "base",
+   "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.14"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

doc/source/user_guide/Examples.rst

@@ -21,6 +21,7 @@ The available examples include:
    ../notebooks/Plot-Ocean-Pole-Tide-Map
    ../notebooks/Plot-Ross-Ice-Shelf-Map
    ../notebooks/Plot-Tide-Forecasts
+   ../notebooks/Plot-Tide-Form-Factor
    ../notebooks/Solve-Synthetic-Tides
 
 .. grid:: 1 2 4 4
@@ -83,6 +84,11 @@ The available examples include:
       :img-top: ../_assets/tide_forecasts.png
       :link: ../notebooks/Plot-Tide-Forecasts.html
 
+    .. grid-item-card::  Plot Tide Form Factor
+      :text-align: center
+      :img-top: ../_assets/tide_form_factor.png
+      :link: ../notebooks/Plot-Tide-Form-Factor.html
+
     .. grid-item-card::  Solve Synthetic Tides
       :text-align: center
       :img-top: ../_assets/solve_synthetic_tides.png