notebooks/Plot Ross Ice Shelf Map.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot Ross Ice Shelf Map\n",
    "=======================\n",
    "\n",
    "Demonstrates plotting hourly tidal displacements for the Ross Ice Shelf\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/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_GOT_model.py`: extract tidal harmonic constants from GSFC GOT 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',\n",
    "    'TPXO9-atlas-v3','TPXO9.1','TPXO8-atlas','TPXO7.2',\n",
    "    'GOT4.7','GOT4.8','GOT4.10','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','h0_CATS02_01')\n",
    "    model_format = 'OTIS'\n",
    "    EPSG = '4326'\n",
    "    TYPE = 'z'\n",
    "elif (MODEL == 'CATS2008'):\n",
    "    grid_file = os.path.join(tide_dir,'CATS2008','grid_CATS2008')\n",
    "    model_file = os.path.join(tide_dir,'CATS2008','hf.CATS2008.out')\n",
    "    model_format = 'OTIS'\n",
    "    EPSG = 'CATS2008'\n",
    "    TYPE = 'z'\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 = ['h_q1_tpxo9_atlas_30.nc.gz','h_o1_tpxo9_atlas_30.nc.gz',\n",
    "        'h_p1_tpxo9_atlas_30.nc.gz','h_k1_tpxo9_atlas_30.nc.gz',\n",
    "        'h_n2_tpxo9_atlas_30.nc.gz','h_m2_tpxo9_atlas_30.nc.gz',\n",
    "        'h_s2_tpxo9_atlas_30.nc.gz','h_k2_tpxo9_atlas_30.nc.gz',\n",
    "        'h_m4_tpxo9_atlas_30.nc.gz','h_ms4_tpxo9_atlas_30.nc.gz',\n",
    "        'h_mn4_tpxo9_atlas_30.nc.gz','h_2n2_tpxo9_atlas_30.nc.gz']\n",
    "    model_format = 'netcdf'\n",
    "    TYPE = 'z'\n",
    "    SCALE = 1.0/1000.0\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 = ['h_q1_tpxo9_atlas_30_v2.nc.gz','h_o1_tpxo9_atlas_30_v2.nc.gz',\n",
    "        'h_p1_tpxo9_atlas_30_v2.nc.gz','h_k1_tpxo9_atlas_30_v2.nc.gz',\n",
    "        'h_n2_tpxo9_atlas_30_v2.nc.gz','h_m2_tpxo9_atlas_30_v2.nc.gz',\n",
    "        'h_s2_tpxo9_atlas_30_v2.nc.gz','h_k2_tpxo9_atlas_30_v2.nc.gz',\n",
    "        'h_m4_tpxo9_atlas_30_v2.nc.gz','h_ms4_tpxo9_atlas_30_v2.nc.gz',\n",
    "        'h_mn4_tpxo9_atlas_30_v2.nc.gz','h_2n2_tpxo9_atlas_30_v2.nc.gz']\n",
    "    model_format = 'netcdf'\n",
    "    TYPE = 'z'\n",
    "    SCALE = 1.0/1000.0\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','h_tpxo9.v1')\n",
    "    model_format = 'OTIS'\n",
    "    EPSG = '4326'\n",
    "    TYPE = 'z'\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','hf.tpxo8_atlas_30_v1')\n",
    "    model_format = 'ATLAS'\n",
    "    EPSG = '4326'\n",
    "    TYPE = 'z'\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','h_tpxo7.2')\n",
    "    model_format = 'OTIS'\n",
    "    EPSG = '4326'\n",
    "    TYPE = 'z'\n",
    "elif (MODEL == 'GOT4.7'):\n",
    "    model_directory = os.path.join(tide_dir,'GOT4.7','grids_oceantide')\n",
    "    model_files = ['q1.d.gz','o1.d.gz','p1.d.gz','k1.d.gz','n2.d.gz',\n",
    "        'm2.d.gz','s2.d.gz','k2.d.gz','s1.d.gz','m4.d.gz']\n",
    "    c = ['q1','o1','p1','k1','n2','m2','s2','k2','s1','m4']\n",
    "    model_format = 'GOT'\n",
    "    SCALE = 1.0/100.0\n",
    "elif (MODEL == 'GOT4.8'):\n",
    "    model_directory = os.path.join(tide_dir,'got4.8','grids_oceantide')\n",
    "    model_files = ['q1.d.gz','o1.d.gz','p1.d.gz','k1.d.gz','n2.d.gz',\n",
    "        'm2.d.gz','s2.d.gz','k2.d.gz','s1.d.gz','m4.d.gz']\n",
    "    c = ['q1','o1','p1','k1','n2','m2','s2','k2','s1','m4']\n",
    "    model_format = 'GOT'\n",
    "    SCALE = 1.0/100.0\n",
    "elif (MODEL == 'GOT4.10'):\n",
    "    model_directory = os.path.join(tide_dir,'GOT4.10c','grids_oceantide')\n",
    "    model_files = ['q1.d.gz','o1.d.gz','p1.d.gz','k1.d.gz','n2.d.gz',\n",
    "        'm2.d.gz','s2.d.gz','k2.d.gz','s1.d.gz','m4.d.gz']\n",
    "    c = ['q1','o1','p1','k1','n2','m2','s2','k2','s1','m4']\n",
    "    model_format = 'GOT'\n",
    "    SCALE = 1.0/100.0\n",
    "elif (MODEL == 'FES2014'):\n",
    "    model_directory = os.path.join(tide_dir,'fes2014','ocean_tide')\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",
    "    TYPE = 'z'\n",
    "    SCALE = 1.0/100.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Setup coordinates for calculating tides"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#-- create an image around the Ross Ice Shelf\n",
    "xlimits = [-740000,520000]\n",
    "ylimits = [-1430000,-300000]\n",
    "spacing = [5e3,-5e3]\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",
    "#-- 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, METHOD='spline', SCALE=SCALE)\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, model_files,\n",
    "        TYPE=TYPE, VERSION=MODEL, 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",
    "    \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 map calculated every hour\n",
    "tide_cm = 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",
    "    #-- convert from meters to centimeters\n",
    "    tide_cm[:,:,hour] = 100.0*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 Ross Ice Shelf Tide Animation\n",
    "projection = ccrs.Stereographic(central_longitude=0.0,\n",
    "    central_latitude=-90.0,true_scale_latitude=-71.0)\n",
    "fig, ax = plt.subplots(num=1, figsize=(9,8),\n",
    "    subplot_kw=dict(projection=projection))\n",
    "#-- plot tide height\n",
    "vmin,vmax = (np.min(tide_cm), np.max(tide_cm))\n",
    "extent = (xlimits[0],xlimits[1],ylimits[0],ylimits[1])\n",
    "im = ax.imshow(np.zeros((ny,nx)), interpolation='nearest',\n",
    "    vmin=vmin, vmax=vmax, transform=projection,\n",
    "    extent=extent, origin='upper', animated=True)\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, pad=0.025, extend='both',\n",
    "    extendfrac=0.0375, shrink=0.85, aspect=22.5, drawedges=False)\n",
    "#-- rasterized colorbar to remove lines\n",
    "cbar.solids.set_rasterized(True)\n",
    "#-- Add label to the colorbar\n",
    "cbar.ax.set_ylabel('{0} Tide Height'.format(MODEL), fontsize=13)\n",
    "cbar.ax.set_xlabel('cm', fontsize=13)\n",
    "cbar.ax.xaxis.set_label_coords(0.50, 1.04)\n",
    "#-- ticks lines all the way across\n",
    "cbar.ax.tick_params(which='both', width=1, length=18,\n",
    "    labelsize=13, direction='in')\n",
    "#-- add title (date and time)\n",
    "ttl = ax.set_title(None, fontsize=13)\n",
    "#-- set x and y limits\n",
    "ax.set_xlim(xlimits)\n",
    "ax.set_ylim(ylimits)\n",
    "\n",
    "# stronger linewidth on frame\n",
    "ax.spines['geo'].set_linewidth(2.0)\n",
    "ax.spines['geo'].set_capstyle('projecting')\n",
    "# adjust subplot within figure\n",
    "fig.subplots_adjust(left=0.02,right=0.98,bottom=0.05,top=0.98)\n",
    "           \n",
    "#-- animate each map\n",
    "def animate_maps(hour):\n",
    "    #-- set map data\n",
    "    im.set_data(tide_cm[:,:,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
}