AF_compare_GIF.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import gsw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.61019015, 0.11271557, -0.36168176, -0.11810858, 0.056160484, -0.0016579428)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sd = 9 # Sill depth\n",
    "file1 = xr.open_dataset('/Volumes/Extreme SSD/MITgcm outputs/strat_qf700/state_' + str(format(sd,'03d')) + '.nc')\n",
    "state1 = file1.isel(X=range(120), Xp1=range(121), Y=range(30,39), Yp1=range(30,40)) # Remove the boundary walls and the shelf region\n",
    "file2 = xr.open_dataset('/Volumes/Extreme SSD/MITgcm outputs/strat_tide003/state_' + str(format(sd,'03d')) + '.nc')\n",
    "state2 = file2.isel(X=range(120), Xp1=range(121), Y=range(30,39), Yp1=range(30,40)) # Remove the boundary walls and the shelf region\n",
    "\n",
    "time = state1.T.data    \n",
    "x = state1.X.data / 1000 # x coornidate in km\n",
    "z = state1.Z.data\n",
    "pres = gsw.p_from_z(z, 0)\n",
    "p = np.tile(pres,[len(x),1]).T\n",
    "\n",
    "\n",
    "# NoTide data\n",
    "pt1 = state1.Temp.data.mean(2) # Dimesions T Z X\n",
    "s1 = state1.S.data.mean(2)\n",
    "rho1 = gsw.rho(s1, pt1, p)\n",
    "u1_0 = 0.5 * (state1.U.data[:, :, :, 1:] + state1.U.data[:, :, :, :-1])\n",
    "u1 = u1_0.mean(2)\n",
    "w1 = state1.W.data.mean(2)\n",
    "mask1 = (pt1 == 0)\n",
    "pt1[mask1] = np.nan\n",
    "s1[mask1] = np.nan\n",
    "u1[mask1] = np.nan\n",
    "w1[mask1] = np.nan\n",
    "rho1[mask1] = np.nan\n",
    "\n",
    "\n",
    "# qf-tide data\n",
    "pt2 = state2.Temp.data.mean(2) # Dimesions T Z X\n",
    "s2 = state2.S.data.mean(2)\n",
    "rho2 = gsw.rho(s2, pt2, p)\n",
    "u2_0 = 0.5 * (state2.U.data[:, :, :, 1:] + state2.U.data[:, :, :, :-1])\n",
    "u2 = u2_0.mean(2)\n",
    "w2 = state2.W.data.mean(2)\n",
    "mask2 = (pt2 == 0)\n",
    "pt2[mask2] = np.nan\n",
    "s2[mask2] = np.nan\n",
    "u2[mask2] = np.nan\n",
    "w2[mask2] = np.nan\n",
    "rho2[mask2] = np.nan\n",
    "\n",
    "np.nanmax(u1), np.nanmax(u2), np.nanmin(u1), np.nanmin(u2), np.nanmean(u1), np.nanmean(u2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/weiyangbao/miniconda3/envs/py37/lib/python3.7/site-packages/ipykernel_launcher.py:50: UserWarning: No contour levels were found within the data range.\n"
     ]
    }
   ],
   "source": [
    "fig_path = \"/Users/weiyangbao/Documents/Results_Figs/Af_compare/Sill009/\"\n",
    "if not os.path.exists(fig_path):\n",
    "    os.makedirs(fig_path)\n",
    "\n",
    "\n",
    "    \n",
    "ct = np.linspace(6,9,10) # Contour range for temperature\n",
    "crho = np.linspace(1016,1020,10)\n",
    "cmap = 'RdYlBu_r'\n",
    "\n",
    "for t in range(120,145):   \n",
    "    fig = plt.figure(figsize=(16,8))\n",
    "    fig.tight_layout()\n",
    "\n",
    "    # NoTide Plot\n",
    "    ax1 = fig.add_subplot(211)\n",
    "    im1 = ax1.contourf(x,z[:24],pt1[t,:24,:], ct, cmap=cmap)\n",
    "    ax1.xaxis.set_visible(False)\n",
    "    ax1.set_xlim([x[1],x[-1]])\n",
    "    ax1.set_facecolor('gray')\n",
    "\n",
    "    cbaxes = fig.add_axes([0.74, 0.92, 0.15, 0.01])\n",
    "    cbar = fig.colorbar(im1,cax=cbaxes, orientation='horizontal')\n",
    "    cbar.set_ticks([6, 7, 8, 9])\n",
    "    cbar.ax.set_title('Potential Temperature (\\N{DEGREE SIGN}C)', size=10, weight='bold')\n",
    "    # Plot isopycnals\n",
    "    CS1 = ax1.contour(x,z[:24],rho1[t,:24,:],crho, colors = 'b')\n",
    "    plt.clabel(CS1, fmt='%1.1f', fontsize=15, inline=1)\n",
    "    ax1.set_ylabel('Depth (m)', size=15)\n",
    "    ax1.tick_params(labelsize=12)\n",
    "    ax1.set_title('Qf-Only', loc='left', size=15, weight='bold')\n",
    "    ax1.set_title('Hour ' + str(t), loc='center', size=20, weight='bold')\n",
    "    \n",
    "    ut1 = u1[t,:,:]\n",
    "    wt1 = w1[t,:,:]\n",
    "    ut1[[1,2,4,5,7,8,10,12,14,16,18,20,22,24,26,28], :] = np.nan\n",
    "    wt1[[1,2,4,5,7,8,10,12,14,16,18,20,22,24,26,28], :] = np.nan\n",
    "    ut1[::2,:] = np.nan\n",
    "    wt1[::2,:] = np.nan   \n",
    "    Q = ax1.quiver(x[::5],z[:24],0.1*ut1[:24,::5], 0.1*wt1[:24,::5], width=0.001, scale=1)\n",
    "    ax1.quiverkey(Q, 0.7, 0.89, 0.02,r'0.2 $m/s$', fontproperties={'weight': 'bold', 'size' : 10}, coordinates='figure')\n",
    "\n",
    "\n",
    "\n",
    "    # qf-tide Plot\n",
    "    ax2 = fig.add_subplot(212)\n",
    "    im2 = ax2.contourf(x,z[:24],pt2[t,:24,:], ct, cmap=cmap)\n",
    "    ax2.set_xlim([x[1],x[-1]])\n",
    "    ax2.set_facecolor('gray')\n",
    "    CS2 = ax2.contour(x,z[:24],rho2[t,:24,:],crho, colors = 'b')\n",
    "    plt.clabel(CS2, fmt='%1.1f', fontsize=15, inline=1)\n",
    "    ax2.set_ylabel('Depth (m)', size=15)\n",
    "    ax2.tick_params(labelsize=12)\n",
    "    ax2.set_title('Tide-Only', loc='left', size=15, weight='bold')\n",
    "    \n",
    "    ut2 = u2[t,:,:]\n",
    "    wt2 = w2[t,:,:]\n",
    "    ut2[[1,2,4,5,7,8,10,12,14,16,18,20,22,24,26,28], :] = np.nan\n",
    "    wt2[[1,2,4,5,7,8,10,12,14,16,18,20,22,24,26,28], :] = np.nan\n",
    "    ut2[::2,:] = np.nan\n",
    "    wt2[::2,:] = np.nan   \n",
    "    ax2.quiver(x[::5],z[:24],0.1*ut2[:24,::5], 0.1*wt2[:24,::5], width=0.001, scale=1)\n",
    "    ax2.set_xlabel('Distance from glacier (km)', size=15)\n",
    "    ax2.tick_params(labelsize=12)\n",
    "    plt.close()\n",
    "    fig.savefig(fig_path + '/time_r' + format(t, '03d') + '.png', dpi = 150, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate GIF from saved output figures\n",
    "import imageio\n",
    "import glob\n",
    "\n",
    "images = []\n",
    "imagenames = sorted(glob.glob(fig_path + '/time_r*.png'))\n",
    "\n",
    "for filename in imagenames:\n",
    "    images.append(imageio.imread(filename))\n",
    "imageio.mimsave(fig_path + 'Comp009_r.gif', images, duration=.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}