CEE_Figure5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "44d5d6a3-1dc1-4f93-a5eb-5384ad05a0da",
   "metadata": {},
   "source": [
    "# Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b919588b-06a9-4259-bcab-1788c9202de0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import cv2\n",
    "import cmaps\n",
    "import cmocean\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import scipy.io as sio\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import matplotlib.gridspec as gridspec\n",
    "from netCDF4 import Dataset\n",
    "from datetime import date\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "from matplotlib.pyplot import Polygon\n",
    "from matplotlib import rcParams\n",
    "from matplotlib.backends.backend_pdf import PdfPages\n",
    "version = mpl.__version__\n",
    "rcParams['font.family'] = 'sans-serif'\n",
    "directory   = '/srv/scratch/z3533156'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e591f60b-a019-46f6-a78d-b1e8a24b1e2e",
   "metadata": {},
   "source": [
    "# Import time functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5da018aa-6ff2-42cf-9540-691b133fcce7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def datestring_to_serial_day(datestring,epochY=1990,epochm=1,epochd=1,epochH=0,epochM=0):\n",
    "    import pandas as pd\n",
    "    import datetime\n",
    "    serial_day_timedelta = pd.to_datetime(datestring) - datetime.datetime(epochY,epochm,epochd,epochH,epochM)\n",
    "    corrected_serial_day_number = serial_day_timedelta.days + serial_day_timedelta.seconds/86400\n",
    "    return corrected_serial_day_number\n",
    "def serial_day_to_datestring(day,epochY=1990,epochm=1,epochd=1,epochH=0,epochM=0):\n",
    "    import datetime\n",
    "    corrected_date = datetime.datetime(epochY,epochm,epochd,epochH,epochM) + datetime.timedelta(day)\n",
    "    return corrected_date.strftime(\"%Y-%m-%d\") \n",
    "def serial_day_to_datestring1(day,epochY=1990,epochm=1,epochd=1,epochH=0,epochM=0):\n",
    "    import datetime\n",
    "    corrected_date = datetime.datetime(epochY,epochm,epochd,epochH,epochM) + datetime.timedelta(day)\n",
    "    return corrected_date.strftime(\"%Y-\") \n",
    "def serial_day_to_datestring2(day,epochY=1990,epochm=1,epochd=1,epochH=0,epochM=0):\n",
    "    import datetime\n",
    "    corrected_date = datetime.datetime(epochY,epochm,epochd,epochH,epochM) + datetime.timedelta(day)\n",
    "    return corrected_date.strftime(\"%m-%d\")  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b971c07-d943-4f62-83f5-10abb2fee7c0",
   "metadata": {},
   "source": [
    "# Read data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "16691bfd-10c6-4278-b560-18908ded636e",
   "metadata": {},
   "outputs": [],
   "source": [
    "AVISO_SLA   = Dataset(directory+'/MHW/AVISO/AVISO_EAC_MHW.nc','r')\n",
    "SLA_daily   = np.transpose(AVISO_SLA.variables['sla'],[2,1,0])\n",
    "SSH_daily   = np.transpose(AVISO_SLA.variables['adt'],[2,1,0])\n",
    "lon0        = AVISO_SLA.variables['lon'][:]\n",
    "lat0        = AVISO_SLA.variables['lat'][:]\n",
    "lon         = np.tile(lon0,[np.size(lat0,0),1]).transpose()\n",
    "lat         = np.tile(lat0,[np.size(lon0,0),1])\n",
    "dataset     = sio.loadmat(directory+'/MHW/Figure_data/Figure1_OSTIA_MHW.mat')\n",
    "time_daily  = dataset['time_daily'][:,0]/86400\n",
    "time        = datestring_to_serial_day(\"1981-01-01\",epochY=1990,epochm=1,epochd=1,epochH=0,epochM=0)\n",
    "daily_t     = np.arange(date(1993,1,1).toordinal(),date(2022,12,31).toordinal()+1,1)\n",
    "daily_t1    = np.arange(date(2021,12,1).toordinal(),date(2022,2,28).toordinal()+1,1)\n",
    "daily_dates = [date.fromordinal(tt.astype(int)) for tt in daily_t1]\n",
    "SLA         = SLA_daily[:,:,(daily_t>=date(2021,12,1).toordinal()) & (daily_t<=date(2022,2,28).toordinal())]\n",
    "SSH         = SSH_daily[:,:,(daily_t>=date(2021,12,1).toordinal()) & (daily_t<=date(2022,2,28).toordinal())]\n",
    "dataset1       = sio.loadmat(directory+'/MHW/Figure_data/Figure_etopo_AVISO.mat')\n",
    "depth          = dataset1['depth'][:,:]\n",
    "depth[lon>155] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f44f7a0c-f904-49b4-ab01-714bedc4872d",
   "metadata": {},
   "source": [
    "# Plot the spatial distribution of weekly sea level anomalies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02ff5b38-2021-47b0-aa93-756f2df21fdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(80, 60))\n",
    "gs  = gridspec.GridSpec(3,4) \n",
    "m   = Basemap(projection='merc',llcrnrlat=-38.0-0.001,urcrnrlat=-28+0.001,llcrnrlon=150.0-0.001,urcrnrlon=160-0.001,resolution='i') \n",
    "x2, y2 = m(lon, lat) \n",
    "padspacescale = 50 \n",
    "labelpadscale = 15 \n",
    "linefont  = 12 \n",
    "labelfont = 70 \n",
    "scale1    = 1.35 \n",
    "scale2    = 1.2 \n",
    "levels1   = np.linspace(-0.5,0.5,50) \n",
    "levels2   = np.linspace(-100,100,50) \n",
    "tick_marks1 = np.linspace(-0.5,0.5,11) \n",
    "tick_marks2 = np.linspace(-100,100,9) \n",
    "sla_levels1 = np.linspace(0,0.5,3) \n",
    "sla_levels2 = np.linspace(-0.5,0,3) \n",
    "sla_levels3 = np.linspace(-0.001,0.001,2) \n",
    "ssh_levels1 = np.linspace(0.949,0.951,2) \n",
    "cmaps1 = cmaps.cmocean_balance \n",
    "labels = ['a','b','c','d','e','f','g','h','i','j','k','l'] \n",
    "for i in range(12): \n",
    "    k = i*7\n",
    "    if i<4: \n",
    "        j = np.mod(i,4) \n",
    "    elif i<8: \n",
    "        j = 4 + np.mod(i,4) \n",
    "    else: \n",
    "        j = 8 + np.mod(i,4)\n",
    "    ax = fig.add_subplot(gs[j]) \n",
    "    l, b, w, h = ax.get_position().bounds \n",
    "    m = Basemap(projection='merc',llcrnrlat=-38.0-0.001,urcrnrlat=-28+0.001,llcrnrlon=150.0-0.001,urcrnrlon=160-0.001,resolution='i') \n",
    "    m.drawcoastlines(color='0.1', linewidth=0.5*linefont)\n",
    "    m.drawmapboundary(color='0.1', linewidth=0.5*linefont) \n",
    "    m.fillcontinents(color='0.95', lake_color='white') \n",
    "    m.drawmeridians(np.arange(150,165,2),labels=[0,0,0,0],linewidth=0.4*linefont,dashes=[2,2],color='gray',fontsize=1.1*labelfont)\n",
    "    CB1 = m.contourf(x2, y2, SLA[:,:,k], \n",
    "    cmap=cmaps1,levels=levels1,origin='lower',extend='both')\n",
    "    CS1 = m.contour(x2, y2, SSH[:,:,k], ssh_levels1,linewidths=1.2*linefont,linestyles='solid',colors='black') \n",
    "    CS2 = m.contour(x2, y2, SLA[:,:,k], sla_levels3,linewidths=1.2*linefont,linestyles='solid',colors='gray') \n",
    "    CS3 = m.contour(x2, y2, depth, np.arange(-201,-199,10), linewidths=1.2*linefont,linestyles='dashed',colors='xkcd:purple')\n",
    "    lon2 = np.array([lon[22,20],lon[60,20],lon[60,33],lon[30,33],lon[22,20]]) \n",
    "    lat2 = np.array([lat[22,20],lat[60,20],lat[60,33],lat[30,33],lon[22,20]]) \n",
    "    x22,y22 = m(lon2,lat2) \n",
    "    lon3 = np.array([lon[22,20],lon[25,20],lon[33,33],lon[30,33],lon[22,20]]) \n",
    "    lat3 = np.array([lat[22,20],lat[25,20],lat[33,33],lat[30,33],lat[22,20]]) \n",
    "    x33,y33 = m(lon3,lat3) \n",
    "    m.plot(x33, y33, linewidth=labelpadscale, linestyle='dashed', color='green') \n",
    "    if np.mod(i,4)==0: \n",
    "        plt.ylabel('Latitude',fontsize=labelfont,labelpad=3.5*padspacescale,family='sans-serif') \n",
    "        m.drawparallels(np.arange(-40,-10, 2),labels=[1,0,0,0],linewidth=0.4*linefont,dashes=[2,2],color='gray',fontsize=1.1*labelfont) \n",
    "    elif np.mod(i,4)==1: \n",
    "        l=l-0.006 \n",
    "        m.drawparallels(np.arange(-40,-10, 2),labels=[0,0,0,0],linewidth=0.4*linefont,dashes=[2,2],color='gray',fontsize=1.1*labelfont) \n",
    "    elif np.mod(i,4)==2: \n",
    "        l=l-0.012 \n",
    "        m.drawparallels(np.arange(-40,-10, 2),labels=[0,0,0,0],linewidth=0.4*linefont,dashes=[2,2],color='gray',fontsize=1.1*labelfont) \n",
    "    else: \n",
    "        l=l-0.018 \n",
    "        m.drawparallels(np.arange(-40,-10, 2),labels=[0,0,0,0],linewidth=0.4*linefont,dashes=[2,2],color='gray',fontsize=1.1*labelfont) \n",
    "    ax.spines['bottom'].set_linewidth(labelpadscale) \n",
    "    ax.spines['left'].set_linewidth(labelpadscale) \n",
    "    ax.spines['top'].set_linewidth(labelpadscale) \n",
    "    ax.spines['right'].set_linewidth(labelpadscale) \n",
    "    cx1,cy1 = m(150.2, -28.8) \n",
    "    cx2,cy2 = m(150.2, -29.6) \n",
    "    plt.text(cx1, cy1,serial_day_to_datestring1(time + time_daily[10561+k]), color='dodgerblue', fontsize=1.2*labelfont, weight='bold',family='sans-serif') \n",
    "    plt.text(cx2, cy2,serial_day_to_datestring2(time + time_daily[10561+k]), color='dodgerblue', fontsize=1.2*labelfont, weight='bold',family='sans-serif') \n",
    "    if i<4: \n",
    "        b = b \n",
    "    elif i<8: \n",
    "        b = b-0.08 \n",
    "    else:\n",
    "        b = b-0.16 \n",
    "    plt.xlabel('Longitude',fontsize=labelfont,labelpad=1.5*padspacescale,family='sans-serif') \n",
    "    m.drawmeridians(np.arange(150,165, 4),labels=[0,0,0,1],linewidth=0.4*linefont,dashes=[2,2],color='.7',fontsize=1.1*labelfont)\n",
    "    ax.set_position([l, b, scale1*w, scale1*h])\n",
    "    plt.tick_params(axis='both',which='major',bottom='on',left='on',length=40.0,width=15,colors='black',direction='out') \n",
    "    plt.title(labels[i], fontsize=1.2*labelfont,loc='left',pad=0.6*padspacescale, weight='bold',family='sans-serif') \n",
    "    cbaxes1 = fig.add_axes([0.15, -0.12, 0.75, 0.02]) \n",
    "    cb1 = plt.colorbar(CB1,orientation='horizontal',cax = cbaxes1) \n",
    "    cb1.set_ticks(tick_marks1) \n",
    "    cb1.set_label(r'Sea level anomalies (m)', fontsize=labelfont,labelpad=0) \n",
    "    cb1.ax.tick_params(labelsize=labelfont) \n",
    "#############################################################################################\n",
    "fig.savefig(directory+'/MHW/Figure_plots/Figure5_weekly_SLA.pdf',dpi=300,bbox_inches = 'tight')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.8 (ipykernel)",
   "language": "python",
   "name": "python-3.10.8"
  },
  "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}