{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Plot Antarctic Tidal Currents\n", "=============================\n", "\n", "Demonstrates plotting hourly tidal currents 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", "Finite Element Solution (FES) provided by AVISO \n", "- https://www.aviso.altimetry.fr/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", "- `calc_astrol_longitudes.py`: computes the basic astronomical mean longitudes \n", "- `calc_delta_time.py`: calculates difference between universal and dynamic time \n", "- `convert_ll_xy.py`: convert lat/lon points to and from projected coordinates \n", "- `load_constituent.py`: loads parameters for a given tidal constituent \n", "- `load_nodal_corrections.py`: load the nodal corrections for tidal constituents \n", "- `infer_minor_corrections.py`: return corrections for minor constituents \n", "- `read_tide_model.py`: extract tidal harmonic constants from OTIS tide models \n", "- `read_netcdf_model.py`: extract tidal harmonic constants from netcdf models \n", "- `read_FES_model.py`: extract tidal harmonic constants from FES tide models \n", "- `predict_tide.py`: predict tidal elevation at a single time using harmonic constants \n", "\n", "This notebook uses Jupyter widgets to set parameters for calculating the tidal maps. \n", "The widgets can be installed as described below. \n", "```\n", "pip3 install --user ipywidgets\n", "jupyter nbextension install --user --py widgetsnbextension\n", "jupyter nbextension enable --user --py widgetsnbextension\n", "jupyter-notebook\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Load modules" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "import pyproj\n", "import datetime\n", "import numpy as np\n", "import matplotlib\n", "matplotlib.rcParams['axes.linewidth'] = 2.0\n", "matplotlib.rcParams[\"animation.html\"] = \"jshtml\"\n", "import matplotlib.pyplot as plt\n", "import matplotlib.animation as animation\n", "import cartopy.crs as ccrs\n", "import ipywidgets as widgets\n", "from IPython.display import HTML\n", "\n", "import pyTMD.time\n", "from pyTMD.calc_delta_time import calc_delta_time\n", "from pyTMD.read_tide_model import extract_tidal_constants\n", "from pyTMD.read_netcdf_model import extract_netcdf_constants\n", "from pyTMD.read_GOT_model import extract_GOT_constants\n", "from pyTMD.read_FES_model import extract_FES_constants\n", "from pyTMD.infer_minor_corrections import infer_minor_corrections\n", "from pyTMD.predict_tide import predict_tide\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 \n", "- Date to run " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#-- set the directory with tide models\n", "dirText = widgets.Text(\n", " value=os.getcwd(),\n", " description='Directory:',\n", " disabled=False\n", ")\n", "\n", "#-- dropdown menu for setting tide model\n", "model_list = ['CATS0201','CATS2008','TPXO9-atlas','TPXO9-atlas-v2','TPXO9.1',\n", " 'TPXO8-atlas','TPXO7.2','FES2014']\n", "modelDropdown = widgets.Dropdown(\n", " options=model_list,\n", " value='CATS2008',\n", " description='Model:',\n", " disabled=False,\n", ")\n", "\n", "#-- date picker widget for setting time\n", "datepick = widgets.DatePicker(\n", " description='Date:',\n", " value = datetime.date.today(),\n", " disabled=False\n", ")\n", "\n", "#-- display widgets for setting directory, model and date\n", "widgets.VBox([dirText,modelDropdown,datepick])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Setup tide model parameters" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#-- directory with tide models\n", "tide_dir = os.path.expanduser(dirText.value)\n", "MODEL = modelDropdown.value\n", "if (MODEL == 'CATS0201'):\n", " grid_file = os.path.join(tide_dir,'cats0201_tmd','grid_CATS')\n", " model_file = os.path.join(tide_dir,'cats0201_tmd','UV0_CATS02_01')\n", " model_format = 'OTIS'\n", " EPSG = '4326'\n", " TYPES = ['u','v']\n", "elif (MODEL == 'CATS2008'):\n", " grid_file = os.path.join(tide_dir,'CATS2008','grid_CATS2008')\n", " model_file = os.path.join(tide_dir,'CATS2008','uv.CATS2008.out')\n", " model_format = 'OTIS'\n", " EPSG = 'CATS2008'\n", " TYPES = ['u','v']\n", "elif (MODEL == 'TPXO9-atlas'):\n", " model_directory = os.path.join(tide_dir,'TPXO9_atlas')\n", " grid_file = 'grid_tpxo9_atlas.nc.gz'\n", " model_files = {}\n", " model_files['u'] = ['u_q1_tpxo9_atlas_30.nc.gz','u_o1_tpxo9_atlas_30.nc.gz',\n", " 'u_p1_tpxo9_atlas_30.nc.gz','u_k1_tpxo9_atlas_30.nc.gz',\n", " 'u_n2_tpxo9_atlas_30.nc.gz','u_m2_tpxo9_atlas_30.nc.gz',\n", " 'u_s2_tpxo9_atlas_30.nc.gz','u_k2_tpxo9_atlas_30.nc.gz',\n", " 'u_m4_tpxo9_atlas_30.nc.gz','u_ms4_tpxo9_atlas_30.nc.gz',\n", " 'u_mn4_tpxo9_atlas_30.nc.gz','u_2n2_tpxo9_atlas_30.nc.gz']\n", " model_files['v'] = ['v_q1_tpxo9_atlas_30.nc.gz','v_o1_tpxo9_atlas_30.nc.gz',\n", " 'v_p1_tpxo9_atlas_30.nc.gz','v_k1_tpxo9_atlas_30.nc.gz',\n", " 'v_n2_tpxo9_atlas_30.nc.gz','v_m2_tpxo9_atlas_30.nc.gz',\n", " 'v_s2_tpxo9_atlas_30.nc.gz','v_k2_tpxo9_atlas_30.nc.gz',\n", " 'v_m4_tpxo9_atlas_30.nc.gz','v_ms4_tpxo9_atlas_30.nc.gz',\n", " 'v_mn4_tpxo9_atlas_30.nc.gz','v_2n2_tpxo9_atlas_30.nc.gz']\n", " model_format = 'netcdf'\n", " TYPES = ['u','v']\n", " SCALE = 1.0/100.0\n", " GZIP = True\n", "elif (MODEL == 'TPXO9-atlas-v2'):\n", " model_directory = os.path.join(tide_dir,'TPXO9_atlas_v2')\n", " grid_file = 'grid_tpxo9_atlas_30_v2.nc.gz'\n", " model_files = {}\n", " model_files['u'] = ['u_q1_tpxo9_atlas_30_v2.nc.gz','u_o1_tpxo9_atlas_30_v2.nc.gz',\n", " 'u_p1_tpxo9_atlas_30_v2.nc.gz','u_k1_tpxo9_atlas_30_v2.nc.gz',\n", " 'u_n2_tpxo9_atlas_30_v2.nc.gz','u_m2_tpxo9_atlas_30_v2.nc.gz',\n", " 'u_s2_tpxo9_atlas_30_v2.nc.gz','u_k2_tpxo9_atlas_30_v2.nc.gz',\n", " 'u_m4_tpxo9_atlas_30_v2.nc.gz','u_ms4_tpxo9_atlas_30_v2.nc.gz',\n", " 'u_mn4_tpxo9_atlas_30_v2.nc.gz','u_2n2_tpxo9_atlas_30_v2.nc.gz']\n", " model_files['v'] = ['v_q1_tpxo9_atlas_30_v2.nc.gz','v_o1_tpxo9_atlas_30_v2.nc.gz',\n", " 'v_p1_tpxo9_atlas_30_v2.nc.gz','v_k1_tpxo9_atlas_30_v2.nc.gz',\n", " 'v_n2_tpxo9_atlas_30_v2.nc.gz','v_m2_tpxo9_atlas_30_v2.nc.gz',\n", " 'v_s2_tpxo9_atlas_30_v2.nc.gz','v_k2_tpxo9_atlas_30_v2.nc.gz',\n", " 'v_m4_tpxo9_atlas_30_v2.nc.gz','v_ms4_tpxo9_atlas_30_v2.nc.gz',\n", " 'v_mn4_tpxo9_atlas_30_v2.nc.gz','v_2n2_tpxo9_atlas_30_v2.nc.gz']\n", " model_format = 'netcdf'\n", " TYPES = ['u','v']\n", " SCALE = 1.0/100.0\n", " GZIP = True\n", "elif (MODEL == 'TPXO9.1'):\n", " grid_file = os.path.join(tide_dir,'TPXO9.1','DATA','grid_tpxo9')\n", " model_file = os.path.join(tide_dir,'TPXO9.1','DATA','u_tpxo9.v1')\n", " model_format = 'OTIS'\n", " EPSG = '4326'\n", " TYPES = ['u','v']\n", "elif (MODEL == 'TPXO8-atlas'):\n", " grid_file = os.path.join(tide_dir,'tpxo8_atlas','grid_tpxo8atlas_30_v1')\n", " model_file = os.path.join(tide_dir,'tpxo8_atlas','uv.tpxo8_atlas_30_v1')\n", " model_format = 'ATLAS'\n", " EPSG = '4326'\n", " TYPES = ['u','v']\n", "elif (MODEL == 'TPXO7.2'):\n", " grid_file = os.path.join(tide_dir,'TPXO7.2_tmd','grid_tpxo7.2')\n", " model_file = os.path.join(tide_dir,'TPXO7.2_tmd','u_tpxo7.2')\n", " model_format = 'OTIS'\n", " EPSG = '4326'\n", " TYPES = ['u','v']\n", "elif (MODEL == 'FES2014'):\n", " model_directory = {}\n", " model_directory['u'] = os.path.join(tide_dir,'fes2014','eastward_velocity')\n", " model_directory['v'] = os.path.join(tide_dir,'fes2014','northward_velocity')\n", " model_files = ['2n2.nc.gz','eps2.nc.gz','j1.nc.gz','k1.nc.gz',\n", " 'k2.nc.gz','l2.nc.gz','la2.nc.gz','m2.nc.gz','m3.nc.gz','m4.nc.gz',\n", " 'm6.nc.gz','m8.nc.gz','mf.nc.gz','mks2.nc.gz','mm.nc.gz',\n", " 'mn4.nc.gz','ms4.nc.gz','msf.nc.gz','msqm.nc.gz','mtm.nc.gz',\n", " 'mu2.nc.gz','n2.nc.gz','n4.nc.gz','nu2.nc.gz','o1.nc.gz','p1.nc.gz',\n", " 'q1.nc.gz','r2.nc.gz','s1.nc.gz','s2.nc.gz','s4.nc.gz','sa.nc.gz',\n", " 'ssa.nc.gz','t2.nc.gz']\n", " c = ['2n2','eps2','j1','k1','k2','l2','lambda2','m2','m3','m4','m6',\n", " 'm8','mf','mks2','mm','mn4','ms4','msf','msqm','mtm','mu2','n2',\n", " 'n4','nu2','o1','p1','q1','r2','s1','s2','s4','sa','ssa','t2']\n", " model_format = 'FES'\n", " TYPES = ['u','v']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Setup coordinates for calculating tidal currents" ] }, { "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 = np.int((xlimits[1]-xlimits[0])/spacing[0])+1\n", "ny = np.int((ylimits[0]-ylimits[1])/spacing[1])+1\n", "#-- convert image coordinates from polar stereographic to latitude/longitude\n", "crs1 = pyproj.CRS.from_string(\"epsg:{0:d}\".format(3031))\n", "crs2 = pyproj.CRS.from_string(\"epsg:{0:d}\".format(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": [ "#### Calculate tide map" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#-- convert from calendar date to days relative to Jan 1, 1992 (48622 MJD)\n", "YMD = datepick.value\n", "tide_time = pyTMD.time.convert_calendar_dates(YMD.year, YMD.month,\n", " YMD.day, hour=np.arange(24))\n", "\n", "#-- save tide currents\n", "tide = {}\n", "#-- iterate over u and v currents\n", "for TYPE in TYPES:\n", " #-- read tidal constants and interpolate to grid points\n", " if model_format in ('OTIS','ATLAS'):\n", " amp,ph,D,c = extract_tidal_constants(lon, lat, grid_file, model_file,\n", " EPSG, TYPE=TYPE, METHOD='spline', GRID=model_format)\n", " DELTAT = np.zeros_like(tide_time)\n", " elif (model_format == 'netcdf'):\n", " amp,ph,D,c = extract_netcdf_constants(lon, lat, model_directory,\n", " grid_file, model_files[TYPE], TYPE=TYPE, METHOD='spline',\n", " SCALE=SCALE, GZIP=GZIP)\n", " DELTAT = np.zeros_like(tide_time)\n", " elif (model_format == 'GOT'):\n", " amp,ph = extract_GOT_constants(lon, lat, model_directory, model_files,\n", " METHOD='spline', SCALE=SCALE)\n", " #-- interpolate delta times from calendar dates to tide time\n", " delta_file = pyTMD.utilities.get_data_path(['data','merged_deltat.data'])\n", " DELTAT = calc_delta_time(delta_file, tide_time)\n", " elif (model_format == 'FES'):\n", " amp,ph = extract_FES_constants(lon, lat, model_directory[TYPE], model_files,\n", " TYPE=TYPE, VERSION=MODEL, METHOD=METHOD, SCALE=SCALE)\n", " #-- interpolate delta times from calendar dates to tide time\n", " delta_file = pyTMD.utilities.get_data_path(['data','merged_deltat.data'])\n", " DELTAT = calc_delta_time(delta_file, tide_time)\n", " \n", " #-- calculate complex phase in radians for Euler's\n", " cph = -1j*ph*np.pi/180.0\n", " #-- calculate constituent oscillation\n", " hc = amp*np.exp(cph)\n", "\n", " #-- allocate for tide current map calculated every hour\n", " tide[TYPE] = np.ma.zeros((ny,nx,24))\n", " for hour in range(24):\n", " #-- predict tidal elevations at time and infer minor corrections\n", " TIDE = predict_tide(tide_time[hour], hc, c, DELTAT=DELTAT[hour],\n", " CORRECTIONS=model_format)\n", " MINOR = infer_minor_corrections(tide_time[hour], hc, c,\n", " DELTAT=DELTAT[hour], CORRECTIONS=model_format)\n", " #-- add major and minor components and reform grid\n", " tide[TYPE][:,:,hour] = np.reshape((TIDE+MINOR),(ny,nx))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Create animation of hourly tidal oscillation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#-- output Antarctic Tidal Current Animation\n", "projection = ccrs.Stereographic(central_longitude=0.0,\n", " central_latitude=-90.0,true_scale_latitude=-71.0)\n", "#-- figure axis and image objects\n", "ax1,im = ({},{})\n", "fig, (ax1['u'],ax1['v']) = plt.subplots(num=1, ncols=2,\n", " figsize=(11.5,7), subplot_kw=dict(projection=projection))\n", "vmin = np.min([tide['u'].min(),tide['v'].min()])\n", "vmax = np.max([tide['u'].max(),tide['v'].max()])\n", "extent = (xlimits[0],xlimits[1],ylimits[0],ylimits[1])\n", "for TYPE,ax in ax1.items():\n", " #-- plot tidal currents\n", " im[TYPE] = ax.imshow(np.zeros((ny,nx)),\n", " interpolation='nearest', vmin=vmin, vmax=vmax,\n", " transform=projection, extent=extent, origin='upper',\n", " animated=True)\n", " #-- add high resolution cartopy coastlines\n", " ax.coastlines('10m')\n", " #-- set x and y limits\n", " ax.set_xlim(xlimits)\n", " ax.set_ylim(ylimits)\n", " # stronger linewidth on frame\n", " ax.spines['geo'].set_linewidth(2.0)\n", " ax.spines['geo'].set_capstyle('projecting')\n", "\n", "#-- Add colorbar with a colorbar axis\n", "#-- Add an axes at position rect [left, bottom, width, height]\n", "cbar_ax = fig.add_axes([0.085, 0.075, 0.83, 0.035])\n", "#-- extend = add extension triangles to upper and lower bounds\n", "#-- options: neither, both, min, max\n", "cbar = fig.colorbar(im['u'], cax=cbar_ax, extend='both',\n", " extendfrac=0.0375, drawedges=False, orientation='horizontal')\n", "#-- rasterized colorbar to remove lines\n", "cbar.solids.set_rasterized(True)\n", "#-- Add label to the colorbar\n", "cbar.ax.set_xlabel('{0} Tidal Velocity'.format(MODEL), fontsize=13)\n", "cbar.ax.set_ylabel('cm/s', fontsize=13, rotation=0)\n", "cbar.ax.yaxis.set_label_coords(1.04, 0.15)\n", "#-- ticks lines all the way across\n", "cbar.ax.tick_params(which='both', width=1, length=18,\n", " labelsize=13, direction='in')\n", "\n", "#-- add title (date and time)\n", "ttl = fig.suptitle(None, fontsize=13)\n", "# adjust subplot within figure\n", "fig.subplots_adjust(left=0.02,right=0.98,bottom=0.1,top=0.98,wspace=0.04)\n", " \n", "#-- animate each map\n", "def animate_maps(hour):\n", " #-- set map data iterating over u and v currents\n", " for TYPE in TYPES:\n", " im[TYPE].set_data(tide[TYPE][:,:,hour])\n", " #-- set title\n", " args = (YMD.year,YMD.month,YMD.day,hour)\n", " ttl.set_text('{0:4d}-{1:02d}-{2:02d}T{3:02d}:00:00'.format(*args))\n", "\n", "#-- set animation\n", "anim = animation.FuncAnimation(fig, animate_maps, frames=24)\n", "%matplotlib inline\n", "HTML(anim.to_jshtml())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "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.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }