{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6836fb72",
   "metadata": {},
   "source": [
    "# Eumetnet postprocessing benchmark dataset demo"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d98b3f5",
   "metadata": {},
   "source": [
    "Installation with [pip](https://pypi.org/) (if needed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53a60658",
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip install matplotlib\n",
    "#!pip install climetlab\n",
    "#!pip install climetlab-eumetnet-postprocessing-benchmark"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4c72959",
   "metadata": {},
   "source": [
    "Loading climetlab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "arabic-sally",
   "metadata": {},
   "outputs": [],
   "source": [
    "import climetlab as cml"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e2b0b17",
   "metadata": {},
   "source": [
    "and matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9919adef",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2b449d3",
   "metadata": {},
   "source": [
    "## Example with a 2 metre temperature field deterministic forecast"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f28463c1",
   "metadata": {},
   "source": [
    "Download of a deterministic high-resolution forecast (the 2 metre temperature)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "affiliated-binding",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = cml.load_dataset(\n",
    "    \"eumetnet-postprocessing-benchmark-training-data-gridded-forecasts-surface\",\n",
    "    date=\"2017-12-02\",\n",
    "    parameter=\"2t\",\n",
    "    kind=\"highres\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a534b8a0",
   "metadata": {},
   "source": [
    "The `ds` object returned is a climetlab datasource that can be converted to various format.\n",
    "Note that the data are cached temporarily on your disk so next time you ask for this forecast I might be faster."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e68721a",
   "metadata": {},
   "source": [
    "For example, conversion to a [xarray](http://xarray.pydata.org/en/stable/index.html) :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "boring-bryan",
   "metadata": {},
   "outputs": [],
   "source": [
    "fcs = ds.to_xarray()\n",
    "fcs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f125e024",
   "metadata": {},
   "source": [
    "Plotting using climetlab and [magics](https://github.com/ecmwf/magics-python) (still in development)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "416c2ff7",
   "metadata": {},
   "outputs": [],
   "source": [
    "cml.plot_map(ds[5], # plot of the 5th time step\n",
    "             foreground=dict(\n",
    "                 map_grid=False,\n",
    "                 map_label=False,\n",
    "                 map_grid_frame=True,\n",
    "                 map_grid_frame_thickness=5,\n",
    "                 map_boundaries=True,\n",
    "                ),\n",
    "            )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b856799d",
   "metadata": {},
   "source": [
    "Retrieving the observations corresponding to the forecasts (in xarray format)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b24131d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = ds.get_observations_as_xarray()\n",
    "obs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afdd5ee9",
   "metadata": {},
   "source": [
    "Plotting the observations with magics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73a56068",
   "metadata": {},
   "outputs": [],
   "source": [
    "cml.plot_map(obs.t2m[0,0,5,0], # plot of the 5th time step\n",
    "             foreground=dict(\n",
    "                 map_grid=False,\n",
    "                 map_label=False,\n",
    "                 map_grid_frame=True,\n",
    "                 map_grid_frame_thickness=5,\n",
    "                 map_boundaries=True,\n",
    "                ),\n",
    "            )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b6113aa",
   "metadata": {},
   "source": [
    "Compute the difference between the forecast and the observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1a89682",
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = fcs - obs\n",
    "diff"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f6a00ec",
   "metadata": {},
   "source": [
    "and plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72923337",
   "metadata": {},
   "outputs": [],
   "source": [
    "cml.plot_map(diff.t2m[0,0,5,0], # plot of the 5th time step\n",
    "             foreground=dict(\n",
    "                 map_grid=False,\n",
    "                 map_label=False,\n",
    "                 map_grid_frame=True,\n",
    "                 map_grid_frame_thickness=5,\n",
    "                 map_boundaries=True,\n",
    "                ),\n",
    "            )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7f58a5c",
   "metadata": {},
   "source": [
    "Save the forecast and observations on the hard drive as a netcdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abc0dd18",
   "metadata": {},
   "outputs": [],
   "source": [
    "fcs.to_netcdf('fcs_hr_2017-12-02.nc')\n",
    "obs.to_netcdf('obs_2017-12-02.nc')\n",
    "!ls"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bd57285",
   "metadata": {},
   "source": [
    "Plotting the observations and the forecast at 4 different points (forecasts are solid lines, observations are dashed lines)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2b85aa2",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "points = [(20, 20), (20, 70), (70, 20), (70, 70)]\n",
    "\n",
    "plt.figure(figsize=(10, 8))\n",
    "\n",
    "for point in points:\n",
    "    line = fcs.t2m.isel(latitude=point[0], longitude=point[1], time=0, surface=0, number=0).plot()\n",
    "    obs.t2m.isel(latitude=point[0], longitude=point[1], time=0, surface=0, number=0).plot(ls='--', color=line[0].get_color())\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27fd272d",
   "metadata": {},
   "source": [
    "## Example with a total precipitation accumulated over the last 6 hours ensemble forecast"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f126e517",
   "metadata": {},
   "source": [
    "Download of an ensemble forecast (the precipitation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "604c29ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = cml.load_dataset(\n",
    "    \"eumetnet-postprocessing-benchmark-training-data-gridded-forecasts-surface-processed\",\n",
    "    date=\"2017-12-02\",\n",
    "    parameter=\"tp\",\n",
    "    kind=\"ensemble\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b955005",
   "metadata": {},
   "outputs": [],
   "source": [
    "fcs = ds.to_xarray()\n",
    "fcs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6903dbf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = ds.get_observations_as_xarray()\n",
    "obs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be1218d4",
   "metadata": {},
   "source": [
    "Plotting the observations and forecast ensemble mean at 4 different points (forecasts are solid lines, observations are dashed lines)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f09a24c",
   "metadata": {},
   "outputs": [],
   "source": [
    "points = [(20, 20), (20, 70), (70, 20), (70, 70)]\n",
    "\n",
    "plt.figure(figsize=(10, 8))\n",
    "\n",
    "for point in points:\n",
    "    line = fcs.tp6.isel(latitude=point[0], longitude=point[1], time=0, surface=0).mean('number').plot()    \n",
    "    obs.tp6.isel(latitude=point[0], longitude=point[1], time=0, surface=0, number=0).plot(ls='--', color=line[0].get_color())\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3b5e595",
   "metadata": {},
   "source": [
    "## Example with a 2 metre temperature field ensemble reforecast"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a361f7eb",
   "metadata": {},
   "source": [
    "Download of an ensemble reforecast (the 2 metre temperature)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ec13b0a",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "ds = cml.load_dataset(\n",
    "    \"eumetnet-postprocessing-benchmark-training-data-gridded-reforecasts-surface\",\n",
    "    date=\"2017-12-28\",\n",
    "    parameter=\"2t\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5632efd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fcs = ds.to_xarray()\n",
    "fcs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "687a954b",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = ds.get_observations_as_xarray()\n",
    "obs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "665b9acf",
   "metadata": {},
   "source": [
    "Plotting the observations and forecast ensemble mean at 4 different points (forecasts are solid lines, observations are dashed lines) **in 1997 with the 2017 model version** !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e1d2c90",
   "metadata": {},
   "outputs": [],
   "source": [
    "points = [(20, 20), (20, 70), (70, 20), (70, 70)]\n",
    "\n",
    "plt.figure(figsize=(10, 8))\n",
    "\n",
    "for point in points:\n",
    "    line = fcs.t2m.isel(latitude=point[0], longitude=point[1], time=0, surface=0).mean('number').plot()\n",
    "    obs.t2m.isel(latitude=point[0], longitude=point[1], time=0, surface=0, number=0).plot(ls='--', color=line[0].get_color())\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c58bb851",
   "metadata": {},
   "source": [
    "## Example with a total precipitation accumulated over the last 6 hours ensemble reforecast"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2917a81d",
   "metadata": {},
   "source": [
    "Download of an ensemble reforecast (the precipitation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "728e80b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = cml.load_dataset(\n",
    "    \"eumetnet-postprocessing-benchmark-training-data-gridded-reforecasts-surface-processed\",\n",
    "    date=\"2017-12-28\",\n",
    "    parameter=\"tp\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4ae99dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "fcs = ds.to_xarray()\n",
    "fcs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa9e3514",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = ds.get_observations_as_xarray()\n",
    "obs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da0f4d49",
   "metadata": {},
   "source": [
    "Plotting the observations and forecast ensemble mean at 4 different points (forecasts are solid lines, observations are dashed lines) **in 1997 with the 2017 model version** !!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "133a38fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "points = [(20, 20), (20, 70), (70, 20), (70, 70)]\n",
    "\n",
    "plt.figure(figsize=(10, 8))\n",
    "\n",
    "for point in points:\n",
    "    line = fcs.tp6.isel(latitude=point[0], longitude=point[1], time=0, surface=0).mean('number').plot()        \n",
    "    obs.tp6.isel(latitude=point[0], longitude=point[1], time=0, surface=0, number=0).plot(ls='--', color=line[0].get_color())\n",
    "    "
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}