{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Sea ice analysis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Work in Progress

\n", "

This notebook is currently under development. Content and results may change!

\n", "
" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%load_ext autoreload\n", "%autoreload " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Basemap module not found. Some regional plots may not function properly\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "ERROR 1: PROJ: proj_create_from_database: Open of /glade/work/gmarques/conda-envs/mom6-tools/share/proj failed\n" ] } ], "source": [ "%matplotlib inline\n", "import warnings, yaml\n", "warnings.filterwarnings(\"ignore\")\n", "from mom6_tools.MOM6grid import MOM6grid\n", "from mom6_tools.m6toolbox import cime_xmlquery\n", "from ncar_jobqueue import NCARCluster\n", "from dask.distributed import Client\n", "import matplotlib.pyplot as plt\n", "import cartopy.crs as ccrs\n", "import cartopy.feature\n", "import numpy as np\n", "import xarray as xr\n", "from mom6_tools.m6plot import polarplot" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Read in the yaml file\n", "diag_config_yml_path = \"diag_config.yml\"\n", "diag_config_yml = yaml.load(open(diag_config_yml_path,'r'), Loader=yaml.Loader)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Rundir directory is: /glade/derecho/scratch/gmarques/g.e30_a03c.GJRAv4.TL319_t232_wgx3_hycom1_N75.2024.079/run\n", "Ice history files are at: /glade/derecho/scratch/gmarques/archive/g.e30_a03c.GJRAv4.TL319_t232_wgx3_hycom1_N75.2024.079/ice/hist/\n", "Casename is: g.e30_a03c.GJRAv4.TL319_t232_wgx3_hycom1_N75.2024.079\n" ] } ], "source": [ "caseroot = diag_config_yml['Case']['CASEROOT']\n", "casename = cime_xmlquery(caseroot, 'CASE')\n", "DOUT_S = cime_xmlquery(caseroot, 'DOUT_S')\n", "rundir = cime_xmlquery(caseroot, 'RUNDIR')\n", "\n", "if DOUT_S:\n", " OUTDIR = cime_xmlquery(caseroot, 'DOUT_S_ROOT')+'/ocn/hist/'\n", " OUTDIR_ice = cime_xmlquery(caseroot, 'DOUT_S_ROOT')+'/ice/hist/'\n", "else:\n", " OUTDIR = cime_xmlquery(caseroot, 'RUNDIR')\n", " OUTDIR_ice = OUTDIR\n", " \n", "print('Rundir directory is:', rundir)\n", "print('Ice history files are at:', OUTDIR_ice)\n", "print('Casename is:', casename)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# create an empty class object\n", "class args:\n", " pass# create an empty class object\n", " \n", "# load avg dates\n", "avg = diag_config_yml['Avg']\n", "\n", "args.start_date = avg['start_date']\n", "args.end_date = avg['end_date']\n", "args.casename = casename\n", "args.ice = casename+diag_config_yml['Fnames']['ice']\n", "args.geom = casename+diag_config_yml['Fnames']['geom']\n", "args.static = casename+diag_config_yml['Fnames']['static']\n", "args.savefigs = False\n", "args.time_series = True\n", "args.nw = 6 # requesting 6 workers" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 2min 4s, sys: 2min 27s, total: 4min 32s\n", "Wall time: 9min 7s\n" ] } ], "source": [ "variables = ['hi', 'aice']\n", "\n", "def preprocess(ds):\n", " '''Return the dataset with variables'''\n", " return ds[variables] \n", "\n", "%time ds = xr.open_mfdataset(OUTDIR_ice+args.ice, preprocess=preprocess)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset> Size: 2GB\n",
       "Dimensions:  (time: 732, nj: 480, ni: 540)\n",
       "Coordinates:\n",
       "  * time     (time) object 6kB 0001-01-16 12:00:00 ... 0061-12-16 12:00:00\n",
       "    TLON     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>\n",
       "    TLAT     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>\n",
       "    ULON     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>\n",
       "    ULAT     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>\n",
       "    NLON     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>\n",
       "    NLAT     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>\n",
       "    ELON     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>\n",
       "    ELAT     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>\n",
       "Dimensions without coordinates: nj, ni\n",
       "Data variables:\n",
       "    hi       (time, nj, ni) float32 759MB dask.array<chunksize=(1, 480, 540), meta=np.ndarray>\n",
       "    aice     (time, nj, ni) float32 759MB dask.array<chunksize=(1, 480, 540), meta=np.ndarray>\n",
       "Attributes:\n",
       "    title:               g.e30_a03c.GJRAv4.TL319_t232_wgx3_hycom1_N75.2024.079\n",
       "    contents:            Diagnostic and Prognostic Variables\n",
       "    source:              CICE Sea Ice Model, unknown_version_name\n",
       "    comment:             All years have exactly 365 days\n",
       "    comment2:            File written on model date 00010201\n",
       "    comment3:            seconds elapsed into model date:      0\n",
       "    time_period_freq:    month_1\n",
       "    time_axis_position:  middle\n",
       "    conventions:         CF-1.0\n",
       "    history:             This dataset was created on 2024-10-02 at 15:04\n",
       "    io_flavor:           io_pio2 cdf1
" ], "text/plain": [ " Size: 2GB\n", "Dimensions: (time: 732, nj: 480, ni: 540)\n", "Coordinates:\n", " * time (time) object 6kB 0001-01-16 12:00:00 ... 0061-12-16 12:00:00\n", " TLON (nj, ni) float32 1MB dask.array\n", " TLAT (nj, ni) float32 1MB dask.array\n", " ULON (nj, ni) float32 1MB dask.array\n", " ULAT (nj, ni) float32 1MB dask.array\n", " NLON (nj, ni) float32 1MB dask.array\n", " NLAT (nj, ni) float32 1MB dask.array\n", " ELON (nj, ni) float32 1MB dask.array\n", " ELAT (nj, ni) float32 1MB dask.array\n", "Dimensions without coordinates: nj, ni\n", "Data variables:\n", " hi (time, nj, ni) float32 759MB dask.array\n", " aice (time, nj, ni) float32 759MB dask.array\n", "Attributes:\n", " title: g.e30_a03c.GJRAv4.TL319_t232_wgx3_hycom1_N75.2024.079\n", " contents: Diagnostic and Prognostic Variables\n", " source: CICE Sea Ice Model, unknown_version_name\n", " comment: All years have exactly 365 days\n", " comment2: File written on model date 00010201\n", " comment3: seconds elapsed into model date: 0\n", " time_period_freq: month_1\n", " time_axis_position: middle\n", " conventions: CF-1.0\n", " history: This dataset was created on 2024-10-02 at 15:04\n", " io_flavor: io_pio2 cdf1" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " Selecting data between 0031-01-01 and 0062-01-01...\n", "CPU times: user 10.4 ms, sys: 0 ns, total: 10.4 ms\n", "Wall time: 44.3 ms\n" ] } ], "source": [ "print('\\n Selecting data between {} and {}...'.format(args.start_date, args.end_date))\n", "%time ds_sel = ds.sel(time=slice(args.start_date, args.end_date))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " Computing time mean...\n", "CPU times: user 1min 18s, sys: 2min 53s, total: 4min 12s\n", "Wall time: 9min 18s\n" ] } ], "source": [ "print('\\n Computing time mean...')\n", "%time ds_sel_mean = ds_sel.mean('time').compute()" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Computing monthly climatology...\n", "CPU times: user 1min 29s, sys: 3min 34s, total: 5min 4s\n", "Wall time: 10min 54s\n" ] } ], "source": [ "print('Computing monthly climatology...')\n", "%time ds_monthly = ds_sel.groupby(\"time.month\").mean('time').compute()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# load obs\n", "obs_path = '/glade/work/gmarques/cesm/datasets/Seaice/'\n", "obs_NH = xr.open_dataset(obs_path+'OBS_NSIDC_sat_NH_T2Ms_sic.nc')\n", "obs_SH = xr.open_dataset(obs_path+'OBS_NSIDC_sat_SH_T2Ms_sic.nc')" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "obs_NH_mean = obs_NH.mean('time').compute()\n", "obs_SH_mean = obs_SH.mean('time').compute()" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "obs_NH_monthly = obs_NH.groupby(\"time.month\").mean('time').compute()\n", "obs_SH_monthly = obs_SH.groupby(\"time.month\").mean('time').compute()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sea ice thickness" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Annual mean" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "ename": "TypeError", "evalue": "gca() got an unexpected keyword argument 'projection'", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn[14], line 9\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[38;5;66;03m# SH\u001b[39;00m\n\u001b[1;32m 8\u001b[0m fig \u001b[38;5;241m=\u001b[39m plt\u001b[38;5;241m.\u001b[39mfigure(figsize\u001b[38;5;241m=\u001b[39m(\u001b[38;5;241m10\u001b[39m,\u001b[38;5;241m8\u001b[39m))\n\u001b[0;32m----> 9\u001b[0m ax \u001b[38;5;241m=\u001b[39m \u001b[43mplt\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mgca\u001b[49m\u001b[43m(\u001b[49m\u001b[43mprojection\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mccrs\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mSouthPolarStereo\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 10\u001b[0m polarplot(hi_mean, grd, title\u001b[38;5;241m=\u001b[39mdcase\u001b[38;5;241m.\u001b[39mcasename, debug\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m, colormap\u001b[38;5;241m=\u001b[39mplt\u001b[38;5;241m.\u001b[39mcm\u001b[38;5;241m.\u001b[39mgist_ncar, clim\u001b[38;5;241m=\u001b[39m(\u001b[38;5;241m0\u001b[39m,\u001b[38;5;241m2\u001b[39m), axis\u001b[38;5;241m=\u001b[39max)\n\u001b[1;32m 11\u001b[0m plt\u001b[38;5;241m.\u001b[39msuptitle(suptitle)\n", "\u001b[0;31mTypeError\u001b[0m: gca() got an unexpected keyword argument 'projection'" ] }, { "data": { "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Sea ice thickness\n", "aice_mean = np.ma.masked_invalid(ds_monthly['aice'].mean('month').values) * 100.\n", "hi_mean = np.ma.masked_invalid(ds_monthly['hi'].mean('month').values)\n", "hi_mean = np.ma.masked_where(aice_mean < 15.0, hi_mean)\n", "suptitle = ('ANN mean sea ice thickness, ' + str(args.start_date) + ' to ' + str(args.end_date))\n", "%matplotlib inline\n", "# SH\n", "fig = plt.figure(figsize=(10,8))\n", "ax = plt.gca(projection=ccrs.SouthPolarStereo())\n", "polarplot(hi_mean, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,2), axis=ax)\n", "plt.suptitle(suptitle)\n", "\n", "# NH\n", "fig = plt.figure(figsize=(10,8))\n", "ax = plt.gca(projection=ccrs.NorthPolarStereo())\n", "polarplot(hi_mean, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,4), axis=ax, proj='NP')\n", "plt.suptitle(suptitle)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### JFM" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "aice_JFM = np.ma.masked_invalid(ds_monthly['aice'].sel(month=[1,2,3]).mean('month').values) * 100.\n", "hi_JFM = np.ma.masked_invalid(ds_monthly['hi'].sel(month=[1,2,3]).mean('month').values)\n", "hi_JFM = np.ma.masked_where(aice_JFM < 15.0, hi_JFM)\n", "suptitle = ('JFM mean sea ice thickness, ' + str(args.start_date) + ' to ' + str(args.end_date))\n", "%matplotlib inline\n", "# SH\n", "fig = plt.figure(figsize=(10,8))\n", "ax = plt.gca(projection=ccrs.SouthPolarStereo())\n", "polarplot(hi_JFM, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,2), axis=ax)\n", "plt.suptitle(suptitle)\n", "\n", "# NH\n", "fig = plt.figure(figsize=(10,8))\n", "ax = plt.gca(projection=ccrs.NorthPolarStereo())\n", "polarplot(hi_JFM, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,4), axis=ax, proj='NP')\n", "plt.suptitle(suptitle)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### AMJ" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "aice_AMJ = np.ma.masked_invalid(ds_monthly['aice'].sel(month=[4,5,6]).mean('month').values) * 100.\n", "hi_AMJ = np.ma.masked_invalid(ds_monthly['hi'].sel(month=[4,5,6]).mean('month').values)\n", "hi_AMJ = np.ma.masked_where(aice_AMJ < 15.0, hi_AMJ)\n", "suptitle = ('AMJ mean sea ice thickness, ' + str(args.start_date) + ' to ' + str(args.end_date))\n", "%matplotlib inline \n", "# SH\n", "fig = plt.figure(figsize=(10,8))\n", "ax = plt.gca(projection=ccrs.SouthPolarStereo())\n", "polarplot(hi_AMJ, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,2), axis=ax)\n", "plt.suptitle(suptitle)\n", "\n", "# NH\n", "fig = plt.figure(figsize=(10,8))\n", "ax = plt.gca(projection=ccrs.NorthPolarStereo())\n", "polarplot(hi_AMJ, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,4), axis=ax, proj='NP')\n", "plt.suptitle(suptitle)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### JAS" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "aice_JAS = np.ma.masked_invalid(ds_monthly['aice'].sel(month=[7,8,9]).mean('month').values) * 100.\n", "hi_JAS = np.ma.masked_invalid(ds_monthly['hi'].sel(month=[7,8,9]).mean('month').values)\n", "hi_JAS = np.ma.masked_where(aice_JAS < 15.0, hi_JAS)\n", "suptitle = ('JAS mean sea ice thickness, ' + str(args.start_date) + ' to ' + str(args.end_date))\n", "%matplotlib inline \n", "# SH\n", "fig = plt.figure(figsize=(10,8))\n", "ax = plt.gca(projection=ccrs.SouthPolarStereo())\n", "polarplot(hi_JAS, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,2), axis=ax)\n", "plt.suptitle(suptitle)\n", "\n", "# NH\n", "fig = plt.figure(figsize=(10,8))\n", "ax = plt.gca(projection=ccrs.NorthPolarStereo())\n", "polarplot(hi_JAS, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,4), axis=ax, proj='NP')\n", "plt.suptitle(suptitle)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### OND" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "aice_OND = np.ma.masked_invalid(ds_monthly['aice'].sel(month=[10,11,12]).mean('month').values) * 100.\n", "hi_OND = np.ma.masked_invalid(ds_monthly['hi'].sel(month=[10,11,12]).mean('month').values)\n", "hi_OND = np.ma.masked_where(aice_OND < 15.0, hi_OND)\n", "suptitle = ('OND mean sea ice thickness, ' + str(args.start_date) + ' to ' + str(args.end_date))\n", "%matplotlib inline \n", "# SH\n", "fig = plt.figure(figsize=(10,8))\n", "ax = plt.gca(projection=ccrs.SouthPolarStereo())\n", "polarplot(hi_OND, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,2), axis=ax)\n", "plt.suptitle(suptitle)\n", "\n", "# NH\n", "fig = plt.figure(figsize=(10,8))\n", "ax = plt.gca(projection=ccrs.NorthPolarStereo())\n", "polarplot(hi_OND, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,4), axis=ax, proj='NP')\n", "plt.suptitle(suptitle)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "(obs_SH_mean.sic * 100.).plot()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tmp = obs_SH_mean['sic'].values * 100.\n", "#tmp1 = np.ma.masked_where(tmp < 15.0, tmp)\n", "fig = plt.figure(figsize=(10,8))\n", "ax = plt.gca(projection=ccrs.SouthPolarStereo())\n", "ax.contour(obs_SH_mean.lon.values, obs_SH_mean.lat.values, tmp, levels=[16.5], transform=ccrs.PlateCarree(), colors='black')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sea ice concentration" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "aice_mean = np.ma.masked_where(aice_mean <= 15.0, aice_mean)\n", "\n", "suptitle = ('ANN mean ice concentration, ' + str(args.start_date) + ' to ' + str(args.end_date))\n", "%matplotlib inline\n", "# SH\n", "fig = plt.figure(figsize=(10,8))\n", "ax = plt.gca(projection=ccrs.SouthPolarStereo())\n", "polarplot(aice_mean, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,100), axis=ax)\n", "ax.contour(obs_SH_mean.lon.values, obs_SH_mean.lat.values, tmp, levels=[16.5], transform=ccrs.PlateCarree(), colors='black', \n", " linestyles='solid')\n", "plt.suptitle(suptitle)\n", "\n", "# NH\n", "fig = plt.figure(figsize=(10,8))\n", "ax = plt.gca(projection=ccrs.NorthPolarStereo())\n", "polarplot(aice_mean, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,100), axis=ax, proj='NP')\n", "plt.suptitle(suptitle)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Northern Hemisphere" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Time series \n", "\n", "### Maximum sea ice thickness " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from mom6_tools.latlon_analysis import create_xarray_dataset, plot_area_ave_stats\n", "from mom6_tools.m6plot import myStats" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dtime = seaice.time.values\n", "# variable name\n", "var = 'hi'\n", "# create datasets\n", "ds_sh = create_xarray_dataset(var,'m',dtime)\n", "ds_nh = create_xarray_dataset(var,'m',dtime)\n", "\n", "# loop in time\n", "for t in range(0,len(dtime)):\n", " # northern hemisphere\n", " tmp = np.ma.masked_invalid(seaice[var].sel(time=dtime[t]).values)\n", " tmp_nh = np.ma.masked_where(grd.geolat < 0, tmp)\n", " # get stats\n", " sMin, sMax, mean, std, rms = myStats(tmp_nh, area=None)\n", " # update Dataset\n", " ds_nh[var][1,t] = sMax\n", " # southern hemisphere\n", " tmp_sh = np.ma.masked_where(grd.geolat > 0, tmp)\n", " # get stats\n", " sMin, sMax, mean, std, rms = myStats(tmp_sh, area=None)\n", " # update Dataset\n", " ds_sh[var][1,t] = sMax" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Southern Hemisphere" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sh = ds_sh.sel(stats='max')\n", "fig = plt.figure(figsize=(12, 6))\n", "ax = fig.add_subplot(111)\n", "sh.hi.plot(ax=ax)\n", "plt.title('Maximum sea ice thickness, southern hemisphere', fontsize=14)\n", "plt.xlabel('Time [years]')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Northern Hemisphere" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nh = ds_nh.sel(stats='max')\n", "fig = plt.figure(figsize=(12, 6))\n", "ax = fig.add_subplot(111)\n", "nh.hi.plot(ax=ax)\n", "plt.title('Maximum sea ice thickness, northern hemisphere', fontsize=14)\n", "plt.xlabel('Time [years]')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sea ice extent" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dtime = seaice.time.values\n", "# variable name\n", "var = 'aice'\n", "# \n", "area_sh = []\n", "area_nh = []\n", "\n", "# loop in time\n", "for t in range(len(dtime)):\n", " tmp = np.ma.masked_invalid(seaice[var].sel(time=dtime[t]).values)\n", " # mask southern hemisphere\n", " tmp_nh = np.ma.masked_where(grd.geolat < 0, grd.area_t)\n", " # mask northern hemisphere\n", " tmp_sh = np.ma.masked_where(grd.geolat > 0, grd.area_t)\n", " # mask if aggregrated concentration <= 0.15\n", " tmp_nh = np.ma.masked_where(tmp <= 0.15, tmp_nh*tmp)\n", " tmp_sh = np.ma.masked_where(tmp <= 0.15, tmp_sh*tmp)\n", " # append data\n", " area_sh.append(np.sum(tmp_sh)*1.0e-12) # 1.0e6 km^2\n", " area_nh.append(np.sum(tmp_nh)*1.0e-12) # 1.0e6 km^2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Southern hemisphere" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(12, 6))\n", "plt.plot(dtime,area_sh)\n", "plt.title('Sea ice area, southern hemisphere', fontsize=14)\n", "plt.xlabel('Time [years]', fontsize=14); plt.ylabel('x 1.0e6 km$^2$', fontsize=14)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Northern hemisphere" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(12, 6))\n", "plt.plot(dtime,area_nh)\n", "plt.title('Sea ice area, northern hemisphere', fontsize=14)\n", "plt.xlabel('Time [years]', fontsize=14); plt.ylabel('x 1.0e6 km$^2$', fontsize=14)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Under development" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.path as mpath\n", "fig = plt.figure(figsize=[10, 8])\n", "ice = (seaice.aice[-1,:,:].data)\n", "ax1 = plt.subplot(1, 1, 1, projection=ccrs.NorthPolarStereo())\n", "ax1.set_extent([-180, 180, 50, 90], ccrs.PlateCarree())\n", "ax1.add_feature(cartopy.feature.LAND)\n", "ax1.gridlines()\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", "ax1.set_boundary(circle, transform=ax1.transAxes)\n", "#colormap = 'rainbow'\n", "cs = ax1.pcolormesh(grd.geolon,grd.geolat,ice,transform=ccrs.PlateCarree(),cmap='Blues_r', shading='flat')\n", "fig.colorbar(cs)\n", "# Add Land\n", "ax1.add_feature( cartopy.feature.LAND, zorder=1, edgecolor='none', facecolor='#fae5c9') #fae5c9')\n", "# add Ocean\n", "ax1.add_feature(cartopy.feature.OCEAN)\n", "# Add coastline\n", "ax1.coastlines(color='black')\n", "# Add lat lon rings\n", "ax1.gridlines(alpha='0.1',color='black')\n", "im1 = ax1.contour(grd.geolon,grd.geolat,ice,[15],colors='red', transform=ccrs.PlateCarree())\n", "plt.title('Model, Years ' + str(args.year_start) + ' to ' + str(args.year_end))" ] } ], "metadata": { "kernelspec": { "display_name": "Python (mom6-tools)", "language": "python", "name": "mom6-tools" }, "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.11.9" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "position": { "height": "588.85px", "left": "1073px", "right": "260px", "top": "172px", "width": "587px" }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }