Refactor nmse_PSL to work with CESM2 f19 output

1. To work with CESM2 output, needed to include a fix_time_dim() function that sets time to the averaging midpoint rather than the endpoint. 2. Output that has been been regridded in post-processing is in proc/tseries/ rather than proc/tseries/regrid/, so user needs to control that from config.yml 3. Runs on the fv1.9x2.5 grid need to compare against coarser observational data (so the location of validation data is set in config.yml)
mnlevy1981 · Aug 7, 2024 · a38b0b0 · a38b0b0
1 parent 98fd753
commit a38b0b0
Show file tree

Hide file tree

Showing 4 changed files with 118 additions and 111 deletions.
diff --git a/examples/key_metrics/config.yml b/examples/key_metrics/config.yml
@@ -99,7 +99,7 @@ compute_notebooks:
   # parameters are specified. Several examples of different
   # types of notebooks are provided.
 
-  # The first key (here simple_no_params_nb) is the name of the
+  # The first key (here infrastructure) is the name of the
   # notebook from nb_path_root, minus the .ipynb
 
     infrastructure:
@@ -111,6 +111,8 @@ compute_notebooks:
       nmse_PSL:
         parameter_groups:
           none:
+            regridded_output: True
+            validation_path: '/glade/campaign/cesm/development/cross-wg/diagnostic_framework/nmse_validation/fv0.9x1.25'
             start_date: '0001-01-01'
             end_date: '0101-01-01'
 

diff --git a/examples/nblibrary/atm/nmse_PSL.ipynb b/examples/nblibrary/atm/nmse_PSL.ipynb
@@ -1,5 +1,22 @@
 {
  "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "3f230d52-dca7-4ce4-98cc-6267fc04893d",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "source": [
+    "# Normalized Mean Square Error\n",
+    "\n",
+    "This notebook computes the normalized mean square error of atmospheric surface pressure.\n",
+    "It is compared to ERA5 observations, as well as the CESM2 large ensemble and CMIP6 model output."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -13,16 +30,14 @@
    },
    "outputs": [],
    "source": [
-    "import xarray as xr\n",
-    "import numpy as np\n",
-    "import matplotlib.pyplot as plt\n",
-    "from nmse_utils import *\n",
-    "from averaging_utils import *\n",
     "import glob\n",
     "\n",
-    "import warnings\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import xarray as xr\n",
     "\n",
-    "warnings.filterwarnings(\"ignore\")"
+    "from nmse_utils import nmse\n",
+    "from averaging_utils import seasonal_climatology_weighted"
    ]
   },
   {
@@ -36,7 +51,9 @@
     "tags": []
    },
    "source": [
-    "### Parameters"
+    "## Parameters\n",
+    "\n",
+    "These variables are set in `config.yml`"
    ]
   },
   {
@@ -54,41 +71,17 @@
    },
    "outputs": [],
    "source": [
-    "sname = \"\"\n",
     "CESM_output_dir = \"\"\n",
     "case_name = \"\"\n",
     "start_date = \"\"\n",
-    "end_date = \"\""
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "14b554f3-40e1-4c12-9f26-c1684ac03fbf",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "print(start_date)\n",
-    "print(end_date)\n",
-    "print(CESM_output_dir)\n",
-    "print(case_name)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "80edad94-426a-4304-9003-78e12960613a",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "dat = xr.open_mfdataset(\n",
-    "    CESM_output_dir + \"/\" + case_name + \"/atm/proc/tseries/regrid/*PSL*.nc\"\n",
-    ")"
+    "end_date = \"\"\n",
+    "validation_path = \"\"\n",
+    "regridded_output = False"
    ]
   },
   {
    "cell_type": "markdown",
-   "id": "e0527e3e-cd26-46b5-8c1e-08882109e12e",
+   "id": "74c7803f-a8c5-445d-9233-0aa2663c58bd",
    "metadata": {
     "editable": true,
     "slideshow": {
@@ -97,13 +90,13 @@
     "tags": []
    },
    "source": [
-    "### Read in validation data and other CMIP models for comparison (precomputed)"
+    "## Read in the current case"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "126e65b3-2b8c-400c-af02-2ad0b0f82e6e",
+   "id": "ccca8e3a-a52f-4202-9704-9d4470eda984",
    "metadata": {
     "editable": true,
     "slideshow": {
@@ -113,92 +106,116 @@
    },
    "outputs": [],
    "source": [
-    "# This will need to be changed if we have a common repo\n",
-    "validation_path = (\n",
-    "    \"/glade/campaign/cgd/cas/islas/python_savs/CUPID/NMSE/validation_data_example/\"\n",
-    ")\n",
+    "def fix_time_dim(dat):\n",
+    "    \"\"\"CESM2 output sets time as the end of the averaging interval (e.g. January average is midnight on February 1st);\n",
+    "    This function sets the time dimension to the midpoint of the averaging interval.\n",
+    "    Note that CESM3 output sets time to the midpoint already, so this function should not change CESM3 data.\"\"\"\n",
+    "    if \"time\" not in dat.dims:\n",
+    "        return dat\n",
+    "    if \"bounds\" not in dat.time.attrs:\n",
+    "        return dat\n",
+    "    time_bounds_avg = dat[dat.time.attrs[\"bounds\"]].mean(\"nbnd\")\n",
+    "    time_bounds_avg.attrs = dat.time.attrs\n",
+    "    dat = dat.assign_coords({\"time\": time_bounds_avg})\n",
+    "    return xr.decode_cf(dat)\n",
     "\n",
-    "# ---ERA5\n",
-    "era5 = xr.open_dataset(validation_path + \"PSL_ERA5.nc\")\n",
-    "era5 = era5 / 100.0  # convert to hPa\n",
     "\n",
-    "# ---CESM2\n",
-    "lens2 = xr.open_dataset(validation_path + \"PSL_LENS2.nc\")\n",
-    "lens2 = lens2 / 100.0  # convert to hPa\n",
-    "lens2[\"lon\"] = era5.lon\n",
-    "lens2[\"lat\"] = era5.lat\n",
+    "if regridded_output:\n",
+    "    file_path = f\"{CESM_output_dir}/{case_name}/atm/proc/tseries/regrid\"\n",
+    "else:\n",
+    "    file_path = f\"{CESM_output_dir}/{case_name}/atm/proc/tseries\"\n",
     "\n",
-    "# ---CMIP6\n",
-    "modelfiles = sorted(glob.glob(validation_path + \"CMIP6/*.nc\"))\n",
-    "datcmip6 = [xr.open_dataset(ifile).mean(\"M\") for ifile in modelfiles]\n",
-    "datcmip6 = xr.concat(datcmip6, dim=\"model\")\n",
-    "datcmip6 = datcmip6 / 100.0\n",
-    "datcmip6[\"lon\"] = era5.lon\n",
-    "datcmip6[\"lat\"] = era5.lat"
+    "dat = (\n",
+    "    fix_time_dim(xr.open_mfdataset(f\"{file_path}/*PSL*.nc\", decode_times=False))\n",
+    "    .sel(time=slice(start_date, end_date))\n",
+    "    .PSL\n",
+    "    / 100.0\n",
+    ")"
    ]
   },
   {
-   "cell_type": "markdown",
-   "id": "ab5f1441-4b6b-4386-b3aa-a0f6fa9109b8",
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "073a2ad0-81e6-4817-9024-4b9b718fabb4",
    "metadata": {
     "editable": true,
     "slideshow": {
      "slide_type": ""
     },
     "tags": []
    },
-   "source": [
-    "### Read in the current case"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "2b4883ef-f3a7-458b-b96e-bdb79d5abdd6",
-   "metadata": {},
    "outputs": [],
    "source": [
-    "dat"
+    "# --Compute seasonal and annual means\n",
+    "dat = seasonal_climatology_weighted(dat).load()"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "09397e30-62b9-474b-a47c-4a1174137bce",
+   "cell_type": "markdown",
+   "id": "e0527e3e-cd26-46b5-8c1e-08882109e12e",
    "metadata": {
     "editable": true,
     "slideshow": {
      "slide_type": ""
     },
     "tags": []
    },
-   "outputs": [],
    "source": [
-    "dat = dat.sel(time=slice(start_date, end_date))\n",
-    "dat = dat.PSL / 100.0"
+    "## Read in validation data and other CMIP models for comparison (precomputed)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "b11ed849-c8b2-4079-bde0-598f60ab4463",
-   "metadata": {},
+   "id": "1ff152b1-2168-4a0d-826b-cf8d11f66ab7",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
    "outputs": [],
    "source": [
-    "# this requires the same grid\n",
-    "dat[\"lon\"] = era5.lon\n",
-    "dat[\"lat\"] = era5.lat"
+    "# Ensure all validation datasets have the same coordinates as the ERA5 data\n",
+    "# (Avoid round-off level differences since all data should be on the same grid)\n",
+    "lon = dat.lon.data\n",
+    "lat = dat.lat.data"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "073a2ad0-81e6-4817-9024-4b9b718fabb4",
-   "metadata": {},
+   "id": "126e65b3-2b8c-400c-af02-2ad0b0f82e6e",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
    "outputs": [],
    "source": [
-    "# --Compute seasonal and annual means\n",
-    "dat = seasonal_climatology_weighted(dat).load()"
+    "# ---ERA5\n",
+    "era5 = xr.open_dataset(f\"{validation_path}/PSL_ERA5.nc\").assign_coords(\n",
+    "    {\"lon\": lon, \"lat\": lat}\n",
+    ")\n",
+    "era5 = era5 / 100.0  # convert to hPa\n",
+    "\n",
+    "# ---CESM2\n",
+    "lens2 = xr.open_dataset(f\"{validation_path}/PSL_LENS2.nc\").assign_coords(\n",
+    "    {\"lon\": lon, \"lat\": lat}\n",
+    ")\n",
+    "lens2 = lens2 / 100.0  # convert to hPa\n",
+    "\n",
+    "# ---CMIP6\n",
+    "modelfiles = sorted(glob.glob(f\"{validation_path}/CMIP6/*.nc\"))\n",
+    "datcmip6 = [\n",
+    "    xr.open_dataset(ifile).assign_coords({\"lon\": lon, \"lat\": lat}).mean(\"M\")\n",
+    "    for ifile in modelfiles\n",
+    "]\n",
+    "datcmip6 = xr.concat(datcmip6, dim=\"model\")\n",
+    "datcmip6 = datcmip6 / 100.0"
    ]
   },
   {
@@ -212,7 +229,7 @@
     "tags": []
    },
    "source": [
-    "### Compute the NMSE"
+    "## Compute the NMSE"
    ]
   },
   {
@@ -400,20 +417,6 @@
     "    fontsize=14,\n",
     ")"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "d01c7860-b4ca-4ab4-a05f-3d0f248a56c2",
-   "metadata": {
-    "editable": true,
-    "slideshow": {
-     "slide_type": ""
-    },
-    "tags": []
-   },
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {

diff --git a/examples/nblibrary/atm/nmse_utils.py b/examples/nblibrary/atm/nmse_utils.py
@@ -11,10 +11,6 @@ def nmse(obs, mod):
     from that
     """
 
-    # make sure the model and obs have the same lons and lats
-    mod["lon"] = obs.lon
-    mod["lat"] = obs.lat
-
     # get the weights and weight by zero if the model or obs is nan
     w = np.cos(np.deg2rad(obs.lat))
     w = w.expand_dims({"lon": obs.lon}, axis=1)

diff --git a/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb b/examples/nblibrary/glc/LIWG_SMB_diagnostic.ipynb
@@ -3,7 +3,13 @@
   {
    "cell_type": "markdown",
    "id": "e45c723f-aa6f-4db4-a197-492132cc8156",
-   "metadata": {},
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
    "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",
@@ -62,17 +68,17 @@
    "source": [
     "# Parameter Defaults\n",
     "\n",
-    "CESM_output_dir = \"/glade/campaign/cesm/development/cross-wg/diagnostic_framework/CESM_output_for_testing\"\n",
-    "case_name = \"b.e23_alpha17f.BLT1850.ne30_t232.092\"  # case name\n",
-    "climo_nyears = 40  # number of years to compute the climatology\n",
-    "last_year = 101\n",
+    "CESM_output_dir = \"\"\n",
+    "case_name = \"\"  # case name\n",
+    "climo_nyears = 0  # number of years to compute the climatology\n",
+    "last_year = 0\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\""
+    "obs_path = \"\"  # directory containing observed dataset\n",
+    "obs_name = \"\"  # file name for observed dataset"
    ]
   },
   {