plot_ice.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e699727f",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import glob\n",
    "import netCDF4 as nc\n",
    "import datetime as dt\n",
    "import sys\n",
    "import gsw as sw\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import matplotlib\n",
    "matplotlib.use('agg')\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.path as mpath\n",
    "import cftime\n",
    "import coast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c57df360",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "\n",
    "var = 'siconc' # thetao, so, uo, vo, siconc, siage, sivol, sithick, siu, siv, \n",
    "i_o = 'SI' # SI or O for sea ice or ocean\n",
    "freq = 'mon' # mon or day\n",
    "time_s = 'highres-future' # 'highres-future' or 'hist-1950'\n",
    "\n",
    "def make_path(var, i_o, freq, time_s):\n",
    "    if 'future' in time_s:\n",
    "        ddir = 'MOHC'\n",
    "    else:\n",
    "        ddir = 'NERC'\n",
    "    root = '/badc/cmip6/data/CMIP6/HighResMIP/' + ddir + '/HadGEM3-GC31-HH/' + time_s + '/r1i1p1f1/'\n",
    "    return root + i_o + freq + '/' + var + '/gn/latest/' + var + '_' + i_o + freq + '_HadGEM3-GC31-HH_' + time_s + '_r1i1p1f1_gn_*.nc'\n",
    "\n",
    "fn_nemo_dat_sic = make_path(var, i_o, freq, time_s)\n",
    "fn_nemo_dat_sit = make_path('sithick', i_o, freq, time_s)\n",
    "fn_nemo_dat_t = make_path('thetao', 'O', freq, time_s)\n",
    "\n",
    "domain_root = '/gws/nopw/j04/nemo_vol5/acc/eORCA12-N512/domain/'\n",
    "fn_nemo_dom1 = domain_root + 'eORCA12_coordinates.nc'\n",
    "fn_nemo_dom = domain_root + 'mesh_mask_eORCA12_v2.4.nc'\n",
    "fn_nemo_bathy = domain_root + 'eORCA12_bathymetry_v2.4.nc'\n",
    "fn_config_t_grid = './config/gc31_nemo_grid_t.json'\n",
    "\n",
    "out_file = './Processed/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d47c6799",
   "metadata": {},
   "outputs": [],
   "source": [
    "flist_sic = sorted(glob.glob(fn_nemo_dat_sic))\n",
    "flist_sit = sorted(glob.glob(fn_nemo_dat_sit))\n",
    "flist_t = sorted(glob.glob(fn_nemo_dat_t))\n",
    "\n",
    "use_xarray = False\n",
    "if use_xarray:\n",
    "    nemo_t = coast.Gridded(fn_data = flist_si[-1], fn_domain = fn_nemo_dom, config=fn_config_t_grid)\n",
    "\n",
    "    print(nemo_t.dataset.longitude.shape)\n",
    "    with nc.Dataset(fn_nemo_dom, 'r') as nc_fid:\n",
    "        tmask = nc_fid.variables['tmask'][0, :, 1:-1, 1:-1]\n",
    "else:\n",
    "    v_map = {}\n",
    "    v_map['e1t'] = 'e1t'\n",
    "    v_map['e2t'] = 'e2t'\n",
    "    v_map['e3t_0'] = 'e3t_0'\n",
    "    v_map['tmask'] = 'tmask'\n",
    "    v_map['lat'] = 'latitude'\n",
    "    v_map['lon'] = 'longitude'\n",
    "    v_map['depth'] = 'lev'\n",
    "    v_map['time'] = 'time'\n",
    "    v_map['temp'] = 'thetao'\n",
    "    v_map['sal'] = 'so' \n",
    "    v_map['siconc'] = 'siconc'\n",
    "    v_map['sithick'] = 'sithick'\n",
    "    \n",
    "    with nc.Dataset(fn_nemo_dom, 'r') as nc_fid:\n",
    "        tmask = nc_fid.variables[v_map['tmask']][0, :, 1:-1, 1:-1]\n",
    "        \n",
    "    with nc.Dataset(flist_sic[0], 'r') as nc_fid:\n",
    "        lat = nc_fid.variables[v_map['lat']][:]\n",
    "        lon = nc_fid.variables[v_map['lon']][:]\n",
    "        siconc = nc_fid.variables[v_map['siconc']][0, ...]\n",
    "    with nc.Dataset(flist_sit[0], 'r') as nc_fid:\n",
    "        sithick = nc_fid.variables[v_map['sithick']][0, ...]\n",
    "        \n",
    "    with nc.Dataset(flist_t[0], 'r') as nc_fid:\n",
    "        lat = nc_fid.variables[v_map['lat']][:]\n",
    "        lon = nc_fid.variables[v_map['lon']][:]\n",
    "        \n",
    "    #siconc = np.ma.masked_where((siconc==1e20) | (siconc==0), siconc)\n",
    "    #mask = siconc.mask\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "debf5187-857b-47bd-90c4-22c1d3cd7cf8",
   "metadata": {},
   "source": [
    "Subset data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dec2799",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "#print(nemo_t.dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8c5164b",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "yi1 = 2800\n",
    "    \n",
    "if use_xarray:\n",
    "    yi2 = nemo_t.dataset.longitude.shape[0]\n",
    "    ind_y = np.arange(yi1, yi2).astype(int)\n",
    "    #print(yi2, ind_y)\n",
    "    nemo_t_subset = nemo_t.isel(y_dim=ind_y)\n",
    "\n",
    "    print(nemo_t_subset.dataset)\n",
    "else:\n",
    "    lat = lat[yi1:, :]\n",
    "    lon = lon[yi1:, :]\n",
    "    siconc = siconc[yi1:, :]\n",
    "    sithick = sithick[yi1:, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90c543ee-0efc-48cc-86a5-7139ab51579c",
   "metadata": {},
   "source": [
    "Time slice. Notes dates are 360 day years so use cftime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8967cfd1-5c17-4c46-bb59-b187eb9ad5ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Does not work\n",
    "#st_date = cftime.datetime(2050, 12, 1)\n",
    "#en_date = cftime.datetime(2051, 1, 1) \n",
    "#nemo_t_subset = nemo_t_subset.time_slice(st_date, en_date)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "481b58fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# change this to decrease resolution but decrease run time\n",
    "sub = 5\n",
    "if use_xarray:\n",
    "    ind_y = np.arange(0, nemo_t_subset.dataset.longitude.shape[0], sub).astype(int)\n",
    "    ind_x = np.arange(0, nemo_t_subset.dataset.longitude.shape[1], sub).astype(int)\n",
    "    nemo_t_subset = nemo_t_subset.isel(y_dim=ind_y, x_dim=ind_x)\n",
    "        \n",
    "    lat = nemo_t_subset.dataset.latitude.values\n",
    "    lon = nemo_t_subset.dataset.longitude.values\n",
    "    lat = np.ma.masked_where((lon==0), lat)\n",
    "    lon = np.ma.masked_where((lon==0), lon)\n",
    "    siconc = nemo_t_subset.dataset.siconc.to_masked_array() # time, lev, j, i\n",
    "else:\n",
    "    lat = lat[::sub, ::sub]\n",
    "    lon = lon[::sub, ::sub]\n",
    "    siconc = siconc[::sub, ::sub]\n",
    "    sithick = sithick[::sub, ::sub]\n",
    "print(siconc.shape)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee167fa5-d9b1-46e5-831a-418c4666e2b7",
   "metadata": {},
   "source": [
    "Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61d0616d",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "def set_circle(ax):\n",
    "  # Compute a circle in axes coordinates, which we can use as a boundary\n",
    "  # for the map. We can pan/zoom as much as we like - the boundary will be\n",
    "  # permanently circular.\n",
    "  theta = np.linspace(0, 2 * np.pi, 100)\n",
    "  center, radius = [0.5, 0.5], 0.5\n",
    "  verts = np.vstack([np.sin(theta), np.cos(theta)]).T\n",
    "  circle = mpath.Path(verts * radius + center)\n",
    "  ax.set_boundary(circle, transform=ax.transAxes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0230b4a",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "data_crs = ccrs.PlateCarree()\n",
    "mrc = ccrs.NorthPolarStereo(central_longitude=0.0)\n",
    "\n",
    "my_cm = plt.cm.plasma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90d7f202",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "fig1 = plt.figure(figsize=(12, 6))\n",
    "ax1 = fig1.add_axes([0.03, 0.03, 0.3, 0.85], projection=mrc)\n",
    "ax2 = fig1.add_axes([0.36, 0.03, 0.3, 0.85], projection=mrc)\n",
    "ax3 = fig1.add_axes([0.69, 0.03, 0.3, 0.85], projection=mrc)\n",
    "cax1 = fig1.add_axes([0.03, 0.94, 0.3, 0.02])\n",
    "cax2 = fig1.add_axes([0.36, 0.94, 0.3, 0.02])\n",
    "cax3 = fig1.add_axes([0.69, 0.94, 0.3, 0.02])\n",
    "\n",
    "cs1 = ax1.pcolormesh(lon, lat, siconc[:, :], transform=data_crs, cmap=my_cm)\n",
    "cs2 = ax2.pcolormesh(lon, lat, sithick[:, :], transform=data_crs, cmap=my_cm)\n",
    "\n",
    "ax1.add_feature(cfeature.LAND, zorder=100)\n",
    "ax2.add_feature(cfeature.LAND, zorder=100)\n",
    "#ax3.add_feature(cfeature.LAND, zorder=100)\n",
    "ax1.gridlines()\n",
    "#ax2.gridlines()\n",
    "#ax3.gridlines()\n",
    "\n",
    "ax1.set_extent([-180, 180, 60, 90], crs=data_crs)\n",
    "ax2.set_extent([-180, 180, 60, 90], crs=data_crs)\n",
    "#ax3.set_extent([-180, 180, 60, 90], crs=data_crs)\n",
    "set_circle(ax1)\n",
    "set_circle(ax2)\n",
    "#set_circle(ax3)\n",
    "\n",
    "\n",
    "ax1.annotate('(a)', (0.05, 0.95), xycoords='axes fraction', bbox=dict(boxstyle=\"round\", fc=\"w\"), zorder=105)\n",
    "ax2.annotate('(b)', (0.05, 0.95), xycoords='axes fraction', bbox=dict(boxstyle=\"round\", fc=\"w\"), zorder=105)\n",
    "ax3.annotate('(c)', (0.05, 0.95), xycoords='axes fraction', bbox=dict(boxstyle=\"round\", fc=\"w\"), zorder=105)\n",
    "\n",
    "fig1.colorbar(cs1, cax=cax1, orientation='horizontal')\n",
    "fig1.colorbar(cs2, cax=cax2, orientation='horizontal')\n",
    "#fig1.colorbar(cs3, cax=cax3, orientation='horizontal')\n",
    "\n",
    "cax1.set_xlabel('Ice Concentration (%)')\n",
    "cax2.set_xlabel('Ice Thickness (m)')\n",
    "#cax3.set_xlabel('Density (kg/m$^{3}$)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3278c989",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "fig1.savefig('./Figures/ice_conc.png')\n"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "encoding": "# coding: utf-8",
   "executable": "/usr/bin/env python",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  },
  "kernelspec": {
   "display_name": "Python 3 + Jaspy",
   "language": "python",
   "name": "jaspy"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}