From a10eea378f69882bbb03669232f32b1422a6b5b8 Mon Sep 17 00:00:00 2001 From: Mattia Almansi Date: Mon, 17 Jul 2023 10:06:49 +0200 Subject: [PATCH] add time consistency --- notebooks/wp5/era5_nino.ipynb | 6 +- notebooks/wp5/era5_time_consistency.ipynb | 317 ++++++++++++++++++++++ 2 files changed, 320 insertions(+), 3 deletions(-) create mode 100644 notebooks/wp5/era5_time_consistency.ipynb diff --git a/notebooks/wp5/era5_nino.ipynb b/notebooks/wp5/era5_nino.ipynb index 0012776..e440570 100644 --- a/notebooks/wp5/era5_nino.ipynb +++ b/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" ] }, { diff --git a/notebooks/wp5/era5_time_consistency.ipynb b/notebooks/wp5/era5_time_consistency.ipynb new file mode 100644 index 0000000..2168de1 --- /dev/null +++ b/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 +}