notebook/3_forecast_preprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5fd5f453-13c1-47c1-ad27-8af80b47ee08",
   "metadata": {},
   "source": [
    "# Forecast preprocessing\n",
    "***\n",
    "\n",
    "**Author**: Chus Casado<br>\n",
    "**Date**: 27-07-2023<br>\n",
    "\n",
    "**Introduction**:<br>\n",
    "This code computes the probability of the forecasted discharge of exceeding a threshold associated with a specific return period.\n",
    "\n",
    "The input data are:\n",
    "* The discharge forecast for the complete set of reporting points, numerical weather prediction (NWP) model and the complete study period. This data is saved in NetCDF format in a hard disk; due to its size it cannot be included in the GitHub repository.\n",
    "* The discharge thresholds associated to each reporting point, from which a specific return period (by default 5 years) will be used.\n",
    "\n",
    "The output is the probability of exceeding the specified return period. This is a matrix of multiple dimensions: station, NWP model, date-time, and lead time. Actually, 2 matrixes are computed, one related to the probability of exceeding the specified threshold ($Q_{rp}$), and another with the probability of exceeding a slighly lower threshold ($0.95\\cdot Q_{rp}$). This lower threshold is used to avoid false positives or false negatives in forecasts very close to the observation, but on opposite sides of the discharge threshold.\n",
    "\n",
    "The results are saved inside the repository as NetCDF, saving one file per reporting point.\n",
    "\n",
    "**Questions**:<br>\n",
    "\n",
    "\n",
    "**Tasks to do**:<br>\n",
    "* [ ] Probably, the `discharge2exceedance` function could be improved so the two discharge thresholds are computed at the same time. Since the function needs to load all the forecast discharge records, it is a waste of time the current procedure, in which all the records need to be loaded twice (once for each discharge threshold).\n",
    "\n",
    "**Interesting links**<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a4ee7a0f-f62b-4a8b-8147-c46a2e2471ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "path_root = os.getcwd()\n",
    "import glob\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "from datetime import datetime, timedelta\n",
    "import time\n",
    "from tqdm.notebook import tqdm\n",
    "import yaml\n",
    "from pathlib import Path\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "os.chdir('../py/')\n",
    "from config import Config\n",
    "from compute import *\n",
    "os.chdir(path_root)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c93f9bd-6b3d-4ee8-a815-4c5048c13ac3",
   "metadata": {},
   "source": [
    "## 1 Configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "de2a6451-65a1-4750-b10c-38d3aab506a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "config_path = Path('../conf')\n",
    "config = Config.load_from_yaml(config_path / 'config_NWP_all_leadtimes.yml')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce48789e-5e03-40f8-8e4f-44f44cf21679",
   "metadata": {},
   "source": [
    "### 1.1 Reporting points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0f5de8cd-69dd-4632-9a33-66368ee90c11",
   "metadata": {},
   "outputs": [],
   "source": [
    "# area threshold\n",
    "area_threshold = config.reporting_points['area']\n",
    "\n",
    "# reporting points\n",
    "path_stations = config.reporting_points['output']\n",
    "file_stations = path_stations / f'reporting_points_over_{area_threshold}km2.parquet'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b7e82e9-4934-4409-b30d-6eea222000d4",
   "metadata": {},
   "source": [
    "### 1.2 Discharge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "83f90b1a-a354-456c-877c-d34230d34b06",
   "metadata": {},
   "outputs": [],
   "source": [
    "# local directory where I have saved the raw discharge data\n",
    "path_in = config.discharge['input']['forecast']\n",
    "\n",
    "# start and end of the study period\n",
    "start = config.discharge['study_period']['start']\n",
    "end = config.discharge['study_period']['end']\n",
    "    \n",
    "# return period\n",
    "rp = config.discharge['return_period']['threshold']\n",
    "# percent buffer over the discharge threshold\n",
    "reducing_factor = config.discharge['return_period']['reducing_factor']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d51149b7-3593-439a-976d-0ad24789c3ed",
   "metadata": {},
   "source": [
    "### 1.3 Exceedance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "94be1247-0e98-477a-bd5b-2e36e5e364db",
   "metadata": {},
   "outputs": [],
   "source": [
    "# path where the output exceedance datasets will be saved\n",
    "path_out = config.exceedance['output']['forecast'] / f'{rp}'\n",
    "path_out.mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b1ba771-7ccf-4837-a6c0-32e084881f1b",
   "metadata": {},
   "source": [
    "## 2 Data\n",
    "\n",
    "### 2.1 Reporting points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4802b22f-39eb-4bea-84fb-307639be5091",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no. stations:\t2371\n"
     ]
    }
   ],
   "source": [
    "# load table of fixed reporing points\n",
    "stations = pd.read_parquet(file_stations)\n",
    "\n",
    "print('no. stations:\\t{0}'.format(stations.shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5de69ffb-357e-4629-bf20-63ffa0cebb14",
   "metadata": {},
   "source": [
    "### 2.2 Discharge forecast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e2a7d7fe-2b33-41db-9d66-344913ea1a12",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "COS:\t2006 files\n",
      "DWD:\t2006 files\n",
      "EUD:\t2006 files\n",
      "EUE:\t2006 files\n"
     ]
    }
   ],
   "source": [
    "# list files\n",
    "fore_files = {model: [] for model in list(models)}\n",
    "for year in np.arange(start.year, end.year + 1):\n",
    "    for month in range(1, 13):    \n",
    "        # list files\n",
    "        for model in models:\n",
    "            fore_files[model] += glob.glob(f'{path_in}/{model}/{year}/{month:02d}/*.nc')\n",
    "\n",
    "# count files and check if all are avaible\n",
    "n_files = pd.Series(data=[len(fore_files[model]) for model in models], index=models)\n",
    "\n",
    "# list of forecast from the beginning to the end of the data\n",
    "start, end = datetime(1900, 1, 1), datetime(2100, 1, 1)\n",
    "for model in models:\n",
    "    st, en = [datetime.strptime(fore_files[model][step][-13:-3], '%Y%m%d%H') for step in [0, -1]]\n",
    "    start = max(st, start)\n",
    "    end = min(en, end)\n",
    "dates = pd.date_range(start, end, freq='12h')\n",
    "\n",
    "# find missing files\n",
    "if any(n_files != len(dates)):\n",
    "    missing = {}\n",
    "    for model in models:\n",
    "        filedates = [datetime.strptime(file[-13:-3], '%Y%m%d%H') for file in fore_files[model]]    \n",
    "        missing[model] = [date for date in dates if date not in filedates]\n",
    "    print('mising files:', missing)\n",
    "\n",
    "# trim files to the period where all models are available\n",
    "for model in models:\n",
    "    fore_files[model] = [file for file in fore_files[model] if start <= datetime.strptime(file[-13:-3], '%Y%m%d%H') <= end]\n",
    "    print('{0}:\\t{1} files'.format(model, len(fore_files[model])))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a6337e9-e2d2-4c4d-87fd-51786396de62",
   "metadata": {},
   "source": [
    "## 3 Analysis: exceedance probability\n",
    "\n",
    "This section will iteratively (station by station) load all the available forecast and compute the probability of exceeding the discharge threshold for each of the meteorological forcings. The result will be a NetCDF file for each station that contains the exceedance probability. These files will be later used in the skill assessment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "cff2b77d-9afc-4920-921c-adf2df95880d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no. new stations:\t\t\t0\n"
     ]
    }
   ],
   "source": [
    "# select stations that haven't been processed before\n",
    "files = glob.glob(f'{path_out}*.nc')\n",
    "if len(files) > 0:\n",
    "    old_stations = [int(file.split('\\\\')[-1].split('.')[0]) for file in files]\n",
    "    new_stations = set(stations.index).difference(old_stations)\n",
    "    stations = stations.loc[new_stations]\n",
    "print('no. new stations:\\t\\t\\t{0}'.format(stations.shape[0]))\n",
    "\n",
    "# generate a DataArray with the discharge threshold of the stations in the catchment\n",
    "thresholds = xr.DataArray(stations[f'rl{rp}'], dims='id', coords={'id': stations.index.astype(str).tolist()})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "556a40aa-1d60-462f-a217-0153b4f8465c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8a00baca84e7404e99c507b8c9704b14",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/4 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d2a3cefac2f34224adf950af908db1f2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2006 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "COS\tE:/casadje/Documents/EFASv4/forecast/COS/2023/06\\COS2023063012.nc\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ad0268516b9d4414af0c43a003e8c420",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2006 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DWD\tE:/casadje/Documents/EFASv4/forecast/DWD/2023/06\\DWD2023063012.nc\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7ea7fee8f8b0482b8bf159c4cfa599f1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2006 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EUD\tE:/casadje/Documents/EFASv4/forecast/EUD/2023/06\\EUD2023063012.nc\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "eed20552f9ef40efb42ceb59274cecac",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2006 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EUE\tE:/casadje/Documents/EFASv4/forecast/EUE/2023/06\\EUE2023063012.nc\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b5c041b4227546e1bab981ea398c5540",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/4 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "680c94b298474251abe86c456ac6494a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2006 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "COS\tE:/casadje/Documents/EFASv4/forecast/COS/2023/06\\COS2023063012.nc\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "85bb061caa3c4a6399cf21fbc1f3dac9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2006 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DWD\tE:/casadje/Documents/EFASv4/forecast/DWD/2023/06\\DWD2023063012.nc\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e94d424178dd43eaa58da999c3a9eb0c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2006 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EUD\tE:/casadje/Documents/EFASv4/forecast/EUD/2023/06\\EUD2023063012.nc\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f3f833e34fe542ee91139b69cf09ab3d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2006 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EUE\tE:/casadje/Documents/EFASv4/forecast/EUE/2023/06\\EUE2023063012.nc\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2e7d5cbb66784aa3b8ab6e8185b6d228",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2371 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "excecution time: 1254.4 s\n"
     ]
    }
   ],
   "source": [
    "start = time.perf_counter()\n",
    "\n",
    "# exceedances over Q5\n",
    "exceedance_q5 = discharge2exceedance(fore_files, thresholds) # * (1 + reducing_factor))\n",
    "# exceedances over 95% of the Q5\n",
    "exceedance_buffer = discharge2exceedance(fore_files, thresholds * (1 - reducing_factor))\n",
    "# join both DataArrays in a single Dataset\n",
    "exceedance = xr.Dataset({'high': exceedance_q5, 'low': exceedance_buffer})\n",
    "\n",
    "# export files station by station\n",
    "for stn in tqdm(exceedance.id.data):\n",
    "    file = f'{stn:>04}.nc'\n",
    "    if file in os.listdir(path_out):\n",
    "        print(f'File {file} already exists')\n",
    "        continue\n",
    "    else:\n",
    "        exceedance.sel(id=stn).to_netcdf(path_out / file)\n",
    "        \n",
    "end = time.perf_counter()\n",
    "\n",
    "print('excecution time: {0:.1f} s'.format(end - start))"
   ]
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}