spatial_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4d047299",
   "metadata": {
    "papermill": {
     "duration": 0.027065,
     "end_time": "2023-03-06T02:49:49.549053",
     "exception": false,
     "start_time": "2023-03-06T02:49:49.521988",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Spatial UNSEEN analysis\n",
    "\n",
    "This notebook performs a spatial analysis of a climate hazard using the UNprecedented Simulated Extremes using ENsembles (UNSEEN) method for the Australian Climate Service (ACS). \n",
    "\n",
    "The UNSEEN method uses model ensembles to generate large samples of plausible events to estimate the likelihood of extreme events. \n",
    "\n",
    "\n",
    "## Notes\n",
    "\n",
    "- The notebook displays the output of a single model and a single climate metric.\n",
    "- This notebook is setup to be run using papermill, which allow the user to specify parameters that will be imported into the notebook (in a new cell after tagged parameters).\n",
    "- The cell tagged as parameters contains default parameters that can be overwritten by the user.\n",
    "- Multiple versions of most plots are generated: \n",
    "    - with and without stippling in regions where the model forecast is significantly different from the observations (i.e., stippling shows where the model and obs distributions are significantly different). However, only the masked versions of each plot will be displayed in the notebook.\n",
    "    - original and/or additive/multiplicative bias corrected \n",
    "- Some of the output plot files names are defined as input arguments (namely the independence and similarity plots), while the majority are generated in the spatial_plots functions (the model forecast filename stem is appended to the plot filename and saved in the user defined figure directory).\n",
    "\n",
    "## Prerequisites\n",
    "\n",
    "- There are a number of prerequisite data files including the climate metric (observations and model), independence and similarity test results and nonstationary GEV parameters. \n",
    "- Requires local copies of [spatial_plots.py](https://github.com/AusClimateService/unseen-projects/blob/master/spatial_plots.py) and [plotting_maps](https://github.com/AusClimateService/plotting_maps) (updated to allow kw norm). \n",
    "    - To avoid errors importing the local modules, either:\n",
    "        - Add the paths to your conda environment using: conda develop /path/to/plotting_maps and conda develop /path/to/unseen-projects/code\n",
    "        - Run (or add to your ~/.bashrc): export PYTHONPATH=\"/path/to/module/:$PYTHONPATH\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bafc012",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Parameters\n",
    "----------\n",
    "metric : str\n",
    "    Metric to analyse (e.g. \"txx\", \"tnn\", \"rx1day\")\n",
    "var : str\n",
    "    Variable name of the metric in the dataset\n",
    "model_name : str\n",
    "    Name of the model\n",
    "obs_name : str\n",
    "    Name of the observational dataset for plot labels\n",
    "time_agg : {\"maximum\", \"minimum\"}, default=\"maximum\"\n",
    "    The method for aggregating the time dimension of the metric to examine extremes (e.g., \"maximum\" for TXx or Rx5day and \"minimum\" for TNn)\n",
    "plot_dict : dict or str\n",
    "    Dictionary of metric plot labels, ticks and colormaps with keys:\n",
    "        'metric' (e.g., 'TXx'), 'var', 'var_name', 'units', 'units_label', 'freq', 'obs_name', 'cmap', 'cmap_anom', 'ticks', 'ticks_anom', 'ticks_param_trend', 'cbar_extend', 'agcd_mask'\n",
    "fig_dir : str\n",
    "    Directory to save the figures\n",
    "reference_time_period : str, optional\n",
    "    Reference time period for the metric in observations (e.g., \"1961-06-30 2023-06-30\")\n",
    "min_lead_spatial_agg : {None, \"mean\", \"median\", \"max\"}\n",
    "    Spatial aggregation for determining the minimum independent lead time\n",
    "shapefile : str, optional\n",
    "    Path to shapefile for used if the the minimum independent lead time is based on spatial aggregation, rather than grid point specific\n",
    "shape_overlap : float, optional\n",
    "    Minimum overlap between shapefile and grid cell (see `shapefile` parameter)\n",
    "similarity_test : {\"ks\", \"ad\"}, default=\"ks\"\n",
    "    Statistical test testing the similarity between model and observations\n",
    "alpha : float, default=0.05\n",
    "    Significance level for similarity test\n",
    "covariate_base : int\n",
    "    Base year to substitute in the non-stationary GEV model (i.e., the year for which AEP is relative to)\n",
    "gev_trend_period : str or list\n",
    "    Start and end year for determining the trend in the AEP (e.g., '[1961, 2025]')\n",
    "plot_additive_bc : bool, default=True\n",
    "    Plot the bias-corrected model data using an additive bias correction. If True, then then all associated prerequisite files must be provided. Note that the similarity test plots will still be plotted if the data file is available\n",
    "plot_multiplicative_bc : bool, default=False\n",
    "    Plot the multiplicative bias-corrected model data\n",
    "obs_file : str\n",
    "    Filename of the metric based on observations\n",
    "model_file : str\n",
    "    Filename of metric from the model\n",
    "model_mulc_bc_file : str, optional\n",
    "    Filename of the additive bias-corrected metric\n",
    "model_mulc_bc_file : str, optional\n",
    "    Filename of the multiplicative bias-corrected metric\n",
    "independence_file : str\n",
    "    PFilename of the independence test file\n",
    "independence_plot : str\n",
    "    Filename for the minimum independent lead time plot\n",
    "similarity_raw_file : str\n",
    "    Filename of the similarity test file based on uncorrected model data\n",
    "similarity_mulc_bc_file : str, optional\n",
    "    Filename of the similarity test file based on additive bias-corrected model data\n",
    "similarity_mulc_bc_file : str, optional\n",
    "    Filename of the similarity test file based on multiplicative bias-corrected model data\n",
    "similarity_raw_plot : str\n",
    "    Filename of the similarity test plot based on uncorrected model data\n",
    "similarity_add_bc_plot : str, optional\n",
    "    Filename of the similarity test plot based on additive bias-corrected model data\n",
    "similarity_mulc_bc_plot : str, optional\n",
    "    Filename of the similarity test plot based on multiplicative bias-corrected model data\n",
    "gev_params_nonstationary_file : str\n",
    "    Filename of the model non-stationary GEV parameters\n",
    "gev_params_nonstationary_add_bc_file : str, optional\n",
    "    Filename of the additive bias-corrected model non-stationary GEV parameters\n",
    "gev_params_nonstationary_mulc_bc_file : str, optional\n",
    "    Filename of the multiplicative bias-corrected model non-stationary GEV parameters\n",
    "\"\"\"\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import cartopy.crs as ccrs\n",
    "from cartopy.mpl.gridliner import LatitudeFormatter, LongitudeFormatter\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import AutoMinorLocator\n",
    "import numpy as np\n",
    "import os\n",
    "from pathlib import Path\n",
    "import xarray as xr\n",
    "\n",
    "from unseen import fileio, independence, similarity, eva, time_utils, spatial_selection\n",
    "from acs_plotting_maps import cmap_dict, tick_dict, regions_dict  # NOQA\n",
    "import spatial_plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56e9e0ab",
   "metadata": {
    "papermill": {
     "duration": 0.024169,
     "end_time": "2023-03-06T02:49:59.159559",
     "exception": false,
     "start_time": "2023-03-06T02:49:59.135390",
     "status": "completed"
    },
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "# Optional parameters\n",
    "alpha = 0.05\n",
    "ensemble_dim = \"ensemble\"\n",
    "init_dim = \"init_date\"\n",
    "lead_dim = \"lead_time\"\n",
    "time_dim = \"time\"\n",
    "lat_dim = \"lat\"\n",
    "lon_dim = \"lon\"\n",
    "similarity_test = \"ks\"\n",
    "time_agg = \"maximum\"\n",
    "plot_additive_bc = True\n",
    "plot_multiplicative_bc = False\n",
    "reference_time_period = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9787631",
   "metadata": {
    "papermill": {
     "duration": 0.042552,
     "end_time": "2023-03-06T02:49:59.256849",
     "exception": false,
     "start_time": "2023-03-06T02:49:59.214297",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Required parameters\n",
    "kwargs = locals()\n",
    "assert \"metric\" in kwargs, \"Must provide a metric name\"\n",
    "assert \"var\" in kwargs, \"Must provide a variable name\"\n",
    "assert \"model_name\" in kwargs, \"Must provide a model name\"\n",
    "\n",
    "if isinstance(reference_time_period, str):\n",
    "    reference_time_period = list(reference_time_period.split(\" \"))\n",
    "assert \"covariate_base\" in kwargs, \"Must provide a nonstationary GEV covariate base year\"\n",
    "assert (\n",
    "    time_agg in spatial_plots.func_dict\n",
    "), f\"Invalid time aggregation method: {time_agg} (options: {list(spatial_plots.func_dict.keys())})\"\n",
    "\n",
    "# Format parameters passed as strings\n",
    "assert \"gev_trend_period\" in kwargs, \"Must provide a GEV trend period (e.g., '[1981, 2010]')\"\n",
    "if isinstance(gev_trend_period, str):\n",
    "    # Convert string to list\n",
    "    gev_trend_period = eval(gev_trend_period)\n",
    "    print(f\"gev_trend_period: {gev_trend_period}\")\n",
    "\n",
    "assert \"plot_dict\" in kwargs, \"Must provide spatial plot dictionary of labels, ticks and colormaps\"\n",
    "if isinstance(plot_dict, str):\n",
    "    # Convert string to dictionary and check for required keys\n",
    "    plot_dict = eval(plot_dict)\n",
    "\n",
    "    for key in [\n",
    "        \"metric\",\n",
    "        \"var\",\n",
    "        \"var_name\",\n",
    "        \"units\",\n",
    "        \"units_label\",\n",
    "        \"freq\",\n",
    "        \"cmap\",\n",
    "        \"cmap_anom\",\n",
    "        \"ticks\",\n",
    "        \"ticks_anom\",\n",
    "        \"ticks_param_trend\",\n",
    "        \"cbar_extend\",\n",
    "        \"agcd_mask\",\n",
    "    ]:\n",
    "        assert key in plot_dict, f\"Missing key: {key} in plot_dict\"\n",
    "    plot_dict[\"agcd_mask\"] = False  # AGCD data/grid only\n",
    "\n",
    "assert os.path.isfile(obs_file), f\"Must provide an observations data file (papermill option -p obs_file [filepath])\"\n",
    "assert os.path.isfile(model_file), \"Must provide an model data file (papermill option -p model_file [filepath])\"\n",
    "assert os.path.isfile(\n",
    "    independence_file\n",
    "), \"Must provide an independence min lead file (papermill option -p independence_file [filepath])\"\n",
    "assert os.path.isfile(\n",
    "    similarity_raw_file\n",
    "), \"Must provide an raw data similarity test file (papermill option -p similarity_raw_file [filepath])\"\n",
    "assert os.path.isfile(\n",
    "    gev_params_nonstationary_file\n",
    "), \"Must provide a non-stationary GEV parameter file (papermill option -p gev_params_nonstationary_file [filepath])\"\n",
    "\n",
    "# Bias corrected files\n",
    "if plot_additive_bc:\n",
    "    assert os.path.isfile(\n",
    "        similarity_add_bc_file\n",
    "    ), \"Must provide an additive bias corrected similarity test file (papermill option -p similarity_add_bias_file [filepath])\"\n",
    "    assert os.path.isfile(\n",
    "        model_add_bc_file\n",
    "    ), \"Must provide a model additive bias corrected data file (papermill option -p model_add_bc_file [filepath])\"\n",
    "    assert os.path.isfile(gev_params_nonstationary_add_bc_file)\n",
    "\n",
    "if plot_multiplicative_bc:\n",
    "    assert os.path.isfile(\n",
    "        similarity_mulc_bc_file\n",
    "    ), \"Must provide an multiplicative bias corrected similarity test file (papermill option -p similarity_mulc_bias_file [filepath])\"\n",
    "    assert os.path.isfile(\n",
    "        model_mulc_bc_file\n",
    "    ), \"Must provide a model multiplicative bias corrected data file (papermill option -p model_mulc_bc_file [filepath])\"\n",
    "    assert os.path.isfile(gev_params_nonstationary_mulc_bc_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d082e6af",
   "metadata": {
    "papermill": {
     "duration": 0.014695,
     "end_time": "2023-03-06T02:50:00.583605",
     "exception": false,
     "start_time": "2023-03-06T02:50:00.568910",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Model ensemble\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "146854fe",
   "metadata": {
    "papermill": {
     "duration": 0.036872,
     "end_time": "2023-03-06T02:50:12.012408",
     "exception": false,
     "start_time": "2023-03-06T02:50:11.975536",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### Model data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1620b927",
   "metadata": {
    "papermill": {
     "duration": 0.093555,
     "end_time": "2023-03-06T02:50:12.140348",
     "exception": false,
     "start_time": "2023-03-06T02:50:12.046793",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "model_ds = fileio.open_dataset(model_file)\n",
    "model_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58d7d91d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert event_time date strings to datetime or cftime\n",
    "try:\n",
    "    event_times = np.vectorize(time_utils.str_to_cftime)(model_ds.event_time, model_ds.time.dt.calendar)\n",
    "    model_ds[\"event_time\"] = (model_ds.event_time.dims, event_times)\n",
    "except ValueError:\n",
    "    model_ds[\"event_time\"] = model_ds.event_time.astype(dtype=\"datetime64[ns]\").compute()\n",
    "model_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32fa28e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_additive_bc:\n",
    "    model_add_bc_ds = fileio.open_dataset(model_add_bc_file)\n",
    "\n",
    "    # Convert event_time date strings to cftime\n",
    "    try:\n",
    "        event_times = np.vectorize(time_utils.str_to_cftime)(\n",
    "            model_add_bc_ds.event_time, model_add_bc_ds.time.dt.calendar\n",
    "        )\n",
    "        model_add_bc_ds[\"event_time\"] = (model_add_bc_ds.event_time.dims, event_times)\n",
    "    except ValueError:\n",
    "        model_add_bc_ds[\"event_time\"] = model_add_bc_ds.event_time.astype(dtype=\"datetime64[ns]\").compute()\n",
    "    model_add_bc_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "197a8470",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_multiplicative_bc:\n",
    "    model_mulc_bc_ds = fileio.open_dataset(model_mulc_bc_file)\n",
    "    # Convert event_time date strings to cftime\n",
    "    try:\n",
    "        event_times = np.vectorize(time_utils.str_to_cftime)(\n",
    "            model_mulc_bc_ds.event_time, model_mulc_bc_ds.time.dt.calendar\n",
    "        )\n",
    "        model_mulc_bc_ds[\"event_time\"] = (model_mulc_bc_ds.event_time.dims, event_times)\n",
    "    except ValueError:\n",
    "        model_mulc_bc_ds[\"event_time\"] = model_mulc_bc_ds.event_time.astype(dtype=\"datetime64[ns]\").compute()\n",
    "    model_mulc_bc_ds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b25bd815",
   "metadata": {
    "papermill": {
     "duration": 0.034589,
     "end_time": "2023-03-06T02:50:12.345116",
     "exception": false,
     "start_time": "2023-03-06T02:50:12.310527",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### Independence tests\n",
    "Perform independence testing to identify which lead times can be considered independent samples. We plot spatial maps of the ensemble mean correlation coefficients for each lead time, the confidence interval bounds and the first lead time in which the correlation is within the bounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bb73ae9",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_independence = xr.open_dataset(independence_file, use_cftime=True)\n",
    "ds_independence"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6febf422",
   "metadata": {},
   "source": [
    "Correlation coefficients for each initialisation month and lead time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f2e209e",
   "metadata": {},
   "outputs": [],
   "source": [
    "cm = ds_independence.r.plot(\n",
    "    col=lead_dim,\n",
    "    row=\"month\",\n",
    "    subplot_kws=dict(projection=ccrs.PlateCarree()),\n",
    "    transform=ccrs.PlateCarree(),\n",
    "    add_colorbar=False,\n",
    "    cmap=plt.cm.seismic,\n",
    ")\n",
    "\n",
    "# Fix hidden axis ticks and labels\n",
    "for i, ax in enumerate(cm.axs.flat):\n",
    "    ax.coastlines()\n",
    "    ax.xaxis.set_major_formatter(LongitudeFormatter())\n",
    "    ax.yaxis.set_major_formatter(LatitudeFormatter())\n",
    "    ax.xaxis.set_minor_locator(AutoMinorLocator())\n",
    "    ax.yaxis.set_minor_locator(AutoMinorLocator())\n",
    "    ax.set_xlabel(None)\n",
    "    ax.set_ylabel(None)\n",
    "\n",
    "    subplotspec = ax.get_subplotspec()\n",
    "    ax.xaxis.set_visible(True)\n",
    "    if subplotspec.is_first_col():\n",
    "        ax.yaxis.set_visible(True)\n",
    "\n",
    "cm.fig.set_constrained_layout(True)\n",
    "cm.fig.get_layout_engine().set(h_pad=0.2)\n",
    "cm.add_colorbar(pad=0.02)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4284747",
   "metadata": {},
   "source": [
    "The lower and upper null correlation bounds at each grid point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93367a2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "cm = ds_independence.ci.plot(\n",
    "    col=\"quantile\",\n",
    "    row=\"month\",\n",
    "    subplot_kws=dict(projection=ccrs.PlateCarree()),\n",
    "    transform=ccrs.PlateCarree(),\n",
    "    add_colorbar=False,\n",
    "    cmap=plt.cm.seismic,\n",
    ")\n",
    "# Fix hidden axis ticks\n",
    "for i, ax in enumerate(cm.axs.flat):\n",
    "    ax.coastlines()\n",
    "    ax.xaxis.set_major_formatter(LongitudeFormatter())\n",
    "    ax.yaxis.set_major_formatter(LatitudeFormatter())\n",
    "    ax.xaxis.set_minor_locator(AutoMinorLocator())\n",
    "    ax.yaxis.set_minor_locator(AutoMinorLocator())\n",
    "    ax.set_xlabel(None)\n",
    "    ax.set_ylabel(None)\n",
    "    subplotspec = ax.get_subplotspec()\n",
    "    ax.xaxis.set_visible(True)\n",
    "    if subplotspec.is_first_col():\n",
    "        ax.yaxis.set_visible(True)\n",
    "\n",
    "cm.fig.set_constrained_layout(True)\n",
    "cm.fig.get_layout_engine().set(h_pad=0.2)\n",
    "cm.add_colorbar(pad=0.02)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f91d305",
   "metadata": {},
   "source": [
    "Plot map of first independent lead time (first lead time where the correlation coefficient is within the null correlation bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b819d698",
   "metadata": {
    "papermill": {
     "duration": 0.058675,
     "end_time": "2023-03-06T02:50:12.448582",
     "exception": false,
     "start_time": "2023-03-06T02:50:12.389907",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "independence.spatial_plot(\n",
    "    ds_independence,\n",
    "    dataset_name=f\"{model_name} {plot_dict['metric']}\",\n",
    "    outfile=independence_plot,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8fa63e3",
   "metadata": {},
   "source": [
    "### Drop dependent lead times based on the minimum independent lead time (may be grid point specific or a spatially averaged value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "536925a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_lead_ds = fileio.open_dataset(\n",
    "    independence_file,\n",
    "    variables=\"min_lead\",\n",
    "    shapefile=shapefile,\n",
    "    shape_overlap=shape_overlap,\n",
    "    spatial_agg=min_lead_spatial_agg,\n",
    ")\n",
    "min_lead = min_lead_ds[\"min_lead\"].load()\n",
    "\n",
    "print(min_lead)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b642bce8",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_ds = model_ds.groupby(f\"{init_dim}.month\").where(model_ds[lead_dim] >= min_lead)\n",
    "model_ds = model_ds.dropna(lead_dim, how=\"all\")\n",
    "model_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6b9abd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_additive_bc:\n",
    "    model_add_bc_ds = model_add_bc_ds.groupby(f\"{init_dim}.month\").where(model_add_bc_ds[lead_dim] >= min_lead)\n",
    "    model_add_bc_ds = model_add_bc_ds.dropna(lead_dim, how=\"all\")\n",
    "    model_add_bc_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1102b779",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_multiplicative_bc:\n",
    "    model_mulc_bc_ds = model_mulc_bc_ds.groupby(f\"{init_dim}.month\").where(model_mulc_bc_ds[lead_dim] >= min_lead)\n",
    "    model_mulc_bc_ds = model_mulc_bc_ds.dropna(lead_dim, how=\"all\")\n",
    "    model_mulc_bc_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f669987",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the model start year after dropping any dependent lead times\n",
    "start_year = model_ds.time.dt.year.min().load().item()\n",
    "start_year"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9b89883",
   "metadata": {},
   "source": [
    "## Similarity tests\n",
    "We can look at p-values for the KS-test and Anderson-Darling test for each lead time.\n",
    "\n",
    "A p value greater than alpha indicates the null hypothesis (that the two samples are from the same population) can't be rejected (i.e., blue indicates where the model and obs are significantly different).\n",
    "\n",
    "We add version of each plot with stippling over regions where the model data distribution is significantly different from the observations (p<=alpha; specific test based on user input)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6386cb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\n",
    "    f\"The map stippling is based on the {similarity_test} test significance (i.e., stippling where p<={alpha} - indicating where the model is significantly different from the observations)\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a80dcfd",
   "metadata": {},
   "outputs": [],
   "source": [
    "similarity_ds = fileio.open_dataset(str(similarity_raw_file))\n",
    "similarity_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a8277a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the similarity test results to the dataset for convenience\n",
    "model_ds[\"pval_mask\"] = similarity_ds[f\"{similarity_test}_pval\"] <= alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d3dc44e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that the colour bars for the KS and AD test p-values are the same, but the AD test p-value maximum is 0.25\n",
    "similarity.similarity_spatial_plot(\n",
    "    similarity_ds,\n",
    "    dataset_name=f\"{model_name} {plot_dict['metric']}\",\n",
    "    outfile=similarity_raw_plot,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "637f17af",
   "metadata": {},
   "source": [
    "### Additive bias corrected similarity tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84a9b1c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "if Path(similarity_add_bc_file).exists():\n",
    "    similarity_add_bc_ds = fileio.open_dataset(str(similarity_add_bc_file))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "641d5dcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_additive_bc:\n",
    "    # Add the similarity test results to the dataset\n",
    "    model_add_bc_ds[\"pval_mask\"] = similarity_add_bc_ds[f\"{similarity_test}_pval\"] <= alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "177c0a4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "if Path(similarity_add_bc_file).exists():\n",
    "    # Plot the similarity test results\n",
    "    similarity.similarity_spatial_plot(\n",
    "        similarity_add_bc_ds,\n",
    "        dataset_name=f\"{model_name} {plot_dict['metric']} (additive bias corrected)\",\n",
    "        outfile=similarity_add_bc_plot,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97e68603",
   "metadata": {},
   "source": [
    "### Multiplicative bias corrected similarity test results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "163a87e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "if Path(similarity_mulc_bc_file).exists():\n",
    "    similarity_mulc_bc_ds = fileio.open_dataset(str(similarity_mulc_bc_file))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3878738",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_multiplicative_bc:\n",
    "    model_mulc_bc_ds[\"pval_mask\"] = similarity_mulc_bc_ds[f\"{similarity_test}_pval\"] <= alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7e5235c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the similarity test results\n",
    "if Path(similarity_mulc_bc_file).exists():\n",
    "    similarity.similarity_spatial_plot(\n",
    "        similarity_mulc_bc_ds,\n",
    "        dataset_name=f\"{model_name} {plot_dict['metric']} (multiplicative bias corrected)\",\n",
    "        outfile=similarity_mulc_bc_plot,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e039e61",
   "metadata": {
    "papermill": {
     "duration": 0.013056,
     "end_time": "2023-03-06T02:50:18.132496",
     "exception": false,
     "start_time": "2023-03-06T02:50:18.119440",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Observational data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09d1baf5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open the observational dataset\n",
    "obs_ds = fileio.open_dataset(obs_file)\n",
    "obs_ds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03ecf4c5",
   "metadata": {},
   "source": [
    "Drop observations before model data is available for plots of model-observational anomalies.\n",
    "\n",
    "A note on model and observational data time periods:\n",
    "\n",
    "For most model-obs comparisions, we use the same start time for the model and observational data. Therefore, we use the start year that is updated after dropping dependent lead times, which may vary for different models or grid points (unless the minimum lead is based on a spatial aggregate) even if they were initialised in the same year. To maximise the number of available samples, we do not constrain the end of the time period. This means there may be a mismatch in the the last years of model and obs samples. \n",
    "\n",
    "\n",
    "Although, we do use the same time periods for models and obs when bias correcting the data: the climatological base period is based on the model initialisation start and end date (defined in the model dataset config file).\n",
    "\n",
    "\n",
    "There are some exceptions to this. To aid multi-model comparison of the GEV-based analysis, we use the same model start year (i.e., before dropping dependent lead times) and observational end year (for calculating trends) and the same covariate_base base year to be substituted in the non-stationary GEV model (which doesn't necessarily have to match the input data time period)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58d908d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop observations before model initialisation start year\n",
    "if reference_time_period is not None:\n",
    "    obs_ds = time_utils.select_time_period(obs_ds, reference_time_period)\n",
    "obs_ds = obs_ds.where(obs_ds.time.dt.year >= start_year, drop=True)\n",
    "obs_ds = obs_ds.dropna(\"time\", how=\"all\")\n",
    "obs_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "366a1173",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_max_event = obs_ds[var].max().load().item()\n",
    "obs_max_event_loc = obs_ds[var].where(obs_ds[var].load() == obs_max_event, drop=True).squeeze()\n",
    "obs_max_event_loc.load()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa3d6416",
   "metadata": {},
   "source": [
    "Combine plot kwargs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf650636",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Store plot related variables using the InfoSet class\n",
    "info = spatial_plots.InfoSet(\n",
    "    name=model_name,\n",
    "    fig_dir=fig_dir,\n",
    "    file=model_file,\n",
    "    obs_name=obs_name,\n",
    "    ds=model_ds,\n",
    "    ds_obs=obs_ds,\n",
    "    bias_correction=None,\n",
    "    **plot_dict,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "597577f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_additive_bc:\n",
    "    info_add_bc = spatial_plots.InfoSet(\n",
    "        name=model_name,\n",
    "        fig_dir=fig_dir,\n",
    "        file=model_add_bc_file,\n",
    "        ds=model_add_bc_ds,\n",
    "        obs_name=obs_name,\n",
    "        ds_obs=obs_ds,\n",
    "        bias_correction=\"additive\",\n",
    "        **plot_dict,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3cd891d",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_multiplicative_bc:\n",
    "    info_mulc_bc = spatial_plots.InfoSet(\n",
    "        name=model_name,\n",
    "        fig_dir=fig_dir,\n",
    "        file=model_mulc_bc_file,\n",
    "        obs_name=obs_name,\n",
    "        ds=model_mulc_bc_ds,\n",
    "        ds_obs=obs_ds,\n",
    "        bias_correction=\"multiplicative\",\n",
    "        **plot_dict,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d9beb65",
   "metadata": {},
   "source": [
    "### Stack dataset dimensions\n",
    "Combine the ensemble, lead and initialisation time dimensions into a single \"sample\" dimension. The stacked data will have shape (n_samples, n_lat, n_lon) and will be used for the remaining analysis. \n",
    "\n",
    "The dataset is also rechunked into a single array chunk along the sample dimension (avoids multiple chunks along core dimension error when using apply_ufunc with dask='parallelized' with CAFE output)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73a2cfd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_ds_stacked = model_ds.stack({\"sample\": [ensemble_dim, init_dim, lead_dim]}, create_index=False)\n",
    "model_ds_stacked = model_ds_stacked.transpose(\"sample\", ...)\n",
    "model_ds_stacked = model_ds_stacked.dropna(\"sample\", how=\"all\")\n",
    "model_ds_stacked = model_ds_stacked.chunk(dict(sample=-1))\n",
    "\n",
    "model_ds_stacked"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5adb4dce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Stack additive bias corrected data\n",
    "if plot_additive_bc:\n",
    "    model_add_bc_ds_stacked = model_add_bc_ds.stack({\"sample\": [ensemble_dim, init_dim, lead_dim]}, create_index=False)\n",
    "    model_add_bc_ds_stacked = model_add_bc_ds_stacked.transpose(\"sample\", ...)\n",
    "    model_add_bc_ds_stacked = model_add_bc_ds_stacked.dropna(\"sample\", how=\"all\")\n",
    "    model_add_bc_ds_stacked = model_add_bc_ds_stacked.chunk(dict(sample=-1))\n",
    "    assert model_ds_stacked[var].shape == model_add_bc_ds_stacked[var].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5341a6c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Stack multiplicative bias corrected data\n",
    "if plot_multiplicative_bc:\n",
    "    model_mulc_bc_ds_stacked = model_mulc_bc_ds.stack(\n",
    "        {\"sample\": [ensemble_dim, init_dim, lead_dim]}, create_index=False\n",
    "    )\n",
    "    model_mulc_bc_ds_stacked = model_mulc_bc_ds_stacked.transpose(\"sample\", ...)\n",
    "    model_mulc_bc_ds_stacked = model_mulc_bc_ds_stacked.dropna(\"sample\", how=\"all\")\n",
    "    model_mulc_bc_ds_stacked = model_mulc_bc_ds_stacked.chunk(dict(sample=-1))\n",
    "    assert model_ds_stacked[var].shape == model_mulc_bc_ds_stacked[var].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e6425ec",
   "metadata": {},
   "source": [
    "Number of samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c3c0d48",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = model_ds_stacked[var].dropna(\"sample\", how=\"any\").sample.size\n",
    "print(f\"Number of samples: {n_samples}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5386183",
   "metadata": {},
   "source": [
    "## Plot histograms of the model and observations at each capital city"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73c6f8d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lat and lon of australian capital cities\n",
    "cities = {\n",
    "    \"Canberra\": (-35.282, 149.128),\n",
    "    \"Sydney\": (-33.8688, 151.2093),\n",
    "    \"Melbourne\": (-37.8136, 144.9631),\n",
    "    \"Brisbane\": (-27.4698, 153.0251),\n",
    "    \"Perth\": (-31.9505, 115.8605),\n",
    "    \"Adelaide\": (-34.9285, 138.6007),\n",
    "    \"Hobart\": (-42.8821, 147.3272),\n",
    "    \"Darwin\": (-12.4634, 130.8456),\n",
    "}\n",
    "# Plot histograms of the model and observations at each capital city\n",
    "fig = plt.figure(figsize=(15, 12))\n",
    "\n",
    "for i, (city, (lat, lon)) in enumerate(cities.items()):\n",
    "    ax = fig.add_subplot(3, 3, i + 1)\n",
    "    obs_ds[var].sel(lat=lat, lon=lon, method=\"nearest\").plot.hist(\n",
    "        bins=20, alpha=0.6, color=\"k\", label=\"Observations\", ax=ax, density=True\n",
    "    )\n",
    "    model_ds_stacked[var].sel(lat=lat, lon=lon, method=\"nearest\").plot.hist(\n",
    "        bins=\"auto\", alpha=0.3, color=\"b\", label=model_name, ax=ax, density=True\n",
    "    )\n",
    "    if plot_additive_bc:\n",
    "        model_add_bc_ds_stacked[var].sel(lat=lat, lon=lon, method=\"nearest\").plot.hist(\n",
    "            bins=\"auto\", alpha=0.3, color=\"g\", label=\"Additive BC\", ax=ax, density=True\n",
    "        )\n",
    "    if plot_multiplicative_bc:\n",
    "        model_mulc_bc_ds_stacked[var].sel(lat=lat, lon=lon, method=\"nearest\").plot.hist(\n",
    "            bins=\"auto\", alpha=0.3, color=\"darkorange\", label=\"Mltiplicative BC\", ax=ax, density=True\n",
    "        )\n",
    "    ax.set_title(f\"{city} {info.metric}\")\n",
    "    # Plot y-axis label on the first column\n",
    "    if i % 3 == 0:\n",
    "        ax.set_ylabel(\"Density\")\n",
    "    ax.set_xlabel(f\"{plot_dict['var_name']} ({plot_dict['units']})\")\n",
    "    ax.legend(fontsize=\"small\")\n",
    "\n",
    "# Add model to suptitle\n",
    "fig.suptitle(model_name)\n",
    "plt.tight_layout()\n",
    "plt.savefig(f\"{fig_dir}/histograms_aus_capitals_{model_name}_{info.metric}.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1710e130",
   "metadata": {},
   "source": [
    "### Plot histograms where the observations and the bias corrected model are significantly different.\n",
    "This is a diagnostic plot to help determine why bias-corrected data may still have failed the similarity test (i.e., does the model have a much longer tail than observations & did this skew the mean bias used to correct the model).\n",
    "To select a few representative locations where the test failed, we search for the worst performing grid point in each NCRA region. We leave a blank subplot for any NCRA regions where there are no grid points that failed the test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e9841c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_region_min_coords(shapes, da):\n",
    "    \"\"\"Get the lat and lon with the smallest similarity pvalue in a region.\"\"\"\n",
    "    # Create a mask of grid points within the region\n",
    "    mask = spatial_selection.centre_mask(shapes, da.lon, da.lat, output=\"2D\")\n",
    "    mask = spatial_selection._squeeze_and_drop_region(mask)\n",
    "    masked_da = da.where(~np.isnan(mask), drop=True)\n",
    "    # Find the coords of da.min() in the region\n",
    "    da_min = masked_da.where((masked_da == masked_da.min()).load(), drop=True)\n",
    "    lat, lon = da_min.lat.item(), da_min.lon.item()\n",
    "    return lat, lon\n",
    "\n",
    "\n",
    "regions = regions_dict[\"ncra_regions\"]\n",
    "print(regions)\n",
    "\n",
    "region_coords = {}\n",
    "# Iterate through the NCRA regions in the shapefile\n",
    "for r in list(regions.label):\n",
    "    # Subset the NCR\n",
    "    shapes = regions[regions.label == r]\n",
    "\n",
    "    # Use bias-corrected or un-corrected model similarity test\n",
    "    if plot_additive_bc:\n",
    "        similarity_da = similarity_add_bc_ds[f\"{similarity_test}_pval\"]\n",
    "    elif plot_multiplicative_bc:\n",
    "        similarity_da = similarity_mulc_bc_ds[f\"{similarity_test}_pval\"]\n",
    "    else:\n",
    "        similarity_da = similarity_ds[f\"{similarity_test}_pval\"]\n",
    "\n",
    "    # Find the grid point with the lowest similarity test p-value in region\n",
    "    lat, lon = get_region_min_coords(shapes, similarity_da)\n",
    "\n",
    "    # Add to the dict to be plotted if the similarity if the null hypothesis is rejected\n",
    "    if similarity_da.sel(lat=lat, lon=lon, method=\"nearest\") <= alpha:\n",
    "        region_coords[r] = (lat, lon)\n",
    "region_coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0c9aebc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot histograms locations in each NCRA region where the similarity test null hypothesis is rejected\n",
    "if plot_additive_bc:\n",
    "    title_model_name = f\"{model_name} (additive bc)\"\n",
    "elif plot_multiplicative_bc:\n",
    "    title_model_name = f\"{model_name} (multiplicative bc)\"\n",
    "else:\n",
    "    title_model_name = model_name\n",
    "\n",
    "fig = plt.figure(figsize=(15, 12))\n",
    "i = 0\n",
    "for i, region in enumerate(list(regions.label)):\n",
    "    if region not in region_coords:\n",
    "        continue\n",
    "    (lat, lon) = region_coords[region]\n",
    "    ax = fig.add_subplot(3, 3, i + 1)\n",
    "\n",
    "    obs_ds[var].sel(lat=lat, lon=lon, method=\"nearest\").plot.hist(\n",
    "        bins=20, alpha=0.6, color=\"k\", label=\"Observations\", ax=ax, density=True\n",
    "    )\n",
    "    model_ds_stacked[var].sel(lat=lat, lon=lon, method=\"nearest\").plot.hist(\n",
    "        bins=\"auto\", alpha=0.4, color=\"b\", label=model_name, ax=ax, density=True\n",
    "    )\n",
    "    if plot_additive_bc:\n",
    "        model_add_bc_ds_stacked[var].sel(lat=lat, lon=lon, method=\"nearest\").plot.hist(\n",
    "            bins=\"auto\", alpha=0.4, color=\"g\", label=\"Additive BC\", ax=ax, density=True\n",
    "        )\n",
    "    elif plot_multiplicative_bc:\n",
    "        model_mulc_bc_ds_stacked[var].sel(lat=lat, lon=lon, method=\"nearest\").plot.hist(\n",
    "            bins=\"auto\", alpha=0.4, color=\"darkorange\", label=\"Mltiplicative BC\", ax=ax, density=True\n",
    "        )\n",
    "    ax.set_title(f\"{region} ({lat:.1f}, {lon:.1f})\")\n",
    "    # Plot y-axis label on the first column\n",
    "    if i % 3 == 0:\n",
    "        ax.set_ylabel(\"Density\")\n",
    "    ax.set_xlabel(f\"{plot_dict['var_name']} ({plot_dict['units']})\")\n",
    "    ax.legend(fontsize=\"small\")\n",
    "\n",
    "\n",
    "# Add model to suptitle\n",
    "fig.suptitle(\n",
    "    f\"{title_model_name} {info.metric}: exemplar locations with significant obs-model differences in each NCRA region\"\n",
    ")\n",
    "plt.tight_layout()\n",
    "plt.savefig(f\"{fig_dir}/histograms_bad_fit_{model_name}_{info.metric}.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e24a842e",
   "metadata": {},
   "source": [
    "\n",
    "## GEV distribution parameters\n",
    "\n",
    "The stationary GEV distribution is described by three parameters: shape ($\\xi$), location ($\\mu$) and scale ($\\sigma$) parameters. \n",
    "\n",
    "The location parameter is qualitatively similar to the median, the scale parameter to the standard deviation and the shape parameter to the skewness of the distribution.\n",
    "\n",
    "The non-stationary GEV model allows the location and scale parameters to vary linearly with time $t$ as follows:\n",
    "\n",
    "$\\mu = \\mu_0 + \\mu_1 \\times t$  \n",
    "$\\sigma = \\sigma_0 + \\sigma_1 \\times t$.\n",
    "\n",
    "The location trend term $\\mu_1$ can be useful for understanding trends in the mean, while the scale trend term $\\sigma_1$ captures the linear trend in the variability of the data (for example changes in the tail behaviours).\n",
    "\n",
    "The shape parameter is assumed to be constant over time, as it tends to be sensitive and difficult to estimate accurately (Coles, 2001).\n",
    "\n",
    "Coles, Stuart. An Introduction to Statistical Modeling of Extreme Values. Springer Series in Statistics, London, Springer London, 2001, [doi.org/10.1007/978-1-4471-3675-0](https://doi.org/10.1007/978-1-4471-3675-0).\n",
    "\n",
    "\n",
    "For any value of $t$, the model provides a specific location ($\\mu$) and scale ($\\sigma$) parameters, resulting in a unique probability density function (PDF) and return value. It is important to note that the model estimates the distribution for a specific year based on the overall trends and relationships in the data. That is, the PDF for any given year is optimised to fit the entire dataset, not just the data from a specific year (like linear regression).\n",
    "\n",
    "\n",
    "Unless otherwise specified (i.e., for spatial_plots.plot_aep and spatial_plots.new_record_probability), the non-stationary GEV is evaluated using the `covariate_base`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6b0e47e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Non-stationary covariate (i.e., a variable that co-varies with the target variable)\n",
    "covariate = model_ds_stacked[time_dim].dt.year\n",
    "dparams_ns = xr.open_dataset(gev_params_nonstationary_file)[var]\n",
    "dparams_ns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3992f1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ensure model ds and GEV parameters have the same lat/lon grid (it is hard to fit a\n",
    "# GEV at ocean points, so we usually apply a mask that may drop some lat/lons)\n",
    "model_ds_stacked = model_ds_stacked.sel(lat=dparams_ns.lat, lon=dparams_ns.lon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7f4fd45",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_additive_bc:\n",
    "    dparams_add_bc_ns = xr.open_dataset(gev_params_nonstationary_add_bc_file)[var]\n",
    "    model_add_bc_ds_stacked = model_add_bc_ds_stacked.sel(lat=dparams_add_bc_ns.lat, lon=dparams_add_bc_ns.lon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8ace33b",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_multiplicative_bc:\n",
    "    dparams_mulc_bc_ns = xr.open_dataset(gev_params_nonstationary_mulc_bc_file)[var]\n",
    "    model_mulc_bc_ds_stacked = model_mulc_bc_ds_stacked.sel(lat=dparams_mulc_bc_ns.lat, lon=dparams_mulc_bc_ns.lon)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87b82982",
   "metadata": {},
   "source": [
    "### Plots of GEV distribution parameters\n",
    "\n",
    "Diagnostic plot of all time-dependent GEV parameters: shape, location (intercept and trend) and scale (intercept and trend). Note that the stationary parameters (with no trend) will be shown if if nonstationary GEV fit is not significantly better than the stationary fit (depending on the user input)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "037c9ff2",
   "metadata": {},
   "outputs": [],
   "source": [
    "eva.spatial_plot_gev_parameters(\n",
    "    dparams_ns,\n",
    "    dataset_name=f\"{model_name} {plot_dict['metric']}\",\n",
    "    outfile=f\"{info.fig_dir}/gev_parameters_{info.filestem()}.png\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "256c71c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_additive_bc:\n",
    "    eva.spatial_plot_gev_parameters(\n",
    "        dparams_add_bc_ns,\n",
    "        dataset_name=f\"{model_name} {plot_dict['metric']} (additive bias corrected)\",\n",
    "        outfile=f\"{info_add_bc.fig_dir}/gev_parameters_{info_add_bc.filestem()}.png\",\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7267c645",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_multiplicative_bc:\n",
    "    eva.spatial_plot_gev_parameters(\n",
    "        dparams_mulc_bc_ns,\n",
    "        dataset_name=f\"{model_name} {plot_dict['metric']} (multiplicative bias corrected)\",\n",
    "        outfile=f\"{info_mulc_bc.fig_dir}/gev_parameters_{info_mulc_bc.filestem()}.png\",\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17ea70af",
   "metadata": {},
   "source": [
    "### Plot the trends in the GEV location parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "352a067b",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_gev_param_trend(info, dparams_ns, param=\"location\", mask=mask)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d7ba075",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()  # Hides the non-masked plot\n",
    "        spatial_plots.plot_gev_param_trend(\n",
    "            info_add_bc,\n",
    "            dparams_add_bc_ns,\n",
    "            param=\"location\",\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4279feda",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()  # Hides the non-masked plot\n",
    "        spatial_plots.plot_gev_param_trend(\n",
    "            info_mulc_bc,\n",
    "            dparams_mulc_bc_ns,\n",
    "            param=\"location\",\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83db614e",
   "metadata": {},
   "source": [
    "Diagnostic plot:\n",
    " the difference between the nonstationary GEV parameters and the best of stationary/nonstationary GEV parameters (stationary GEV parameters in regions where it is a better fit than the nonstationary parameters). By default, the analysis is based on the nonstationary GEV parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98e36235",
   "metadata": {},
   "outputs": [],
   "source": [
    "if Path(gev_params_best_file).exists():\n",
    "    dparams_best = xr.open_dataset(gev_params_best_file)[var]\n",
    "    dparams_diff = dparams_ns - dparams_best\n",
    "    eva.spatial_plot_gev_parameters(\n",
    "        dparams_diff,\n",
    "        dataset_name=f\"{info.name} (non-stationary - best)\",\n",
    "        outfile=f\"{info.fig_dir}/gev_parameters_best_diff_{info.filestem()}.png\",\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df6e3773",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_additive_bc and Path(gev_params_best_add_bc_file).exists():\n",
    "    dparams_add_bc_best = xr.open_dataset(gev_params_best_add_bc_file)[var]\n",
    "    dparams_add_bc_diff = dparams_add_bc_ns - dparams_add_bc_best\n",
    "    eva.spatial_plot_gev_parameters(\n",
    "        dparams_add_bc_diff,\n",
    "        dataset_name=f\"{info_add_bc.name} (non-stationary - best)\",\n",
    "        outfile=f\"{info.fig_dir}/gev_parameters_best_diff_{info_add_bc.filestem()}.png\",\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0548d4e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_multiplicative_bc and Path(gev_params_best_mulc_bc_file).exists():\n",
    "    dparams_mulc_bc_best = xr.open_dataset(gev_params_best_mulc_bc_file)[var]\n",
    "    dparams_mulc_bc_diff = dparams_mulc_bc_ns - dparams_mulc_bc_best\n",
    "    eva.spatial_plot_gev_parameters(\n",
    "        dparams_mulc_bc_diff,\n",
    "        dataset_name=f\"{info_mulc_bc.name} (non-stationary - best)\",\n",
    "        outfile=f\"{info.fig_dir}/gev_parameters_best_diff_{info_mulc_bc.filestem()}.png\",\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bc57686",
   "metadata": {},
   "source": [
    "### Plot the trends in the GEV scale parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "379c2725",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_gev_param_trend(info, dparams_ns, param=\"scale\", mask=mask)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51fe2ad2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bias corrected data\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_gev_param_trend(\n",
    "            info_add_bc,\n",
    "            dparams_add_bc_ns,\n",
    "            param=\"scale\",\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46ef2445",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_gev_param_trend(\n",
    "            info_mulc_bc,\n",
    "            dparams_mulc_bc_ns,\n",
    "            param=\"scale\",\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "165a2e53",
   "metadata": {},
   "source": [
    "## Record event year\n",
    "Plot the year when the record event occurred for each grid point. \n",
    "\n",
    "Lighter colours indicate when more recent record breaking events occurred. We generally expect the year of the record event to be random. If there is a trend in the record event year, it may indicate a trend in the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16fd082e",
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_plots.plot_event_year(info, model_ds_stacked, time_agg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "876e9a19",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_additive_bc:\n",
    "    spatial_plots.plot_event_year(info_add_bc, model_add_bc_ds_stacked, time_agg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9292c213",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_multiplicative_bc:\n",
    "    spatial_plots.plot_event_year(info_mulc_bc, model_mulc_bc_ds_stacked, time_agg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac9770e2",
   "metadata": {},
   "source": [
    "## Most common month when event occurs\n",
    "Plot of the most common month when the event occurred for each grid point. \n",
    "\n",
    "While the most common month is shown, there may be other months when events when events commonly occur."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "454e3d2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_event_month_mode(info, model_ds_stacked, mask=mask)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e0d87c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_event_month_mode(info_add_bc, model_add_bc_ds_stacked, mask=mask)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41a82283",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_event_month_mode(info_mulc_bc, model_mulc_bc_ds_stacked, mask=mask)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98d99237",
   "metadata": {},
   "source": [
    "\n",
    "## Map of metric median\n",
    "Plot the median value of the metric for each grid point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c4ec4de",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_time_agg(info, model_ds_stacked, \"median\", mask=mask)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e00a1e2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_time_agg(info_add_bc, model_add_bc_ds_stacked, \"median\", mask=mask)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2de9b089",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_time_agg(info_mulc_bc, model_mulc_bc_ds_stacked, \"median\", mask=mask)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4da16db",
   "metadata": {},
   "source": [
    "### Model minus observations (median anomaly)\n",
    "Plot the difference between the model and observational median "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95207dd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_obs_anom(info, model_ds_stacked, obs_ds, \"median\", \"anom\", mask=mask)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "349b9259",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_obs_anom(info_add_bc, model_add_bc_ds_stacked, obs_ds, \"median\", \"anom\", mask=mask)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8408d48f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_obs_anom(info_mulc_bc, model_mulc_bc_ds_stacked, obs_ds, \"median\", \"anom\", mask=mask)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17576f05",
   "metadata": {},
   "source": [
    "# Map of metric maximum or minimum\n",
    "Plot the maximum or minimum value of the metric for each grid point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c206d2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_time_agg(info, model_ds_stacked, time_agg, mask=mask)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "946f3783",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_time_agg(info_add_bc, model_add_bc_ds_stacked, time_agg, mask=mask)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98aad3bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_time_agg(info_mulc_bc, model_mulc_bc_ds_stacked, time_agg, mask=mask)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfe778eb",
   "metadata": {},
   "source": [
    "## Model and observation record event difference\n",
    "\n",
    "Plot the difference between the model and observed record event using the following metrics:\n",
    "- absolute differenc (model - obs)\n",
    "- anomaly percentage (model - obs / obs)\n",
    "- standardised anomaly (model - obs / obs standard deviation)\n",
    "- ratio of the theoretical maximum to the observed maximum (Model GEV-estimated 2000 year return level / obs)\n",
    "\n",
    "Note that the metric in observations is conservatively regridded to the coarser model grid in order to plot the difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a792a6b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "anom_metrics = [\"anom\", \"anom_pct\", \"anom_std\", \"anom_2000yr\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f44345c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "for anom_metric in anom_metrics:\n",
    "    for mask in [None, model_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_obs_anom(\n",
    "            info,\n",
    "            model_ds_stacked,\n",
    "            obs_ds,\n",
    "            time_agg,\n",
    "            anom_metric,\n",
    "            dparams_ns,\n",
    "            covariate_base,\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77137d47",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for anom_metric in anom_metrics:\n",
    "        for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "            plt.close()\n",
    "            spatial_plots.plot_obs_anom(\n",
    "                info_add_bc,\n",
    "                model_add_bc_ds_stacked,\n",
    "                obs_ds,\n",
    "                time_agg,\n",
    "                anom_metric,\n",
    "                dparams_add_bc_ns,\n",
    "                covariate_base,\n",
    "                mask=mask,\n",
    "            )\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "563c79b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for anom_metric in anom_metrics:\n",
    "        for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "            plt.close()\n",
    "            spatial_plots.plot_obs_anom(\n",
    "                info_mulc_bc,\n",
    "                model_mulc_bc_ds_stacked,\n",
    "                obs_ds,\n",
    "                time_agg,\n",
    "                anom_metric,\n",
    "                dparams_mulc_bc_ns,\n",
    "                covariate_base,\n",
    "                mask=mask,\n",
    "            )\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a2a015e",
   "metadata": {},
   "source": [
    "# Subsampling analysis\n",
    "Here we assess how similar the model data is to the observations with comparable sample sizes. We randomly subsample the model data to match the number of observations and calculate the maximum/minimum values of the metric for each subsample. We repeat this process 1000 times and plot the median of the maximum/minimum values. Note that the subsampling function will plot the masked and unmasked versions of the plot (to avoid re-calculating the subsample, which is quite slow).\n",
    "\n",
    "First we'll plot the observed maximum/minimum values for comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c46c06fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the metric based on observations\n",
    "info_obs = spatial_plots.InfoSet(\n",
    "    name=info.obs_name,\n",
    "    fig_dir=fig_dir,\n",
    "    file=obs_file,\n",
    "    ds=obs_ds,\n",
    "    bias_correction=None,\n",
    "    **plot_dict,\n",
    ")\n",
    "spatial_plots.plot_time_agg(info_obs, obs_ds, time_agg, savefig=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25c841f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_plots.plot_time_agg_subsampled(info, model_ds_stacked, obs_ds, time_agg, resamples=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b19f7962",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    spatial_plots.plot_time_agg_subsampled(\n",
    "        info_add_bc,\n",
    "        model_add_bc_ds_stacked,\n",
    "        obs_ds,\n",
    "        time_agg,\n",
    "        resamples=1000,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f04c8110",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    spatial_plots.plot_time_agg_subsampled(\n",
    "        info_mulc_bc,\n",
    "        model_mulc_bc_ds_stacked,\n",
    "        obs_ds,\n",
    "        time_agg,\n",
    "        resamples=1000,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbfb04db",
   "metadata": {},
   "source": [
    "## Annual exceedance probability \n",
    "### Empirical exceedance probability\n",
    "Plot the empirical annual exceedance probability for different thresholds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "280d952a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ARI: 10, 100, 1000 years (i.e., 10% AEP is equiv to a 1-in-10-year event)\n",
    "aep_list = [10, 1, 0.1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d0671aa",
   "metadata": {},
   "source": [
    "### Plot of the empirical 10% AEP (1-in-10-year event)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c716c71",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_aep_empirical(\n",
    "        info,\n",
    "        model_ds_stacked,\n",
    "        aep=aep_list[0],\n",
    "        mask=mask,\n",
    "    )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fa9bb21",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_aep_empirical(\n",
    "            info_add_bc,\n",
    "            model_add_bc_ds_stacked,\n",
    "            aep=aep_list[0],\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce0cf57a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_aep_empirical(\n",
    "            info_mulc_bc,\n",
    "            model_mulc_bc_ds_stacked,\n",
    "            aep=aep_list[0],\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fa3815c",
   "metadata": {},
   "source": [
    "### Plot of the empirical 1% AEP (1-in-100-year event)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4cdfda2",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_aep_empirical(\n",
    "        info,\n",
    "        model_ds_stacked,\n",
    "        aep=aep_list[1],\n",
    "        mask=mask,\n",
    "    )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da891ce2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_aep_empirical(\n",
    "            info_add_bc,\n",
    "            model_add_bc_ds_stacked,\n",
    "            aep=aep_list[1],\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94cf60bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_aep_empirical(\n",
    "            info_mulc_bc,\n",
    "            model_mulc_bc_ds_stacked,\n",
    "            aep=aep_list[1],\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4a9acb6",
   "metadata": {},
   "source": [
    "### Plot of the empirical 0.1% AEP (1-in-1000-year event)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b70c71c",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_aep_empirical(\n",
    "        info,\n",
    "        model_ds_stacked,\n",
    "        aep=aep_list[2],\n",
    "        mask=mask,\n",
    "    )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f587da4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_aep_empirical(\n",
    "            info_add_bc,\n",
    "            model_add_bc_ds_stacked,\n",
    "            aep=aep_list[2],\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "946d33c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_aep_empirical(\n",
    "            info_mulc_bc,\n",
    "            model_mulc_bc_ds_stacked,\n",
    "            aep=aep_list[2],\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d47a0f29",
   "metadata": {},
   "source": [
    "### GEV-based exceedance probability\n",
    "Plot the GEV-based annual exceedance probability for different thresholds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b284332",
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is the start and end year we use generate trends in the AEP\n",
    "times = xr.DataArray(gev_trend_period, dims=info.time_dim)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1e05d09",
   "metadata": {},
   "source": [
    "### Plot of 10% AEP (1-in-10-year event) using the non-stationary GEV (past year, current year and the change per decade)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b9cbe60",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_aep(\n",
    "        info,\n",
    "        dparams_ns,\n",
    "        times,\n",
    "        aep=aep_list[0],\n",
    "        mask=mask,\n",
    "    )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99161dbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_aep(\n",
    "            info_add_bc,\n",
    "            dparams_add_bc_ns,\n",
    "            times,\n",
    "            aep=aep_list[0],\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d899f414",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_aep(\n",
    "            info_mulc_bc,\n",
    "            dparams_mulc_bc_ns,\n",
    "            times,\n",
    "            aep=aep_list[0],\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ce911af",
   "metadata": {},
   "source": [
    "### Plot of 1% AEP (1-in-100-year event) using the non-stationary GEV (past year, current year and the change per decade)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0099e09b",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_aep(\n",
    "        info,\n",
    "        dparams_ns,\n",
    "        times,\n",
    "        aep=aep_list[1],\n",
    "        mask=mask,\n",
    "    )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88d2aaaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_aep(\n",
    "            info_add_bc,\n",
    "            dparams_add_bc_ns,\n",
    "            times,\n",
    "            aep=aep_list[1],\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ee98ec3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_aep(\n",
    "            info_mulc_bc,\n",
    "            dparams_mulc_bc_ns,\n",
    "            times,\n",
    "            aep=aep_list[1],\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5f3bd86",
   "metadata": {},
   "source": [
    "### Plot of 0.1% AEP (1-in-1000-year event) using the non-stationary GEV (past year, current year and the change per decade)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54a9f94a",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_aep(\n",
    "        info,\n",
    "        dparams_ns,\n",
    "        times,\n",
    "        aep=aep_list[2],\n",
    "        mask=mask,\n",
    "    )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4ac6ade",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_aep(\n",
    "            info_add_bc,\n",
    "            dparams_add_bc_ns,\n",
    "            times,\n",
    "            aep=aep_list[2],\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cef94e71",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_aep(\n",
    "            info_mulc_bc,\n",
    "            dparams_mulc_bc_ns,\n",
    "            times,\n",
    "            aep=aep_list[2],\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94be63cd",
   "metadata": {},
   "source": [
    "## Annual reccurence of observed max/min event\n",
    "\n",
    "### Empirical return period\n",
    "Plot the empirical return period for the observed maximum or minimum event"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db710915",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_obs_ari_empirical(\n",
    "        info,\n",
    "        obs_ds,\n",
    "        model_ds_stacked,\n",
    "        time_agg=time_agg,\n",
    "        mask=mask,\n",
    "    )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e9a1708",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_obs_ari_empirical(\n",
    "            info_add_bc,\n",
    "            obs_ds,\n",
    "            model_add_bc_ds_stacked,\n",
    "            time_agg=time_agg,\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "034c0763",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_obs_ari_empirical(\n",
    "            info_mulc_bc,\n",
    "            obs_ds,\n",
    "            model_mulc_bc_ds_stacked,\n",
    "            time_agg=time_agg,\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66277f19",
   "metadata": {},
   "source": [
    "### GEV-based return period\n",
    "Plot the GEV-based return period for the observed maximum or minimum event"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8766d372",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_obs_ari(\n",
    "        info,\n",
    "        obs_ds,\n",
    "        model_ds_stacked,\n",
    "        dparams_ns,\n",
    "        covariate_base,\n",
    "        time_agg=time_agg,\n",
    "        mask=mask,\n",
    "    )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eea76556",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_obs_ari(\n",
    "            info_add_bc,\n",
    "            obs_ds,\n",
    "            model_add_bc_ds_stacked,\n",
    "            dparams_add_bc_ns,\n",
    "            covariate_base,\n",
    "            time_agg=time_agg,\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63c3e983",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_obs_ari(\n",
    "            info_mulc_bc,\n",
    "            obs_ds,\n",
    "            model_mulc_bc_ds_stacked,\n",
    "            dparams_mulc_bc_ns,\n",
    "            covariate_base,\n",
    "            time_agg=time_agg,\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7efefd9b",
   "metadata": {},
   "source": [
    "# Probability of breaking the observed record\n",
    "Plot the probability of breaking the observed record in the next 10 years.\n",
    "\n",
    "### Empirical-based probability of breaking the record"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cc5d46d",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_new_record_probability_empirical(\n",
    "        info,\n",
    "        obs_ds,\n",
    "        model_ds_stacked,\n",
    "        time_agg,\n",
    "        ari=10,\n",
    "        mask=mask,\n",
    "    )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "359f9fed",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_new_record_probability_empirical(\n",
    "            info_add_bc,\n",
    "            obs_ds,\n",
    "            model_add_bc_ds_stacked,\n",
    "            time_agg,\n",
    "            ari=10,\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10604e3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_new_record_probability_empirical(\n",
    "            info_mulc_bc,\n",
    "            obs_ds,\n",
    "            model_mulc_bc_ds_stacked,\n",
    "            time_agg,\n",
    "            ari=10,\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37367750",
   "metadata": {},
   "source": [
    "### GEV-based probability of breaking the observed record"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5802dbf5",
   "metadata": {},
   "outputs": [],
   "source": [
    "for mask in [None, model_ds_stacked.pval_mask]:\n",
    "    plt.close()\n",
    "    spatial_plots.plot_new_record_probability(\n",
    "        info,\n",
    "        obs_ds,\n",
    "        model_ds_stacked,\n",
    "        dparams_ns,\n",
    "        covariate_base,\n",
    "        time_agg,\n",
    "        ari=10,\n",
    "        mask=mask,\n",
    "    )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6204de2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additive bias corrected\n",
    "if plot_additive_bc:\n",
    "    for mask in [None, model_add_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_new_record_probability(\n",
    "            info_add_bc,\n",
    "            obs_ds,\n",
    "            model_add_bc_ds_stacked,\n",
    "            dparams_add_bc_ns,\n",
    "            covariate_base,\n",
    "            time_agg,\n",
    "            ari=10,\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c34de27c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiplicative bias corrected\n",
    "if plot_multiplicative_bc:\n",
    "    for mask in [None, model_mulc_bc_ds_stacked.pval_mask]:\n",
    "        plt.close()\n",
    "        spatial_plots.plot_new_record_probability(\n",
    "            info_mulc_bc,\n",
    "            obs_ds,\n",
    "            model_mulc_bc_ds_stacked,\n",
    "            dparams_mulc_bc_ns,\n",
    "            covariate_base,\n",
    "            time_agg,\n",
    "            ari=10,\n",
    "            mask=mask,\n",
    "        )\n",
    "    plt.show()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python (unseen)",
   "language": "python",
   "name": "unseen"
  },
  "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.11.9"
  },
  "papermill": {
   "default_parameters": {},
   "duration": 208.255213,
   "end_time": "2023-03-06T02:53:13.009470",
   "environment_variables": {},
   "exception": true,
   "input_path": "analysis.ipynb",
   "output_path": "analysis_CanESM5.ipynb",
   "parameters": {
    "agcd_file": "/g/data/xv83/dbi599/east-coast-rain/data/Rx15day_AGCD-CSIRO_r005_1900-2022_annual-aug-to-sep_east-coast-flood-region.zarr.zip",
    "file_list": "/home/599/dbi599/east-coast-rain/file_lists/CanESM5_dcppA-hindcast_pr_files.txt",
    "independence_plot": "/g/data/xv83/dbi599/east-coast-rain/figures/independence-test_Rx15day_CanESM5-dcppA-hindcast_196101-201701_annual-aug-to-sep_east-coast-flood-region.png",
    "min_lead": 1,
    "model_add_bc_file": "/g/data/xv83/dbi599/east-coast-rain/data/Rx15day_CanESM5-dcppA-hindcast_196101-201701_annual-aug-to-sep_east-coast-flood-region_bias-corrected-AGCD-CSIRO-additive.zarr.zip",
    "model_file": "/g/data/xv83/dbi599/east-coast-rain/data/Rx15day_CanESM5-dcppA-hindcast_196101-201701_annual-aug-to-sep_east-coast-flood-region.zarr.zip",
    "model_mulc_bc_file": "/g/data/xv83/dbi599/east-coast-rain/data/Rx15day_CanESM5-dcppA-hindcast_196101-201701_annual-aug-to-sep_east-coast-flood-region_bias-corrected-AGCD-CSIRO-multiplicative.zarr.zip",
    "model_name": "CanESM5",
    "nino_file": "/g/data/xv83/dbi599/east-coast-rain/data/nino34-anomaly_CanESM5-dcppA-hindcast_196101-201701_base-196101-201701.nc",
    "region_name": "east-coast-flood-region",
    "shape_file": "/g/data/xv83/dbi599/east-coast-rain/shapefiles/east-coast-flood.shp",
    "similarity_add_bc_file": "/g/data/xv83/dbi599/east-coast-rain/data/similarity-test_Rx15day_CanESM5-dcppA-hindcast_196101-201701_annual-aug-to-sep_east-coast-flood-region_bias-corrected-AGCD-CSIRO-additive.zarr.zip",
    "similarity_mulc_bc_file": "/g/data/xv83/dbi599/east-coast-rain/data/similarity-test_Rx15day_CanESM5-dcppA-hindcast_196101-201701_annual-aug-to-sep_east-coast-flood-region_bias-corrected-AGCD-CSIRO-multiplicative.zarr.zip",
    "similarity_raw_file": "/g/data/xv83/dbi599/east-coast-rain/data/similarity-test_Rx15day_CanESM5-dcppA-hindcast_196101-201701_annual-aug-to-sep_east-coast-flood-region_AGCD-CSIRO.zarr.zip"
   },
   "start_time": "2023-03-06T02:49:44.754257",
   "version": "2.3.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}