Fig7-8_seasonality_amplitude_phase_panel.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "75484b8f-0f71-4c9a-b8c0-40a771b50438",
   "metadata": {},
   "source": [
    "# Seasonal Amplitude and Phase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9eb32fb-da42-41d9-9450-67cb353a73fb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Import modules\n",
    "%matplotlib inline\n",
    "%run /g/data/w40/ri9247/code/aus_precip_benchmarking/master_functions_bmf.ipynb\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import fnmatch\n",
    "import pandas as pd\n",
    "from astropy.stats import circcorrcoef\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import matplotlib.ticker as mticker\n",
    "import matplotlib.gridspec as gridspec\n",
    "from matplotlib.patches import Rectangle\n",
    "import cartopy as cartopy\n",
    "import cartopy.crs as ccrs\n",
    "from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e18ee31-be75-4a05-98fc-b43f257fb72e",
   "metadata": {},
   "source": [
    "## Define Spatiotemporal Boundaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3118deb-2000-40ec-9089-d38c178ee275",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Define region of interest (lat and lon boundaries for all of Australia)\n",
    "lat_slice = slice(-44.5,-10)\n",
    "lon_slice = slice(112,156.25)\n",
    "time_slice = slice(\"1976-01-01\", \"2005-12-30\")\n",
    "\n",
    "# Define temporal scale for SPI variable (1 for 3-months, 2 for 6-months, 3 for 12-months); use \"None\" if not using the SPI variable\n",
    "iscale = None\n",
    "\n",
    "# Define season as a list of month numbers. If not wanting a seasonal breakdown, use None.\n",
    "season = None \n",
    "season_name = 'Annual' \n",
    "\n",
    "# Paths to mask datasets\n",
    "qc_mask_ds = xr.open_dataset('/g/data/w40/ri9247/CORDEX-Australasia/data/obs/AUS-44i_grid/no_indices/mask/AUS-44i_combined_quality_mask_no_oceans.nc')\n",
    "\n",
    "# Extract masks over Australia\n",
    "qc_mask = qc_mask_ds.station_mask.sel(lat=lat_slice, lon=lon_slice)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2243121a-4e16-4fa4-b906-37a5e87c7c46",
   "metadata": {},
   "source": [
    "## Define Benchmarking Thresholds for Metrics for Amplitude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "879cf33a-0a03-4a15-b45b-741a84d11367",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "amplitude_threshold = 0.6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "101d03af-b8a0-423d-beb3-8b981fdc42ec",
   "metadata": {},
   "source": [
    "## Define Keywords for Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d1ff48d-e4fe-43f4-9739-1aa5f3bf3203",
   "metadata": {
    "tags": []
   },
   "source": [
    "Keywords needed are bracketed {}. <br><br>\n",
    "**RCM** database is organized following: <br>\n",
    "**Climpact Indices:** <br>\n",
    "parent_directory/{time_period}/{index-keyword}/{variable}_{time_average} _dataset_file.nc \n",
    "<br><br>\n",
    "**Observations** database is organized as follows: <br>\n",
    "obs_parent_directory/grid_type}/climpact/{variable}/{variable} _{time_average}_agcd_historical_v1_1950-2020.nc <br><br>\n",
    "**Keyword Options:** <br>\n",
    "grid_type: 'AUS-44i_grid', 'native_grid', 'one_degree_grid' <br>\n",
    "time_period: 'historical', 'rcp85' <br>\n",
    "variable: See list of Climpact Indices at: https://climpact-sci.org/indices; must be all-lowercase <br>\n",
    "time_average: 'ANN', 'MON'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d040045-cf0a-4d50-b117-1320995aecb0",
   "metadata": {},
   "source": [
    "## Define list for subset of models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13333b00-3209-44b4-b194-e7bfe3f4e0d6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# List names of model subsets\n",
    "subset_names = [\n",
    "    \"ACCESS1-0   CCAM-1704\"\n",
    "    , \"ACCESS1-0   CCAM-2008\"\n",
    "    , \"ACCESS1-0   WRF360J\"\n",
    "    , \"ACCESS1-0   WRF360K\"\n",
    "    , \"CanESM2   CCAM-2008\"\n",
    "    , \"CanESM2   WRF360J\"\n",
    "    , \"CNRM-CM5   CCAM-1704\"\n",
    "    , \"GFDL-ESM2M   CCAM-1704\"\n",
    "    , \"GFDL-ESM2M   CCAM-2008\"\n",
    "    , \"HadGEM2-CC   CCAM-1704\"\n",
    "    , \"HadGEM2-ES   CCLM5-0-15\"\n",
    "    , \"HadGEM2-ES   REMO2015\"\n",
    "    , \"MIROC5   CCAM-1704\"\n",
    "    , \"MIROC5   CCAM-2008\"\n",
    "    , \"MPI-ESM-LR   CCLM5-0-15\"\n",
    "    , \"MPI-ESM-LR   REMO2015\"\n",
    "    , \"MPI-ESM-MR   RegCM4-7\"\n",
    "    , \"NorESM1-M   CCAM-1704\"\n",
    "    , \"NorESM1-M   RegCM4-7\"\n",
    "    , \"NorESM1-M   REMO2015\"\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2257425-3fad-4704-be3e-bf6a98c0340a",
   "metadata": {},
   "source": [
    "### Using Climpact Indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be34e777-1bcb-49c7-8a86-cb169eda9756",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Define keywords to be used with f-strings to build the path to datasets\n",
    "time_period = 'historical'\n",
    "variable = 'prcptot'\n",
    "time_average = 'MON'\n",
    "grid_type = 'AUS-44i_grid'\n",
    "\n",
    "# Define paths to data\n",
    "model_master_path = '/g/data/ks32/CLEX_Data/CORDEX_Australasia_Indices/v1-0/'\n",
    "obs_master_path = '/g/data/w40/ri9247/CORDEX-Australasia/data/obs/'\n",
    "\n",
    "model_data_path = model_master_path + f'{time_period}/{variable}/'\n",
    "agcd_data_path = obs_master_path + f'{grid_type}/climpact/{variable}/{variable}_{time_average}_agcd_historical_v1_1950-2020.nc'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e44e3cf-a9ff-46a0-b20e-a1029531235c",
   "metadata": {},
   "source": [
    "## Get Pandas Dataframe with list of dataset file paths for subset of model simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c51e456f-0c78-4824-9ef0-d4cef286c2d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get all model paths in full ensemble (this will be stored in a Pandas DataFrame)\n",
    "model_paths = get_model_files(model_data_path, time_average)\n",
    "\n",
    "# Get file paths for subset of models\n",
    "model_paths_subset = get_model_files_subset(model_paths, subset_names)\n",
    "pd.set_option('display.max_colwidth', None)\n",
    "\n",
    "# Option to print paths to confirm we get the correct files\n",
    "model_paths_subset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "015f8542-04d0-46fa-a57b-2b12fbf752b5",
   "metadata": {},
   "source": [
    "## Calculate Amplitude and Phase of the Seasonal Cycle at each Grid Point"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21ece2e2-9c33-4df9-8cc7-02b37a9b859b",
   "metadata": {},
   "source": [
    "## Amplitude"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52c1af5f-ad62-4cb8-a603-98c616465b5a",
   "metadata": {},
   "source": [
    "### Observations (AGCD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31c0b822-ade8-4c17-a0a3-4b0bac410c73",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define data type for function\n",
    "data_type = \"obs\"\n",
    "\n",
    "# Calculate the amplitude for the Observational dataset\n",
    "obs_amplitude = get_amplitude_of_annual_cycle(agcd_data_path, variable, time_slice, lat_slice, lon_slice, data_type, iscale, qc_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62c5b5ee-a176-491e-bb25-97927fef1430",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Option to quickly plot the amplitude result\n",
    "obs_amplitude.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3aabb5a4-a62f-482d-a4c0-ef114b5f42fe",
   "metadata": {},
   "source": [
    "### Models (stored in a dictionary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2187fdf9-60cf-4578-9844-d74499e0717d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note data_type is model so units are properly converted if daily or monthly mean data\n",
    "data_type = 'model_climpact'\n",
    "\n",
    "# Intialize an empty dictionary to store model climatologies as DataArrays\n",
    "amplitude_model_dict = {}\n",
    "\n",
    "# Loop through RCM simulations and calculate the model bias; store in a dictionary\n",
    "for i, row in model_paths_subset.iterrows():\n",
    "    amplitude_xr = get_amplitude_of_annual_cycle(model_paths_subset.iloc[i,1], variable, time_slice, lat_slice, lon_slice, data_type, iscale, qc_mask)\n",
    "    amplitude_model_dict[model_paths_subset.iloc[i,0]] = amplitude_xr\n",
    "\n",
    "# Option to print the dictionary keys after completion\n",
    "print(amplitude_model_dict.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a723e36-4849-4a57-b1d8-fdfed6d72212",
   "metadata": {},
   "source": [
    "### Calculate NRMSE of Amplitude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a897e15-2639-490e-ad46-360abfba9f42",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Initialize Pandas DataFrame to store circular correlation coefficient values\n",
    "model_master_nrmse_df = pd.DataFrame(columns=['dataset_name', 'amp_nrmse'])\n",
    "\n",
    "# Loop through the dictionaries to calculate the NRMSE for amplitude and phase for each simulation against obs\n",
    "for model_name in amplitude_model_dict.keys():\n",
    "    \n",
    "    # Calculate circular correlation coefficient\n",
    "    amp_nrmse = get_nrmse_mean_ppdata(obs_amplitude, amplitude_model_dict[model_name])\n",
    "    \n",
    "    # Create Pandas DataFrame for each model\n",
    "    model_amp_nrmse_df = pd.DataFrame({'dataset_name': f'{model_name}', 'amp_nrmse': amp_nrmse.item(0)}, index=[0])\n",
    "    \n",
    "    # Add model dataframe to master DataFrame\n",
    "    model_master_nrmse_df = pd.concat([model_master_nrmse_df, model_amp_nrmse_df], ignore_index=True)\n",
    "\n",
    "# Sort values for plotting\n",
    "model_master_nrmse_df = model_master_nrmse_df.sort_values('amp_nrmse', ascending=True)\n",
    "\n",
    "# Option to print the DataFrame for a sanity check\n",
    "model_master_nrmse_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b454cae9-1757-4797-a8d1-279e737517b1",
   "metadata": {},
   "source": [
    "## Phase"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a47fdd01-8a89-4e7d-af3e-bfe491124268",
   "metadata": {},
   "source": [
    "### Observations (AGCD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c14dc33a-b9f5-4636-a641-3c0f4f7c1687",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Calculate the phase of the observational product\n",
    "obs_phase = get_phase_of_annual_cycle(agcd_data_path, variable, time_slice, lat_slice, lon_slice, iscale, qc_mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "838875e3-9e9d-4f3b-a31c-cf53554efa13",
   "metadata": {},
   "source": [
    "### Models (stored in dictionary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "331642c8-fd3d-4aad-8a65-4b554634d6df",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Intialize an empty dictionary to store model climatologies as DataArrays\n",
    "phase_model_dict = {}\n",
    "\n",
    "# Loop through RCM simulations and calculate the model bias; store in a dictionary\n",
    "for i, row in model_paths_subset.iterrows():\n",
    "    phase_xr = get_phase_of_annual_cycle(model_paths_subset.iloc[i,1], variable, time_slice, lat_slice, lon_slice, iscale, qc_mask)\n",
    "    phase_model_dict[model_paths_subset.iloc[i,0]] = phase_xr\n",
    "\n",
    "# Option to print the dictionary keys after completion\n",
    "print(phase_model_dict.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6def94e-aec7-439c-a1af-92797680addb",
   "metadata": {},
   "source": [
    "### Convert Phase to Radians"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fc82e68-0f71-4e28-9b71-76a707b71989",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Convert input data to angles around a circle (divide by 12 to get values 0-1, then multiply by 2pi)\n",
    "obs_radians = ((obs_phase + 1) / 12) * 2 * np.pi\n",
    "\n",
    "# Create new dictionary for model phase as radians\n",
    "model_radians_dict = {}\n",
    "\n",
    "# Loop through phase of models dictionary and convert phase to radians; store in new dictionary\n",
    "for model_name in phase_model_dict.keys():\n",
    "    \n",
    "    # Convert phase to radians\n",
    "    model_radians = ((phase_model_dict[model_name] + 1) / 12) * 2 * np.pi\n",
    "    \n",
    "    # Store in Dictionary\n",
    "    model_radians_dict[model_name] = model_radians\n",
    "\n",
    "# Option to print the dictionary keys after completion\n",
    "print(model_radians_dict.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de61e880-1c65-4bcf-9765-1c3798653dec",
   "metadata": {},
   "source": [
    "### Calculate Circular Correlation Coefficient and Bootstrapped Confidence Interval (95th Percentile) Against AGCD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b4ae8f2-199e-4659-af8b-87fa1aa40200",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Define resampling percentage\n",
    "resample_percentage = 0.6\n",
    "\n",
    "# Convert Obs to Numpy Array and remove NaN values\n",
    "obs_alpha_i = np.asarray(obs_radians)\n",
    "obs_alpha = obs_alpha_i[~np.isnan(obs_alpha_i)]\n",
    "\n",
    "# Initialize Pandas DataFrame to store circular correlation coefficient values\n",
    "model_master_circular_cor_df = pd.DataFrame(columns=['dataset_name', 'phase_circ_corr', 'lower_ci', 'upper_ci'])\n",
    "\n",
    "# Loop through the dictionaries to calculate the spatial correlation for amplitude and phase for each simulation against obs\n",
    "for model_name in model_radians_dict.keys():\n",
    "    \n",
    "    # Convert dataset to numpy array and remove NaN values\n",
    "    model_alpha_i = np.asarray(model_radians_dict[model_name])\n",
    "    model_alpha = model_alpha_i[~np.isnan(model_alpha_i)]\n",
    "    \n",
    "    # Calculate circular correlation coefficient\n",
    "    phase_circ_corr = circcorrcoef(obs_alpha, model_alpha)\n",
    "    \n",
    "    # Calculate 95th Percentile Confidence Interval using Bootstrapping\n",
    "    bootstrap_results = []\n",
    "    \n",
    "    for i in range(5000):\n",
    "        bootstrap_results.append(paired_bootstrap(obs_alpha, model_alpha, resample_percentage))\n",
    "\n",
    "    results = np.array(bootstrap_results)\n",
    "    lower_ci,upper_ci = np.percentile(results,[2.5, 97.5])\n",
    "    \n",
    "    # Create Pandas DataFrame for each model\n",
    "    model_phase_circ_corr_df = pd.DataFrame({'dataset_name': f'{model_name}', 'phase_circ_corr': phase_circ_corr.item(0), 'lower_ci': lower_ci, 'upper_ci': upper_ci}, index=[0])\n",
    "    \n",
    "    # Add model dataframe to master DataFrame\n",
    "    model_master_circular_cor_df = pd.concat([model_master_circular_cor_df, model_phase_circ_corr_df], ignore_index=True)\n",
    "\n",
    "# Sort values for plotting\n",
    "model_master_circular_cor_df = model_master_circular_cor_df.sort_values('phase_circ_corr', ascending=False)\n",
    "\n",
    "# Option to print the DataFrame showing the confidence intervals\n",
    "model_master_circular_cor_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8106df5e-35d3-4f4a-8fb6-4ab279cf4091",
   "metadata": {},
   "source": [
    "## Set Up Panel Plot Sorted by Spatial Area Average"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cb429d1-c8f0-4f9b-8883-f8d25d672856",
   "metadata": {},
   "source": [
    "### Set up mapping constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb64750b-1634-42b0-ac70-a2a1735344e8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Define map constants and specifications\n",
    "extent = [112, 156.25, -44.5, -10]\n",
    "proj = ccrs.PlateCarree()\n",
    "\n",
    "# Axes indices where I want lat and/or lon tickmarks on the map\n",
    "# 5X4 Plot\n",
    "lat_label = [0,1,5,9,13]\n",
    "lon_label = [18,19,20]\n",
    "both_label = [17]\n",
    "\n",
    "# Lat/Lon values to show for axis labels\n",
    "lat_ticks = [-10,-20,-30,-40]\n",
    "lon_ticks = [115,125,135,145,155]\n",
    "\n",
    "my_colors_four_seasons =['#4b92c3', '#a5c8e1', '#ce660a', '#ff8b26', '#ffbf86', '#c02324', '#e26768', '#eea8a9', '#238023', '#56b356', '#bfe2bf', '#185f90']\n",
    "\n",
    "cmap_phase = mpl.colors.ListedColormap(my_colors_four_seasons, name='colorblind_friendly_4seasons')\n",
    "\n",
    "# Set Min/Max values for phase (12 months)\n",
    "mmin = 0\n",
    "mmax = 11\n",
    "\n",
    "# Print the colormap\n",
    "cmap_phase\n",
    "\n",
    "# colorblind-friendly colors available at: https://gist.github.com/thriveth/8560036"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfe26847-7d3c-49f4-977e-1dbf16c45fa1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Define the color for the masked area\n",
    "cmap_mask_i = plt.get_cmap('gray',5)\n",
    "cmap_mask_i = cmap_mask_i.reversed()\n",
    "cmap_mask = truncate_colormap(cmap_mask_i, 0.5, 1)\n",
    "\n",
    "# Option to print the colormap\n",
    "cmap_mask"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96d63804-dc12-4c63-86dc-7e1c85ff5d12",
   "metadata": {},
   "source": [
    "## Create Sorted Panel Plots - Portrait (for Paper)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d60336ae-85cc-4ec6-aad2-c0d079b1e6c2",
   "metadata": {
    "tags": []
   },
   "source": [
    "## PHASE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b6f45a4-2a54-4ecb-abfc-da84e4e16152",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Define figure size\n",
    "fig = plt.figure(figsize=(12,16))\n",
    "fig.suptitle(f\"Climatological Rainfall Phase (1976-2005)\", fontsize=16, y=0.99)\n",
    "\n",
    "# Setup axes for all subplots\n",
    "gs = gridspec.GridSpec(7,4,width_ratios=[1,1,1,1], height_ratios=[1,1,1,1,1,1,0.1])\n",
    "\n",
    "row_max = 5\n",
    "col_max = 3\n",
    "row = 1\n",
    "col = 0\n",
    "\n",
    "axs = {}\n",
    "\n",
    "# Define AGCD as the only dataset on the top row\n",
    "axs['AGCD'] = fig.add_subplot(gs[0,0], projection=proj)\n",
    "\n",
    "# Set up Axes labels (this loops through the sorted Pandas DF to assign axes positions based on wetness)\n",
    "for i, df_row_mod in model_master_circular_cor_df.iterrows():\n",
    "    axs[f'{df_row_mod[0]}'] = fig.add_subplot(gs[row,col], projection=proj)\n",
    "    \n",
    "    if col == col_max:\n",
    "        row = row + 1\n",
    "        col = 0 \n",
    "    else:\n",
    "        col = col + 1\n",
    "\n",
    "for name, ax in axs.items():\n",
    "    ax.set_title(name, fontsize=13)\n",
    "    \n",
    "# Plot Phase Maps \n",
    "# Define empty dictionary of plots\n",
    "plots = {}\n",
    "\n",
    "#Plot Phase Map for AGCD data\n",
    "plots['AGCD'] = obs_phase.plot(ax=axs['AGCD'], add_labels=False, add_colorbar=False, vmin=mmin, vmax=mmax, cmap=cmap_phase)\n",
    "(mask.where(mask==0.)).plot(ax=axs['AGCD'], add_labels=False, add_colorbar=False, cmap=cmap_mask)\n",
    "\n",
    "# Plot Phase maps and add mask   \n",
    "for model_name in phase_model_dict.keys():\n",
    "    \n",
    "    # Define Circular Correlation Coefficient Score and Confidence Interval for each model\n",
    "    phase_cir_corr_score = model_master_circular_cor_df[model_master_circular_cor_df['dataset_name'] == f'{model_name}'].iloc[0,1].round(2)\n",
    "    phase_lower_ci = model_master_circular_cor_df[model_master_circular_cor_df['dataset_name'] == f'{model_name}'].iloc[0,2].round(2)\n",
    "    phase_upper_ci = model_master_circular_cor_df[model_master_circular_cor_df['dataset_name'] == f'{model_name}'].iloc[0,3].round(2)\n",
    "    \n",
    "    # Plot Data and add grey mask \n",
    "    plots[model_name] = phase_model_dict[model_name].plot(ax=axs[model_name], add_labels=False, add_colorbar=False, vmin=mmin, vmax=mmax, cmap=cmap_phase)\n",
    "    (mask.where(mask==0.)).plot(ax=axs[model_name], add_labels=False, add_colorbar=False, cmap=cmap_mask)\n",
    "    \n",
    "    # Check if models meet Phase benchmark and add highlight passing metrics\n",
    "    if ((phase_cir_corr_score >= phase_lower_ci and phase_cir_corr_score <= phase_upper_ci) and (phase_cir_corr_score >0)):\n",
    "        axs[model_name].text(0.06, 0.06,'CCor:\\n%.2f' % (phase_cir_corr_score), fontsize=12.0, transform=axs[model_name].transAxes,\n",
    "                            bbox=dict(facecolor='lavender', edgecolor='blueviolet',boxstyle='round'))\n",
    "            \n",
    "    else:\n",
    "        axs[model_name].text(0.06, 0.06,'CCor:\\n%.2f' % (phase_cir_corr_score), fontsize=12.0, transform=axs[model_name].transAxes)\n",
    "\n",
    "# Add coastlines and set map extent\n",
    "for ax in axs.values():\n",
    "    ax.coastlines()\n",
    "    #ax.set_extent(extent)\n",
    "    \n",
    "\n",
    "# Add lat and long values on first column and bottom row; add empty tick marks to the rest of the plots    \n",
    "for i, ax in enumerate(fig.axes):\n",
    "    ax.set_xticks(lon_ticks, crs=ccrs.PlateCarree())\n",
    "    ax.set_yticks(lat_ticks, crs=ccrs.PlateCarree())\n",
    "    ax.xaxis.set_ticklabels([])\n",
    "    ax.yaxis.set_ticklabels([])\n",
    "    if i in lat_label:\n",
    "        ax.set_yticks(lat_ticks, crs=ccrs.PlateCarree())\n",
    "        lat_formatter = LatitudeFormatter()\n",
    "        ax.yaxis.set_major_formatter(lat_formatter)\n",
    "    elif i in lon_label:\n",
    "        ax.set_xticks(lon_ticks, crs=ccrs.PlateCarree())\n",
    "        lon_formatter = LongitudeFormatter(zero_direction_label=True)\n",
    "        ax.xaxis.set_major_formatter(lon_formatter)\n",
    "    elif i in both_label:\n",
    "        ax.set_xticks(lon_ticks, crs=ccrs.PlateCarree())\n",
    "        ax.set_yticks(lat_ticks, crs=ccrs.PlateCarree())\n",
    "        lon_formatter = LongitudeFormatter(zero_direction_label=True)\n",
    "        lat_formatter = LatitudeFormatter()\n",
    "        ax.xaxis.set_major_formatter(lon_formatter)\n",
    "        ax.yaxis.set_major_formatter(lat_formatter)\n",
    "    else:\n",
    "        continue\n",
    "\n",
    "\n",
    "# Plot colorbar\n",
    "cbar_gs = gridspec.GridSpecFromSubplotSpec(1,1, subplot_spec=gs[6,:], hspace=0.5)\n",
    "cbar_ax = fig.add_subplot(cbar_gs[0,0])\n",
    "tick_pos = np.linspace(0.5, 10.5, 12)\n",
    "plt.colorbar(plots['CanESM2   WRF360J'], cbar_ax, ticks=tick_pos, orientation='horizontal')\n",
    "\n",
    "# Specify colorbar labels as months and set colorbar title\n",
    "cbar_ax.set_xticklabels(['J','F','M','A','M','J','J','A','S','O','N','D'])\n",
    "cbar_ax.set_xlabel(\"Month of Maximum Rainfall\",fontsize=12.0)\n",
    "\n",
    "# Plot with tight layout\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "353b0147-c210-4172-8420-fef983801484",
   "metadata": {},
   "source": [
    "## AMPLITUDE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae92485d-7d71-419f-851a-2a2899eb4998",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Update colorbar for Amplitude\n",
    "cmap = plt.get_cmap('YlGnBu',16)\n",
    "\n",
    "# Set Min/Max values for Amplitude\n",
    "mmin = 0\n",
    "mmax = 200\n",
    "\n",
    "# Print the colormap\n",
    "cmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8525a2c4-3296-45cb-a3f7-88605a584022",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Define figure size\n",
    "fig = plt.figure(figsize=(12,16))\n",
    "fig.suptitle(f\"Climatological Rainfall Amplitude (1976-2005)\", fontsize=16, y=0.99)\n",
    "\n",
    "# Setup axes for all subplots\n",
    "gs = gridspec.GridSpec(7,4,width_ratios=[1,1,1,1], height_ratios=[1,1,1,1,1,1,0.1])\n",
    "\n",
    "row_max = 5\n",
    "col_max = 3\n",
    "row = 1\n",
    "col = 0\n",
    "\n",
    "axs = {}\n",
    "\n",
    "# Define AGCD as only dataset on the top row\n",
    "axs['AGCD'] = fig.add_subplot(gs[0,0], projection=proj)\n",
    "\n",
    "# Set up Axes labels (this loops through the sorted Pandas DF to assign axes positions based on wetness)\n",
    "for i, df_row_mod in model_master_nrmse_df.iterrows():\n",
    "    axs[f'{df_row_mod[0]}'] = fig.add_subplot(gs[row,col], projection=proj)\n",
    "    \n",
    "    if col == col_max:\n",
    "        row = row + 1\n",
    "        col = 0 \n",
    "    else:\n",
    "        col = col + 1\n",
    "\n",
    "for name, ax in axs.items():\n",
    "    ax.set_title(name, fontsize=13)\n",
    "    \n",
    "# Plot Amplitude Maps \n",
    "# Define empty dictionary of plots\n",
    "plots = {}\n",
    "\n",
    "#Plot Phase Map for AGCD data\n",
    "plots['AGCD'] = obs_amplitude.plot(ax=axs['AGCD'], add_labels=False, add_colorbar=False, vmin=mmin, vmax=mmax, cmap=cmap)\n",
    "(mask.where(mask==0.)).plot(ax=axs['AGCD'], add_labels=False, add_colorbar=False, cmap=cmap_mask)\n",
    "\n",
    "\n",
    "# Plot Phase maps and add mask   \n",
    "for model_name in amplitude_model_dict.keys():\n",
    "    \n",
    "    # Define NRMSE for each model\n",
    "    amplitude_score = model_master_nrmse_df[model_master_nrmse_df['dataset_name'] == f'{model_name}'].iloc[0,1]\n",
    "    \n",
    "    # Plot Data and add grey mask \n",
    "    plots[model_name] = amplitude_model_dict[model_name].plot(ax=axs[model_name], add_labels=False, add_colorbar=False, vmin=mmin, vmax=mmax, cmap=cmap)\n",
    "    (mask.where(mask==0.)).plot(ax=axs[model_name], add_labels=False, add_colorbar=False, cmap=cmap_mask)\n",
    "    \n",
    "    # Check if models meet Amplitude benchmark and highlight passing metrics\n",
    "    if amplitude_score <= amplitude_threshold:\n",
    "        axs[model_name].text(0.06, 0.06,'NRMSE:\\n%.2f' % (amplitude_score), fontsize=12.0, transform=axs[model_name].transAxes,\n",
    "                            bbox=dict(facecolor='lavender', edgecolor='blueviolet',boxstyle='round'))\n",
    "            \n",
    "    elif amplitude_score >= amplitude_threshold:\n",
    "        axs[model_name].text(0.06, 0.06,'NRMSE:\\n%.2f' % (amplitude_score), fontsize=12.0, transform=axs[model_name].transAxes)\n",
    "\n",
    "# Add coastlines and set map extent\n",
    "for ax in axs.values():\n",
    "    ax.coastlines()\n",
    "    #ax.set_extent(extent)\n",
    "    #ax.add_feature(cartopy.feature.OCEAN)\n",
    "    \n",
    "\n",
    "# Add lat and long values on first column and bottom row; add empty tick marks to the rest of the plots    \n",
    "for i, ax in enumerate(fig.axes):\n",
    "    ax.set_xticks(lon_ticks, crs=ccrs.PlateCarree())\n",
    "    ax.set_yticks(lat_ticks, crs=ccrs.PlateCarree())\n",
    "    ax.xaxis.set_ticklabels([])\n",
    "    ax.yaxis.set_ticklabels([])\n",
    "    if i in lat_label:\n",
    "        ax.set_yticks(lat_ticks, crs=ccrs.PlateCarree())\n",
    "        lat_formatter = LatitudeFormatter()\n",
    "        ax.yaxis.set_major_formatter(lat_formatter)\n",
    "    elif i in lon_label:\n",
    "        ax.set_xticks(lon_ticks, crs=ccrs.PlateCarree())\n",
    "        lon_formatter = LongitudeFormatter(zero_direction_label=True)\n",
    "        ax.xaxis.set_major_formatter(lon_formatter)\n",
    "    elif i in both_label:\n",
    "        ax.set_xticks(lon_ticks, crs=ccrs.PlateCarree())\n",
    "        ax.set_yticks(lat_ticks, crs=ccrs.PlateCarree())\n",
    "        lon_formatter = LongitudeFormatter(zero_direction_label=True)\n",
    "        lat_formatter = LatitudeFormatter()\n",
    "        ax.xaxis.set_major_formatter(lon_formatter)\n",
    "        ax.yaxis.set_major_formatter(lat_formatter)\n",
    "    else:\n",
    "        continue\n",
    "\n",
    "\n",
    "# Plot colorbars\n",
    "cbar_gs = gridspec.GridSpecFromSubplotSpec(1,1, subplot_spec=gs[6,:], hspace=0.5)\n",
    "cbar_ax = fig.add_subplot(cbar_gs[0,0])\n",
    "plt.colorbar(plots['CanESM2   WRF360J'], cbar_ax, orientation='horizontal', extend='max')\n",
    "cbar_ax.set_xlabel(\"Amplitude (mm)\",fontsize=12.0)\n",
    "\n",
    "# Plot with tight layout\n",
    "plt.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}