ProbabilityRatio_calcu_E3SM.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3c23fa9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import scipy.io as scio\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import cartopy.feature as cfeature\n",
    "import cartopy.crs as ccrs\n",
    "import netCDF4 as nc\n",
    "import numpy as np\n",
    "from Load_data import Data_from_nc\n",
    "import xarray as xr\n",
    "import math\n",
    "import gc\n",
    "\n",
    "data_land_mask = Data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/surfdata_0.9x1.25_hist_78pfts_CMIP6_simyr1850_c201126.nc')\n",
    "data_mask = data_land_mask.load_variable('PFTDATA_MASK')\n",
    "data_mask = data_mask[29:, :]\n",
    "\n",
    "#create func to compute PR\n",
    "def comp_PR(ctl,scen,q):\n",
    "    #first get quantile from ctl\n",
    "    qctl=ctl.quantile(q,dim='time')\n",
    "    \n",
    "    del ctl\n",
    "    gc.collect()\n",
    "    \n",
    "\n",
    "    #get subsample scenarios\n",
    "    qscen=xr.where(scen>qctl,1,0)\n",
    "    pscen=qscen.sum('time')/len(qscen.time)\n",
    "    \n",
    "    del qscen\n",
    "    gc.collect()\n",
    "\n",
    "    #compute PR\n",
    "    PR_scen=pscen/(1-q)\n",
    "    PR_scen = np.where(data_mask==0, np.nan, PR_scen)\n",
    "    PR_scen = np.where(PR_scen==0, np.nan, PR_scen)\n",
    "    \n",
    "    return PR_scen\n",
    "\n",
    "def get_data_from_mat(file, variable):\n",
    "    var_dict = scio.loadmat(file)\n",
    "    var = var_dict[variable]\n",
    "    var = var[:, 29:, :]\n",
    "    return var\n",
    "\n",
    "data_surface = Data_from_nc(\\\n",
    "                            '/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/' + \\\n",
    "                            'surfdata_irrigation_method.nc'\\\n",
    "                           )   #load the data\n",
    "\n",
    "data_irrigation_method = data_surface.load_variable('irrigation_method')\n",
    "\n",
    "data_lon = data_surface.load_variable('LONGXY')\n",
    "data_lon = data_lon[0,:]\n",
    "data_lat = data_surface.load_variable('LATIXY')\n",
    "data_lat = data_lat[29:,0]\n",
    "data_lat = data_lat[:]\n",
    "data_time = np.array(range(262800))\n",
    "\n",
    "def get_data_for_year(metric, period, exp, ens):\n",
    "    data_out = get_data_from_mat('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/' +\\\n",
    "                                             metric + '_' + period + '_3_hourly_' + exp + ens + '_1.mat', metric)\n",
    "    data_out = np.array(data_out)\n",
    "    return data_out\n",
    "\n",
    "\n",
    "def nparray_to_xarray_for_this_study(np_array):\n",
    "    data_xrray = xr.DataArray(np_array, coords={'time': data_time, 'y': data_lat, 'x': data_lon}, dims=[\"time\", \"y\", \"x\"])\n",
    "    return data_xrray\n",
    "    \n",
    "def calcu_PR_between_two_array(data_ref, data_new, perc):\n",
    "    data_xr_ref = xr.DataArray(data_ref, coords={'time': data_time, 'y': data_lat, 'x': data_lon}, dims=[\"time\", \"y\", \"x\"])\n",
    "    data_xr_new = xr.DataArray(data_new, coords={'time': data_time, 'y': data_lat, 'x': data_lon}, dims=[\"time\", \"y\", \"x\"])\n",
    "    PR_result = comp_PR(data_xr_ref,data_xr_new, perc)\n",
    "    return PR_result\n",
    "\n",
    "\n",
    "def calcu_3_PRs(metric, perc):\n",
    "    \n",
    "    print('test')\n",
    "    data_IRR01_1901_1930 = get_data_for_year(metric, '1901', 'IRR', '01')\n",
    "    data_NOI01_1901_1930 = get_data_for_year(metric, '1901', 'NOI', '01')\n",
    "    \n",
    "    data_IRR02_1901_1930 = get_data_for_year(metric, '1901', 'IRR', '02')\n",
    "    data_NOI02_1901_1930 = get_data_for_year(metric, '1901', 'NOI', '02')\n",
    "    \n",
    "    print('test1')\n",
    "    \n",
    "    for period in range(1902, 1931):\n",
    "        data_IRR01_year = get_data_for_year(metric, str(period), 'IRR', '01')\n",
    "        data_NOI01_year = get_data_for_year(metric, str(period), 'NOI', '01')\n",
    "        \n",
    "        data_IRR02_year = get_data_for_year(metric, str(period), 'IRR', '02')\n",
    "        data_NOI02_year = get_data_for_year(metric, str(period), 'NOI', '02')\n",
    "        \n",
    "\n",
    "        data_IRR01_1901_1930 = np.concatenate((data_IRR01_1901_1930, data_IRR01_year), axis = 0)\n",
    "        data_NOI01_1901_1930 = np.concatenate((data_NOI01_1901_1930, data_NOI01_year), axis = 0)\n",
    "        \n",
    "        data_IRR02_1901_1930 = np.concatenate((data_IRR02_1901_1930, data_IRR02_year), axis = 0)\n",
    "        data_NOI02_1901_1930 = np.concatenate((data_NOI02_1901_1930, data_NOI02_year), axis = 0)\n",
    "        \n",
    "    print('test2')\n",
    "    \n",
    "    data_IRR01_1985_2014 = get_data_for_year(metric, '1985', 'IRR', '01')\n",
    "    data_NOI01_1985_2014 = get_data_for_year(metric, '1985', 'NOI', '01')\n",
    "    \n",
    "    data_IRR02_1985_2014 = get_data_for_year(metric, '1985', 'IRR', '02')\n",
    "    data_NOI02_1985_2014 = get_data_for_year(metric, '1985', 'NOI', '02')\n",
    "    \n",
    "    print('test3')\n",
    "    \n",
    "    for period in range(1986, 2015):\n",
    "        data_IRR01_year = get_data_for_year(metric, str(period), 'IRR', '01')\n",
    "        data_NOI01_year = get_data_for_year(metric, str(period), 'NOI', '01')\n",
    "        \n",
    "        data_IRR02_year = get_data_for_year(metric, str(period), 'IRR', '02')\n",
    "        data_NOI02_year = get_data_for_year(metric, str(period), 'NOI', '02')\n",
    "        \n",
    "\n",
    "        \n",
    "        data_IRR01_1985_2014 = np.concatenate((data_IRR01_1985_2014, data_IRR01_year), axis = 0)\n",
    "        data_NOI01_1985_2014 = np.concatenate((data_NOI01_1985_2014, data_NOI01_year), axis = 0)\n",
    "        \n",
    "        data_IRR02_1985_2014 = np.concatenate((data_IRR02_1985_2014, data_IRR02_year), axis = 0)\n",
    "        data_NOI02_1985_2014 = np.concatenate((data_NOI02_1985_2014, data_NOI02_year), axis = 0)\n",
    "        \n",
    "    data_IRR_1901_1930 = np.concatenate((data_IRR01_1901_1930,\n",
    "                                         data_IRR02_1901_1930,\n",
    "                                         ), axis = 0)\n",
    "    \n",
    "    print('test4')\n",
    "          \n",
    "    del data_IRR01_1901_1930, data_IRR02_1901_1930\n",
    "    gc.collect()\n",
    "    data_time = np.array(range(data_IRR_1901_1930.shape[0]))\n",
    "    data_xr_ref = xr.DataArray(data_IRR_1901_1930, coords={'time': data_time, 'y': data_lat, 'x': data_lon}, dims=[\"time\", \"y\", \"x\"])\n",
    "    \n",
    "    del data_IRR_1901_1930\n",
    "    gc.collect()\n",
    "    \n",
    "    qctl=data_xr_ref.quantile(perc,dim='time')\n",
    "    print('test5')\n",
    "    del data_xr_ref\n",
    "    gc.collect()\n",
    "    \n",
    "\n",
    "    \n",
    "    data_IRR_1985_2014 = np.concatenate((data_IRR01_1985_2014,\n",
    "                                         data_IRR02_1985_2014,\n",
    "                                         ), axis = 0)\n",
    "    \n",
    "    del data_IRR01_1985_2014, data_IRR02_1985_2014\n",
    "    gc.collect()\n",
    "    \n",
    "    data_time = np.array(range(data_IRR_1985_2014.shape[0]))\n",
    "    data_xr_new = xr.DataArray(data_IRR_1985_2014, coords={'time': data_time, 'y': data_lat, 'x': data_lon}, dims=[\"time\", \"y\", \"x\"])\n",
    "    \n",
    "    \n",
    "    #get subsample scenarios\n",
    "    qscen=xr.where(data_xr_new>qctl,1,0)\n",
    "    \n",
    "    del data_xr_new\n",
    "    gc.collect()\n",
    "    \n",
    "    pscen=qscen.sum('time')/len(qscen.time)\n",
    "    \n",
    "    del qscen\n",
    "    gc.collect()\n",
    "\n",
    "    #compute PR\n",
    "    PR_scen=pscen/(1-perc)\n",
    "    #PR_scen = np.where(data_mask==0, np.nan, PR_scen)\n",
    "    PR_scen = np.where(PR_scen==0, np.nan, PR_scen)\n",
    "    \n",
    "    PR_1 = PR_scen\n",
    "    \n",
    "    del PR_scen\n",
    "    gc.collect()\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    print('test6')\n",
    "    data_NOI_1901_1930 = np.concatenate((data_NOI01_1901_1930,\n",
    "                                         data_NOI02_1901_1930,\n",
    "                                         ), axis = 0)\n",
    "    \n",
    "    del data_NOI01_1901_1930, data_NOI02_1901_1930\n",
    "    gc.collect()\n",
    "    data_time = np.array(range(data_NOI_1901_1930.shape[0]))\n",
    "    data_xr_ref = xr.DataArray(data_NOI_1901_1930, coords={'time': data_time, 'y': data_lat, 'x': data_lon}, dims=[\"time\", \"y\", \"x\"])\n",
    "    \n",
    "    del data_NOI_1901_1930\n",
    "    gc.collect()\n",
    "    \n",
    "    qctl=data_xr_ref.quantile(perc,dim='time')\n",
    "    \n",
    "    del data_xr_ref\n",
    "    gc.collect()\n",
    "    \n",
    "\n",
    "    \n",
    "    \n",
    "    data_NOI_1985_2014 = np.concatenate((data_NOI01_1985_2014,\n",
    "                                         data_NOI02_1985_2014,\n",
    "                                         ), axis = 0)\n",
    "    \n",
    "    del data_NOI01_1985_2014, data_NOI02_1985_2014\n",
    "    gc.collect()\n",
    "    data_time = np.array(range(data_NOI_1985_2014.shape[0]))\n",
    "    data_xr_new = xr.DataArray(data_NOI_1985_2014, coords={'time': data_time, 'y': data_lat, 'x': data_lon}, dims=[\"time\", \"y\", \"x\"])\n",
    "    \n",
    "    \n",
    "    \n",
    "    #get subsample scenarios\n",
    "    qscen=xr.where(data_xr_new>qctl,1,0)\n",
    "    \n",
    "    del data_xr_new\n",
    "    gc.collect()\n",
    "    \n",
    "    pscen=qscen.sum('time')/len(qscen.time)\n",
    "    \n",
    "    del qscen\n",
    "    gc.collect()\n",
    "\n",
    "    #compute PR\n",
    "    PR_scen=pscen/(1-perc)\n",
    "    #PR_scen = np.where(data_mask==0, np.nan, PR_scen)\n",
    "    PR_scen = np.where(PR_scen==0, np.nan, PR_scen)\n",
    "    \n",
    "    PR_2 = PR_scen\n",
    "    \n",
    "    del PR_scen\n",
    "    gc.collect()\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    data_time = np.array(range(data_NOI_1985_2014.shape[0]))\n",
    "    data_xr_ref = xr.DataArray(data_NOI_1985_2014, coords={'time': data_time, 'y': data_lat, 'x': data_lon}, dims=[\"time\", \"y\", \"x\"])\n",
    "    \n",
    "    del data_NOI_1985_2014\n",
    "    gc.collect()\n",
    "    \n",
    "    qctl=data_xr_ref.quantile(perc,dim='time')\n",
    "    \n",
    "    del data_xr_ref\n",
    "    gc.collect()\n",
    "    \n",
    "\n",
    "    data_time = np.array(range(data_IRR_1985_2014.shape[0]))\n",
    "    data_xr_new = xr.DataArray(data_IRR_1985_2014, coords={'time': data_time, 'y': data_lat, 'x': data_lon}, dims=[\"time\", \"y\", \"x\"])\n",
    "    \n",
    "    \n",
    "    \n",
    "    #get subsample scenarios\n",
    "    qscen=xr.where(data_xr_new>qctl,1,0)\n",
    "    \n",
    "    del data_xr_new\n",
    "    gc.collect()\n",
    "    \n",
    "    pscen=qscen.sum('time')/len(qscen.time)\n",
    "    \n",
    "    del qscen\n",
    "    gc.collect()\n",
    "\n",
    "    #compute PR\n",
    "    PR_scen=pscen/(1-perc)\n",
    "    #PR_scen = np.where(data_mask==0, np.nan, PR_scen)\n",
    "    PR_scen = np.where(PR_scen==0, np.nan, PR_scen)\n",
    "    \n",
    "    PR_3 = PR_scen\n",
    "    \n",
    "    del PR_scen\n",
    "    gc.collect()\n",
    "    \n",
    "    return PR_1, PR_2, PR_3\n",
    "\n",
    "def plot_PR_single(ax, data_xarray, title, metric, perc):\n",
    "    cmap='bwr'\n",
    "    vmin=1/4\n",
    "    vmax=4\n",
    "    alpha=0.05\n",
    "    levels = [1/4, 1/3, 1/2, 2/3, 1.5, 2, 3, 4]\n",
    "    im = data_xarray.plot(ax=ax, cmap=cmap,vmin=vmin,vmax=vmax,levels=levels, extend='both', add_colorbar=False, add_labels=False)\n",
    "    cb = plt.colorbar(im,fraction= 0.3, pad= 0.04, extend='both', orientation='horizontal')\n",
    "    cb.ax.tick_params(labelsize=12)\n",
    "    cb.set_ticks([0.25, 0.33333333, 0.5, 0.66666667, 1.5, 2, 3, 4])\n",
    "    cb.set_ticklabels([\"1/4\", \"1/3\", \"1/2\", \"2/3\", \"3/2\", \"2\", \"3\", \"4\"])\n",
    "    ax.coastlines(color='dimgray', linewidth=0.5)\n",
    "    ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "    ax.set_title(title, loc='right', fontsize=14)\n",
    "    ax.set_title(metric + '_' + perc, loc='left', fontsize=14)\n",
    "    ax.add_feature(cfeature.BORDERS, linestyle='-', linewidth=0.1)\n",
    "    \n",
    "def plot_PR_single_no_cb(ax, data_xarray, title, metric, perc):\n",
    "    cmap='bwr'\n",
    "    vmin=1/4\n",
    "    vmax=4\n",
    "    alpha=0.05\n",
    "    levels = [1/4, 1/3, 1/2, 2/3, 1, 1.5, 2, 3, 4]\n",
    "    im = data_xarray.plot(ax=ax, cmap=cmap,vmin=vmin,vmax=vmax,levels=levels, extend='both', add_colorbar=False, add_labels=False)\n",
    "    #cb = plt.colorbar(im,fraction= 0.02, pad= 0.04, extend='both')\n",
    "    #cb.set_ticks([0.25, 0.33333344, 0.5, 0.66666667, 1, 1.5, 2, 3, 4])\n",
    "    #cb.set_ticklabels([\"1/4\", \"1/3\", \"1/2\", \"2/3\", \"1\", \"3/2\", \"2\", \"3\", \"4\"])\n",
    "    ax.coastlines(color='dimgray', linewidth=0.5)\n",
    "    ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "    ax.set_title(title, loc='right', fontsize=10)\n",
    "    ax.set_title(metric + '_' + perc, loc='left', fontsize=10)\n",
    "    ax.add_feature(cfeature.BORDERS, linestyle='-', linewidth=0.1)\n",
    "    \n",
    "def plot_all(PR_1, PR_2, PR_3, metric, perc):\n",
    "    fig = plt.figure(figsize=(15,6),dpi=300)\n",
    "    proj=ccrs.PlateCarree()\n",
    "    data_xr_1 = xr.DataArray(PR_1, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(311, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_1, 'all forcings', metric, perc)\n",
    "    \n",
    "\n",
    "    data_xr_2 = xr.DataArray(PR_2, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(312, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_2, 'other forcings', metric, perc)\n",
    "\n",
    "    data_xr_3 = xr.DataArray(PR_3, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(313, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_3, 'irrigation expansion', metric, perc)\n",
    "\n",
    "\n",
    "    \n",
    "def just_do_all(metric, perc, perc_str):\n",
    "    PR_1, PR_2, PR_3 = calcu_PR_6_ensembles(metric, perc)\n",
    "    plot_all(PR_1, PR_2, PR_3, metric, perc_str)\n",
    "    \n",
    "    \n",
    "def just_do_all_even_more(metric, varname):\n",
    "    PR_1, PR_2, PR_3 = calcu_3_PRs(metric, 0.99)\n",
    "\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/PR_all_\" + \\\n",
    "                     metric + \"_99.0%_0.9x1.25.mat\", {'PR_all': PR_1})\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/PR_oth_\" + \\\n",
    "                     metric + \"_99.0%_0.9x1.25.mat\", {'PR_oth': PR_2})\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/PR_irr_\" + \\\n",
    "                     metric + \"_99.0%_0.9x1.25.mat\", {'PR_irr': PR_3})\n",
    "    \n",
    "    \n",
    "    PR_4, PR_5, PR_6 = calcu_3_PRs(metric, 0.995)\n",
    "\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/PR_all_\" + \\\n",
    "                     metric + \"_99.5%_0.9x1.25.mat\", {'PR_all': PR_4})\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/PR_oth_\" + \\\n",
    "                     metric + \"_99.5%_0.9x1.25.mat\", {'PR_oth': PR_5})\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/PR_irr_\" + \\\n",
    "                     metric + \"_99.5%_0.9x1.25.mat\", {'PR_irr': PR_6})\n",
    "    \n",
    "    \n",
    "    PR_7, PR_8, PR_9 = calcu_3_PRs(metric, 0.999)\n",
    "\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/PR_all_\" + \\\n",
    "                     metric + \"_99.9%_0.9x1.25.mat\", {'PR_all': PR_7})\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/PR_oth_\" + \\\n",
    "                     metric + \"_99.9%_0.9x1.25.mat\", {'PR_oth': PR_8})\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/PR_irr_\" + \\\n",
    "                     metric + \"_99.9%_0.9x1.25.mat\", {'PR_irr': PR_9})\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    fig = plt.figure(figsize=(18,12),dpi=300)\n",
    "    proj=ccrs.PlateCarree()\n",
    "    \n",
    "    \n",
    "    \n",
    "    data_xr_1 = xr.DataArray(PR_1, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(331, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_1, 'all forcings', varname, '99.0%')\n",
    "    ax.text(0.01, 0.92, 'a', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')\n",
    "\n",
    "    \n",
    "    \n",
    "    data_xr_4 = xr.DataArray(PR_4, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(332, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_4, 'all forcings', varname, '99.5%')\n",
    "    ax.text(0.01, 0.92, 'b', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "\n",
    "    \n",
    "    data_xr_7 = xr.DataArray(PR_7, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(333, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_7, 'all forcings', varname, '99.9%')\n",
    "    ax.text(0.01, 0.92, 'c', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "    \n",
    "    \n",
    "    data_xr_2 = xr.DataArray(PR_2, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(334, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_2, 'all forcings', varname, '99.0%')\n",
    "    ax.text(0.01, 0.92, 'a', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')\n",
    "\n",
    "    \n",
    "    \n",
    "    data_xr_5 = xr.DataArray(PR_5, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(335, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_5, 'all forcings', varname, '99.5%')\n",
    "    ax.text(0.01, 0.92, 'b', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "\n",
    "    \n",
    "    data_xr_8 = xr.DataArray(PR_8, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(336, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_8, 'all forcings', varname, '99.9%')\n",
    "    ax.text(0.01, 0.92, 'c', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    data_xr_3 = xr.DataArray(PR_3, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(337, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_3, 'all forcings', varname, '99.0%')\n",
    "    ax.text(0.01, 0.92, 'a', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')\n",
    "\n",
    "    \n",
    "    \n",
    "    data_xr_6 = xr.DataArray(PR_6, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(338, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_6, 'all forcings', varname, '99.5%')\n",
    "    ax.text(0.01, 0.92, 'b', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "\n",
    "    \n",
    "    data_xr_9 = xr.DataArray(PR_9, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(339, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_9, 'all forcings', varname, '99.9%')\n",
    "    ax.text(0.01, 0.92, 'c', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50a4306d",
   "metadata": {},
   "outputs": [],
   "source": [
    "just_do_all_even_more('TSA', 'MIROC T_2m')\n",
    "just_do_all_even_more('AT', 'CESM2 AT')\n",
    "just_do_all_even_more('DI', 'CESM2 DI')\n",
    "just_do_all_even_more('HUMIDEX', 'CESM2 HU')\n",
    "just_do_all_even_more('THIC', 'CESM2 THIC')\n",
    "just_do_all_even_more('THIP', 'CESM2 THIP')\n",
    "just_do_all_even_more('T_s_65', 'CESM2 T_s_65')\n",
    "just_do_all_even_more('T_s_80', 'CESM2 T_s_80')\n",
    "just_do_all_even_more('T_w', 'CESM2 T_w')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}