diff --git a/doc/source/_assets/tide_form_factor.png b/doc/source/_assets/tide_form_factor.png new file mode 100644 index 00000000..3aa18a5e Binary files /dev/null and b/doc/source/_assets/tide_form_factor.png differ diff --git a/doc/source/notebooks/Plot-Tide-Form-Factor.ipynb b/doc/source/notebooks/Plot-Tide-Form-Factor.ipynb new file mode 100644 index 00000000..4f784280 --- /dev/null +++ b/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 `) 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 +} diff --git a/doc/source/user_guide/Examples.rst b/doc/source/user_guide/Examples.rst index ea5d8601..75d00eae 100644 --- a/doc/source/user_guide/Examples.rst +++ b/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 diff --git a/pyTMD/arguments.py b/pyTMD/arguments.py index 52843f21..93c3d5cf 100755 --- a/pyTMD/arguments.py +++ b/pyTMD/arguments.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" arguments.py -Written by Tyler Sutterley (10/2024) +Written by Tyler Sutterley (11/2024) Calculates the nodal corrections for tidal constituents Modification of ARGUMENTS fortran subroutine by Richard Ray 03/1999 @@ -38,6 +38,7 @@ Ocean Tides", Journal of Atmospheric and Oceanic Technology, (2002). UPDATE HISTORY: + Updated 11/2024: allow variable case for Doodson number formalisms Updated 10/2024: can convert Doodson numbers formatted as strings update Doodson number conversions to follow Cartwright X=10 convention add function to parse Cartwright/Tayler/Edden tables @@ -478,13 +479,13 @@ def doodson_number( else: return None # extract identifier in formalism - if (kwargs['formalism'] == 'Cartwright'): + if (kwargs['formalism'].title() == 'Cartwright'): # extract Cartwright number numbers = np.array(coefficients[:6,0]) - elif (kwargs['formalism'] == 'Doodson'): + elif (kwargs['formalism'].title() == 'Doodson'): # convert from coefficients to Doodson number numbers = _to_doodson_number(coefficients[:,0], **kwargs) - elif (kwargs['formalism'] == 'Extended'): + elif (kwargs['formalism'].title() == 'Extended'): # convert to extended Doodson number in UKHO format numbers = _to_extended_doodson(coefficients[:,0], **kwargs) else: @@ -502,13 +503,13 @@ def doodson_number( numbers[c] = None continue # convert from coefficients to Doodson number - if (kwargs['formalism'] == 'Cartwright'): + if (kwargs['formalism'].title() == 'Cartwright'): # extract Cartwright number numbers[c] = np.array(coefficients[:6,0]) - elif (kwargs['formalism'] == 'Doodson'): + elif (kwargs['formalism'].title() == 'Doodson'): # convert from coefficients to Doodson number numbers[c] = _to_doodson_number(coefficients[:,0], **kwargs) - elif (kwargs['formalism'] == 'Extended'): + elif (kwargs['formalism'].title() == 'Extended'): # convert to extended Doodson number in UKHO format numbers[c] = _to_extended_doodson(coefficients[:,0], **kwargs) # return the Doodson or Cartwright number diff --git a/pyTMD/io/constituents.py b/pyTMD/io/constituents.py index 6102ce7b..517269c8 100644 --- a/pyTMD/io/constituents.py +++ b/pyTMD/io/constituents.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" constituents.py -Written by Tyler Sutterley (10/2024) +Written by Tyler Sutterley (11/2024) Basic tide model constituent class PYTHON DEPENDENCIES: @@ -10,6 +10,7 @@ https://numpy.org/doc/stable/user/numpy-for-matlab-users.html UPDATE HISTORY: + Updated 11/2024: added property for Extended Doodson numbers Updated 10/2024: added property for the shape of constituent fields Updated 09/2024: add more known constituents to string parser function Updated 08/2024: add GOT prime nomenclature for 3rd degree constituents @@ -166,7 +167,7 @@ def phase(self, field: str): @property def doodson_number(self): - """constituent Doodson number + """Doodson number for constituents """ doodson_numbers = [] # for each constituent ID @@ -183,7 +184,7 @@ def doodson_number(self): @property def cartwright_number(self): - """constituent Cartwright numbers + """Cartwright numbers for constituents """ cartwright_numbers = [] # for each constituent ID @@ -198,6 +199,23 @@ def cartwright_number(self): # return the list of Cartwright numbers return cartwright_numbers + @property + def extended_doodson(self): + """Extended Doodson numbers for constituents + """ + extended_numbers = [] + # for each constituent ID + for f in self.fields: + try: + # try to get the Extended Doodson number + XDO = pyTMD.arguments.doodson_number(f, formalism='Extended') + except (AssertionError, ValueError) as exc: + XDO = None + # add Extended Doodson number to the combined list + extended_numbers.append(XDO) + # return the list of Extended Doodson numbers + return extended_numbers + @property def shape(self): """Shape of constituent fields diff --git a/pyTMD/io/model.py b/pyTMD/io/model.py index 14cdcd9e..dace713b 100644 --- a/pyTMD/io/model.py +++ b/pyTMD/io/model.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" model.py -Written by Tyler Sutterley (10/2024) +Written by Tyler Sutterley (11/2024) Retrieves tide model parameters for named tide models and from model definition files @@ -11,6 +11,7 @@ https://numpy.org/doc/stable/user/numpy-for-matlab-users.html UPDATE HISTORY: + Updated 11/2024: use Love numbers for long-period tides in node equilibrium Updated 10/2024: add wrapper functions to read and interpolate constants add functions to append node tide equilibrium values to amplitudes remove redundant default keyword arguments to readers and interpolators @@ -421,7 +422,7 @@ def long_name(self) -> str: elif (self.type == 'z') and (self.variable == 'tide_load'): return 'load_tide_elevation' elif (self.type == 'z') and (self.variable == 'tide_lpe'): - return 'Equilibrium_Tide' + return 'equilibrium_tide_elevation' elif (self.type == ['u','v']): return dict(u='zonal_tidal_current', v='meridional_tidal_current') else: @@ -1203,9 +1204,10 @@ def node_equilibrium(lat: np.ndarray): """ # Cartwright and Edden potential amplitude amajor = 0.027929# node - # love numbers - k2 = 0.302 - h2 = 0.609 + # Love numbers for long-period tides (Wahr, 1981) + k2 = 0.299 + h2 = 0.606 + # tilt factor: response with respect to the solid earth gamma_2 = (1.0 + k2 - h2) # 2nd degree Legendre polynomials P20 = 0.5*(3.0*np.sin(lat*np.pi/180.0)**2 - 1.0) diff --git a/pyTMD/predict.py b/pyTMD/predict.py index 14212edd..abde477b 100644 --- a/pyTMD/predict.py +++ b/pyTMD/predict.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" predict.py -Written by Tyler Sutterley (10/2024) +Written by Tyler Sutterley (11/2024) Prediction routines for ocean, load, equilibrium and solid earth tides REFERENCES: @@ -20,6 +20,7 @@ spatial.py: utilities for working with geospatial data UPDATE HISTORY: + Updated 11/2024: use Love numbers for long-period tides when inferring Updated 10/2024: use PREM as the default Earth model for Love numbers more descriptive error message if cannot infer minor constituents updated calculation of long-period equilibrium tides @@ -771,6 +772,7 @@ def _infer_diurnal( j1, = j # Love numbers of degree 2 for constituent h2, k2, l2 = _body_tide_love_numbers(omajor[i]) + # tilt factor: response with respect to the solid earth gamma_2 = (1.0 + k2 - h2) # "normalize" tide values z[:,i] = zmajor[:,j1]/(amajor[i]*gamma_2) @@ -833,6 +835,7 @@ def _infer_diurnal( for k in minor_indices: # Love numbers of degree 2 for constituent h2, k2, l2 = _body_tide_love_numbers(omega[k]) + # tilt factor: response with respect to the solid earth gamma_2 = (1.0 + k2 - h2) # interpolate from major constituents if (kwargs['method'].lower() == 'linear'): @@ -1228,12 +1231,15 @@ def equilibrium_tide( G = np.dot(fargs, coef) Z = np.inner(np.cos(G*np.pi/180.0), CTE) - # Multiply by gamma_2 * normalization * P20(lat) - k2 = 0.302 - h2 = 0.609 + # Love numbers for long-period tides (Wahr, 1981) + k2 = 0.299 + h2 = 0.606 + # tilt factor: response with respect to the solid earth gamma_2 = (1.0 + k2 - h2) + # 2nd degree Legendre polynomials P20 = 0.5*(3.0*np.sin(lat*np.pi/180.0)**2 - 1.0) # calculate long-period equilibrium tide and convert to meters + # Multiply by gamma_2 * normalization * P20(lat) if (nlat != nt): lpet = gamma_2*np.sqrt((4.0+1.0)/(4.0*np.pi))*np.outer(P20,Z/100.0) else: