add time consistency

notebooks/wp5/era5_nino.ipynb

@@ -34,10 +34,10 @@
   },
   {
    "cell_type": "markdown",
-   "id": "f271ebf8",
+   "id": "6337fe6b",
    "metadata": {},
    "source": [
-    "## Import packages"
+    "## Define Parameters"
    ]
   },
   {
@@ -185,7 +185,7 @@
     "        )\n",
     "    )\n",
     "df = xr.concat(datasets, \"reduction\").to_pandas()\n",
-    "df"
+    "df.T"
    ]
   },
   {

notebooks/wp5/era5_time_consistency.ipynb

@@ -0,0 +1,317 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "e4dfa17b",
+   "metadata": {},
+   "source": [
+    "# Time consistency in ERA5"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d91979f4",
+   "metadata": {},
+   "source": [
+    "## Import packages"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "745d15bc-3f3e-4d69-8925-963a86f2ebe0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import xarray as xr\n",
+    "from c3s_eqc_automatic_quality_control import diagnostics, download\n",
+    "\n",
+    "plt.style.use(\"seaborn-v0_8-notebook\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a19f5d8a",
+   "metadata": {},
+   "source": [
+    "## Define Parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3a1c2535-0087-4eb9-8ab9-b0e9c17e01d7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Time period\n",
+    "start = \"1940-01\"\n",
+    "stop = None  # None: present"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1997f010",
+   "metadata": {},
+   "source": [
+    "## Define request"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a30542c2-07d2-44bf-9755-1e07bccee3c7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "collection_id = \"reanalysis-era5-pressure-levels-monthly-means\"\n",
+    "variable = \"temperature\"\n",
+    "request = {\n",
+    "    \"format\": \"grib\",\n",
+    "    \"product_type\": \"monthly_averaged_reanalysis\",\n",
+    "    \"variable\": [\n",
+    "        \"temperature\",\n",
+    "        \"u_component_of_wind\",\n",
+    "        \"vertical_velocity\",\n",
+    "        \"geopotential\",\n",
+    "        \"ozone_mass_mixing_ratio\",\n",
+    "        \"relative_humidity\",\n",
+    "        \"fraction_of_cloud_cover\",\n",
+    "    ],\n",
+    "    \"pressure_level\": [\n",
+    "        \"1\",\n",
+    "        \"2\",\n",
+    "        \"3\",\n",
+    "        \"5\",\n",
+    "        \"7\",\n",
+    "        \"10\",\n",
+    "        \"20\",\n",
+    "        \"30\",\n",
+    "        \"50\",\n",
+    "        \"70\",\n",
+    "        \"100\",\n",
+    "        \"125\",\n",
+    "        \"150\",\n",
+    "        \"175\",\n",
+    "        \"200\",\n",
+    "        \"225\",\n",
+    "        \"250\",\n",
+    "        \"300\",\n",
+    "        \"350\",\n",
+    "        \"400\",\n",
+    "        \"450\",\n",
+    "        \"500\",\n",
+    "        \"550\",\n",
+    "        \"600\",\n",
+    "        \"650\",\n",
+    "        \"700\",\n",
+    "        \"750\",\n",
+    "        \"775\",\n",
+    "        \"800\",\n",
+    "        \"825\",\n",
+    "        \"850\",\n",
+    "        \"875\",\n",
+    "        \"900\",\n",
+    "        \"925\",\n",
+    "        \"950\",\n",
+    "        \"975\",\n",
+    "        \"1000\",\n",
+    "    ],\n",
+    "    \"time\": \"00:00\",\n",
+    "}\n",
+    "requests = download.update_request_date(request, start=start, stop=stop)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "dba2c821",
+   "metadata": {},
+   "source": [
+    "## Download and transform"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1a587d84-5c3f-4e14-9ad5-204cc12c8ef7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds = download.download_and_transform(\n",
+    "    collection_id,\n",
+    "    requests,\n",
+    "    transform_func=diagnostics.spatial_weighted_mean,\n",
+    "    chunks={\"year\": 1, \"variable\": 1},\n",
+    ")\n",
+    "# Convert plev to hPa\n",
+    "with xr.set_options(keep_attrs=True):\n",
+    "    ds[\"plev\"] = ds[\"plev\"] / 100\n",
+    "    ds[\"z\"] /= 9.8\n",
+    "ds[\"plev\"].attrs[\"units\"] = \"hPa\"\n",
+    "ds[\"z\"].attrs[\"units\"] = \"m\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "458d21a2",
+   "metadata": {},
+   "source": [
+    "## Compute anomaly"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0c8f9ea5-95c4-4daf-8efd-a896f2d23d0a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "group = \"forecast_reference_time.month\"\n",
+    "with xr.set_options(keep_attrs=True):\n",
+    "    ds_anoma = ds.groupby(group) - ds.groupby(group).mean()\n",
+    "for varname, da in ds_anoma.data_vars.items():\n",
+    "    da.attrs[\"long_name\"] += \" anomaly\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "983c8b3b",
+   "metadata": {},
+   "source": [
+    "## Show min and max values"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "59a35bd1-c25e-40f6-8ff2-ed57f04c0bcf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "datasets = []\n",
+    "for reduction in (\"min\", \"max\"):\n",
+    "    datasets.append(\n",
+    "        getattr(ds_anoma.sel(plev=slice(1000, 10)), reduction)().expand_dims(\n",
+    "            reduction=[reduction]\n",
+    "        )\n",
+    "    )\n",
+    "df = xr.concat(datasets, \"reduction\").to_pandas()\n",
+    "df.T"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1ca34cc2",
+   "metadata": {},
+   "source": [
+    "## Plot Hovmöller diagrams"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f2cf69ed-fd1a-4e82-8fc8-162ac04cf36a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_dict = {\n",
+    "    \"t\": {\"levels\": np.arange(-5, 5.5, 0.5), \"cmap\": \"RdBu_r\", \"yscale\": \"log\"},\n",
+    "    \"u\": {\"levels\": np.arange(-20, 20 + 2, 2), \"cmap\": \"PuOr\", \"yscale\": \"log\"},\n",
+    "    \"w\": {\n",
+    "        \"levels\": np.arange(-0.8e-3, 0.8e-3 + 0.8e-4, 0.8e-4),\n",
+    "        \"cmap\": \"PuOr\",\n",
+    "        \"yscale\": \"linear\",\n",
+    "    },\n",
+    "    \"z\": {\"levels\": np.arange(-300, 300 + 30, 30), \"cmap\": \"seismic\", \"yscale\": \"log\"},\n",
+    "    \"o3\": {\n",
+    "        \"levels\": np.arange(-1.0e-6, 1.0e-6 + 1.0e-7, 1.0e-7),\n",
+    "        \"cmap\": \"RdGy_r\",\n",
+    "        \"yscale\": \"log\",\n",
+    "    },\n",
+    "    \"r\": {\"levels\": np.arange(-5, 5 + 0.5, 0.5), \"cmap\": \"BrBG\", \"yscale\": \"linear\"},\n",
+    "    \"cc\": {\n",
+    "        \"levels\": np.arange(-0.01, 0.01 + 0.001, 0.001),\n",
+    "        \"cmap\": \"PRGn\",\n",
+    "        \"yscale\": \"linear\",\n",
+    "    },\n",
+    "}\n",
+    "\n",
+    "eruptions = [\n",
+    "    {\n",
+    "        \"volcano\": \"Pinatubo\",\n",
+    "        \"date\": \"1991-06-15\",\n",
+    "    },\n",
+    "    {\n",
+    "        \"volcano\": \"El Chichón, Mount St. Helens\",\n",
+    "        \"date\": \"1981-01-01\",\n",
+    "        \"linestyle\": \"--\",\n",
+    "    },\n",
+    "    {\n",
+    "        \"volcano\": \"Agung\",\n",
+    "        \"date\": \"1963-02-24\",\n",
+    "    },\n",
+    "    {\n",
+    "        \"volcano\": \"Bezymianny\",\n",
+    "        \"date\": \"1956-01-01\",\n",
+    "    },\n",
+    "]\n",
+    "\n",
+    "zooms = {\n",
+    "    \"Pinatubo\": slice(\"1988-01-01\", \"1995-12-31\"),\n",
+    "    \"Agung\": slice(\"1962-01-01\", \"1967-12-31\"),\n",
+    "    \"El Chicon\": slice(\"1980-01-01\", \"1985-12-31\"),\n",
+    "    \"Entire period\": slice(\"1940-01-01\", \"2022-12-31\"),\n",
+    "}\n",
+    "\n",
+    "for title, zoom in zooms.items():\n",
+    "    ds_to_plot = ds_anoma.sel(forecast_reference_time=zoom)\n",
+    "    fig, axs = plt.subplots(len(ds_anoma.data_vars), figsize=(15, 24))\n",
+    "    for (varname, kwargs), ax in zip(plot_dict.items(), axs):\n",
+    "        ds_to_plot[varname].plot.contourf(\n",
+    "            y=\"plev\", yincrease=False, extend=\"both\", ax=ax, **kwargs\n",
+    "        )\n",
+    "        for eruption in eruptions:\n",
+    "            date = np.datetime64(eruption[\"date\"])\n",
+    "            if (\n",
+    "                not ds_to_plot[\"forecast_reference_time\"].min()\n",
+    "                <= date\n",
+    "                <= ds_to_plot[\"forecast_reference_time\"].max()\n",
+    "            ):\n",
+    "                continue\n",
+    "            ax.axvline(date, color=\"black\", linestyle=eruption.get(\"linestyle\"), lw=0.5)\n",
+    "            ax.text(\n",
+    "                date,\n",
+    "                ax.get_ylim()[0],\n",
+    "                \" \" + eruption[\"volcano\"],\n",
+    "                rotation=\"vertical\",\n",
+    "                ha=\"left\",\n",
+    "                va=\"bottom\",\n",
+    "            )\n",
+    "    fig.suptitle(title, fontsize=16, x=0.45, y=0.92)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}