ProbabilityValue_calcu_CESM2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3c23fa9b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "node5447.dodrio.os:rank0.python: Failed to modify UD QP to INIT on mlx5_0: Operation not permitted\n",
      "node5447.dodrio.os:rank0.python: Failed to modify UD QP to INIT on mlx5_0: Operation not permitted\n",
      "node5447.dodrio.os:rank0.python: Failed to modify UD QP to INIT on mlx5_0: Operation not permitted\n",
      "node5447.dodrio.os:rank0.python: Failed to modify UD QP to INIT on mlx5_0: Operation not permitted\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "--------------------------------------------------------------------------\n",
      "Open MPI failed an OFI Libfabric library call (fi_endpoint).  This is highly\n",
      "unusual; your job may behave unpredictably (and/or abort) after this.\n",
      "\n",
      "  Local host: node5447\n",
      "  Location: mtl_ofi_component.c:513\n",
      "  Error: Invalid argument (22)\n",
      "--------------------------------------------------------------------------\n"
     ]
    }
   ],
   "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_this_study_1(metric, period, exp, ens):\n",
    "    data_out = get_data_from_mat('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/post_process/3_hourly/' +\\\n",
    "                                             metric + '_' + period + '_3_hourly_' + exp + ens + '_1.mat', metric)\n",
    "    data_out = np.array(data_out)\n",
    "    return data_out\n",
    "\n",
    "def get_data_for_this_study_2(metric, period, exp, ens):\n",
    "    data_out = get_data_from_mat('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/post_process/3_hourly/' +\\\n",
    "                                             metric + '_' + period + '_3_hourly_' + exp + ens + '_2.mat', metric)\n",
    "    data_out = np.array(data_out)\n",
    "    return data_out\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_PR_6_ensembles(metric):\n",
    "    \n",
    "    data_IRR01_1901_1910_1 = get_data_for_this_study_1(metric, '1901_1910', 'IRR', '01')\n",
    "    data_IRR01_1901_1910_2 = get_data_for_this_study_2(metric, '1901_1910', 'IRR', '01')\n",
    "    \n",
    "    data_IRR01_1911_1920_1 = get_data_for_this_study_1(metric, '1911_1920', 'IRR', '01')\n",
    "    data_IRR01_1911_1920_2 = get_data_for_this_study_2(metric, '1911_1920', 'IRR', '01')\n",
    "    \n",
    "    data_IRR01_1921_1930_1 = get_data_for_this_study_1(metric, '1921_1930', 'IRR', '01')\n",
    "    data_IRR01_1921_1930_2 = get_data_for_this_study_2(metric, '1921_1930', 'IRR', '01')\n",
    "    \n",
    "    data_IRR02_1901_1910_1 = get_data_for_this_study_1(metric, '1901_1910', 'IRR', '02')\n",
    "    data_IRR02_1901_1910_2 = get_data_for_this_study_2(metric, '1901_1910', 'IRR', '02')\n",
    "    \n",
    "    data_IRR02_1911_1920_1 = get_data_for_this_study_1(metric, '1911_1920', 'IRR', '02')\n",
    "    data_IRR02_1911_1920_2 = get_data_for_this_study_2(metric, '1911_1920', 'IRR', '02')\n",
    "    \n",
    "    data_IRR02_1921_1930_1 = get_data_for_this_study_1(metric, '1921_1930', 'IRR', '02')\n",
    "    data_IRR02_1921_1930_2 = get_data_for_this_study_2(metric, '1921_1930', 'IRR', '02')\n",
    "    \n",
    "    data_IRR03_1901_1910_1 = get_data_for_this_study_1(metric, '1901_1910', 'IRR', '03')\n",
    "    data_IRR03_1901_1910_2 = get_data_for_this_study_2(metric, '1901_1910', 'IRR', '03')\n",
    "    \n",
    "    data_IRR03_1911_1920_1 = get_data_for_this_study_1(metric, '1911_1920', 'IRR', '03')\n",
    "    data_IRR03_1911_1920_2 = get_data_for_this_study_2(metric, '1911_1920', 'IRR', '03')\n",
    "    \n",
    "    data_IRR03_1921_1930_1 = get_data_for_this_study_1(metric, '1921_1930', 'IRR', '03')\n",
    "    data_IRR03_1921_1930_2 = get_data_for_this_study_2(metric, '1921_1930', 'IRR', '03')\n",
    "    \n",
    "    data_IRR_1901_1930 = np.concatenate((data_IRR01_1901_1910_1,\n",
    "                                         data_IRR01_1901_1910_2,\n",
    "                                         data_IRR01_1911_1920_1,\n",
    "                                         data_IRR01_1911_1920_2,\n",
    "                                         data_IRR01_1921_1930_1,\n",
    "                                         data_IRR01_1921_1930_2,\n",
    "                                         \n",
    "                                         data_IRR02_1901_1910_1,\n",
    "                                         data_IRR02_1901_1910_2,\n",
    "                                         data_IRR02_1911_1920_1,\n",
    "                                         data_IRR02_1911_1920_2,\n",
    "                                         data_IRR02_1921_1930_1,\n",
    "                                         data_IRR02_1921_1930_2,\n",
    "                                         \n",
    "                                         data_IRR03_1901_1910_1,\n",
    "                                         data_IRR03_1901_1910_2,\n",
    "                                         data_IRR03_1911_1920_1,\n",
    "                                         data_IRR03_1911_1920_2,\n",
    "                                         data_IRR03_1921_1930_1,\n",
    "                                         data_IRR03_1921_1930_2), axis = 0)\n",
    "    \n",
    "    del data_IRR01_1901_1910_1, data_IRR01_1901_1910_2, data_IRR02_1901_1910_1, data_IRR02_1901_1910_2, data_IRR03_1901_1910_1, data_IRR03_1901_1910_2\n",
    "    gc.collect()\n",
    "    \n",
    "    del data_IRR01_1911_1920_1, data_IRR01_1911_1920_2, data_IRR02_1911_1920_1, data_IRR02_1911_1920_2, data_IRR03_1911_1920_1, data_IRR03_1911_1920_2\n",
    "    gc.collect()\n",
    "    \n",
    "    del data_IRR01_1921_1930_1, data_IRR01_1921_1930_2, data_IRR02_1921_1930_1, data_IRR02_1921_1930_2, data_IRR03_1921_1930_1, data_IRR03_1921_1930_2\n",
    "    gc.collect()\n",
    "    \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",
    "    qctl_IRR_1901_1930_99 = data_xr_ref.quantile(0.99,dim='time')\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/PercentileValue_IRR_1901_1930_\" + \\\n",
    "                     metric + \"_99.0%.mat\", {'qctl_IRR_1901_1930_99': qctl_IRR_1901_1930_99})\n",
    "    \n",
    "    qctl_IRR_1901_1930_995 = data_xr_ref.quantile(0.995,dim='time')\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/PercentileValue_IRR_1901_1930_\" + \\\n",
    "                     metric + \"_99.5%.mat\", {'qctl_IRR_1901_1930_995': qctl_IRR_1901_1930_995})\n",
    "    \n",
    "    qctl_IRR_1901_1930_999 = data_xr_ref.quantile(0.999,dim='time')\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/PercentileValue_IRR_1901_1930_\" + \\\n",
    "                     metric + \"_99.9%.mat\", {'qctl_IRR_1901_1930_999': qctl_IRR_1901_1930_999})\n",
    "    \n",
    "    del data_IRR_1901_1930\n",
    "    gc.collect()\n",
    "    \n",
    "    data_IRR01_1985_1994_1 = get_data_for_this_study_1(metric, '1985_1994', 'IRR', '01')\n",
    "    data_IRR01_1985_1994_2 = get_data_for_this_study_2(metric, '1985_1994', 'IRR', '01')\n",
    "    \n",
    "    data_IRR01_1995_2004_1 = get_data_for_this_study_1(metric, '1995_2004', 'IRR', '01')\n",
    "    data_IRR01_1995_2004_2 = get_data_for_this_study_2(metric, '1995_2004', 'IRR', '01')\n",
    "    \n",
    "    data_IRR01_2005_2014_1 = get_data_for_this_study_1(metric, '2005_2014', 'IRR', '01')\n",
    "    data_IRR01_2005_2014_2 = get_data_for_this_study_2(metric, '2005_2014', 'IRR', '01')\n",
    "    \n",
    "    data_IRR02_1985_1994_1 = get_data_for_this_study_1(metric, '1985_1994', 'IRR', '02')\n",
    "    data_IRR02_1985_1994_2 = get_data_for_this_study_2(metric, '1985_1994', 'IRR', '02')\n",
    "    \n",
    "    data_IRR02_1995_2004_1 = get_data_for_this_study_1(metric, '1995_2004', 'IRR', '02')\n",
    "    data_IRR02_1995_2004_2 = get_data_for_this_study_2(metric, '1995_2004', 'IRR', '02')\n",
    "    \n",
    "    data_IRR02_2005_2014_1 = get_data_for_this_study_1(metric, '2005_2014', 'IRR', '02')\n",
    "    data_IRR02_2005_2014_2 = get_data_for_this_study_2(metric, '2005_2014', 'IRR', '02')\n",
    "    \n",
    "    data_IRR03_1985_1994_1 = get_data_for_this_study_1(metric, '1985_1994', 'IRR', '03')\n",
    "    data_IRR03_1985_1994_2 = get_data_for_this_study_2(metric, '1985_1994', 'IRR', '03')\n",
    "    \n",
    "    data_IRR03_1995_2004_1 = get_data_for_this_study_1(metric, '1995_2004', 'IRR', '03')\n",
    "    data_IRR03_1995_2004_2 = get_data_for_this_study_2(metric, '1995_2004', 'IRR', '03')\n",
    "    \n",
    "    data_IRR03_2005_2014_1 = get_data_for_this_study_1(metric, '2005_2014', 'IRR', '03')\n",
    "    data_IRR03_2005_2014_2 = get_data_for_this_study_2(metric, '2005_2014', 'IRR', '03')\n",
    "    \n",
    "\n",
    "    \n",
    "    \n",
    "    data_IRR_1985_2014 = np.concatenate((data_IRR01_1985_1994_1,\n",
    "                                         data_IRR01_1985_1994_2,\n",
    "                                         data_IRR01_1995_2004_1,\n",
    "                                         data_IRR01_1995_2004_2,\n",
    "                                         data_IRR01_2005_2014_1,\n",
    "                                         data_IRR01_2005_2014_2,\n",
    "                                         \n",
    "                                         data_IRR02_1985_1994_1,\n",
    "                                         data_IRR02_1985_1994_2,\n",
    "                                         data_IRR02_1995_2004_1,\n",
    "                                         data_IRR02_1995_2004_2,\n",
    "                                         data_IRR02_2005_2014_1,\n",
    "                                         data_IRR02_2005_2014_2,\n",
    "                                         \n",
    "                                         data_IRR03_1985_1994_1,\n",
    "                                         data_IRR03_1985_1994_2,\n",
    "                                         data_IRR03_1995_2004_1,\n",
    "                                         data_IRR03_1995_2004_2,\n",
    "                                         data_IRR03_2005_2014_1,\n",
    "                                         data_IRR03_2005_2014_2,), axis = 0)\n",
    "    \n",
    "    del data_IRR01_1985_1994_1, data_IRR01_1985_1994_2, data_IRR02_1985_1994_1, data_IRR02_1985_1994_2, data_IRR03_1985_1994_1, data_IRR03_1985_1994_2\n",
    "    gc.collect()\n",
    "    \n",
    "    del data_IRR01_1995_2004_1, data_IRR01_1995_2004_2, data_IRR02_1995_2004_1, data_IRR02_1995_2004_2, data_IRR03_1995_2004_1, data_IRR03_1995_2004_2\n",
    "    gc.collect()\n",
    "\n",
    "    del data_IRR01_2005_2014_1, data_IRR01_2005_2014_2, data_IRR02_2005_2014_1, data_IRR02_2005_2014_2, data_IRR03_2005_2014_1, data_IRR03_2005_2014_2\n",
    "    gc.collect()\n",
    "    \n",
    "    data_xr_ref = xr.DataArray(data_IRR_1985_2014, coords={'time': data_time, 'y': data_lat, 'x': data_lon}, dims=[\"time\", \"y\", \"x\"])\n",
    "    \n",
    "    qctl_IRR_1985_2014_99 = data_xr_ref.quantile(0.99,dim='time')\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/PercentileValue_IRR_1985_2014_\" + \\\n",
    "                     metric + \"_99.0%.mat\", {'qctl_IRR_1985_2014_99': qctl_IRR_1985_2014_99})\n",
    "    qctl_IRR_1985_2014_995 = data_xr_ref.quantile(0.995,dim='time')\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/PercentileValue_IRR_1985_2014_\" + \\\n",
    "                     metric + \"_99.5%.mat\", {'qctl_IRR_1985_2014_995': qctl_IRR_1985_2014_995})\n",
    "    qctl_IRR_1985_2014_999 = data_xr_ref.quantile(0.999,dim='time')\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/PercentileValue_IRR_1985_2014_\" + \\\n",
    "                     metric + \"_99.9%.mat\", {'qctl_IRR_1985_2014_999': qctl_IRR_1985_2014_999})\n",
    "    \n",
    "    del data_IRR_1985_2014\n",
    "    gc.collect()\n",
    "    \n",
    "    \n",
    "    data_NOI01_1985_1994_1 = get_data_for_this_study_1(metric, '1985_1994', 'NOI', '01')\n",
    "    data_NOI01_1985_1994_2 = get_data_for_this_study_2(metric, '1985_1994', 'NOI', '01')\n",
    "    \n",
    "    data_NOI01_1995_2004_1 = get_data_for_this_study_1(metric, '1995_2004', 'NOI', '01')\n",
    "    data_NOI01_1995_2004_2 = get_data_for_this_study_2(metric, '1995_2004', 'NOI', '01')\n",
    "    \n",
    "    data_NOI01_2005_2014_1 = get_data_for_this_study_1(metric, '2005_2014', 'NOI', '01')\n",
    "    data_NOI01_2005_2014_2 = get_data_for_this_study_2(metric, '2005_2014', 'NOI', '01')\n",
    "\n",
    "    data_NOI02_1985_1994_1 = get_data_for_this_study_1(metric, '1985_1994', 'NOI', '02')\n",
    "    data_NOI02_1985_1994_2 = get_data_for_this_study_2(metric, '1985_1994', 'NOI', '02')\n",
    "    \n",
    "    data_NOI02_1995_2004_1 = get_data_for_this_study_1(metric, '1995_2004', 'NOI', '02')\n",
    "    data_NOI02_1995_2004_2 = get_data_for_this_study_2(metric, '1995_2004', 'NOI', '02')\n",
    "    \n",
    "    data_NOI02_2005_2014_1 = get_data_for_this_study_1(metric, '2005_2014', 'NOI', '02')\n",
    "    data_NOI02_2005_2014_2 = get_data_for_this_study_2(metric, '2005_2014', 'NOI', '02')\n",
    "\n",
    "    data_NOI03_1985_1994_1 = get_data_for_this_study_1(metric, '1985_1994', 'NOI', '03')\n",
    "    data_NOI03_1985_1994_2 = get_data_for_this_study_2(metric, '1985_1994', 'NOI', '03')\n",
    "    \n",
    "    data_NOI03_1995_2004_1 = get_data_for_this_study_1(metric, '1995_2004', 'NOI', '03')\n",
    "    data_NOI03_1995_2004_2 = get_data_for_this_study_2(metric, '1995_2004', 'NOI', '03')\n",
    "    \n",
    "    data_NOI03_2005_2014_1 = get_data_for_this_study_1(metric, '2005_2014', 'NOI', '03')\n",
    "    data_NOI03_2005_2014_2 = get_data_for_this_study_2(metric, '2005_2014', 'NOI', '03')\n",
    "    \n",
    "    \n",
    "    \n",
    "    data_NOI_1985_2014 = np.concatenate((data_NOI01_1985_1994_1,\n",
    "                                         data_NOI01_1985_1994_2,\n",
    "                                         data_NOI01_1995_2004_1,\n",
    "                                         data_NOI01_1995_2004_2,\n",
    "                                         data_NOI01_2005_2014_1,\n",
    "                                         data_NOI01_2005_2014_2,\n",
    "                                         \n",
    "                                         data_NOI02_1985_1994_1,\n",
    "                                         data_NOI02_1985_1994_2,\n",
    "                                         data_NOI02_1995_2004_1,\n",
    "                                         data_NOI02_1995_2004_2,\n",
    "                                         data_NOI02_2005_2014_1,\n",
    "                                         data_NOI02_2005_2014_2,\n",
    "                                         \n",
    "                                         data_NOI03_1985_1994_1,\n",
    "                                         data_NOI03_1985_1994_2,\n",
    "                                         data_NOI03_1995_2004_1,\n",
    "                                         data_NOI03_1995_2004_2,\n",
    "                                         data_NOI03_2005_2014_1,\n",
    "                                         data_NOI03_2005_2014_2,), axis = 0)\n",
    "    \n",
    "    del data_NOI01_1985_1994_1, data_NOI01_1985_1994_2, data_NOI02_1985_1994_1, data_NOI02_1985_1994_2, data_NOI03_1985_1994_1, data_NOI03_1985_1994_2\n",
    "    gc.collect()\n",
    "    \n",
    "    del data_NOI01_1995_2004_1, data_NOI01_1995_2004_2, data_NOI02_1995_2004_1, data_NOI02_1995_2004_2, data_NOI03_1995_2004_1, data_NOI03_1995_2004_2\n",
    "    gc.collect()\n",
    "\n",
    "    del data_NOI01_2005_2014_1, data_NOI01_2005_2014_2, data_NOI02_2005_2014_1, data_NOI02_2005_2014_2, data_NOI03_2005_2014_1, data_NOI03_2005_2014_2\n",
    "    gc.collect()\n",
    "    \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",
    "    qctl_NOI_1985_2014_99 = data_xr_ref.quantile(0.99,dim='time')\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/PercentileValue_NOI_1985_2014_\" + \\\n",
    "                     metric + \"_99.0%.mat\", {'qctl_NOI_1985_2014_99': qctl_NOI_1985_2014_99})\n",
    "    qctl_NOI_1985_2014_995 = data_xr_ref.quantile(0.995,dim='time')\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/PercentileValue_NOI_1985_2014_\" + \\\n",
    "                     metric + \"_99.5%.mat\", {'qctl_NOI_1985_2014_995': qctl_NOI_1985_2014_995})\n",
    "    qctl_NOI_1985_2014_999 = data_xr_ref.quantile(0.999,dim='time')\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/PercentileValue_NOI_1985_2014_\" + \\\n",
    "                     metric + \"_99.9%.mat\", {'qctl_NOI_1985_2014_999': qctl_NOI_1985_2014_999})\n",
    "    del data_NOI_1985_2014\n",
    "    gc.collect()\n",
    "    \n",
    "    \n",
    "    data_NOI01_1901_1910_1 = get_data_for_this_study_1(metric, '1901_1910', 'NOI', '01')\n",
    "    data_NOI01_1901_1910_2 = get_data_for_this_study_2(metric, '1901_1910', 'NOI', '01')\n",
    "    \n",
    "    data_NOI01_1911_1920_1 = get_data_for_this_study_1(metric, '1911_1920', 'NOI', '01')\n",
    "    data_NOI01_1911_1920_2 = get_data_for_this_study_2(metric, '1911_1920', 'NOI', '01')\n",
    "    \n",
    "    data_NOI01_1921_1930_1 = get_data_for_this_study_1(metric, '1921_1930', 'NOI', '01')\n",
    "    data_NOI01_1921_1930_2 = get_data_for_this_study_2(metric, '1921_1930', 'NOI', '01')\n",
    "    \n",
    "    data_NOI02_1901_1910_1 = get_data_for_this_study_1(metric, '1901_1910', 'NOI', '02')\n",
    "    data_NOI02_1901_1910_2 = get_data_for_this_study_2(metric, '1901_1910', 'NOI', '02')\n",
    "    \n",
    "    data_NOI02_1911_1920_1 = get_data_for_this_study_1(metric, '1911_1920', 'NOI', '02')\n",
    "    data_NOI02_1911_1920_2 = get_data_for_this_study_2(metric, '1911_1920', 'NOI', '02')\n",
    "    \n",
    "    data_NOI02_1921_1930_1 = get_data_for_this_study_1(metric, '1921_1930', 'NOI', '02')\n",
    "    data_NOI02_1921_1930_2 = get_data_for_this_study_2(metric, '1921_1930', 'NOI', '02')\n",
    "    \n",
    "    data_NOI03_1901_1910_1 = get_data_for_this_study_1(metric, '1901_1910', 'NOI', '03')\n",
    "    data_NOI03_1901_1910_2 = get_data_for_this_study_2(metric, '1901_1910', 'NOI', '03')\n",
    "    \n",
    "    data_NOI03_1911_1920_1 = get_data_for_this_study_1(metric, '1911_1920', 'NOI', '03')\n",
    "    data_NOI03_1911_1920_2 = get_data_for_this_study_2(metric, '1911_1920', 'NOI', '03')\n",
    "    \n",
    "    data_NOI03_1921_1930_1 = get_data_for_this_study_1(metric, '1921_1930', 'NOI', '03')\n",
    "    data_NOI03_1921_1930_2 = get_data_for_this_study_2(metric, '1921_1930', 'NOI', '03')\n",
    "    \n",
    "    \n",
    "    \n",
    "    data_NOI_1901_1930 = np.concatenate((data_NOI01_1901_1910_1,\n",
    "                                         data_NOI01_1901_1910_2,\n",
    "                                         data_NOI01_1911_1920_1,\n",
    "                                         data_NOI01_1911_1920_2,\n",
    "                                         data_NOI01_1921_1930_1,\n",
    "                                         data_NOI01_1921_1930_2,\n",
    "                                         \n",
    "                                         data_NOI02_1901_1910_1,\n",
    "                                         data_NOI02_1901_1910_2,\n",
    "                                         data_NOI02_1911_1920_1,\n",
    "                                         data_NOI02_1911_1920_2,\n",
    "                                         data_NOI02_1921_1930_1,\n",
    "                                         data_NOI02_1921_1930_2,\n",
    "                                         \n",
    "                                         data_NOI03_1901_1910_1,\n",
    "                                         data_NOI03_1901_1910_2,\n",
    "                                         data_NOI03_1911_1920_1,\n",
    "                                         data_NOI03_1911_1920_2,\n",
    "                                         data_NOI03_1921_1930_1,\n",
    "                                         data_NOI03_1921_1930_2), axis = 0)\n",
    "    del data_NOI01_1901_1910_1, data_NOI01_1901_1910_2, data_NOI02_1901_1910_1, data_NOI02_1901_1910_2, data_NOI03_1901_1910_1, data_NOI03_1901_1910_2\n",
    "    gc.collect()\n",
    "\n",
    "    del data_NOI01_1911_1920_1, data_NOI01_1911_1920_2, data_NOI02_1911_1920_1, data_NOI02_1911_1920_2, data_NOI03_1911_1920_1, data_NOI03_1911_1920_2\n",
    "    gc.collect()\n",
    "    \n",
    "    del data_NOI01_1921_1930_1, data_NOI01_1921_1930_2, data_NOI02_1921_1930_1, data_NOI02_1921_1930_2, data_NOI03_1921_1930_1, data_NOI03_1921_1930_2\n",
    "    gc.collect()\n",
    "    \n",
    "    \n",
    "    \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",
    "    qctl_NOI_1901_1930_99 = data_xr_ref.quantile(0.99,dim='time')\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/PercentileValue_NOI_1901_1930_\" + \\\n",
    "                     metric + \"_99.0%.mat\", {'qctl_NOI_1901_1930_99': qctl_NOI_1901_1930_99})\n",
    "    qctl_NOI_1901_1930_995 = data_xr_ref.quantile(0.995,dim='time')\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/PercentileValue_NOI_1901_1930_\" + \\\n",
    "                     metric + \"_99.5%.mat\", {'qctl_NOI_1901_1930_995': qctl_NOI_1901_1930_995})\n",
    "    qctl_NOI_1901_1930_999 = data_xr_ref.quantile(0.999,dim='time')\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/PercentileValue_NOI_1901_1930_\" + \\\n",
    "                     metric + \"_99.9%.mat\", {'qctl_NOI_1901_1930_999': qctl_NOI_1901_1930_999})\n",
    "    del data_NOI_1901_1930\n",
    "    gc.collect()\n",
    "    \n",
    "    qctl_ALL_impact_99 = qctl_IRR_1985_2014_99 - qctl_IRR_1901_1930_99\n",
    "    qctl_OTH_impact_99 = qctl_NOI_1985_2014_99 - qctl_NOI_1901_1930_99\n",
    "    qctl_IRR_impact_99 = qctl_IRR_1985_2014_99 - qctl_NOI_1985_2014_99\n",
    "    \n",
    "    qctl_ALL_impact_995 = qctl_IRR_1985_2014_995 - qctl_IRR_1901_1930_995\n",
    "    qctl_OTH_impact_995 = qctl_NOI_1985_2014_995 - qctl_NOI_1901_1930_995\n",
    "    qctl_IRR_impact_995 = qctl_IRR_1985_2014_995 - qctl_NOI_1985_2014_995\n",
    "    \n",
    "    qctl_ALL_impact_999 = qctl_IRR_1985_2014_999 - qctl_IRR_1901_1930_999\n",
    "    qctl_OTH_impact_999 = qctl_NOI_1985_2014_999 - qctl_NOI_1901_1930_999\n",
    "    qctl_IRR_impact_999 = qctl_IRR_1985_2014_999 - qctl_NOI_1985_2014_999\n",
    "    \n",
    "    return qctl_NOI_1901_1930_99, qctl_NOI_1901_1930_995, qctl_NOI_1901_1930_999, qctl_IRR_1901_1930_99, qctl_IRR_1901_1930_995, qctl_IRR_1901_1930_999, qctl_NOI_1985_2014_99, qctl_NOI_1985_2014_995, qctl_NOI_1985_2014_999, qctl_IRR_1985_2014_99, qctl_IRR_1985_2014_995, qctl_IRR_1985_2014_999\n",
    "    \n",
    "    \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 = [-4, -2, -1, -0.2, 0.2, 1, 2, 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([-4, -2, -1, -0.2, 0.2, 1, 2, 4])\n",
    "    cb.set_ticklabels([\"-4\", \"-2\", \"-1\", \"-0.2\", \"+0.2\", \"+1\", \"+2\", \"+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 just_do_all_even_more(metric, varname):\n",
    "    qctl_NOI_1901_1930_99, qctl_NOI_1901_1930_995, qctl_NOI_1901_1930_999, qctl_IRR_1901_1930_99, qctl_IRR_1901_1930_995, qctl_IRR_1901_1930_999, qctl_NOI_1985_2014_99, qctl_NOI_1985_2014_995, qctl_NOI_1985_2014_999, qctl_IRR_1985_2014_99, qctl_IRR_1985_2014_995, qctl_IRR_1985_2014_999 = calcu_PR_6_ensembles(metric)\n",
    "    \n",
    "    fig = plt.figure(figsize=(18,12),dpi=300)\n",
    "    proj=ccrs.PlateCarree()\n",
    "    \n",
    "    \n",
    "    data_xr_1 = xr.DataArray(qctl_IRR_1985_2014_99-qctl_IRR_1901_1930_99, 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",
    "    data_xr_2 = xr.DataArray(qctl_NOI_1985_2014_99-qctl_NOI_1901_1930_99, 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, 'other forcings', varname, '99.0%')\n",
    "    ax.text(0.01, 0.92, 'd', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "    data_xr_3 = xr.DataArray(qctl_IRR_1985_2014_99-qctl_NOI_1985_2014_99, 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, 'irrigation expansion', varname, '99.0%')\n",
    "    ax.text(0.01, 0.92, 'g', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "    \n",
    "    data_xr_4 = xr.DataArray(qctl_IRR_1985_2014_995-qctl_IRR_1901_1930_995, 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",
    "    data_xr_5 = xr.DataArray(qctl_NOI_1985_2014_995-qctl_NOI_1901_1930_995, 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, 'other forcings', varname, '99.5%')\n",
    "    ax.text(0.01, 0.92, 'e', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')\n",
    "\n",
    "    data_xr_6 = xr.DataArray(qctl_IRR_1985_2014_995-qctl_NOI_1985_2014_995, 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, 'irrigation expansion', varname, '99.5%')\n",
    "    ax.text(0.01, 0.92, 'h', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "    \n",
    "    \n",
    "    data_xr_7 = xr.DataArray(qctl_IRR_1985_2014_999-qctl_IRR_1901_1930_999, 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",
    "    data_xr_8 = xr.DataArray(qctl_NOI_1985_2014_999-qctl_NOI_1901_1930_999, 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, 'other forcings', varname, '99.9%')\n",
    "    ax.text(0.01, 0.92, 'f', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')\n",
    "\n",
    "    data_xr_9 = xr.DataArray(qctl_IRR_1985_2014_999-qctl_NOI_1985_2014_999, 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, 'irrigation expansion', varname, '99.9%')\n",
    "    ax.text(0.01, 0.92, 'i', color='dimgrey', fontsize=10, transform=ax.transAxes, weight='bold')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a7a75e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "just_do_all_even_more('TSA', 'CESM2 T_2m')\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
}