From 814ea687356d0a00752424d84682aba9e385ff82 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Thu, 18 Jul 2024 14:55:02 -0600 Subject: [PATCH 01/11] Add glacier notebook to nblibrary This is the first notebook that runs under the key_metric example --- examples/key_metrics/config.yml | 53 +- .../nblibrary/glc/LIWG_SMB_diagnostic.ipynb | 818 ++++++++++++++++++ 2 files changed, 850 insertions(+), 21 deletions(-) create mode 100644 examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb diff --git a/examples/key_metrics/config.yml b/examples/key_metrics/config.yml index cc397e1..db3066a 100644 --- a/examples/key_metrics/config.yml +++ b/examples/key_metrics/config.yml @@ -115,19 +115,29 @@ compute_notebooks: # config_path: . # config_fil_str: "config_f.cam6_3_119.FLTHIST_ne30.r328_gamma0.33_soae.001.yaml" -# ocn: -# ocean_surface: + glc: + LIWG_SMB_diagnostic: + parameter_groups: + none: + case_name2: 'f.e23_alpha17f.FLTHIST_ne30.roughtopo.099' + CESM2_output_dir: '/glade/campaign/cesm/collections/collections/CESM2-LE/timeseries/glc/proc/tseries/year_1/smb' + CESM2_case_name: 'b.e21.BHISTsmbb.f09_g17.LE2-1251.019.cism.h.smb.1850-2015.nc' + path_obs: '/glade/u/home/gunterl/obs_diagnostic_cesm/' + obs_name: 'GrIS_MARv3.12_climo_1960_1999.nc' + last_year: 2005 + +# ice: +# seaice: # parameter_groups: # none: -# Case: b.e23_alpha16b.BLT1850.ne30_t232.054 -# savefigs: False -# mom6_tools_config: -# start_date: '0091-01-01' -# end_date: '0101-01-01' -# Fnames: -# native: 'mom6.h.native.????-??.nc' -# static: 'mom6.h.static.nc' -# oce_cat: /glade/u/home/gmarques/libs/oce-catalogs/reference-datasets.yml +# cases: +# - g.e23_a16g.GJRAv4.TL319_t232_hycom1_N75.2024.005 +# - g.e23_a16g.GJRAv4.TL319_t232_zstar_N65.2024.004 +# begyr1: 245 +# endyr1: 305 +# begyr2: 245 +# endyr2: 305 +# nyears: 25 # lnd: # land_comparison: @@ -140,18 +150,19 @@ compute_notebooks: # - 1850pAD # - 1850pSASU -# ice: -# seaice: +# ocn: +# ocean_surface: # parameter_groups: # none: -# cases: -# - g.e23_a16g.GJRAv4.TL319_t232_hycom1_N75.2024.005 -# - g.e23_a16g.GJRAv4.TL319_t232_zstar_N65.2024.004 -# begyr1: 245 -# endyr1: 305 -# begyr2: 245 -# endyr2: 305 -# nyears: 25 +# Case: b.e23_alpha16b.BLT1850.ne30_t232.054 +# savefigs: False +# mom6_tools_config: +# start_date: '0091-01-01' +# end_date: '0101-01-01' +# Fnames: +# native: 'mom6.h.native.????-??.nc' +# static: 'mom6.h.static.nc' +# oce_cat: /glade/u/home/gmarques/libs/oce-catalogs/reference-datasets.yml ########### JUPYTER BOOK CONFIG ########### diff --git a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb new file mode 100644 index 0000000..65fc2a0 --- /dev/null +++ b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb @@ -0,0 +1,818 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e45c723f-aa6f-4db4-a197-492132cc8156", + "metadata": {}, + "source": [ + "# Land ice SMB model comparison\n", + "This notebook compares the downscaled output of surface mass balance (SMB) over the Greenland ice sheet (GrIS) to the regional model MAR. In what follows, we interchangeably call the MAR data \"observation\".\n", + "\\\n", + "Note1: the MAR data are processed as a climatology spanning 1960-1999.\\\n", + "Note2: the MAR data are available at a uniform resolution of 1km using the same projection as the CISM grid. This notebook requires the interpolation of the MAR data on the CISM grid. The interpolation is done in this notebook (for now) to allow for the eventuality of the CISM grid or the MAR grid to change in the future. \\\n", + "creation: 05-26-24 \\\n", + "contact: Gunter Leguy (gunterl@ucar.edu)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "81bedf90-288c-4cfa-add5-b199ca9bcf72", + "metadata": {}, + "outputs": [], + "source": [ + "# Import packages\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.cm as mcm\n", + "from netCDF4 import Dataset\n", + "import os\n", + "from scipy.interpolate import RegularGridInterpolator\n", + "\n", + "# to display figures in notebook after executing the code.\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e02f1300-e0e1-448c-b534-68146555d660", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "parameters" + ] + }, + "outputs": [], + "source": [ + "# Parameter Defaults\n", + "\n", + "CESM_output_dir = \"/glade/campaign/cesm/development/cross-wg/diagnostic_framework/CESM_output_for_testing\"\n", + "CESM2_output_dir = \"/glade/campaign/cesm/collections/collections/CESM2-LE/timeseries/glc/proc/tseries/year_1/smb\"\n", + "\n", + "case_name = \"b.e23_alpha17f.BLT1850.ne30_t232.092\" # case name\n", + "case_name2 = \"f.e23_alpha17f.FLTHIST_ne30.roughtopo.099\" # case name\n", + "CESM2_case_name = \"b.e21.BHISTsmbb.f09_g17.LE2-1251.019.cism.h.smb.1850-2015.nc\"\n", + "path_obs = \"/glade/u/home/gunterl/obs_diagnostic_cesm\" # path to observed dataset\n", + "obs_name = \"GrIS_MARv3.12_climo_1960_1999.nc\"\n", + "last_year = 2005" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cef60ddb-9ff4-4a14-a8ea-0d5740b6c18a", + "metadata": {}, + "outputs": [], + "source": [ + "path_glc = f\"{CESM_output_dir}/{case_name}/glc/hist\" # path to glc output\n", + "file_glc = f\"{path_glc}/{case_name}.cism.gris.initial_hist.0001-01-01-00000.nc\" # name of glc file output\n", + "\n", + "path_cpl = f\"{CESM_output_dir}/{case_name2}/cpl/hist\" # path to cpl output\n", + "file_cpl = f\"{path_cpl}/{case_name2}.cpl.hx.1yr2glc.1996-01-01-00000.nc\" # name of last cpl simulation output\n", + "\n", + "file_obs = f\"{path_obs}/{obs_name}\" # name of observed dataset file\n", + "file_cesm2 = f\"{CESM2_output_dir}/{CESM2_case_name}\" # name of CESM2 dataset" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "52ba1301-0c4c-481c-a57b-a9614160807b", + "metadata": {}, + "outputs": [], + "source": [ + "## CISM grid information and loading a field used for filtering\n", + "\n", + "# Reading the information we need from the glc file\n", + "nid = Dataset(file_glc, \"r\")\n", + "x_cism = nid.variables[\"x1\"][:]\n", + "y_cism = nid.variables[\"y1\"][:]\n", + "thk_cism = np.squeeze(nid.variables[\"thk\"][0, :, :])\n", + "nid.close()\n", + "\n", + "# Defining the grid dimensions\n", + "## For the CISM grid\n", + "nx_cism = len(x_cism)\n", + "ny_cism = len(y_cism)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "69891e0a-0cc6-48bd-a9bf-b61e19da13b4", + "metadata": {}, + "outputs": [], + "source": [ + "## The observed dataset\n", + "nid = Dataset(file_obs, \"r\")\n", + "x_obs = nid.variables[\"x\"][:]\n", + "y_obs = nid.variables[\"y\"][:]\n", + "smb_obs_src = np.squeeze(nid.variables[\"SMB\"][0, :, :])\n", + "nid.close()\n", + "\n", + "## For the observed grid\n", + "nx_obs = len(x_obs)\n", + "ny_obs = len(y_obs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05eee529-d42f-4872-8cf2-f484ca44bf3f", + "metadata": {}, + "outputs": [], + "source": [ + "# Constants\n", + "ny_climo = 10 # number of years to compute the climatology\n", + "res = np.abs(x_cism[1] - x_cism[0]) # CISM output resolution\n", + "\n", + "rhoi = 917 # ice density kg/m3\n", + "rhow = 1000 # water density kg/m3\n", + "sec_in_yr = 60 * 60 * 24 * 365 # seconds in a year\n", + "\n", + "smb_convert = sec_in_yr / rhoi * 1000 # converting kg m-2 s-1 ice to mm y-1 w.e.\n", + "kg_to_Gt = 1e-12 # Converting kg to Gt\n", + "mm_to_Gt = rhow * 1e-3 * res**2 * kg_to_Gt # converting mm/yr to Gt/yr" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4cf8782d-1500-4707-bd3b-88aa2e2c9120", + "metadata": {}, + "outputs": [], + "source": [ + "# Functions used in this notebook\n", + "\n", + "\n", + "def set_plot_prop_clean(ax):\n", + " \"\"\"\n", + " This function cleans up the figures from unnecessary default figure properties.\n", + "\n", + " \"\"\"\n", + " ax.invert_yaxis()\n", + " ax.set_xlabel(\"\")\n", + " ax.set_ylabel(\"\")\n", + " ax.set_xticklabels(\"\")\n", + " ax.set_yticklabels(\"\")\n", + " ax.set_xticks([])\n", + " ax.set_yticks([])\n", + "\n", + "\n", + "def rmse(prediction, target):\n", + " \"\"\"\n", + " This function returns the root mean square error for the SMB.\n", + " Input:\n", + " prediction = field to predict\n", + " target = field to compare with the prediction\n", + " \"\"\"\n", + " return np.sqrt(((prediction - target) ** 2).mean())\n", + "\n", + "\n", + "def net_avrg(data):\n", + " \"\"\"\n", + " This function returns the net average of a data field\n", + " \"\"\"\n", + " return np.sum(np.sum(data, axis=0), axis=0)\n", + "\n", + "\n", + "def read_smb(file):\n", + " \"\"\"\n", + " This function reads the CISM SMB dataset from a CESM simulation output\n", + " in the cpl directory. The output is adjusted to be converted to mm/yr w.e unit.\n", + "\n", + " Input:\n", + " file: name of the file to extract the SMB\n", + " \"\"\"\n", + " nid = Dataset(file, \"r\")\n", + " smb_cism = np.squeeze(nid.variables[\"glc1Exp_Flgl_qice\"][0, :, :]) * smb_convert\n", + " nid.close()\n", + " return smb_cism" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fb320812-c7aa-408b-af7c-b252bcb1acb9", + "metadata": {}, + "outputs": [], + "source": [ + "# Loading the data\n", + "nid = Dataset(file_cpl, \"r\")\n", + "smb_cism = np.squeeze(nid.variables[\"glc1Exp_Flgl_qice\"][0, :, :]) * smb_convert\n", + "nid.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bb81be67-98d6-4924-a90e-930d9b2caed8", + "metadata": {}, + "outputs": [], + "source": [ + "# creating the SMB climatology for comparison with observation\n", + "\n", + "# Initializing a field for the climatology\n", + "smb_cism_climo = np.zeros((ny_cism, nx_cism))\n", + "\n", + "# Counter for available year (only needed if the number of years available is smaller\n", + "# than the number of years requested to create the climatology.\n", + "count_yr = 0\n", + "\n", + "for k in range(ny_climo):\n", + "\n", + " year_to_read = last_year - k\n", + " str_yr = str(year_to_read).zfill(4)\n", + " file_cpl = f\"{path_cpl}/{case_name2}.cpl.hx.1yr2glc.{str_yr}-01-01-00000.nc\"\n", + "\n", + " if not os.path.isfile(file_cpl):\n", + " print(\"The couple file for time\", str_yr, \"does not exist.\")\n", + " print(\n", + " \"We will only use the files that existed until now to create the SMB climatology.\"\n", + " )\n", + " break\n", + "\n", + " smb_cism_climo = smb_cism_climo + read_smb(file_cpl)\n", + " count_yr = count_yr + 1\n", + "\n", + "# Averaging the climo data\n", + "smb_cism_climo = smb_cism_climo / count_yr\n", + "\n", + "print(\"number of years used in climatology = \", count_yr)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ac96bb16-7bd8-4d7b-b00b-d315feeb1a5d", + "metadata": {}, + "outputs": [], + "source": [ + "# Interpolating the observed data onto the CISM grid\n", + "\n", + "# Defining the interpolation functions\n", + "myInterpFunction_smb_obs = RegularGridInterpolator(\n", + " (x_obs, y_obs),\n", + " smb_obs_src.transpose(),\n", + " method=\"linear\",\n", + " bounds_error=False,\n", + " fill_value=None,\n", + ")\n", + "\n", + "# Initializing the glacier ID variable\n", + "smb_obs_cism = np.zeros((ny_cism, nx_cism))\n", + "\n", + "# Performing the interpolation\n", + "for j in range(ny_cism):\n", + " point_y = np.zeros(nx_cism)\n", + " point_y[:] = y_cism[j]\n", + " pts = (x_cism[:], point_y[:])\n", + " smb_obs_cism[j, :] = myInterpFunction_smb_obs(pts)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "288f34cf-808f-4de1-8372-276fdbcfdc46", + "metadata": {}, + "outputs": [], + "source": [ + "# Filtering out fill values\n", + "smb_obs_cism = np.where(smb_obs_cism > 1e20, 0, smb_obs_cism)\n", + "mask = thk_cism[:, :] == 0\n", + "smb_obs_cism = np.where(mask, 0, smb_obs_cism)\n", + "smb_cism = np.where(mask, 0, smb_cism)\n", + "smb_cism_climo = np.where(mask, 0, smb_cism_climo)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c7973cfe-64e0-47d4-a1b6-73cd9e62fdb2", + "metadata": {}, + "outputs": [], + "source": [ + "# Comparing the SMB\n", + "\n", + "# Computing the SMB difference\n", + "smb_diff_cesm3 = smb_cism_climo - smb_obs_cism\n", + "\n", + "# Computing the net averages SMB\n", + "avg_smb_cism_climo = np.round(net_avrg(smb_cism_climo) * mm_to_Gt, 2)\n", + "avg_smb_obs_cism = np.round(net_avrg(smb_obs_cism) * mm_to_Gt, 2)\n", + "avg_smb_diff_cesm3 = np.round(net_avrg(smb_diff_cesm3) * mm_to_Gt, 2)\n", + "\n", + "# Colormap choice\n", + "my_cmap = mcm.get_cmap(\"Spectral\")\n", + "# my_cmap = mcm.get_cmap('RdYlBu_r')\n", + "my_cmap_diff = mcm.get_cmap(\"bwr_r\")\n", + "\n", + "\n", + "# Colorbar bounds\n", + "vmin = -2000\n", + "vmax = 2000\n", + "\n", + "# Figure\n", + "fig, ax = plt.subplots(1, 3, sharey=True, figsize=[22, 9])\n", + "\n", + "## Left panel\n", + "last_panel0 = ax[0].imshow(smb_cism_climo[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap)\n", + "ax[0].set_title(\"SMB CESM3 (mm/y w.e.)\", fontsize=16)\n", + "set_plot_prop_clean(ax[0])\n", + "ax[0].annotate(\"net avg =\" + str(avg_smb_cism_climo) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", + "\n", + "pos = ax[0].get_position()\n", + "cax = fig.add_axes([0.35, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "\n", + "cbar = fig.colorbar(last_panel0, cax=cax)\n", + "cbar.ax.tick_params(labelsize=16)\n", + "\n", + "\n", + "## Center panel\n", + "last_panel1 = ax[1].imshow(smb_obs_cism[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap)\n", + "\n", + "ax[1].set_title(\"SMB Obs (mm/y w.e.)\", fontsize=16)\n", + "set_plot_prop_clean(ax[1])\n", + "ax[1].annotate(\"net avg =\" + str(avg_smb_obs_cism) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", + "\n", + "\n", + "pos = ax[1].get_position()\n", + "cax = fig.add_axes([0.62, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "\n", + "cbar = fig.colorbar(last_panel1, cax=cax)\n", + "cbar.ax.tick_params(labelsize=16)\n", + "\n", + "\n", + "## Right panel\n", + "\n", + "vmin = -2000\n", + "vmax = 2000\n", + "\n", + "last_panel2 = ax[2].imshow(smb_diff_cesm3, vmin=vmin, vmax=vmax, cmap=my_cmap_diff)\n", + "\n", + "ax[2].set_title(\"SMB diff (CESM3 - Obs, mm/yr w.e.)\", fontsize=16)\n", + "set_plot_prop_clean(ax[2])\n", + "\n", + "ax[2].annotate(\"net avg =\" + str(avg_smb_diff_cesm3) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", + "\n", + "pos = ax[2].get_position()\n", + "cax = fig.add_axes([0.89, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "\n", + "cbar = fig.colorbar(last_panel2, cax=cax)\n", + "cbar.ax.tick_params(labelsize=16)\n", + "\n", + "fig.subplots_adjust(right=0.89)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b7bbb2a4-c7ce-4327-9610-2f6ae0d7f3ce", + "metadata": {}, + "outputs": [], + "source": [ + "# Integrated SMB time series\n", + "\n", + "# Initializing a field for the climatology\n", + "avg_smb_cism_timeseries = np.zeros(last_year)\n", + "\n", + "# Counter for available year (only needed if the number of years available is smaller\n", + "# than the number of years requested to create the climatology.\n", + "count_yr = 0\n", + "\n", + "for k in range(last_year):\n", + "\n", + " year_to_read = last_year - k\n", + " str_yr = str(year_to_read).zfill(4)\n", + " file_cpl = f\"{path_cpl}/{case_name2}.cpl.hx.1yr2glc.\" + str_yr + \"-01-01-00000.nc\"\n", + "\n", + " if not os.path.isfile(file_cpl):\n", + " print(\"The couple file for time\", str_yr, \"does not exist.\")\n", + " print(\n", + " \"We will only use the files that existed until now to create the time series.\"\n", + " )\n", + " break\n", + "\n", + " smb_temp = read_smb(file_cpl)\n", + " smb_temp = np.where(mask, 0, smb_temp)\n", + "\n", + " avg_smb_cism_timeseries[year_to_read - 1] = np.round(\n", + " net_avrg(smb_temp) * mm_to_Gt, 2\n", + " )\n", + " count_yr = count_yr + 1\n", + "\n", + " del smb_temp\n", + "\n", + "first_year = year_to_read\n", + "\n", + "print(\"number of years used in climatology = \", count_yr)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "01af1cd1-0351-452c-99e7-125546469f69", + "metadata": {}, + "outputs": [], + "source": [ + "# Plotting the SMB spatially averaged time series\n", + "\n", + "\n", + "time = np.arange(first_year, last_year)\n", + "nt = len(time)\n", + "\n", + "\n", + "avg_smb_climo_timeseries = np.zeros(nt)\n", + "avg_smb_obs_cism_timeseries = np.zeros(nt)\n", + "\n", + "avg_smb_climo_timeseries[:] = avg_smb_cism_climo\n", + "avg_smb_obs_cism_timeseries[:] = avg_smb_obs_cism\n", + "\n", + "\n", + "x_ticks = np.arange(first_year, last_year, 10)\n", + "tickx = x_ticks\n", + "\n", + "ymin = 150\n", + "ymax = 650\n", + "y_step = 50\n", + "y_ticks = np.arange(ymin, ymax + y_step, y_step)\n", + "\n", + "\n", + "plt.figure(figsize=(16, 7))\n", + "\n", + "line = \"-\"\n", + "color = \"blue\"\n", + "sizefont = 16\n", + "linewidth = 2\n", + "\n", + "\n", + "# Plotting NESM3 experiments\n", + "plt.subplot(111)\n", + "\n", + "label = \"CISM\"\n", + "plt.plot(\n", + " time,\n", + " avg_smb_cism_timeseries[first_year::],\n", + " line,\n", + " ms=3,\n", + " mfc=color,\n", + " color=color,\n", + " label=label,\n", + " linewidth=linewidth,\n", + ")\n", + "\n", + "label = \"CISM climatology\"\n", + "color = \"k\"\n", + "plt.plot(\n", + " time,\n", + " avg_smb_climo_timeseries[:],\n", + " line,\n", + " ms=3,\n", + " mfc=color,\n", + " color=color,\n", + " label=label,\n", + " linewidth=linewidth,\n", + ")\n", + "\n", + "label = \"Observation\"\n", + "color = \"r\"\n", + "plt.plot(\n", + " time,\n", + " avg_smb_obs_cism_timeseries[:],\n", + " line,\n", + " ms=3,\n", + " mfc=color,\n", + " color=color,\n", + " label=label,\n", + " linewidth=linewidth,\n", + ")\n", + "\n", + "\n", + "plt.xlim([first_year, last_year])\n", + "plt.xticks(x_ticks, tickx, fontsize=sizefont)\n", + "plt.xlabel(r\"$Time$ (y)\", fontsize=sizefont)\n", + "plt.ylabel(\"SMB average evolution (Gt/yr)\", multialignment=\"center\", fontsize=sizefont)\n", + "plt.ylim([ymin, ymax])\n", + "plt.yticks(fontsize=sizefont)\n", + "plt.legend(loc=\"upper left\", ncol=1, frameon=True, borderaxespad=0)\n", + "\n", + "plt.title(\"SMB average evolution\", fontsize=sizefont)" + ] + }, + { + "cell_type": "markdown", + "id": "04c05f38-6115-4b54-8965-718baa1448a5", + "metadata": {}, + "source": [ + "# Adding comparison with CESM2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2093fc4-042c-4670-9537-175e3df8dffb", + "metadata": {}, + "outputs": [], + "source": [ + "# Defining the paths and files for CESM2\n", + "\n", + "\n", + "nid = Dataset(file_cesm2, \"r\")\n", + "smb_cesm2 = np.squeeze(nid.variables[\"smb\"][0, :, :])\n", + "nid.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a8fa10a5-ab3f-467b-8614-bdf6cb702fd8", + "metadata": {}, + "outputs": [], + "source": [ + "smb_cesm2 = np.where(mask, 0, smb_cesm2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f43f1382-4927-4de7-bd04-3b59bb225db7", + "metadata": {}, + "outputs": [], + "source": [ + "# Comparing the SMB\n", + "\n", + "# Computing the SMB difference\n", + "# smb_diff = smb_cism_climo - smb_obs_cism\n", + "\n", + "# Computing the net averages SMB\n", + "avg_smb_cesm2 = np.round(net_avrg(smb_cesm2) * mm_to_Gt, 2)\n", + "\n", + "\n", + "# Colormap choice\n", + "my_cmap = mcm.get_cmap(\"Spectral\")\n", + "# my_cmap = mcm.get_cmap('RdYlBu_r')\n", + "my_cmap_diff = mcm.get_cmap(\"bwr_r\")\n", + "\n", + "\n", + "# Colorbar bounds\n", + "vmin = -2000\n", + "vmax = 2000\n", + "\n", + "# Figure\n", + "fig, ax = plt.subplots(1, 3, sharey=True, figsize=[22, 9])\n", + "\n", + "## Left panel\n", + "last_panel0 = ax[0].imshow(smb_cesm2[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap)\n", + "ax[0].set_title(\"SMB CESM2 (mm/y w.e.)\", fontsize=16)\n", + "set_plot_prop_clean(ax[0])\n", + "ax[0].annotate(\"net avg =\" + str(avg_smb_cesm2) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", + "\n", + "pos = ax[0].get_position()\n", + "cax = fig.add_axes([0.35, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "\n", + "cbar = fig.colorbar(last_panel0, cax=cax)\n", + "cbar.ax.tick_params(labelsize=16)\n", + "\n", + "\n", + "## Center panel\n", + "last_panel1 = ax[1].imshow(smb_cism_climo[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap)\n", + "\n", + "ax[1].set_title(\"SMB CESM3 (mm/y w.e.)\", fontsize=16)\n", + "set_plot_prop_clean(ax[1])\n", + "ax[1].annotate(\"net avg =\" + str(avg_smb_cism_climo) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", + "\n", + "\n", + "pos = ax[1].get_position()\n", + "cax = fig.add_axes([0.62, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "\n", + "cbar = fig.colorbar(last_panel1, cax=cax)\n", + "cbar.ax.tick_params(labelsize=16)\n", + "\n", + "\n", + "## Right panel\n", + "\n", + "last_panel2 = ax[2].imshow(smb_obs_cism[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap)\n", + "\n", + "ax[2].set_title(\"SMB Obs (mm/y w.e.)\", fontsize=16)\n", + "set_plot_prop_clean(ax[2])\n", + "ax[2].annotate(\"net avg =\" + str(avg_smb_obs_cism) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", + "\n", + "\n", + "pos = ax[2].get_position()\n", + "cax = fig.add_axes([0.62, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "\n", + "cbar = fig.colorbar(last_panel1, cax=cax)\n", + "cbar.ax.tick_params(labelsize=16)\n", + "\n", + "\n", + "# last_panel2 = ax[2].imshow(smb_diff_climo, vmin=vmin, vmax=vmax, cmap=my_cmap)\n", + "\n", + "# ax[2].set_title('SMB diff (CISM - Obs, mm/yr w.e.)',fontsize=16)\n", + "# set_plot_prop_clean(ax[2])\n", + "\n", + "# ax[2].annotate('net avg =' + str(avg_smb_diff) + ' Gt/yr',xy=(5,5),fontsize=16)\n", + "\n", + "# pos = ax[2].get_position()\n", + "# cax = fig.add_axes([0.89, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "\n", + "# cbar = fig.colorbar(last_panel2,cax=cax)\n", + "# cbar.ax.tick_params(labelsize=16)\n", + "\n", + "fig.subplots_adjust(right=0.89)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "85d4a99e-73b1-46eb-b9e2-fad93c155c96", + "metadata": {}, + "outputs": [], + "source": [ + "# Computing the difference w.r.t MAR obs\n", + "smb_diff_cesm2 = smb_cesm2 - smb_obs_cism\n", + "\n", + "# Computing the net averages SMB\n", + "avg_smb_cesm2 = np.round(net_avrg(smb_cesm2) * mm_to_Gt, 2)\n", + "avg_smb_diff_cesm2 = np.round(net_avrg(smb_diff_cesm2) * mm_to_Gt, 2)\n", + "\n", + "# Masking out ocean cells\n", + "mask2 = smb_cesm2[:, :] == 0\n", + "smb_cesm2_plot = np.where(mask2, np.float64(\"Nan\"), smb_cesm2)\n", + "smb_obs_cism_plot = np.where(mask2, np.float64(\"Nan\"), smb_obs_cism)\n", + "smb_cism_climo_plot = np.where(mask2, np.float64(\"Nan\"), smb_cism_climo)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "86e9b67b-ac4e-4551-96ea-ffff814fa038", + "metadata": {}, + "outputs": [], + "source": [ + "# Comparing the SMB\n", + "\n", + "# Computing the SMB difference\n", + "# smb_diff = smb_cism_climo - smb_obs_cism\n", + "\n", + "\n", + "# Colormap choice\n", + "my_cmap = mcm.get_cmap(\"Spectral\")\n", + "# my_cmap = mcm.get_cmap('RdYlBu_r')\n", + "my_cmap_diff = mcm.get_cmap(\"bwr_r\")\n", + "\n", + "\n", + "# Colorbar bounds\n", + "vmin = -2000\n", + "vmax = 2000\n", + "\n", + "# Figure\n", + "fig, ax = plt.subplots(1, 3, sharey=True, figsize=[16, 9])\n", + "\n", + "## Left panel\n", + "last_panel0 = ax[0].imshow(smb_cesm2_plot[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap)\n", + "# ax[0].set_title('SMB CESM2 (mm/y w.e.)',fontsize=16)\n", + "ax[0].set_title(\"CESM2\", fontsize=16)\n", + "set_plot_prop_clean(ax[0])\n", + "ax[0].annotate(\"net avg =\" + str(avg_smb_cesm2) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", + "ax[0].set_ylabel(\"SMB (mm/y w.e)\", fontsize=16)\n", + "# pos = ax[0].get_position()\n", + "# cax = fig.add_axes([0.35, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "\n", + "# cbar = fig.colorbar(last_panel0, cax=cax)\n", + "# cbar.ax.tick_params(labelsize=16)\n", + "\n", + "\n", + "## Center panel\n", + "last_panel1 = ax[1].imshow(\n", + " smb_cism_climo_plot[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap\n", + ")\n", + "\n", + "ax[1].set_title(\"CESM3\", fontsize=16)\n", + "set_plot_prop_clean(ax[1])\n", + "ax[1].annotate(\"net avg =\" + str(avg_smb_cism_climo) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", + "\n", + "\n", + "# pos = ax[1].get_position()\n", + "# cax = fig.add_axes([0.62, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "\n", + "# cbar = fig.colorbar(last_panel1,cax=cax)\n", + "# cbar.ax.tick_params(labelsize=16)\n", + "\n", + "\n", + "## Right panel\n", + "\n", + "last_panel2 = ax[2].imshow(smb_obs_cism_plot[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap)\n", + "\n", + "ax[2].set_title(\"MARv3.12 1960-1999\", fontsize=16)\n", + "set_plot_prop_clean(ax[2])\n", + "ax[2].annotate(\"net avg =\" + str(avg_smb_obs_cism) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", + "\n", + "\n", + "pos = ax[2].get_position()\n", + "cax = fig.add_axes([0.9, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "\n", + "cbar = fig.colorbar(last_panel1, cax=cax)\n", + "cbar.ax.tick_params(labelsize=16)\n", + "\n", + "\n", + "# last_panel2 = ax[2].imshow(smb_diff_climo, vmin=vmin, vmax=vmax, cmap=my_cmap)\n", + "\n", + "# ax[2].set_title('SMB diff (CISM - Obs, mm/yr w.e.)',fontsize=16)\n", + "# set_plot_prop_clean(ax[2])\n", + "\n", + "# ax[2].annotate('net avg =' + str(avg_smb_diff) + ' Gt/yr',xy=(5,5),fontsize=16)\n", + "\n", + "# pos = ax[2].get_position()\n", + "# cax = fig.add_axes([0.89, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "\n", + "# cbar = fig.colorbar(last_panel2,cax=cax)\n", + "# cbar.ax.tick_params(labelsize=16)\n", + "\n", + "fig.subplots_adjust(right=0.89)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "307999aa-83e6-4dc4-b131-af7261d9b663", + "metadata": {}, + "outputs": [], + "source": [ + "# Comparing the SMB\n", + "\n", + "# Colormap choice\n", + "my_cmap = mcm.get_cmap(\"Spectral\")\n", + "# my_cmap = mcm.get_cmap('RdYlBu_r')\n", + "my_cmap_diff = mcm.get_cmap(\"bwr_r\")\n", + "\n", + "\n", + "# Colorbar bounds\n", + "vmin = -2000\n", + "vmax = 2000\n", + "\n", + "# Figure\n", + "fig, ax = plt.subplots(1, 2, sharey=True, figsize=[12, 9])\n", + "\n", + "## Left panel\n", + "last_panel0 = ax[0].imshow(\n", + " smb_diff_cesm2[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap_diff\n", + ")\n", + "ax[0].set_title(\"CESM2 - MAR\", fontsize=16)\n", + "set_plot_prop_clean(ax[0])\n", + "ax[0].annotate(\"net avg =\" + str(avg_smb_diff_cesm2) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", + "ax[0].set_ylabel(\"SMB difference (mm/yr w.e.)\", fontsize=16)\n", + "\n", + "# pos = ax[0].get_position()\n", + "# cax = fig.add_axes([0.465, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "\n", + "# cbar = fig.colorbar(last_panel0, cax=cax)\n", + "# cbar.ax.tick_params(labelsize=16)\n", + "\n", + "\n", + "## Right panel\n", + "last_panel1 = ax[1].imshow(\n", + " smb_diff_cesm3[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap_diff\n", + ")\n", + "\n", + "ax[1].set_title(\"CESM3 - MAR\", fontsize=16)\n", + "set_plot_prop_clean(ax[1])\n", + "ax[1].annotate(\"net avg =\" + str(avg_smb_diff_cesm3) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", + "\n", + "pos = ax[1].get_position()\n", + "cax = fig.add_axes([0.89, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "\n", + "cbar = fig.colorbar(last_panel1, cax=cax)\n", + "cbar.ax.tick_params(labelsize=16)\n", + "\n", + "\n", + "# fig.tight_layout()\n", + "fig.subplots_adjust(right=0.89)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:cupid-analysis]", + "language": "python", + "name": "conda-env-cupid-analysis-py" + }, + "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.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 4917465b2364af4f1f4cdbb0361cade2fa8613bb Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Fri, 19 Jul 2024 17:42:33 -0600 Subject: [PATCH 02/11] LIWG notebook cleanup Removed the CESM3 vs CESM2 comparison section, refactored code to avoid duplication. Still to do: create utils.py so additional notebooks can share these refactored functions. Probably some other code clean-up as well... --- examples/key_metrics/config.yml | 8 +- .../nblibrary/glc/LIWG_SMB_diagnostic.ipynb | 748 +++++++----------- 2 files changed, 293 insertions(+), 463 deletions(-) diff --git a/examples/key_metrics/config.yml b/examples/key_metrics/config.yml index db3066a..2b3e590 100644 --- a/examples/key_metrics/config.yml +++ b/examples/key_metrics/config.yml @@ -119,12 +119,10 @@ compute_notebooks: LIWG_SMB_diagnostic: parameter_groups: none: - case_name2: 'f.e23_alpha17f.FLTHIST_ne30.roughtopo.099' - CESM2_output_dir: '/glade/campaign/cesm/collections/collections/CESM2-LE/timeseries/glc/proc/tseries/year_1/smb' - CESM2_case_name: 'b.e21.BHISTsmbb.f09_g17.LE2-1251.019.cism.h.smb.1850-2015.nc' - path_obs: '/glade/u/home/gunterl/obs_diagnostic_cesm/' + obs_path: '/glade/u/home/gunterl/obs_diagnostic_cesm/' obs_name: 'GrIS_MARv3.12_climo_1960_1999.nc' - last_year: 2005 + climo_nyears: 40 + last_year: 101 # ice: # seaice: diff --git a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb index 65fc2a0..5420bef 100644 --- a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb +++ b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb @@ -51,14 +51,16 @@ "# Parameter Defaults\n", "\n", "CESM_output_dir = \"/glade/campaign/cesm/development/cross-wg/diagnostic_framework/CESM_output_for_testing\"\n", - "CESM2_output_dir = \"/glade/campaign/cesm/collections/collections/CESM2-LE/timeseries/glc/proc/tseries/year_1/smb\"\n", - "\n", "case_name = \"b.e23_alpha17f.BLT1850.ne30_t232.092\" # case name\n", - "case_name2 = \"f.e23_alpha17f.FLTHIST_ne30.roughtopo.099\" # case name\n", - "CESM2_case_name = \"b.e21.BHISTsmbb.f09_g17.LE2-1251.019.cism.h.smb.1850-2015.nc\"\n", - "path_obs = \"/glade/u/home/gunterl/obs_diagnostic_cesm\" # path to observed dataset\n", - "obs_name = \"GrIS_MARv3.12_climo_1960_1999.nc\"\n", - "last_year = 2005" + "climo_nyears = 40 # number of years to compute the climatology\n", + "last_year = 101\n", + "\n", + "base_case_output_dir = CESM_output_dir\n", + "base_case_name = None\n", + "base_last_year = last_year\n", + "\n", + "obs_path = \"/glade/u/home/gunterl/obs_diagnostic_cesm\" # path to observed dataset\n", + "obs_name = \"GrIS_MARv3.12_climo_1960_1999.nc\"" ] }, { @@ -68,14 +70,17 @@ "metadata": {}, "outputs": [], "source": [ - "path_glc = f\"{CESM_output_dir}/{case_name}/glc/hist\" # path to glc output\n", - "file_glc = f\"{path_glc}/{case_name}.cism.gris.initial_hist.0001-01-01-00000.nc\" # name of glc file output\n", + "case_init_file = f\"{CESM_output_dir}/{case_name}/glc/hist/{case_name}.cism.gris.initial_hist.0001-01-01-00000.nc\" # name of glc file output\n", "\n", - "path_cpl = f\"{CESM_output_dir}/{case_name2}/cpl/hist\" # path to cpl output\n", - "file_cpl = f\"{path_cpl}/{case_name2}.cpl.hx.1yr2glc.1996-01-01-00000.nc\" # name of last cpl simulation output\n", + "case_path = f\"{CESM_output_dir}/{case_name}/cpl/hist\" # path to glc output\n", + "case_file = f\"{case_path}/{case_name}.cpl.hx.1yr2glc.{last_year:04d}-01-01-00000.nc\" # name of glc file output\n", + "obs_file = f\"{obs_path}/{obs_name}\" # name of observed dataset file\n", "\n", - "file_obs = f\"{path_obs}/{obs_name}\" # name of observed dataset file\n", - "file_cesm2 = f\"{CESM2_output_dir}/{CESM2_case_name}\" # name of CESM2 dataset" + "if base_case_name:\n", + " base_case_path = (\n", + " f\"{base_case_output_dir}/{base_case_name}/cpl/hist\" # path to cpl output\n", + " )\n", + " base_file = f\"{base_case_path}/{base_case_name}.cpl.hx.1yr2glc.{base_last_year:04d}-01-01-00000.nc\" # name of last cpl simulation output" ] }, { @@ -88,7 +93,7 @@ "## CISM grid information and loading a field used for filtering\n", "\n", "# Reading the information we need from the glc file\n", - "nid = Dataset(file_glc, \"r\")\n", + "nid = Dataset(case_init_file, \"r\")\n", "x_cism = nid.variables[\"x1\"][:]\n", "y_cism = nid.variables[\"y1\"][:]\n", "thk_cism = np.squeeze(nid.variables[\"thk\"][0, :, :])\n", @@ -108,7 +113,7 @@ "outputs": [], "source": [ "## The observed dataset\n", - "nid = Dataset(file_obs, \"r\")\n", + "nid = Dataset(obs_file, \"r\")\n", "x_obs = nid.variables[\"x\"][:]\n", "y_obs = nid.variables[\"y\"][:]\n", "smb_obs_src = np.squeeze(nid.variables[\"SMB\"][0, :, :])\n", @@ -127,7 +132,6 @@ "outputs": [], "source": [ "# Constants\n", - "ny_climo = 10 # number of years to compute the climatology\n", "res = np.abs(x_cism[1] - x_cism[0]) # CISM output resolution\n", "\n", "rhoi = 917 # ice density kg/m3\n", @@ -147,8 +151,6 @@ "outputs": [], "source": [ "# Functions used in this notebook\n", - "\n", - "\n", "def set_plot_prop_clean(ax):\n", " \"\"\"\n", " This function cleans up the figures from unnecessary default figure properties.\n", @@ -194,6 +196,51 @@ " return smb_cism" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "44eda4c4-3c23-450d-83d8-c76204cefdc4", + "metadata": {}, + "outputs": [], + "source": [ + "# More functions used in this notebook\n", + "params = {\n", + " \"ny_cism\": ny_cism,\n", + " \"nx_cism\": nx_cism,\n", + " \"climo_nyears\": climo_nyears,\n", + "}\n", + "\n", + "\n", + "def create_climo(path, case_name, last_year, params):\n", + " # Initializing a field for the climatology\n", + " climo_out = np.zeros((params[\"ny_cism\"], params[\"nx_cism\"]))\n", + "\n", + " # Counter for available year (only needed if the number of years available is smaller\n", + " # than the number of years requested to create the climatology.\n", + " count_yr = 0\n", + "\n", + " for k in range(params[\"climo_nyears\"]):\n", + "\n", + " year_to_read = last_year - k\n", + " filename = (\n", + " f\"{path}/{case_name}.cpl.hx.1yr2glc.{year_to_read:04d}-01-01-00000.nc\"\n", + " )\n", + "\n", + " if not os.path.isfile(filename):\n", + " print(f\"The couple file for time {year_to_read} does not exist.\")\n", + " print(\n", + " \"We will only use the files that existed until now to create the SMB climatology.\"\n", + " )\n", + " break\n", + "\n", + " climo_out = climo_out + read_smb(filename)\n", + " count_yr = count_yr + 1\n", + "\n", + " print(\"number of years used in climatology = \", count_yr)\n", + " # Averaging the climo data\n", + " return climo_out / count_yr" + ] + }, { "cell_type": "code", "execution_count": null, @@ -202,9 +249,10 @@ "outputs": [], "source": [ "# Loading the data\n", - "nid = Dataset(file_cpl, \"r\")\n", - "smb_cism = np.squeeze(nid.variables[\"glc1Exp_Flgl_qice\"][0, :, :]) * smb_convert\n", - "nid.close()" + "if base_case_name:\n", + " nid = Dataset(base_file, \"r\")\n", + " smb_cism = np.squeeze(nid.variables[\"glc1Exp_Flgl_qice\"][0, :, :]) * smb_convert\n", + " nid.close()" ] }, { @@ -214,35 +262,14 @@ "metadata": {}, "outputs": [], "source": [ - "# creating the SMB climatology for comparison with observation\n", - "\n", - "# Initializing a field for the climatology\n", - "smb_cism_climo = np.zeros((ny_cism, nx_cism))\n", - "\n", - "# Counter for available year (only needed if the number of years available is smaller\n", - "# than the number of years requested to create the climatology.\n", - "count_yr = 0\n", - "\n", - "for k in range(ny_climo):\n", - "\n", - " year_to_read = last_year - k\n", - " str_yr = str(year_to_read).zfill(4)\n", - " file_cpl = f\"{path_cpl}/{case_name2}.cpl.hx.1yr2glc.{str_yr}-01-01-00000.nc\"\n", - "\n", - " if not os.path.isfile(file_cpl):\n", - " print(\"The couple file for time\", str_yr, \"does not exist.\")\n", - " print(\n", - " \"We will only use the files that existed until now to create the SMB climatology.\"\n", - " )\n", - " break\n", - "\n", - " smb_cism_climo = smb_cism_climo + read_smb(file_cpl)\n", - " count_yr = count_yr + 1\n", - "\n", - "# Averaging the climo data\n", - "smb_cism_climo = smb_cism_climo / count_yr\n", - "\n", - "print(\"number of years used in climatology = \", count_yr)" + "# creating the SMB climatology for new case\n", + "smb_case_climo = create_climo(case_path, case_name, last_year, params)\n", + "\n", + "# creating the SMB climatology for base_case\n", + "if base_case_name:\n", + " smb_base_climo = create_climo(\n", + " base_case_path, base_case_name, base_last_year, params\n", + " )" ] }, { @@ -264,14 +291,14 @@ ")\n", "\n", "# Initializing the glacier ID variable\n", - "smb_obs_cism = np.zeros((ny_cism, nx_cism))\n", + "smb_obs_climo = np.zeros((ny_cism, nx_cism))\n", "\n", "# Performing the interpolation\n", "for j in range(ny_cism):\n", " point_y = np.zeros(nx_cism)\n", " point_y[:] = y_cism[j]\n", " pts = (x_cism[:], point_y[:])\n", - " smb_obs_cism[j, :] = myInterpFunction_smb_obs(pts)" + " smb_obs_climo[j, :] = myInterpFunction_smb_obs(pts)" ] }, { @@ -282,33 +309,63 @@ "outputs": [], "source": [ "# Filtering out fill values\n", - "smb_obs_cism = np.where(smb_obs_cism > 1e20, 0, smb_obs_cism)\n", + "smb_obs_climo = np.where(smb_obs_climo > 1e20, 0, smb_obs_climo)\n", "mask = thk_cism[:, :] == 0\n", - "smb_obs_cism = np.where(mask, 0, smb_obs_cism)\n", - "smb_cism = np.where(mask, 0, smb_cism)\n", - "smb_cism_climo = np.where(mask, 0, smb_cism_climo)" + "smb_obs_climo = np.where(mask, 0, smb_obs_climo)\n", + "smb_case_climo = np.where(mask, 0, smb_case_climo)\n", + "if base_case_name:\n", + " smb_cism = np.where(mask, 0, smb_cism)\n", + " smb_base_climo = np.where(mask, 0, smb_base_climo)" ] }, { "cell_type": "code", "execution_count": null, - "id": "c7973cfe-64e0-47d4-a1b6-73cd9e62fdb2", + "id": "5514587c-54b3-4c5b-8444-fa877dcc3ed0", "metadata": {}, "outputs": [], "source": [ - "# Comparing the SMB\n", + "def plot_contour(data, ax, title, vmin, vmax, cmap):\n", + " avg_data = np.round(net_avrg(data) * mm_to_Gt, 2)\n", + " last_panel0 = ax.imshow(smb_case_climo[:, :], vmin=vmin, vmax=vmax, cmap=cmap)\n", + " ax.set_title(title, fontsize=16)\n", + " set_plot_prop_clean(ax)\n", + " ax.annotate(\"net avg =\" + str(avg_data) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", + "\n", + " pos = ax.get_position()\n", + " cax = fig.add_axes([0.35, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "\n", + " cbar = fig.colorbar(last_panel0, cax=cax)\n", + " cbar.ax.tick_params(labelsize=16)\n", + "\n", + "\n", + "def plot_contour_diff(data_new, data_old, ax, title, vmin, vmax, cmap):\n", + " avg_data = np.round(net_avrg(data_new - data_old) * mm_to_Gt, 2)\n", + " last_panel2 = ax.imshow(data_new - data_old, vmin=vmin, vmax=vmax, cmap=cmap)\n", + "\n", + " ax.set_title(title, fontsize=16)\n", + " set_plot_prop_clean(ax)\n", + "\n", + " ax.annotate(\"net avg =\" + str(avg_data) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", "\n", - "# Computing the SMB difference\n", - "smb_diff_cesm3 = smb_cism_climo - smb_obs_cism\n", + " pos = ax.get_position()\n", + " cax = fig.add_axes([0.89, pos.y0, 0.02, pos.y1 - pos.y0])\n", "\n", - "# Computing the net averages SMB\n", - "avg_smb_cism_climo = np.round(net_avrg(smb_cism_climo) * mm_to_Gt, 2)\n", - "avg_smb_obs_cism = np.round(net_avrg(smb_obs_cism) * mm_to_Gt, 2)\n", - "avg_smb_diff_cesm3 = np.round(net_avrg(smb_diff_cesm3) * mm_to_Gt, 2)\n", + " cbar = fig.colorbar(last_panel2, cax=cax)\n", + " cbar.ax.tick_params(labelsize=16)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c7973cfe-64e0-47d4-a1b6-73cd9e62fdb2", + "metadata": {}, + "outputs": [], + "source": [ + "# Comparing SMB new run vs obs\n", "\n", "# Colormap choice\n", "my_cmap = mcm.get_cmap(\"Spectral\")\n", - "# my_cmap = mcm.get_cmap('RdYlBu_r')\n", "my_cmap_diff = mcm.get_cmap(\"bwr_r\")\n", "\n", "\n", @@ -320,477 +377,252 @@ "fig, ax = plt.subplots(1, 3, sharey=True, figsize=[22, 9])\n", "\n", "## Left panel\n", - "last_panel0 = ax[0].imshow(smb_cism_climo[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap)\n", - "ax[0].set_title(\"SMB CESM3 (mm/y w.e.)\", fontsize=16)\n", - "set_plot_prop_clean(ax[0])\n", - "ax[0].annotate(\"net avg =\" + str(avg_smb_cism_climo) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", - "\n", - "pos = ax[0].get_position()\n", - "cax = fig.add_axes([0.35, pos.y0, 0.02, pos.y1 - pos.y0])\n", - "\n", - "cbar = fig.colorbar(last_panel0, cax=cax)\n", - "cbar.ax.tick_params(labelsize=16)\n", - "\n", + "plot_contour(\n", + " smb_case_climo, ax[0], f\"{case_name}\\nSMB (mm/y w.e.)\", vmin, vmax, my_cmap\n", + ")\n", "\n", "## Center panel\n", - "last_panel1 = ax[1].imshow(smb_obs_cism[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap)\n", - "\n", - "ax[1].set_title(\"SMB Obs (mm/y w.e.)\", fontsize=16)\n", - "set_plot_prop_clean(ax[1])\n", - "ax[1].annotate(\"net avg =\" + str(avg_smb_obs_cism) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", - "\n", - "\n", - "pos = ax[1].get_position()\n", - "cax = fig.add_axes([0.62, pos.y0, 0.02, pos.y1 - pos.y0])\n", - "\n", - "cbar = fig.colorbar(last_panel1, cax=cax)\n", - "cbar.ax.tick_params(labelsize=16)\n", - "\n", + "plot_contour(smb_obs_climo, ax[1], \"SMB Obs\\n(mm/y w.e.)\", vmin, vmax, my_cmap)\n", "\n", "## Right panel\n", - "\n", - "vmin = -2000\n", - "vmax = 2000\n", - "\n", - "last_panel2 = ax[2].imshow(smb_diff_cesm3, vmin=vmin, vmax=vmax, cmap=my_cmap_diff)\n", - "\n", - "ax[2].set_title(\"SMB diff (CESM3 - Obs, mm/yr w.e.)\", fontsize=16)\n", - "set_plot_prop_clean(ax[2])\n", - "\n", - "ax[2].annotate(\"net avg =\" + str(avg_smb_diff_cesm3) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", - "\n", - "pos = ax[2].get_position()\n", - "cax = fig.add_axes([0.89, pos.y0, 0.02, pos.y1 - pos.y0])\n", - "\n", - "cbar = fig.colorbar(last_panel2, cax=cax)\n", - "cbar.ax.tick_params(labelsize=16)\n", - "\n", - "fig.subplots_adjust(right=0.89)" + "plot_contour_diff(\n", + " smb_case_climo,\n", + " smb_obs_climo,\n", + " ax[2],\n", + " \"SMB bias (mm/yr w.e.)\",\n", + " vmin,\n", + " vmax,\n", + " my_cmap_diff,\n", + ")" ] }, { "cell_type": "code", "execution_count": null, - "id": "b7bbb2a4-c7ce-4327-9610-2f6ae0d7f3ce", + "id": "176594bc-53a1-4934-8210-7aa7d62f5659", "metadata": {}, "outputs": [], "source": [ - "# Integrated SMB time series\n", + "# Comparing SMB new run vs obs\n", "\n", - "# Initializing a field for the climatology\n", - "avg_smb_cism_timeseries = np.zeros(last_year)\n", + "if base_case_name:\n", + " # Colormap choice\n", + " my_cmap = mcm.get_cmap(\"Spectral\")\n", + " my_cmap_diff = mcm.get_cmap(\"bwr_r\")\n", "\n", - "# Counter for available year (only needed if the number of years available is smaller\n", - "# than the number of years requested to create the climatology.\n", - "count_yr = 0\n", - "\n", - "for k in range(last_year):\n", - "\n", - " year_to_read = last_year - k\n", - " str_yr = str(year_to_read).zfill(4)\n", - " file_cpl = f\"{path_cpl}/{case_name2}.cpl.hx.1yr2glc.\" + str_yr + \"-01-01-00000.nc\"\n", - "\n", - " if not os.path.isfile(file_cpl):\n", - " print(\"The couple file for time\", str_yr, \"does not exist.\")\n", - " print(\n", - " \"We will only use the files that existed until now to create the time series.\"\n", - " )\n", - " break\n", + " # Colorbar bounds\n", + " vmin = -2000\n", + " vmax = 2000\n", "\n", - " smb_temp = read_smb(file_cpl)\n", - " smb_temp = np.where(mask, 0, smb_temp)\n", + " # Figure\n", + " fig, ax = plt.subplots(1, 3, sharey=True, figsize=[22, 9])\n", "\n", - " avg_smb_cism_timeseries[year_to_read - 1] = np.round(\n", - " net_avrg(smb_temp) * mm_to_Gt, 2\n", + " ## Left panel\n", + " plot_contour(\n", + " smb_case_climo, ax[0], f\"{case_name}\\nSMB (mm/y w.e.)\", vmin, vmax, my_cmap\n", " )\n", - " count_yr = count_yr + 1\n", "\n", - " del smb_temp\n", - "\n", - "first_year = year_to_read\n", + " ## Center panel\n", + " plot_contour(\n", + " smb_base_climo, ax[1], f\"{base_case_name}\\nSMB (mm/y w.e.)\", vmin, vmax, my_cmap\n", + " )\n", "\n", - "print(\"number of years used in climatology = \", count_yr)" + " ## Right panel\n", + " plot_contour_diff(\n", + " smb_case_climo,\n", + " smb_base_climo,\n", + " ax[2],\n", + " \"SMB difference (mm/yr w.e.)\",\n", + " vmin,\n", + " vmax,\n", + " my_cmap_diff,\n", + " )" ] }, { "cell_type": "code", "execution_count": null, - "id": "01af1cd1-0351-452c-99e7-125546469f69", + "id": "ad3acc67-de0b-454c-995b-df279c1ee5a8", "metadata": {}, "outputs": [], "source": [ - "# Plotting the SMB spatially averaged time series\n", - "\n", - "\n", - "time = np.arange(first_year, last_year)\n", - "nt = len(time)\n", - "\n", - "\n", - "avg_smb_climo_timeseries = np.zeros(nt)\n", - "avg_smb_obs_cism_timeseries = np.zeros(nt)\n", - "\n", - "avg_smb_climo_timeseries[:] = avg_smb_cism_climo\n", - "avg_smb_obs_cism_timeseries[:] = avg_smb_obs_cism\n", - "\n", - "\n", - "x_ticks = np.arange(first_year, last_year, 10)\n", - "tickx = x_ticks\n", - "\n", - "ymin = 150\n", - "ymax = 650\n", - "y_step = 50\n", - "y_ticks = np.arange(ymin, ymax + y_step, y_step)\n", + "# Integrated SMB time series\n", + "def compute_annual_climo(path, case_name, last_year, params):\n", + " # Initializing a field for the climatology\n", + " avg_smb_timeseries = np.zeros(last_year)\n", "\n", + " # Counter for available year (only needed if the number of years available is smaller\n", + " # than the number of years requested to create the climatology.\n", + " count_yr = 0\n", "\n", - "plt.figure(figsize=(16, 7))\n", + " for k in range(last_year):\n", "\n", - "line = \"-\"\n", - "color = \"blue\"\n", - "sizefont = 16\n", - "linewidth = 2\n", + " year_to_read = last_year - k\n", + " file_name = (\n", + " f\"{path}/{case_name}.cpl.hx.1yr2glc.{year_to_read:04d}-01-01-00000.nc\"\n", + " )\n", "\n", + " if not os.path.isfile(file_name):\n", + " print(\"The couple file for time\", year_to_read, \"does not exist.\")\n", + " print(\n", + " \"We will only use the files that existed until now to create the time series.\"\n", + " )\n", + " break\n", "\n", - "# Plotting NESM3 experiments\n", - "plt.subplot(111)\n", + " smb_temp = read_smb(file_name)\n", + " smb_temp = np.where(params[\"mask\"], 0, smb_temp)\n", "\n", - "label = \"CISM\"\n", - "plt.plot(\n", - " time,\n", - " avg_smb_cism_timeseries[first_year::],\n", - " line,\n", - " ms=3,\n", - " mfc=color,\n", - " color=color,\n", - " label=label,\n", - " linewidth=linewidth,\n", - ")\n", - "\n", - "label = \"CISM climatology\"\n", - "color = \"k\"\n", - "plt.plot(\n", - " time,\n", - " avg_smb_climo_timeseries[:],\n", - " line,\n", - " ms=3,\n", - " mfc=color,\n", - " color=color,\n", - " label=label,\n", - " linewidth=linewidth,\n", - ")\n", + " avg_smb_timeseries[year_to_read - 1] = np.round(\n", + " net_avrg(smb_temp) * mm_to_Gt, 2\n", + " )\n", + " count_yr = count_yr + 1\n", "\n", - "label = \"Observation\"\n", - "color = \"r\"\n", - "plt.plot(\n", - " time,\n", - " avg_smb_obs_cism_timeseries[:],\n", - " line,\n", - " ms=3,\n", - " mfc=color,\n", - " color=color,\n", - " label=label,\n", - " linewidth=linewidth,\n", - ")\n", + " if count_yr == params[\"climo_nyears\"]:\n", + " break\n", "\n", + " del smb_temp\n", "\n", - "plt.xlim([first_year, last_year])\n", - "plt.xticks(x_ticks, tickx, fontsize=sizefont)\n", - "plt.xlabel(r\"$Time$ (y)\", fontsize=sizefont)\n", - "plt.ylabel(\"SMB average evolution (Gt/yr)\", multialignment=\"center\", fontsize=sizefont)\n", - "plt.ylim([ymin, ymax])\n", - "plt.yticks(fontsize=sizefont)\n", - "plt.legend(loc=\"upper left\", ncol=1, frameon=True, borderaxespad=0)\n", + " first_year = year_to_read\n", "\n", - "plt.title(\"SMB average evolution\", fontsize=sizefont)" - ] - }, - { - "cell_type": "markdown", - "id": "04c05f38-6115-4b54-8965-718baa1448a5", - "metadata": {}, - "source": [ - "# Adding comparison with CESM2" + " print(\"number of years used in climatology = \", count_yr)\n", + " return first_year, avg_smb_timeseries" ] }, { "cell_type": "code", "execution_count": null, - "id": "d2093fc4-042c-4670-9537-175e3df8dffb", + "id": "a38682c9-dc87-4d7b-887d-8abbbe8a7265", "metadata": {}, "outputs": [], "source": [ - "# Defining the paths and files for CESM2\n", - "\n", + "# Integrated SMB time series\n", + "params[\"mask\"] = mask\n", + "first_year, avg_smb_case_climo = compute_annual_climo(\n", + " case_path, case_name, last_year, params\n", + ")\n", "\n", - "nid = Dataset(file_cesm2, \"r\")\n", - "smb_cesm2 = np.squeeze(nid.variables[\"smb\"][0, :, :])\n", - "nid.close()" + "if base_case_name:\n", + " base_first_year, avg_smb_base_case_climo = compute_annual_climo(\n", + " base_case_path, base_case_name, base_last_year, params\n", + " )" ] }, { "cell_type": "code", "execution_count": null, - "id": "a8fa10a5-ab3f-467b-8614-bdf6cb702fd8", + "id": "a7510286-eae4-4d4c-be90-3e3f9c732eeb", "metadata": {}, "outputs": [], "source": [ - "smb_cesm2 = np.where(mask, 0, smb_cesm2)" + "def plot_line(data, time, line, color, label, linewidth):\n", + " plt.plot(\n", + " time,\n", + " data,\n", + " line,\n", + " ms=3,\n", + " mfc=color,\n", + " color=color,\n", + " label=label,\n", + " linewidth=linewidth,\n", + " )" ] }, { "cell_type": "code", "execution_count": null, - "id": "f43f1382-4927-4de7-bd04-3b59bb225db7", + "id": "01af1cd1-0351-452c-99e7-125546469f69", "metadata": {}, "outputs": [], "source": [ - "# Comparing the SMB\n", - "\n", - "# Computing the SMB difference\n", - "# smb_diff = smb_cism_climo - smb_obs_cism\n", - "\n", - "# Computing the net averages SMB\n", - "avg_smb_cesm2 = np.round(net_avrg(smb_cesm2) * mm_to_Gt, 2)\n", - "\n", - "\n", - "# Colormap choice\n", - "my_cmap = mcm.get_cmap(\"Spectral\")\n", - "# my_cmap = mcm.get_cmap('RdYlBu_r')\n", - "my_cmap_diff = mcm.get_cmap(\"bwr_r\")\n", - "\n", - "\n", - "# Colorbar bounds\n", - "vmin = -2000\n", - "vmax = 2000\n", - "\n", - "# Figure\n", - "fig, ax = plt.subplots(1, 3, sharey=True, figsize=[22, 9])\n", - "\n", - "## Left panel\n", - "last_panel0 = ax[0].imshow(smb_cesm2[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap)\n", - "ax[0].set_title(\"SMB CESM2 (mm/y w.e.)\", fontsize=16)\n", - "set_plot_prop_clean(ax[0])\n", - "ax[0].annotate(\"net avg =\" + str(avg_smb_cesm2) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", - "\n", - "pos = ax[0].get_position()\n", - "cax = fig.add_axes([0.35, pos.y0, 0.02, pos.y1 - pos.y0])\n", - "\n", - "cbar = fig.colorbar(last_panel0, cax=cax)\n", - "cbar.ax.tick_params(labelsize=16)\n", - "\n", - "\n", - "## Center panel\n", - "last_panel1 = ax[1].imshow(smb_cism_climo[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap)\n", - "\n", - "ax[1].set_title(\"SMB CESM3 (mm/y w.e.)\", fontsize=16)\n", - "set_plot_prop_clean(ax[1])\n", - "ax[1].annotate(\"net avg =\" + str(avg_smb_cism_climo) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", - "\n", - "\n", - "pos = ax[1].get_position()\n", - "cax = fig.add_axes([0.62, pos.y0, 0.02, pos.y1 - pos.y0])\n", - "\n", - "cbar = fig.colorbar(last_panel1, cax=cax)\n", - "cbar.ax.tick_params(labelsize=16)\n", - "\n", - "\n", - "## Right panel\n", - "\n", - "last_panel2 = ax[2].imshow(smb_obs_cism[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap)\n", - "\n", - "ax[2].set_title(\"SMB Obs (mm/y w.e.)\", fontsize=16)\n", - "set_plot_prop_clean(ax[2])\n", - "ax[2].annotate(\"net avg =\" + str(avg_smb_obs_cism) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", - "\n", - "\n", - "pos = ax[2].get_position()\n", - "cax = fig.add_axes([0.62, pos.y0, 0.02, pos.y1 - pos.y0])\n", - "\n", - "cbar = fig.colorbar(last_panel1, cax=cax)\n", - "cbar.ax.tick_params(labelsize=16)\n", - "\n", - "\n", - "# last_panel2 = ax[2].imshow(smb_diff_climo, vmin=vmin, vmax=vmax, cmap=my_cmap)\n", - "\n", - "# ax[2].set_title('SMB diff (CISM - Obs, mm/yr w.e.)',fontsize=16)\n", - "# set_plot_prop_clean(ax[2])\n", - "\n", - "# ax[2].annotate('net avg =' + str(avg_smb_diff) + ' Gt/yr',xy=(5,5),fontsize=16)\n", - "\n", - "# pos = ax[2].get_position()\n", - "# cax = fig.add_axes([0.89, pos.y0, 0.02, pos.y1 - pos.y0])\n", - "\n", - "# cbar = fig.colorbar(last_panel2,cax=cax)\n", - "# cbar.ax.tick_params(labelsize=16)\n", + "# Plotting the SMB spatially averaged time series\n", "\n", - "fig.subplots_adjust(right=0.89)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "85d4a99e-73b1-46eb-b9e2-fad93c155c96", - "metadata": {}, - "outputs": [], - "source": [ - "# Computing the difference w.r.t MAR obs\n", - "smb_diff_cesm2 = smb_cesm2 - smb_obs_cism\n", - "\n", - "# Computing the net averages SMB\n", - "avg_smb_cesm2 = np.round(net_avrg(smb_cesm2) * mm_to_Gt, 2)\n", - "avg_smb_diff_cesm2 = np.round(net_avrg(smb_diff_cesm2) * mm_to_Gt, 2)\n", - "\n", - "# Masking out ocean cells\n", - "mask2 = smb_cesm2[:, :] == 0\n", - "smb_cesm2_plot = np.where(mask2, np.float64(\"Nan\"), smb_cesm2)\n", - "smb_obs_cism_plot = np.where(mask2, np.float64(\"Nan\"), smb_obs_cism)\n", - "smb_cism_climo_plot = np.where(mask2, np.float64(\"Nan\"), smb_cism_climo)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "86e9b67b-ac4e-4551-96ea-ffff814fa038", - "metadata": {}, - "outputs": [], - "source": [ - "# Comparing the SMB\n", + "# TODO: include base case, base case climo (horizontal line), new case, new case climo, and obs climo\n", + "# Note: base case is 10 years of historical, new case is PI.\n", + "# what comparisons make sense when base case is HIST and new case is 1850?\n", "\n", - "# Computing the SMB difference\n", - "# smb_diff = smb_cism_climo - smb_obs_cism\n", "\n", + "time = np.arange(first_year, last_year)\n", + "if base_case_name:\n", + " base_time = np.arange(base_first_year, base_last_year) + last_year - base_last_year\n", + " base_nt = len(base_time)\n", + "nt = len(time)\n", "\n", - "# Colormap choice\n", - "my_cmap = mcm.get_cmap(\"Spectral\")\n", - "# my_cmap = mcm.get_cmap('RdYlBu_r')\n", - "my_cmap_diff = mcm.get_cmap(\"bwr_r\")\n", + "avg_smb_obs_timeseries = np.zeros(nt)\n", + "avg_smb_case_timeseries = np.zeros(nt)\n", + "if base_case_name:\n", + " avg_smb_base_timeseries = np.zeros(base_nt)\n", "\n", + "avg_smb_obs_timeseries[:] = np.round(net_avrg(smb_obs_climo) * mm_to_Gt, 2)\n", + "avg_smb_case_timeseries[:] = np.round(net_avrg(smb_case_climo) * mm_to_Gt, 2)\n", + "if base_case_name:\n", + " avg_smb_base_timeseries[:] = np.round(net_avrg(smb_base_climo) * mm_to_Gt, 2)\n", "\n", - "# Colorbar bounds\n", - "vmin = -2000\n", - "vmax = 2000\n", "\n", - "# Figure\n", - "fig, ax = plt.subplots(1, 3, sharey=True, figsize=[16, 9])\n", + "x_ticks = np.arange(first_year, last_year + 2, 5)\n", + "tickx = x_ticks\n", "\n", - "## Left panel\n", - "last_panel0 = ax[0].imshow(smb_cesm2_plot[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap)\n", - "# ax[0].set_title('SMB CESM2 (mm/y w.e.)',fontsize=16)\n", - "ax[0].set_title(\"CESM2\", fontsize=16)\n", - "set_plot_prop_clean(ax[0])\n", - "ax[0].annotate(\"net avg =\" + str(avg_smb_cesm2) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", - "ax[0].set_ylabel(\"SMB (mm/y w.e)\", fontsize=16)\n", - "# pos = ax[0].get_position()\n", - "# cax = fig.add_axes([0.35, pos.y0, 0.02, pos.y1 - pos.y0])\n", + "ymin = 100\n", + "ymax = 600\n", + "y_step = 50\n", + "y_ticks = np.arange(ymin, ymax + y_step, y_step)\n", "\n", - "# cbar = fig.colorbar(last_panel0, cax=cax)\n", - "# cbar.ax.tick_params(labelsize=16)\n", "\n", + "plt.figure(figsize=(16, 7))\n", "\n", - "## Center panel\n", - "last_panel1 = ax[1].imshow(\n", - " smb_cism_climo_plot[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap\n", + "# Plotting annual / spatial means\n", + "plt.subplot(111)\n", + "plot_line(\n", + " avg_smb_case_climo[first_year::],\n", + " time,\n", + " line=\"-\",\n", + " color=\"blue\",\n", + " label=f\"{case_name}\",\n", + " linewidth=2,\n", ")\n", - "\n", - "ax[1].set_title(\"CESM3\", fontsize=16)\n", - "set_plot_prop_clean(ax[1])\n", - "ax[1].annotate(\"net avg =\" + str(avg_smb_cism_climo) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", - "\n", - "\n", - "# pos = ax[1].get_position()\n", - "# cax = fig.add_axes([0.62, pos.y0, 0.02, pos.y1 - pos.y0])\n", - "\n", - "# cbar = fig.colorbar(last_panel1,cax=cax)\n", - "# cbar.ax.tick_params(labelsize=16)\n", - "\n", - "\n", - "## Right panel\n", - "\n", - "last_panel2 = ax[2].imshow(smb_obs_cism_plot[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap)\n", - "\n", - "ax[2].set_title(\"MARv3.12 1960-1999\", fontsize=16)\n", - "set_plot_prop_clean(ax[2])\n", - "ax[2].annotate(\"net avg =\" + str(avg_smb_obs_cism) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", - "\n", - "\n", - "pos = ax[2].get_position()\n", - "cax = fig.add_axes([0.9, pos.y0, 0.02, pos.y1 - pos.y0])\n", - "\n", - "cbar = fig.colorbar(last_panel1, cax=cax)\n", - "cbar.ax.tick_params(labelsize=16)\n", - "\n", - "\n", - "# last_panel2 = ax[2].imshow(smb_diff_climo, vmin=vmin, vmax=vmax, cmap=my_cmap)\n", - "\n", - "# ax[2].set_title('SMB diff (CISM - Obs, mm/yr w.e.)',fontsize=16)\n", - "# set_plot_prop_clean(ax[2])\n", - "\n", - "# ax[2].annotate('net avg =' + str(avg_smb_diff) + ' Gt/yr',xy=(5,5),fontsize=16)\n", - "\n", - "# pos = ax[2].get_position()\n", - "# cax = fig.add_axes([0.89, pos.y0, 0.02, pos.y1 - pos.y0])\n", - "\n", - "# cbar = fig.colorbar(last_panel2,cax=cax)\n", - "# cbar.ax.tick_params(labelsize=16)\n", - "\n", - "fig.subplots_adjust(right=0.89)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "307999aa-83e6-4dc4-b131-af7261d9b663", - "metadata": {}, - "outputs": [], - "source": [ - "# Comparing the SMB\n", - "\n", - "# Colormap choice\n", - "my_cmap = mcm.get_cmap(\"Spectral\")\n", - "# my_cmap = mcm.get_cmap('RdYlBu_r')\n", - "my_cmap_diff = mcm.get_cmap(\"bwr_r\")\n", - "\n", - "\n", - "# Colorbar bounds\n", - "vmin = -2000\n", - "vmax = 2000\n", - "\n", - "# Figure\n", - "fig, ax = plt.subplots(1, 2, sharey=True, figsize=[12, 9])\n", - "\n", - "## Left panel\n", - "last_panel0 = ax[0].imshow(\n", - " smb_diff_cesm2[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap_diff\n", + "plot_line(\n", + " avg_smb_case_timeseries[:],\n", + " time,\n", + " line=\":\",\n", + " color=\"blue\",\n", + " label=f\"{case_name} (mean)\",\n", + " linewidth=2,\n", ")\n", - "ax[0].set_title(\"CESM2 - MAR\", fontsize=16)\n", - "set_plot_prop_clean(ax[0])\n", - "ax[0].annotate(\"net avg =\" + str(avg_smb_diff_cesm2) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", - "ax[0].set_ylabel(\"SMB difference (mm/yr w.e.)\", fontsize=16)\n", - "\n", - "# pos = ax[0].get_position()\n", - "# cax = fig.add_axes([0.465, pos.y0, 0.02, pos.y1 - pos.y0])\n", - "\n", - "# cbar = fig.colorbar(last_panel0, cax=cax)\n", - "# cbar.ax.tick_params(labelsize=16)\n", - "\n", - "\n", - "## Right panel\n", - "last_panel1 = ax[1].imshow(\n", - " smb_diff_cesm3[:, :], vmin=vmin, vmax=vmax, cmap=my_cmap_diff\n", + "if base_case_name:\n", + " plot_line(\n", + " avg_smb_base_case_climo[base_first_year::],\n", + " base_time,\n", + " line=\"-\",\n", + " color=\"red\",\n", + " label=f\"{base_case_name}\",\n", + " linewidth=2,\n", + " )\n", + " plot_line(\n", + " avg_smb_base_timeseries[:],\n", + " base_time,\n", + " line=\":\",\n", + " color=\"red\",\n", + " label=f\"{base_case_name} (mean)\",\n", + " linewidth=2,\n", + " )\n", + "plot_line(\n", + " avg_smb_obs_timeseries[:],\n", + " time,\n", + " line=\"-\",\n", + " color=\"black\",\n", + " label=\"Observations (mean)\",\n", + " linewidth=2,\n", ")\n", "\n", - "ax[1].set_title(\"CESM3 - MAR\", fontsize=16)\n", - "set_plot_prop_clean(ax[1])\n", - "ax[1].annotate(\"net avg =\" + str(avg_smb_diff_cesm3) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", - "\n", - "pos = ax[1].get_position()\n", - "cax = fig.add_axes([0.89, pos.y0, 0.02, pos.y1 - pos.y0])\n", - "\n", - "cbar = fig.colorbar(last_panel1, cax=cax)\n", - "cbar.ax.tick_params(labelsize=16)\n", - "\n", + "sizefont = 16\n", + "plt.xlim([first_year, last_year])\n", + "plt.xticks(x_ticks, tickx, fontsize=sizefont)\n", + "plt.xlabel(r\"$Time$ (y)\", fontsize=sizefont)\n", + "plt.ylabel(\"SMB average evolution (Gt/yr)\", multialignment=\"center\", fontsize=sizefont)\n", + "plt.ylim([ymin, ymax])\n", + "plt.yticks(fontsize=sizefont)\n", + "plt.legend(loc=\"upper left\", ncol=1, frameon=True, borderaxespad=0)\n", "\n", - "# fig.tight_layout()\n", - "fig.subplots_adjust(right=0.89)" + "plt.title(\"SMB average evolution\", fontsize=sizefont);" ] } ], From b5e1150a88153511d3c2469b876f6092b9f895a2 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Wed, 24 Jul 2024 15:19:13 -0600 Subject: [PATCH 03/11] Pull functions out of glacier notebook Created a single utils module for now --- .../nblibrary/glc/LIWG_SMB_diagnostic.ipynb | 281 +++--------------- examples/nblibrary/glc/utils.py | 171 +++++++++++ 2 files changed, 220 insertions(+), 232 deletions(-) create mode 100644 examples/nblibrary/glc/utils.py diff --git a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb index 5420bef..01ca495 100644 --- a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb +++ b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb @@ -29,6 +29,8 @@ "import os\n", "from scipy.interpolate import RegularGridInterpolator\n", "\n", + "import utils\n", + "\n", "# to display figures in notebook after executing the code.\n", "%matplotlib inline" ] @@ -134,68 +136,11 @@ "# Constants\n", "res = np.abs(x_cism[1] - x_cism[0]) # CISM output resolution\n", "\n", - "rhoi = 917 # ice density kg/m3\n", "rhow = 1000 # water density kg/m3\n", - "sec_in_yr = 60 * 60 * 24 * 365 # seconds in a year\n", - "\n", - "smb_convert = sec_in_yr / rhoi * 1000 # converting kg m-2 s-1 ice to mm y-1 w.e.\n", "kg_to_Gt = 1e-12 # Converting kg to Gt\n", "mm_to_Gt = rhow * 1e-3 * res**2 * kg_to_Gt # converting mm/yr to Gt/yr" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "4cf8782d-1500-4707-bd3b-88aa2e2c9120", - "metadata": {}, - "outputs": [], - "source": [ - "# Functions used in this notebook\n", - "def set_plot_prop_clean(ax):\n", - " \"\"\"\n", - " This function cleans up the figures from unnecessary default figure properties.\n", - "\n", - " \"\"\"\n", - " ax.invert_yaxis()\n", - " ax.set_xlabel(\"\")\n", - " ax.set_ylabel(\"\")\n", - " ax.set_xticklabels(\"\")\n", - " ax.set_yticklabels(\"\")\n", - " ax.set_xticks([])\n", - " ax.set_yticks([])\n", - "\n", - "\n", - "def rmse(prediction, target):\n", - " \"\"\"\n", - " This function returns the root mean square error for the SMB.\n", - " Input:\n", - " prediction = field to predict\n", - " target = field to compare with the prediction\n", - " \"\"\"\n", - " return np.sqrt(((prediction - target) ** 2).mean())\n", - "\n", - "\n", - "def net_avrg(data):\n", - " \"\"\"\n", - " This function returns the net average of a data field\n", - " \"\"\"\n", - " return np.sum(np.sum(data, axis=0), axis=0)\n", - "\n", - "\n", - "def read_smb(file):\n", - " \"\"\"\n", - " This function reads the CISM SMB dataset from a CESM simulation output\n", - " in the cpl directory. The output is adjusted to be converted to mm/yr w.e unit.\n", - "\n", - " Input:\n", - " file: name of the file to extract the SMB\n", - " \"\"\"\n", - " nid = Dataset(file, \"r\")\n", - " smb_cism = np.squeeze(nid.variables[\"glc1Exp_Flgl_qice\"][0, :, :]) * smb_convert\n", - " nid.close()\n", - " return smb_cism" - ] - }, { "cell_type": "code", "execution_count": null, @@ -203,56 +148,12 @@ "metadata": {}, "outputs": [], "source": [ - "# More functions used in this notebook\n", "params = {\n", " \"ny_cism\": ny_cism,\n", " \"nx_cism\": nx_cism,\n", " \"climo_nyears\": climo_nyears,\n", - "}\n", - "\n", - "\n", - "def create_climo(path, case_name, last_year, params):\n", - " # Initializing a field for the climatology\n", - " climo_out = np.zeros((params[\"ny_cism\"], params[\"nx_cism\"]))\n", - "\n", - " # Counter for available year (only needed if the number of years available is smaller\n", - " # than the number of years requested to create the climatology.\n", - " count_yr = 0\n", - "\n", - " for k in range(params[\"climo_nyears\"]):\n", - "\n", - " year_to_read = last_year - k\n", - " filename = (\n", - " f\"{path}/{case_name}.cpl.hx.1yr2glc.{year_to_read:04d}-01-01-00000.nc\"\n", - " )\n", - "\n", - " if not os.path.isfile(filename):\n", - " print(f\"The couple file for time {year_to_read} does not exist.\")\n", - " print(\n", - " \"We will only use the files that existed until now to create the SMB climatology.\"\n", - " )\n", - " break\n", - "\n", - " climo_out = climo_out + read_smb(filename)\n", - " count_yr = count_yr + 1\n", - "\n", - " print(\"number of years used in climatology = \", count_yr)\n", - " # Averaging the climo data\n", - " return climo_out / count_yr" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "fb320812-c7aa-408b-af7c-b252bcb1acb9", - "metadata": {}, - "outputs": [], - "source": [ - "# Loading the data\n", - "if base_case_name:\n", - " nid = Dataset(base_file, \"r\")\n", - " smb_cism = np.squeeze(nid.variables[\"glc1Exp_Flgl_qice\"][0, :, :]) * smb_convert\n", - " nid.close()" + " \"mm_to_Gt\": mm_to_Gt,\n", + "}" ] }, { @@ -263,11 +164,11 @@ "outputs": [], "source": [ "# creating the SMB climatology for new case\n", - "smb_case_climo = create_climo(case_path, case_name, last_year, params)\n", + "smb_case_climo = utils.create_climo(case_path, case_name, last_year, params)\n", "\n", "# creating the SMB climatology for base_case\n", "if base_case_name:\n", - " smb_base_climo = create_climo(\n", + " smb_base_climo = utils.create_climo(\n", " base_case_path, base_case_name, base_last_year, params\n", " )" ] @@ -314,47 +215,9 @@ "smb_obs_climo = np.where(mask, 0, smb_obs_climo)\n", "smb_case_climo = np.where(mask, 0, smb_case_climo)\n", "if base_case_name:\n", - " smb_cism = np.where(mask, 0, smb_cism)\n", " smb_base_climo = np.where(mask, 0, smb_base_climo)" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "5514587c-54b3-4c5b-8444-fa877dcc3ed0", - "metadata": {}, - "outputs": [], - "source": [ - "def plot_contour(data, ax, title, vmin, vmax, cmap):\n", - " avg_data = np.round(net_avrg(data) * mm_to_Gt, 2)\n", - " last_panel0 = ax.imshow(smb_case_climo[:, :], vmin=vmin, vmax=vmax, cmap=cmap)\n", - " ax.set_title(title, fontsize=16)\n", - " set_plot_prop_clean(ax)\n", - " ax.annotate(\"net avg =\" + str(avg_data) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", - "\n", - " pos = ax.get_position()\n", - " cax = fig.add_axes([0.35, pos.y0, 0.02, pos.y1 - pos.y0])\n", - "\n", - " cbar = fig.colorbar(last_panel0, cax=cax)\n", - " cbar.ax.tick_params(labelsize=16)\n", - "\n", - "\n", - "def plot_contour_diff(data_new, data_old, ax, title, vmin, vmax, cmap):\n", - " avg_data = np.round(net_avrg(data_new - data_old) * mm_to_Gt, 2)\n", - " last_panel2 = ax.imshow(data_new - data_old, vmin=vmin, vmax=vmax, cmap=cmap)\n", - "\n", - " ax.set_title(title, fontsize=16)\n", - " set_plot_prop_clean(ax)\n", - "\n", - " ax.annotate(\"net avg =\" + str(avg_data) + \" Gt/yr\", xy=(5, 5), fontsize=16)\n", - "\n", - " pos = ax.get_position()\n", - " cax = fig.add_axes([0.89, pos.y0, 0.02, pos.y1 - pos.y0])\n", - "\n", - " cbar = fig.colorbar(last_panel2, cax=cax)\n", - " cbar.ax.tick_params(labelsize=16)" - ] - }, { "cell_type": "code", "execution_count": null, @@ -377,22 +240,33 @@ "fig, ax = plt.subplots(1, 3, sharey=True, figsize=[22, 9])\n", "\n", "## Left panel\n", - "plot_contour(\n", - " smb_case_climo, ax[0], f\"{case_name}\\nSMB (mm/y w.e.)\", vmin, vmax, my_cmap\n", + "utils.plot_contour(\n", + " smb_case_climo,\n", + " fig,\n", + " ax[0],\n", + " f\"{case_name}\\nSMB (mm/y w.e.)\",\n", + " vmin,\n", + " vmax,\n", + " my_cmap,\n", + " mm_to_Gt,\n", ")\n", "\n", "## Center panel\n", - "plot_contour(smb_obs_climo, ax[1], \"SMB Obs\\n(mm/y w.e.)\", vmin, vmax, my_cmap)\n", + "utils.plot_contour(\n", + " smb_obs_climo, fig, ax[1], \"SMB Obs\\n(mm/y w.e.)\", vmin, vmax, my_cmap, mm_to_Gt\n", + ")\n", "\n", "## Right panel\n", - "plot_contour_diff(\n", + "utils.plot_contour_diff(\n", " smb_case_climo,\n", " smb_obs_climo,\n", + " fig,\n", " ax[2],\n", " \"SMB bias (mm/yr w.e.)\",\n", " vmin,\n", " vmax,\n", " my_cmap_diff,\n", + " mm_to_Gt,\n", ")" ] }, @@ -404,7 +278,6 @@ "outputs": [], "source": [ "# Comparing SMB new run vs obs\n", - "\n", "if base_case_name:\n", " # Colormap choice\n", " my_cmap = mcm.get_cmap(\"Spectral\")\n", @@ -418,17 +291,29 @@ " fig, ax = plt.subplots(1, 3, sharey=True, figsize=[22, 9])\n", "\n", " ## Left panel\n", - " plot_contour(\n", - " smb_case_climo, ax[0], f\"{case_name}\\nSMB (mm/y w.e.)\", vmin, vmax, my_cmap\n", + " utils.plot_contour(\n", + " smb_case_climo,\n", + " ax[0],\n", + " f\"{case_name}\\nSMB (mm/y w.e.)\",\n", + " vmin,\n", + " vmax,\n", + " my_cmap,\n", + " mm_to_Gt,\n", " )\n", "\n", " ## Center panel\n", - " plot_contour(\n", - " smb_base_climo, ax[1], f\"{base_case_name}\\nSMB (mm/y w.e.)\", vmin, vmax, my_cmap\n", + " utils.plot_contour(\n", + " smb_base_climo,\n", + " ax[1],\n", + " f\"{base_case_name}\\nSMB (mm/y w.e.)\",\n", + " vmin,\n", + " vmax,\n", + " my_cmap,\n", + " mm_to_Gt,\n", " )\n", "\n", " ## Right panel\n", - " plot_contour_diff(\n", + " utils.plot_contour_diff(\n", " smb_case_climo,\n", " smb_base_climo,\n", " ax[2],\n", @@ -436,58 +321,10 @@ " vmin,\n", " vmax,\n", " my_cmap_diff,\n", + " mm_to_Gt,\n", " )" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "ad3acc67-de0b-454c-995b-df279c1ee5a8", - "metadata": {}, - "outputs": [], - "source": [ - "# Integrated SMB time series\n", - "def compute_annual_climo(path, case_name, last_year, params):\n", - " # Initializing a field for the climatology\n", - " avg_smb_timeseries = np.zeros(last_year)\n", - "\n", - " # Counter for available year (only needed if the number of years available is smaller\n", - " # than the number of years requested to create the climatology.\n", - " count_yr = 0\n", - "\n", - " for k in range(last_year):\n", - "\n", - " year_to_read = last_year - k\n", - " file_name = (\n", - " f\"{path}/{case_name}.cpl.hx.1yr2glc.{year_to_read:04d}-01-01-00000.nc\"\n", - " )\n", - "\n", - " if not os.path.isfile(file_name):\n", - " print(\"The couple file for time\", year_to_read, \"does not exist.\")\n", - " print(\n", - " \"We will only use the files that existed until now to create the time series.\"\n", - " )\n", - " break\n", - "\n", - " smb_temp = read_smb(file_name)\n", - " smb_temp = np.where(params[\"mask\"], 0, smb_temp)\n", - "\n", - " avg_smb_timeseries[year_to_read - 1] = np.round(\n", - " net_avrg(smb_temp) * mm_to_Gt, 2\n", - " )\n", - " count_yr = count_yr + 1\n", - "\n", - " if count_yr == params[\"climo_nyears\"]:\n", - " break\n", - "\n", - " del smb_temp\n", - "\n", - " first_year = year_to_read\n", - "\n", - " print(\"number of years used in climatology = \", count_yr)\n", - " return first_year, avg_smb_timeseries" - ] - }, { "cell_type": "code", "execution_count": null, @@ -497,36 +334,16 @@ "source": [ "# Integrated SMB time series\n", "params[\"mask\"] = mask\n", - "first_year, avg_smb_case_climo = compute_annual_climo(\n", + "first_year, avg_smb_case_climo = utils.compute_annual_climo(\n", " case_path, case_name, last_year, params\n", ")\n", "\n", "if base_case_name:\n", - " base_first_year, avg_smb_base_case_climo = compute_annual_climo(\n", + " base_first_year, avg_smb_base_case_climo = utils.compute_annual_climo(\n", " base_case_path, base_case_name, base_last_year, params\n", " )" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "a7510286-eae4-4d4c-be90-3e3f9c732eeb", - "metadata": {}, - "outputs": [], - "source": [ - "def plot_line(data, time, line, color, label, linewidth):\n", - " plt.plot(\n", - " time,\n", - " data,\n", - " line,\n", - " ms=3,\n", - " mfc=color,\n", - " color=color,\n", - " label=label,\n", - " linewidth=linewidth,\n", - " )" - ] - }, { "cell_type": "code", "execution_count": null, @@ -552,10 +369,10 @@ "if base_case_name:\n", " avg_smb_base_timeseries = np.zeros(base_nt)\n", "\n", - "avg_smb_obs_timeseries[:] = np.round(net_avrg(smb_obs_climo) * mm_to_Gt, 2)\n", - "avg_smb_case_timeseries[:] = np.round(net_avrg(smb_case_climo) * mm_to_Gt, 2)\n", + "avg_smb_obs_timeseries[:] = np.round(utils.net_avrg(smb_obs_climo) * mm_to_Gt, 2)\n", + "avg_smb_case_timeseries[:] = np.round(utils.net_avrg(smb_case_climo) * mm_to_Gt, 2)\n", "if base_case_name:\n", - " avg_smb_base_timeseries[:] = np.round(net_avrg(smb_base_climo) * mm_to_Gt, 2)\n", + " avg_smb_base_timeseries[:] = np.round(utils.net_avrg(smb_base_climo) * mm_to_Gt, 2)\n", "\n", "\n", "x_ticks = np.arange(first_year, last_year + 2, 5)\n", @@ -571,7 +388,7 @@ "\n", "# Plotting annual / spatial means\n", "plt.subplot(111)\n", - "plot_line(\n", + "utils.plot_line(\n", " avg_smb_case_climo[first_year::],\n", " time,\n", " line=\"-\",\n", @@ -579,7 +396,7 @@ " label=f\"{case_name}\",\n", " linewidth=2,\n", ")\n", - "plot_line(\n", + "utils.plot_line(\n", " avg_smb_case_timeseries[:],\n", " time,\n", " line=\":\",\n", @@ -588,7 +405,7 @@ " linewidth=2,\n", ")\n", "if base_case_name:\n", - " plot_line(\n", + " utils.plot_line(\n", " avg_smb_base_case_climo[base_first_year::],\n", " base_time,\n", " line=\"-\",\n", @@ -596,7 +413,7 @@ " label=f\"{base_case_name}\",\n", " linewidth=2,\n", " )\n", - " plot_line(\n", + " utils.plot_line(\n", " avg_smb_base_timeseries[:],\n", " base_time,\n", " line=\":\",\n", @@ -604,7 +421,7 @@ " label=f\"{base_case_name} (mean)\",\n", " linewidth=2,\n", " )\n", - "plot_line(\n", + "utils.plot_line(\n", " avg_smb_obs_timeseries[:],\n", " time,\n", " line=\"-\",\n", diff --git a/examples/nblibrary/glc/utils.py b/examples/nblibrary/glc/utils.py new file mode 100644 index 0000000..eae4c2d --- /dev/null +++ b/examples/nblibrary/glc/utils.py @@ -0,0 +1,171 @@ +from __future__ import annotations + +import os + +import numpy as np +from matplotlib import pyplot as plt +from netCDF4 import Dataset + + +def set_plot_prop_clean(ax): + """ + This function cleans up the figures from unnecessary default figure properties. + + """ + ax.invert_yaxis() + ax.set_xlabel("") + ax.set_ylabel("") + ax.set_xticklabels("") + ax.set_yticklabels("") + ax.set_xticks([]) + ax.set_yticks([]) + + +def rmse(prediction, target): + """ + This function returns the root mean square error for the SMB. + Input: + prediction = field to predict + target = field to compare with the prediction + """ + return np.sqrt(((prediction - target) ** 2).mean()) + + +def net_avrg(data): + """ + This function returns the net average of a data field + """ + return np.sum(np.sum(data, axis=0), axis=0) + + +def read_smb(file): + """ + This function reads the CISM SMB dataset from a CESM simulation output + in the cpl directory. The output is adjusted to be converted to mm/yr w.e unit. + + Input: + file: name of the file to extract the SMB + """ + rhoi = 917 # ice density kg/m3 + sec_in_yr = 60 * 60 * 24 * 365 # seconds in a year + smb_convert = sec_in_yr / rhoi * 1000 # converting kg m-2 s-1 ice to mm y-1 w.e. + + nid = Dataset(file, "r") + smb_cism = np.squeeze(nid.variables["glc1Exp_Flgl_qice"][0, :, :]) * smb_convert + nid.close() + return smb_cism + + +def create_climo(path, case_name, last_year, params): + # Initializing a field for the climatology + climo_out = np.zeros((params["ny_cism"], params["nx_cism"])) + + # Counter for available year (only needed if the number of years available is smaller + # than the number of years requested to create the climatology. + count_yr = 0 + + for k in range(params["climo_nyears"]): + + year_to_read = last_year - k + filename = ( + f"{path}/{case_name}.cpl.hx.1yr2glc.{year_to_read:04d}-01-01-00000.nc" + ) + + if not os.path.isfile(filename): + print(f"The couple file for time {year_to_read} does not exist.") + print( + "We will only use the files that existed until now to create the SMB climatology.", + ) + break + + climo_out = climo_out + read_smb(filename) + count_yr = count_yr + 1 + + print("number of years used in climatology = ", count_yr) + # Averaging the climo data + return climo_out / count_yr + + +def plot_contour(data, fig, ax, title, vmin, vmax, cmap, mm_to_Gt): + avg_data = np.round(net_avrg(data) * mm_to_Gt, 2) + last_panel0 = ax.imshow(data[:, :], vmin=vmin, vmax=vmax, cmap=cmap) + ax.set_title(title, fontsize=16) + set_plot_prop_clean(ax) + ax.annotate("net avg =" + str(avg_data) + " Gt/yr", xy=(5, 5), fontsize=16) + + pos = ax.get_position() + cax = fig.add_axes([0.35, pos.y0, 0.02, pos.y1 - pos.y0]) + + cbar = fig.colorbar(last_panel0, cax=cax) + cbar.ax.tick_params(labelsize=16) + + +def plot_contour_diff(data_new, data_old, fig, ax, title, vmin, vmax, cmap, mm_to_Gt): + avg_data = np.round(net_avrg(data_new - data_old) * mm_to_Gt, 2) + last_panel2 = ax.imshow(data_new - data_old, vmin=vmin, vmax=vmax, cmap=cmap) + + ax.set_title(title, fontsize=16) + set_plot_prop_clean(ax) + + ax.annotate("net avg =" + str(avg_data) + " Gt/yr", xy=(5, 5), fontsize=16) + + pos = ax.get_position() + cax = fig.add_axes([0.89, pos.y0, 0.02, pos.y1 - pos.y0]) + + cbar = fig.colorbar(last_panel2, cax=cax) + cbar.ax.tick_params(labelsize=16) + + +def compute_annual_climo(path, case_name, last_year, params): + # Initializing a field for the climatology + avg_smb_timeseries = np.zeros(last_year) + + # Counter for available year (only needed if the number of years available is smaller + # than the number of years requested to create the climatology. + count_yr = 0 + + for k in range(last_year): + + year_to_read = last_year - k + file_name = ( + f"{path}/{case_name}.cpl.hx.1yr2glc.{year_to_read:04d}-01-01-00000.nc" + ) + + if not os.path.isfile(file_name): + print("The couple file for time", year_to_read, "does not exist.") + print( + "We will only use the files that existed until now to create the time series.", + ) + break + + smb_temp = read_smb(file_name) + smb_temp = np.where(params["mask"], 0, smb_temp) + + avg_smb_timeseries[year_to_read - 1] = np.round( + net_avrg(smb_temp) * params["mm_to_Gt"], + 2, + ) + count_yr = count_yr + 1 + + if count_yr == params["climo_nyears"]: + break + + del smb_temp + + first_year = year_to_read + + print("number of years used in climatology = ", count_yr) + return first_year, avg_smb_timeseries + + +def plot_line(data, time, line, color, label, linewidth): + plt.plot( + time, + data, + line, + ms=3, + mfc=color, + color=color, + label=label, + linewidth=linewidth, + ) From 3df4d5e7d1e2f976bc151aab0fa6b9c2ad1f2a49 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Wed, 24 Jul 2024 15:27:00 -0600 Subject: [PATCH 04/11] Remove unused import --- examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb index 01ca495..d295987 100644 --- a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb +++ b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb @@ -26,7 +26,6 @@ "import matplotlib.pyplot as plt\n", "import matplotlib.cm as mcm\n", "from netCDF4 import Dataset\n", - "import os\n", "from scipy.interpolate import RegularGridInterpolator\n", "\n", "import utils\n", From f1b015d2a668b49c801792e0182ee8936ce358e2 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Wed, 24 Jul 2024 16:10:39 -0600 Subject: [PATCH 05/11] Cleaned up glacier notebook a little more 1. Reorganized cells and added markdown headers (parameter configuration, set up grid, make datasets, generate plots). 2. Fixed issue with plotting base maps (tested by comparing first 40 years of run to last 40 years of run) 3. Added years included in climatology to the title of plots --- .../nblibrary/glc/LIWG_SMB_diagnostic.ipynb | 138 ++++++++++++------ 1 file changed, 93 insertions(+), 45 deletions(-) diff --git a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb index d295987..84610ea 100644 --- a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb +++ b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb @@ -34,6 +34,17 @@ "%matplotlib inline" ] }, + { + "cell_type": "markdown", + "id": "24029ab9-fe52-4c7e-bb41-847c7f3a0b1d", + "metadata": {}, + "source": [ + "## Parameter configuration\n", + "\n", + "Some parameters are set in CUPiD's `config.yml` file,\n", + "others are derived from these parameters." + ] + }, { "cell_type": "code", "execution_count": null, @@ -84,6 +95,16 @@ " base_file = f\"{base_case_path}/{base_case_name}.cpl.hx.1yr2glc.{base_last_year:04d}-01-01-00000.nc\" # name of last cpl simulation output" ] }, + { + "cell_type": "markdown", + "id": "70ee51c9-1c10-475d-a45b-b97583a3a5a9", + "metadata": {}, + "source": [ + "## Set up grid\n", + "\n", + "Read in the grid data, compute resolution and other grid-specific parameters" + ] + }, { "cell_type": "code", "execution_count": null, @@ -106,25 +127,6 @@ "ny_cism = len(y_cism)" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "69891e0a-0cc6-48bd-a9bf-b61e19da13b4", - "metadata": {}, - "outputs": [], - "source": [ - "## The observed dataset\n", - "nid = Dataset(obs_file, \"r\")\n", - "x_obs = nid.variables[\"x\"][:]\n", - "y_obs = nid.variables[\"y\"][:]\n", - "smb_obs_src = np.squeeze(nid.variables[\"SMB\"][0, :, :])\n", - "nid.close()\n", - "\n", - "## For the observed grid\n", - "nx_obs = len(x_obs)\n", - "ny_obs = len(y_obs)" - ] - }, { "cell_type": "code", "execution_count": null, @@ -155,6 +157,37 @@ "}" ] }, + { + "cell_type": "markdown", + "id": "68fca423-582b-4179-8771-16250a5f1904", + "metadata": {}, + "source": [ + "## Make datasets\n", + "\n", + "Read in observations and CESM output.\n", + "Also do necessary computations\n", + "(global mean for time series, temporal mean for climatology)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "69891e0a-0cc6-48bd-a9bf-b61e19da13b4", + "metadata": {}, + "outputs": [], + "source": [ + "## The observed dataset\n", + "nid = Dataset(obs_file, \"r\")\n", + "x_obs = nid.variables[\"x\"][:]\n", + "y_obs = nid.variables[\"y\"][:]\n", + "smb_obs_src = np.squeeze(nid.variables[\"SMB\"][0, :, :])\n", + "nid.close()\n", + "\n", + "## For the observed grid\n", + "nx_obs = len(x_obs)\n", + "ny_obs = len(y_obs)" + ] + }, { "cell_type": "code", "execution_count": null, @@ -217,6 +250,37 @@ " smb_base_climo = np.where(mask, 0, smb_base_climo)" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "a38682c9-dc87-4d7b-887d-8abbbe8a7265", + "metadata": {}, + "outputs": [], + "source": [ + "# Integrated SMB time series\n", + "params[\"mask\"] = mask\n", + "first_year, avg_smb_case_climo = utils.compute_annual_climo(\n", + " case_path, case_name, last_year, params\n", + ")\n", + "\n", + "if base_case_name:\n", + " base_first_year, avg_smb_base_case_climo = utils.compute_annual_climo(\n", + " base_case_path, base_case_name, base_last_year, params\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "1641747b-4997-45ad-bf70-981ed97688dd", + "metadata": {}, + "source": [ + "## Generate plots\n", + "\n", + "Map comparing CESM to observation,\n", + "possibly map comparing CESM to older case,\n", + "and time series of spatial mean SMB." + ] + }, { "cell_type": "code", "execution_count": null, @@ -243,7 +307,7 @@ " smb_case_climo,\n", " fig,\n", " ax[0],\n", - " f\"{case_name}\\nSMB (mm/y w.e.)\",\n", + " f\"{case_name}\\nSMB (mm/y w.e.)\\nMean from {first_year:04d} - {last_year:04d}\",\n", " vmin,\n", " vmax,\n", " my_cmap,\n", @@ -292,8 +356,9 @@ " ## Left panel\n", " utils.plot_contour(\n", " smb_case_climo,\n", + " fig,\n", " ax[0],\n", - " f\"{case_name}\\nSMB (mm/y w.e.)\",\n", + " f\"{case_name}\\nSMB (mm/y w.e.)\\nMean from {first_year:04d} - {last_year:04d}\",\n", " vmin,\n", " vmax,\n", " my_cmap,\n", @@ -303,8 +368,9 @@ " ## Center panel\n", " utils.plot_contour(\n", " smb_base_climo,\n", + " fig,\n", " ax[1],\n", - " f\"{base_case_name}\\nSMB (mm/y w.e.)\",\n", + " f\"{base_case_name}\\nSMB (mm/y w.e.)\\nMean from {base_first_year:04d} - {base_last_year:04d}\",\n", " vmin,\n", " vmax,\n", " my_cmap,\n", @@ -315,6 +381,7 @@ " utils.plot_contour_diff(\n", " smb_case_climo,\n", " smb_base_climo,\n", + " fig,\n", " ax[2],\n", " \"SMB difference (mm/yr w.e.)\",\n", " vmin,\n", @@ -324,25 +391,6 @@ " )" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "a38682c9-dc87-4d7b-887d-8abbbe8a7265", - "metadata": {}, - "outputs": [], - "source": [ - "# Integrated SMB time series\n", - "params[\"mask\"] = mask\n", - "first_year, avg_smb_case_climo = utils.compute_annual_climo(\n", - " case_path, case_name, last_year, params\n", - ")\n", - "\n", - "if base_case_name:\n", - " base_first_year, avg_smb_base_case_climo = utils.compute_annual_climo(\n", - " base_case_path, base_case_name, base_last_year, params\n", - " )" - ] - }, { "cell_type": "code", "execution_count": null, @@ -392,7 +440,7 @@ " time,\n", " line=\"-\",\n", " color=\"blue\",\n", - " label=f\"{case_name}\",\n", + " label=f\"{case_name} ({first_year:04d} - {last_year:04d})\",\n", " linewidth=2,\n", ")\n", "utils.plot_line(\n", @@ -400,7 +448,7 @@ " time,\n", " line=\":\",\n", " color=\"blue\",\n", - " label=f\"{case_name} (mean)\",\n", + " label=f\"{case_name} (mean from {first_year:04d} - {last_year:04d})\",\n", " linewidth=2,\n", ")\n", "if base_case_name:\n", @@ -409,7 +457,7 @@ " base_time,\n", " line=\"-\",\n", " color=\"red\",\n", - " label=f\"{base_case_name}\",\n", + " label=f\"{base_case_name} ({base_first_year:04d} - {base_last_year:04d})\",\n", " linewidth=2,\n", " )\n", " utils.plot_line(\n", @@ -417,7 +465,7 @@ " base_time,\n", " line=\":\",\n", " color=\"red\",\n", - " label=f\"{base_case_name} (mean)\",\n", + " label=f\"{base_case_name} (mean from {base_first_year:04d} - {base_last_year:04d})\",\n", " linewidth=2,\n", " )\n", "utils.plot_line(\n", From e9ddda307b2c0a3bf570cae59e41b98959471471 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Wed, 24 Jul 2024 16:19:19 -0600 Subject: [PATCH 06/11] Update book_toc section of config.yml Commented out the atmosphere, ocean, land, and sea ice section definitions and added a section for land ice. --- examples/key_metrics/config.yml | 32 ++++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/examples/key_metrics/config.yml b/examples/key_metrics/config.yml index 2b3e590..ceed3e6 100644 --- a/examples/key_metrics/config.yml +++ b/examples/key_metrics/config.yml @@ -182,31 +182,35 @@ book_toc: # Parts group notebooks into different sections in the Jupyter book # table of contents, so you can organize different parts of your project. - - caption: Atmosphere + # - caption: Atmosphere - # Each chapter is the name of one of the notebooks that you executed - # in compute_notebooks above, also without .ipynb - chapters: - - file: atm/adf_quick_run + # # Each chapter is the name of one of the notebooks that you executed + # # in compute_notebooks above, also without .ipynb + # chapters: + # - file: atm/adf_quick_run - - caption: Ocean - chapters: - - file: ocn/ocean_surface + # - caption: Ocean + # chapters: + # - file: ocn/ocean_surface - - caption: Land - chapters: - - file: lnd/land_comparison + # - caption: Land + # chapters: + # - file: lnd/land_comparison + + # - caption: Sea Ice + # chapters: + # - file: ice/seaice - - caption: Sea Ice + - caption: Land Ice chapters: - - file: ice/seaice + - file: glc/LIWG_SMB_diagnostic ##################################### # Keys for Jupyter Book _config.yml # ##################################### book_config_keys: - title: Example project # Title of your jupyter book + title: CESM Key Metrics # Title of your jupyter book # Other keys can be added here, see https://jupyterbook.org/en/stable/customize/config.html ### for many more options From 6b94626e13b61f3abc7a195a712dcb6a7b7fe868 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Fri, 26 Jul 2024 13:06:05 -0600 Subject: [PATCH 07/11] First steps towards using xarray Started replacing netcdf4 calls with xarray equivalents. I think the observational dataset read is still entirely netcdf4-based, so that needs to be updated. Also, I want to try calling xr.open_mfdataset() once per case and then passing datasets rather than doing the read twice (once for temporal mean, once for spatial mean). --- .../nblibrary/glc/LIWG_SMB_diagnostic.ipynb | 94 ++++++------------ examples/nblibrary/glc/utils.py | 96 +++++++++---------- 2 files changed, 73 insertions(+), 117 deletions(-) diff --git a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb index 84610ea..1731701 100644 --- a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb +++ b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb @@ -25,8 +25,8 @@ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import matplotlib.cm as mcm\n", - "from netCDF4 import Dataset\n", "from scipy.interpolate import RegularGridInterpolator\n", + "import xarray as xr\n", "\n", "import utils\n", "\n", @@ -108,23 +108,16 @@ { "cell_type": "code", "execution_count": null, - "id": "52ba1301-0c4c-481c-a57b-a9614160807b", + "id": "27373d08-084b-4c3f-8c3f-d5c8a445b2dc", "metadata": {}, "outputs": [], "source": [ - "## CISM grid information and loading a field used for filtering\n", - "\n", - "# Reading the information we need from the glc file\n", - "nid = Dataset(case_init_file, \"r\")\n", - "x_cism = nid.variables[\"x1\"][:]\n", - "y_cism = nid.variables[\"y1\"][:]\n", - "thk_cism = np.squeeze(nid.variables[\"thk\"][0, :, :])\n", - "nid.close()\n", - "\n", - "# Defining the grid dimensions\n", - "## For the CISM grid\n", - "nx_cism = len(x_cism)\n", - "ny_cism = len(y_cism)" + "## Get grid from initial_hist stream\n", + "thk_init_da = xr.open_dataset(case_init_file).isel(time=0)[\"thk\"]\n", + "mask = thk_init_da.data[:, :] == 0\n", + "\n", + "# Shape of array is (ny, nx)\n", + "grid_dims = thk_init_da.shape" ] }, { @@ -135,7 +128,9 @@ "outputs": [], "source": [ "# Constants\n", - "res = np.abs(x_cism[1] - x_cism[0]) # CISM output resolution\n", + "res = np.abs(\n", + " thk_init_da[\"x1\"].data[1] - thk_init_da[\"x1\"].data[0]\n", + ") # CISM output resolution\n", "\n", "rhow = 1000 # water density kg/m3\n", "kg_to_Gt = 1e-12 # Converting kg to Gt\n", @@ -150,10 +145,10 @@ "outputs": [], "source": [ "params = {\n", - " \"ny_cism\": ny_cism,\n", - " \"nx_cism\": nx_cism,\n", + " \"grid_dims\": grid_dims,\n", " \"climo_nyears\": climo_nyears,\n", " \"mm_to_Gt\": mm_to_Gt,\n", + " \"mask\": mask,\n", "}" ] }, @@ -169,25 +164,6 @@ "(global mean for time series, temporal mean for climatology)." ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "69891e0a-0cc6-48bd-a9bf-b61e19da13b4", - "metadata": {}, - "outputs": [], - "source": [ - "## The observed dataset\n", - "nid = Dataset(obs_file, \"r\")\n", - "x_obs = nid.variables[\"x\"][:]\n", - "y_obs = nid.variables[\"y\"][:]\n", - "smb_obs_src = np.squeeze(nid.variables[\"SMB\"][0, :, :])\n", - "nid.close()\n", - "\n", - "## For the observed grid\n", - "nx_obs = len(x_obs)\n", - "ny_obs = len(y_obs)" - ] - }, { "cell_type": "code", "execution_count": null, @@ -213,41 +189,30 @@ "outputs": [], "source": [ "# Interpolating the observed data onto the CISM grid\n", + "smb_obs_da = xr.open_dataset(obs_file).isel(time=0)[\"SMB\"]\n", "\n", "# Defining the interpolation functions\n", "myInterpFunction_smb_obs = RegularGridInterpolator(\n", - " (x_obs, y_obs),\n", - " smb_obs_src.transpose(),\n", + " (smb_obs_da[\"x\"].data, smb_obs_da[\"y\"].data),\n", + " smb_obs_da.data.transpose(),\n", " method=\"linear\",\n", " bounds_error=False,\n", " fill_value=None,\n", ")\n", "\n", "# Initializing the glacier ID variable\n", - "smb_obs_climo = np.zeros((ny_cism, nx_cism))\n", + "smb_obs_climo = np.zeros(grid_dims)\n", "\n", "# Performing the interpolation\n", - "for j in range(ny_cism):\n", - " point_y = np.zeros(nx_cism)\n", - " point_y[:] = y_cism[j]\n", - " pts = (x_cism[:], point_y[:])\n", - " smb_obs_climo[j, :] = myInterpFunction_smb_obs(pts)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "288f34cf-808f-4de1-8372-276fdbcfdc46", - "metadata": {}, - "outputs": [], - "source": [ + "for j in range(grid_dims[0]):\n", + " point_y = np.zeros(grid_dims[1])\n", + " point_y[:] = thk_init_da[\"y1\"].data[j]\n", + " pts = (thk_init_da[\"x1\"].data[:], point_y[:])\n", + " smb_obs_climo[j, :] = myInterpFunction_smb_obs(pts)\n", + "\n", "# Filtering out fill values\n", "smb_obs_climo = np.where(smb_obs_climo > 1e20, 0, smb_obs_climo)\n", - "mask = thk_cism[:, :] == 0\n", - "smb_obs_climo = np.where(mask, 0, smb_obs_climo)\n", - "smb_case_climo = np.where(mask, 0, smb_case_climo)\n", - "if base_case_name:\n", - " smb_base_climo = np.where(mask, 0, smb_base_climo)" + "smb_obs_climo = np.where(mask, 0, smb_obs_climo)" ] }, { @@ -258,7 +223,6 @@ "outputs": [], "source": [ "# Integrated SMB time series\n", - "params[\"mask\"] = mask\n", "first_year, avg_smb_case_climo = utils.compute_annual_climo(\n", " case_path, case_name, last_year, params\n", ")\n", @@ -405,9 +369,11 @@ "# what comparisons make sense when base case is HIST and new case is 1850?\n", "\n", "\n", - "time = np.arange(first_year, last_year)\n", + "time = np.arange(first_year, last_year + 1)\n", "if base_case_name:\n", - " base_time = np.arange(base_first_year, base_last_year) + last_year - base_last_year\n", + " base_time = (\n", + " np.arange(base_first_year, base_last_year + 1) + last_year - base_last_year\n", + " )\n", " base_nt = len(base_time)\n", "nt = len(time)\n", "\n", @@ -436,7 +402,7 @@ "# Plotting annual / spatial means\n", "plt.subplot(111)\n", "utils.plot_line(\n", - " avg_smb_case_climo[first_year::],\n", + " avg_smb_case_climo[first_year - 1 : :],\n", " time,\n", " line=\"-\",\n", " color=\"blue\",\n", @@ -453,7 +419,7 @@ ")\n", "if base_case_name:\n", " utils.plot_line(\n", - " avg_smb_base_case_climo[base_first_year::],\n", + " avg_smb_base_case_climo[base_first_year - 1 : :],\n", " base_time,\n", " line=\"-\",\n", " color=\"red\",\n", diff --git a/examples/nblibrary/glc/utils.py b/examples/nblibrary/glc/utils.py index eae4c2d..391d643 100644 --- a/examples/nblibrary/glc/utils.py +++ b/examples/nblibrary/glc/utils.py @@ -3,6 +3,7 @@ import os import numpy as np +import xarray as xr from matplotlib import pyplot as plt from netCDF4 import Dataset @@ -56,14 +57,13 @@ def read_smb(file): return smb_cism -def create_climo(path, case_name, last_year, params): - # Initializing a field for the climatology - climo_out = np.zeros((params["ny_cism"], params["nx_cism"])) - - # Counter for available year (only needed if the number of years available is smaller - # than the number of years requested to create the climatology. - count_yr = 0 +def _get_cesm_output(path, case_name, last_year, params): + # Set parameters + rhoi = 917 # ice density kg/m3 + sec_in_yr = 60 * 60 * 24 * 365 # seconds in a year + smb_convert = sec_in_yr / rhoi * 1000 # converting kg m-2 s-1 ice to mm y-1 w.e. + filenames = [] for k in range(params["climo_nyears"]): year_to_read = last_year - k @@ -78,12 +78,44 @@ def create_climo(path, case_name, last_year, params): ) break - climo_out = climo_out + read_smb(filename) - count_yr = count_yr + 1 + filenames.append(filename) + + climo_out = ( + xr.open_mfdataset(filenames)["glc1Exp_Flgl_qice"].compute() * smb_convert + ) + # Mask out data that is 0 in initial condition + for k in range(len(climo_out["time"])): + climo_out.data[k, :, :] = np.where( + params["mask"], + 0, + climo_out.isel(time=k).data, + ) + print("number of years used in climatology = ", len(climo_out["time"])) + return climo_out + + +def create_climo(path, case_name, last_year, params): + + climo_out = _get_cesm_output(path, case_name, last_year, params) - print("number of years used in climatology = ", count_yr) # Averaging the climo data - return climo_out / count_yr + return climo_out.mean("time").data + + +def compute_annual_climo(path, case_name, last_year, params): + # Initializing a field for the climatology + avg_smb_timeseries = np.zeros(last_year) + climo_out = _get_cesm_output(path, case_name, last_year, params) + for k in range(len(climo_out["time"])): + # index into avg_smb_timeseries; want largest k to index (last_year-1) + kk = last_year - (len(climo_out["time"]) - k) + # note that mm_to_Gt has 1 / res**2 factor, so we want a sum rather than mean + avg_smb_timeseries[kk] = np.round( + climo_out.isel(time=k).sum() * params["mm_to_Gt"], + 2, + ).data + + return last_year - len(climo_out["time"]) + 1, avg_smb_timeseries def plot_contour(data, fig, ax, title, vmin, vmax, cmap, mm_to_Gt): @@ -116,48 +148,6 @@ def plot_contour_diff(data_new, data_old, fig, ax, title, vmin, vmax, cmap, mm_t cbar.ax.tick_params(labelsize=16) -def compute_annual_climo(path, case_name, last_year, params): - # Initializing a field for the climatology - avg_smb_timeseries = np.zeros(last_year) - - # Counter for available year (only needed if the number of years available is smaller - # than the number of years requested to create the climatology. - count_yr = 0 - - for k in range(last_year): - - year_to_read = last_year - k - file_name = ( - f"{path}/{case_name}.cpl.hx.1yr2glc.{year_to_read:04d}-01-01-00000.nc" - ) - - if not os.path.isfile(file_name): - print("The couple file for time", year_to_read, "does not exist.") - print( - "We will only use the files that existed until now to create the time series.", - ) - break - - smb_temp = read_smb(file_name) - smb_temp = np.where(params["mask"], 0, smb_temp) - - avg_smb_timeseries[year_to_read - 1] = np.round( - net_avrg(smb_temp) * params["mm_to_Gt"], - 2, - ) - count_yr = count_yr + 1 - - if count_yr == params["climo_nyears"]: - break - - del smb_temp - - first_year = year_to_read - - print("number of years used in climatology = ", count_yr) - return first_year, avg_smb_timeseries - - def plot_line(data, time, line, color, label, linewidth): plt.plot( time, From 5b6e3d565dc910284b29dbc6577b45f840033df1 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Fri, 26 Jul 2024 17:42:22 -0600 Subject: [PATCH 08/11] More xarray to replace netcdf4 Read in 40 years of SMB data with read_cesm_smb(), and then use xarray functions (basically taking means or sums over different dimensions) to compute the different climatologies. This greatly reduces the number of functions needed in utils - it's just read_cesm_smb() and some plotting routines. --- .../nblibrary/glc/LIWG_SMB_diagnostic.ipynb | 47 ++++---- examples/nblibrary/glc/utils.py | 107 ++++-------------- 2 files changed, 44 insertions(+), 110 deletions(-) diff --git a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb index 1731701..d661686 100644 --- a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb +++ b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb @@ -172,13 +172,15 @@ "outputs": [], "source": [ "# creating the SMB climatology for new case\n", - "smb_case_climo = utils.create_climo(case_path, case_name, last_year, params)\n", + "smb_case = utils.read_cesm_smb(case_path, case_name, last_year, params)\n", + "smb_case_climo = smb_case.mean(\"time\")\n", "\n", "# creating the SMB climatology for base_case\n", "if base_case_name:\n", - " smb_base_climo = utils.create_climo(\n", + " smb_base_case = utils.read_cesm_smb(\n", " base_case_path, base_case_name, base_last_year, params\n", - " )" + " )\n", + " smb_base_climo = smb_base_case.mean(\"time\")" ] }, { @@ -201,18 +203,19 @@ ")\n", "\n", "# Initializing the glacier ID variable\n", - "smb_obs_climo = np.zeros(grid_dims)\n", + "smb_obs_climo = xr.DataArray(np.zeros(grid_dims), dims=[\"glc1Exp_ny\", \"glc1Exp_nx\"])\n", "\n", "# Performing the interpolation\n", "for j in range(grid_dims[0]):\n", " point_y = np.zeros(grid_dims[1])\n", " point_y[:] = thk_init_da[\"y1\"].data[j]\n", " pts = (thk_init_da[\"x1\"].data[:], point_y[:])\n", - " smb_obs_climo[j, :] = myInterpFunction_smb_obs(pts)\n", + " smb_obs_climo.data[j, :] = myInterpFunction_smb_obs(pts)\n", "\n", "# Filtering out fill values\n", - "smb_obs_climo = np.where(smb_obs_climo > 1e20, 0, smb_obs_climo)\n", - "smb_obs_climo = np.where(mask, 0, smb_obs_climo)" + "smb_obs_climo.data = np.where(\n", + " np.logical_or(mask, smb_obs_climo > 1e20), 0, smb_obs_climo\n", + ")" ] }, { @@ -223,13 +226,13 @@ "outputs": [], "source": [ "# Integrated SMB time series\n", - "first_year, avg_smb_case_climo = utils.compute_annual_climo(\n", - " case_path, case_name, last_year, params\n", - ")\n", + "first_year = last_year - len(smb_case[\"time\"]) + 1\n", + "avg_smb_case_climo = smb_case.sum([\"glc1Exp_ny\", \"glc1Exp_nx\"]) * params[\"mm_to_Gt\"]\n", "\n", "if base_case_name:\n", - " base_first_year, avg_smb_base_case_climo = utils.compute_annual_climo(\n", - " base_case_path, base_case_name, base_last_year, params\n", + " base_first_year = base_last_year - len(smb_base_case[\"time\"]) + 1\n", + " avg_smb_base_case_climo = (\n", + " smb_base_case.sum([\"glc1Exp_ny\", \"glc1Exp_nx\"]) * params[\"mm_to_Gt\"]\n", " )" ] }, @@ -284,9 +287,8 @@ ")\n", "\n", "## Right panel\n", - "utils.plot_contour_diff(\n", - " smb_case_climo,\n", - " smb_obs_climo,\n", + "utils.plot_contour(\n", + " smb_case_climo - smb_obs_climo,\n", " fig,\n", " ax[2],\n", " \"SMB bias (mm/yr w.e.)\",\n", @@ -342,9 +344,8 @@ " )\n", "\n", " ## Right panel\n", - " utils.plot_contour_diff(\n", - " smb_case_climo,\n", - " smb_base_climo,\n", + " utils.plot_contour(\n", + " smb_case_climo - smb_base_climo,\n", " fig,\n", " ax[2],\n", " \"SMB difference (mm/yr w.e.)\",\n", @@ -382,10 +383,10 @@ "if base_case_name:\n", " avg_smb_base_timeseries = np.zeros(base_nt)\n", "\n", - "avg_smb_obs_timeseries[:] = np.round(utils.net_avrg(smb_obs_climo) * mm_to_Gt, 2)\n", - "avg_smb_case_timeseries[:] = np.round(utils.net_avrg(smb_case_climo) * mm_to_Gt, 2)\n", + "avg_smb_obs_timeseries[:] = np.round(smb_obs_climo.sum() * mm_to_Gt, 2)\n", + "avg_smb_case_timeseries[:] = np.round(smb_case_climo.sum() * mm_to_Gt, 2)\n", "if base_case_name:\n", - " avg_smb_base_timeseries[:] = np.round(utils.net_avrg(smb_base_climo) * mm_to_Gt, 2)\n", + " avg_smb_base_timeseries[:] = np.round(smb_base_climo.sum() * mm_to_Gt, 2)\n", "\n", "\n", "x_ticks = np.arange(first_year, last_year + 2, 5)\n", @@ -402,7 +403,7 @@ "# Plotting annual / spatial means\n", "plt.subplot(111)\n", "utils.plot_line(\n", - " avg_smb_case_climo[first_year - 1 : :],\n", + " avg_smb_case_climo,\n", " time,\n", " line=\"-\",\n", " color=\"blue\",\n", @@ -419,7 +420,7 @@ ")\n", "if base_case_name:\n", " utils.plot_line(\n", - " avg_smb_base_case_climo[base_first_year - 1 : :],\n", + " avg_smb_base_case_climo,\n", " base_time,\n", " line=\"-\",\n", " color=\"red\",\n", diff --git a/examples/nblibrary/glc/utils.py b/examples/nblibrary/glc/utils.py index 391d643..4bc7cbf 100644 --- a/examples/nblibrary/glc/utils.py +++ b/examples/nblibrary/glc/utils.py @@ -5,59 +5,13 @@ import numpy as np import xarray as xr from matplotlib import pyplot as plt -from netCDF4 import Dataset -def set_plot_prop_clean(ax): - """ - This function cleans up the figures from unnecessary default figure properties. - - """ - ax.invert_yaxis() - ax.set_xlabel("") - ax.set_ylabel("") - ax.set_xticklabels("") - ax.set_yticklabels("") - ax.set_xticks([]) - ax.set_yticks([]) - - -def rmse(prediction, target): +def read_cesm_smb(path, case_name, last_year, params): """ - This function returns the root mean square error for the SMB. - Input: - prediction = field to predict - target = field to compare with the prediction + This function reads CESM coupler history files and returns + an xarray DataArray containing surface mass balance in units mm/y """ - return np.sqrt(((prediction - target) ** 2).mean()) - - -def net_avrg(data): - """ - This function returns the net average of a data field - """ - return np.sum(np.sum(data, axis=0), axis=0) - - -def read_smb(file): - """ - This function reads the CISM SMB dataset from a CESM simulation output - in the cpl directory. The output is adjusted to be converted to mm/yr w.e unit. - - Input: - file: name of the file to extract the SMB - """ - rhoi = 917 # ice density kg/m3 - sec_in_yr = 60 * 60 * 24 * 365 # seconds in a year - smb_convert = sec_in_yr / rhoi * 1000 # converting kg m-2 s-1 ice to mm y-1 w.e. - - nid = Dataset(file, "r") - smb_cism = np.squeeze(nid.variables["glc1Exp_Flgl_qice"][0, :, :]) * smb_convert - nid.close() - return smb_cism - - -def _get_cesm_output(path, case_name, last_year, params): # Set parameters rhoi = 917 # ice density kg/m3 sec_in_yr = 60 * 60 * 24 * 365 # seconds in a year @@ -94,33 +48,28 @@ def _get_cesm_output(path, case_name, last_year, params): return climo_out -def create_climo(path, case_name, last_year, params): - - climo_out = _get_cesm_output(path, case_name, last_year, params) - - # Averaging the climo data - return climo_out.mean("time").data +# -------------------- # +# PLOTTING FUNCTIONS # +# -------------------- # -def compute_annual_climo(path, case_name, last_year, params): - # Initializing a field for the climatology - avg_smb_timeseries = np.zeros(last_year) - climo_out = _get_cesm_output(path, case_name, last_year, params) - for k in range(len(climo_out["time"])): - # index into avg_smb_timeseries; want largest k to index (last_year-1) - kk = last_year - (len(climo_out["time"]) - k) - # note that mm_to_Gt has 1 / res**2 factor, so we want a sum rather than mean - avg_smb_timeseries[kk] = np.round( - climo_out.isel(time=k).sum() * params["mm_to_Gt"], - 2, - ).data +def set_plot_prop_clean(ax): + """ + This function cleans up the figures from unnecessary default figure properties. - return last_year - len(climo_out["time"]) + 1, avg_smb_timeseries + """ + ax.invert_yaxis() + ax.set_xlabel("") + ax.set_ylabel("") + ax.set_xticklabels("") + ax.set_yticklabels("") + ax.set_xticks([]) + ax.set_yticks([]) -def plot_contour(data, fig, ax, title, vmin, vmax, cmap, mm_to_Gt): - avg_data = np.round(net_avrg(data) * mm_to_Gt, 2) - last_panel0 = ax.imshow(data[:, :], vmin=vmin, vmax=vmax, cmap=cmap) +def plot_contour(da, fig, ax, title, vmin, vmax, cmap, mm_to_Gt): + avg_data = np.round(da.sum().data * mm_to_Gt, 2) + last_panel0 = ax.imshow(da.data[:, :], vmin=vmin, vmax=vmax, cmap=cmap) ax.set_title(title, fontsize=16) set_plot_prop_clean(ax) ax.annotate("net avg =" + str(avg_data) + " Gt/yr", xy=(5, 5), fontsize=16) @@ -132,22 +81,6 @@ def plot_contour(data, fig, ax, title, vmin, vmax, cmap, mm_to_Gt): cbar.ax.tick_params(labelsize=16) -def plot_contour_diff(data_new, data_old, fig, ax, title, vmin, vmax, cmap, mm_to_Gt): - avg_data = np.round(net_avrg(data_new - data_old) * mm_to_Gt, 2) - last_panel2 = ax.imshow(data_new - data_old, vmin=vmin, vmax=vmax, cmap=cmap) - - ax.set_title(title, fontsize=16) - set_plot_prop_clean(ax) - - ax.annotate("net avg =" + str(avg_data) + " Gt/yr", xy=(5, 5), fontsize=16) - - pos = ax.get_position() - cax = fig.add_axes([0.89, pos.y0, 0.02, pos.y1 - pos.y0]) - - cbar = fig.colorbar(last_panel2, cax=cax) - cbar.ax.tick_params(labelsize=16) - - def plot_line(data, time, line, color, label, linewidth): plt.plot( time, From f714646f21d4ec697957959bd1f9894c533fd6b7 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Tue, 30 Jul 2024 16:36:56 -0600 Subject: [PATCH 09/11] Add docstrings in utils.py Also renamed data -> da in plot_line() to be consistent with plot_contour() --- examples/nblibrary/glc/utils.py | 29 ++++++++++++++++++++++++++--- 1 file changed, 26 insertions(+), 3 deletions(-) diff --git a/examples/nblibrary/glc/utils.py b/examples/nblibrary/glc/utils.py index 4bc7cbf..b93300d 100644 --- a/examples/nblibrary/glc/utils.py +++ b/examples/nblibrary/glc/utils.py @@ -56,7 +56,6 @@ def read_cesm_smb(path, case_name, last_year, params): def set_plot_prop_clean(ax): """ This function cleans up the figures from unnecessary default figure properties. - """ ax.invert_yaxis() ax.set_xlabel("") @@ -68,6 +67,19 @@ def set_plot_prop_clean(ax): def plot_contour(da, fig, ax, title, vmin, vmax, cmap, mm_to_Gt): + """ + Plot a contour map of surface mass balance (assumed to be in da.data). + Also computes global mean, in Gt, and prints average in lower left corner. + Arguments: + da - xr.DataArray containing SMB in units of mm/yr + fig - matplotlib.figure.Figure + ax - matplotlib.axes.Axes + title - string containing title of plot + vmin - minimum value for contours + vmax - maximum value for contours + cmap - matplotlib.colors.Colormap + mm_to_Gt - conversion factor for mm/yr -> Gt/yr + """ avg_data = np.round(da.sum().data * mm_to_Gt, 2) last_panel0 = ax.imshow(da.data[:, :], vmin=vmin, vmax=vmax, cmap=cmap) ax.set_title(title, fontsize=16) @@ -81,10 +93,21 @@ def plot_contour(da, fig, ax, title, vmin, vmax, cmap, mm_to_Gt): cbar.ax.tick_params(labelsize=16) -def plot_line(data, time, line, color, label, linewidth): +def plot_line(da, time, line, color, label, linewidth): + """ + Plot a time series of spatially averaged surface mass balance (assumed to + be in da.data). + Arguments: + da - xr.DataArray containing spatially averaged SMB in units of Gt/yr + time - np.array containing time dimension + line - style of line to use in plot + color - color of line in plot + label - label of line in legend + linewidth - thickness of line in plot + """ plt.plot( time, - data, + da, line, ms=3, mfc=color, From b80a8c5649b94394d257f759cd8fd9af592ed2fc Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Tue, 30 Jul 2024 17:09:24 -0600 Subject: [PATCH 10/11] Fix colorbar in plot_contour() Every panel was overwriting the location of the colorbar, so only the bias contour levels were shown. Now the first two plots (with same vmin and vmax) have a colorbar between them while bias plot colorbar is off to the right --- examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb | 15 ++++++++++++++- examples/nblibrary/glc/utils.py | 9 +++++---- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb index d661686..bc123c9 100644 --- a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb +++ b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb @@ -274,6 +274,7 @@ " smb_case_climo,\n", " fig,\n", " ax[0],\n", + " 0.35,\n", " f\"{case_name}\\nSMB (mm/y w.e.)\\nMean from {first_year:04d} - {last_year:04d}\",\n", " vmin,\n", " vmax,\n", @@ -283,7 +284,15 @@ "\n", "## Center panel\n", "utils.plot_contour(\n", - " smb_obs_climo, fig, ax[1], \"SMB Obs\\n(mm/y w.e.)\", vmin, vmax, my_cmap, mm_to_Gt\n", + " smb_obs_climo,\n", + " fig,\n", + " ax[1],\n", + " 0.35,\n", + " \"SMB Obs\\n(mm/y w.e.)\",\n", + " vmin,\n", + " vmax,\n", + " my_cmap,\n", + " mm_to_Gt,\n", ")\n", "\n", "## Right panel\n", @@ -291,6 +300,7 @@ " smb_case_climo - smb_obs_climo,\n", " fig,\n", " ax[2],\n", + " 0.89,\n", " \"SMB bias (mm/yr w.e.)\",\n", " vmin,\n", " vmax,\n", @@ -324,6 +334,7 @@ " smb_case_climo,\n", " fig,\n", " ax[0],\n", + " 0.35,\n", " f\"{case_name}\\nSMB (mm/y w.e.)\\nMean from {first_year:04d} - {last_year:04d}\",\n", " vmin,\n", " vmax,\n", @@ -336,6 +347,7 @@ " smb_base_climo,\n", " fig,\n", " ax[1],\n", + " 0.35,\n", " f\"{base_case_name}\\nSMB (mm/y w.e.)\\nMean from {base_first_year:04d} - {base_last_year:04d}\",\n", " vmin,\n", " vmax,\n", @@ -348,6 +360,7 @@ " smb_case_climo - smb_base_climo,\n", " fig,\n", " ax[2],\n", + " 0.89,\n", " \"SMB difference (mm/yr w.e.)\",\n", " vmin,\n", " vmax,\n", diff --git a/examples/nblibrary/glc/utils.py b/examples/nblibrary/glc/utils.py index b93300d..2a11bc3 100644 --- a/examples/nblibrary/glc/utils.py +++ b/examples/nblibrary/glc/utils.py @@ -66,7 +66,7 @@ def set_plot_prop_clean(ax): ax.set_yticks([]) -def plot_contour(da, fig, ax, title, vmin, vmax, cmap, mm_to_Gt): +def plot_contour(da, fig, ax, left, title, vmin, vmax, cmap, mm_to_Gt): """ Plot a contour map of surface mass balance (assumed to be in da.data). Also computes global mean, in Gt, and prints average in lower left corner. @@ -74,6 +74,7 @@ def plot_contour(da, fig, ax, title, vmin, vmax, cmap, mm_to_Gt): da - xr.DataArray containing SMB in units of mm/yr fig - matplotlib.figure.Figure ax - matplotlib.axes.Axes + left - left dimension of rect (dimensions for colorbar) title - string containing title of plot vmin - minimum value for contours vmax - maximum value for contours @@ -81,15 +82,15 @@ def plot_contour(da, fig, ax, title, vmin, vmax, cmap, mm_to_Gt): mm_to_Gt - conversion factor for mm/yr -> Gt/yr """ avg_data = np.round(da.sum().data * mm_to_Gt, 2) - last_panel0 = ax.imshow(da.data[:, :], vmin=vmin, vmax=vmax, cmap=cmap) + last_panel = ax.imshow(da.data[:, :], vmin=vmin, vmax=vmax, cmap=cmap) ax.set_title(title, fontsize=16) set_plot_prop_clean(ax) ax.annotate("net avg =" + str(avg_data) + " Gt/yr", xy=(5, 5), fontsize=16) pos = ax.get_position() - cax = fig.add_axes([0.35, pos.y0, 0.02, pos.y1 - pos.y0]) + cax = fig.add_axes([left, pos.y0, 0.02, pos.y1 - pos.y0]) - cbar = fig.colorbar(last_panel0, cax=cax) + cbar = fig.colorbar(last_panel, cax=cax) cbar.ax.tick_params(labelsize=16) From e6246e4d8171cd158e8d265496b0e572d0684738 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Tue, 30 Jul 2024 19:49:11 -0600 Subject: [PATCH 11/11] Fix misleading comment --- examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb index bc123c9..64a356c 100644 --- a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb +++ b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb @@ -316,7 +316,7 @@ "metadata": {}, "outputs": [], "source": [ - "# Comparing SMB new run vs obs\n", + "# Comparing SMB new run vs base case\n", "if base_case_name:\n", " # Colormap choice\n", " my_cmap = mcm.get_cmap(\"Spectral\")\n",