Event_hours_calcu_E3SM.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "395fd365",
   "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",
    "def get_data_from_mat(file, variable):\n",
    "    var_dict = scio.loadmat(file)\n",
    "    var = var_dict[variable]\n",
    "#     var[var > 1000] = np.nan\n",
    "#     var[var < -1000] = np.nan\n",
    "    var = var[:, 29:] \n",
    "    return var\n",
    "\n",
    "def get_data_from_mat_for_calcu(file, variable):\n",
    "    var_dict = scio.loadmat(file)\n",
    "    var = var_dict[variable]\n",
    "    var = var[:, 29:]\n",
    "    return var.T\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",
    "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",
    "ar6_region = get_data_from_mat_for_calcu('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/ar6_region.mat','ar6_region')\n",
    "irr_diff = get_data_from_mat_for_calcu('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/irr_diff_out.mat', 'irr_diff_out')\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",
    "\n",
    "\n",
    "\n",
    "\n",
    "\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",
    "def get_data_for_year_delta_irr(metric, period, exp, ens, irr_low, irr_high, irr_diff):\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",
    "    data_out[:, irr_diff< irr_low] = np.nan\n",
    "    data_out[:, irr_diff>=irr_high] = np.nan\n",
    "    return data_out\n",
    "\n",
    "AREA = get_data_from_mat('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/input_data/AREA.mat', 'AREA') # Area for CLM_data\n",
    "AREA = AREA.T\n",
    "\n",
    "def get_hour_km2_years(data_period, qctl):\n",
    "    data_time_one_year = np.array(range(2920))\n",
    "    \n",
    "    data_period_year = data_period[0 : 2920] \n",
    "    data_time_one_year = np.array(range(data_period_year.shape[0]))\n",
    "    data_xr_new = xr.DataArray(data_period_year, coords={'time': data_time_one_year, 'y': data_lat, 'x': data_lon}, dims=[\"time\", \"y\", \"x\"])\n",
    "    qscen=xr.where(data_xr_new>qctl,3,0) # three hours resolution\n",
    "    sumscen=qscen.sum('time')\n",
    "    hour_km_2_ = sumscen * AREA\n",
    "    hour_km_2 = np.nansum(hour_km_2_)\n",
    "    \n",
    "    return hour_km_2\n",
    "    \n",
    "def plot_line(ax, index, data_mean, data_min, data_max, title, title1, label, color, ylabel, xlabel, marker):\n",
    "    ax.text(-0.05, 1.05, index, color='dimgrey', fontsize=12, transform=ax.transAxes, weight='bold')\n",
    "    ax.plot(range(1901, 2015), data_mean, label = label, color = color, linewidth = 1, marker = marker, markersize=1.2, alpha=0.8)\n",
    "    ax.fill_between(range(1901, 2015), data_min, data_max, color = color, alpha=0.2)\n",
    "    plt.title(title, loc='right')\n",
    "    plt.title(title1, loc='left')\n",
    "    plt.ylabel(ylabel, fontsize = 14)\n",
    "    plt.xlabel(xlabel, fontsize = 14)\n",
    "    plt.yticks(fontsize = 12)\n",
    "    plt.xticks(fontsize=12)\n",
    "    plt.legend()\n",
    "    \n",
    "\n",
    "def all_delta_irr(metric, perc, unit, variable_to_show, irr_low, irr_high):\n",
    "    data_NOI01_1901_1930 = get_data_for_year_delta_irr(metric, '1901', 'NOI', '01', irr_low, irr_high, irr_diff)\n",
    "    data_NOI02_1901_1930 = get_data_for_year_delta_irr(metric, '1901', 'NOI', '02', irr_low, irr_high, irr_diff)\n",
    "\n",
    "    \n",
    "    for period in range(1902, 1931):\n",
    "        data_NOI01_year = get_data_for_year_delta_irr(metric, str(period), 'NOI', '01', irr_low, irr_high, irr_diff)\n",
    "        data_NOI01_1901_1930 = np.concatenate((data_NOI01_1901_1930, data_NOI01_year), axis = 0)\n",
    "        data_NOI02_year = get_data_for_year_delta_irr(metric, str(period), 'NOI', '02', irr_low, irr_high, irr_diff)\n",
    "        data_NOI02_1901_1930 = np.concatenate((data_NOI02_1901_1930, data_NOI02_year), axis = 0)\n",
    "\n",
    "    \n",
    "    data_NOI_1901_1930 = np.concatenate((data_NOI01_1901_1930,\n",
    "                                         data_NOI02_1901_1930\n",
    "                                         ), axis = 0)\n",
    "    del data_NOI01_1901_1930, data_NOI02_1901_1930\n",
    "    gc.collect()\n",
    "    \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",
    "    qctl=data_xr_ref.quantile(perc,dim='time')\n",
    "    \n",
    "    hour_km_2_IRR01 = np.zeros(114) \n",
    "    hour_km_2_NOI01 = np.zeros(114)\n",
    "    \n",
    "    hour_km_2_IRR02 = np.zeros(114) \n",
    "    hour_km_2_NOI02 = np.zeros(114)\n",
    "    \n",
    "\n",
    "    \n",
    "    for period in range(1901, 2015):\n",
    "        data_IRR01_year = get_data_for_year_delta_irr(metric, str(period), 'IRR', '01', irr_low, irr_high, irr_diff)\n",
    "        data_IRR02_year = get_data_for_year_delta_irr(metric, str(period), 'IRR', '02', irr_low, irr_high, irr_diff)\n",
    "        \n",
    "        \n",
    "        data_NOI01_year = get_data_for_year_delta_irr(metric, str(period), 'NOI', '01', irr_low, irr_high, irr_diff)\n",
    "        data_NOI02_year = get_data_for_year_delta_irr(metric, str(period), 'NOI', '02', irr_low, irr_high, irr_diff)\n",
    "        \n",
    "        hour_km_2_IRR01[period - 1901] = get_hour_km2_years(data_IRR01_year, qctl)\n",
    "        hour_km_2_NOI01[period - 1901] = get_hour_km2_years(data_NOI01_year, qctl)\n",
    "        \n",
    "        hour_km_2_IRR02[period - 1901] = get_hour_km2_years(data_IRR02_year, qctl)\n",
    "        hour_km_2_NOI02[period - 1901] = get_hour_km2_years(data_NOI02_year, qctl)\n",
    "\n",
    "        \n",
    "    data_IRR = np.vstack((hour_km_2_IRR01,\n",
    "                               hour_km_2_IRR02))\n",
    "    data_NOI = np.vstack((hour_km_2_NOI01,\n",
    "                               hour_km_2_NOI02))\n",
    "\n",
    "    return data_IRR, data_NOI\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "42699039",
   "metadata": {},
   "outputs": [],
   "source": [
    "def all(var, percentile, unit, var_to_show, irr_low, irr_high, name, index):\n",
    "    data_IRR, data_NOI=all_delta_irr(var, percentile, unit, var_to_show, irr_low, irr_high)\n",
    "    \n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/data_IRR_\" + \\\n",
    "                     var + '_' + str(percentile) + '_' + str(irr_low) + '_' + str(irr_high) + '.mat', {'data_IRR': data_IRR})\n",
    "    scio.savemat(\"/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/data_NOI_\" + \\\n",
    "                     var + '_' + str(percentile) + '_' + str(irr_low) + '_' + str(irr_high) + '.mat', {'data_NOI': data_NOI})\n",
    "    data_IRR_75 = np.percentile(data_IRR, 75, axis = 0)\n",
    "    data_IRR_50 = np.percentile(data_IRR, 50, axis = 0)\n",
    "    data_IRR_25 = np.percentile(data_IRR, 25, axis = 0)\n",
    "\n",
    "    data_NOI_75 = np.percentile(data_NOI, 75, axis = 0)\n",
    "    data_NOI_50 = np.percentile(data_NOI, 50, axis = 0)\n",
    "    data_NOI_25 = np.percentile(data_NOI, 25, axis = 0)\n",
    "    \n",
    "    fig = plt.figure(figsize=(12, 6),dpi=300)\n",
    "    fig.subplots_adjust(hspace=0.4, wspace=0.2, left = 0.05, right = 0.95, top = 0.95, bottom = 0.05)\n",
    "    ax1 = plt.subplot(224, frameon=True)\n",
    "    plt.title(name, loc='right')\n",
    "    ax1.spines[\"top\"].set_visible(False)\n",
    "    ax1.spines[\"right\"].set_visible(False)\n",
    "    plot_line(ax1, index, data_NOI_50, data_NOI_25, data_NOI_75, '', '', '1901irr', 'r', r'$\\Delta$(km^3*hours)', 'year', '')\n",
    "    plot_line(ax1, index, data_IRR_50, data_IRR_25, data_IRR_75, '', '', 'tranirr', 'b', r'$\\Delta$(km^3*hours)', 'year', '')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58f46161",
   "metadata": {},
   "outputs": [],
   "source": [
    "all('T_w', 0.99, r'$\\Delta$$\\mathregular{10^6 km^3}$ x hours', 'T_2m 99%', 0.01, 20, '', 'a ' + 'SAS')\n",
    "all('T_w', 0.99, r'$\\Delta$$\\mathregular{10^6 km^3}$ x hours', 'T_2m 99%', 20, 30, '', 'a ' + 'SAS')\n",
    "all('T_w', 0.99, r'$\\Delta$$\\mathregular{10^6 km^3}$ x hours', 'T_2m 99%', 30, 40, '', 'a ' + 'SAS')\n",
    "all('T_w', 0.99, r'$\\Delta$$\\mathregular{10^6 km^3}$ x hours', 'T_2m 99%', 40, 100, '', 'a ' + 'SAS')\n",
    "all('T_w', 0.995, r'$\\Delta$$\\mathregular{10^6 km^3}$ x hours', 'T_2m 99%', 0.01, 20, '', 'a ' + 'SAS')\n",
    "all('T_w', 0.995, r'$\\Delta$$\\mathregular{10^6 km^3}$ x hours', 'T_2m 99%', 20, 30, '', 'a ' + 'SAS')\n",
    "all('T_w', 0.995, r'$\\Delta$$\\mathregular{10^6 km^3}$ x hours', 'T_2m 99%', 30, 40, '', 'a ' + 'SAS')\n",
    "all('T_w', 0.995, r'$\\Delta$$\\mathregular{10^6 km^3}$ x hours', 'T_2m 99%', 40, 100, '', 'a ' + 'SAS')\n",
    "all('T_w', 0.999, r'$\\Delta$$\\mathregular{10^6 km^3}$ x hours', 'T_2m 99%', 0.01, 20, '', 'a ' + 'SAS')\n",
    "all('T_w', 0.999, r'$\\Delta$$\\mathregular{10^6 km^3}$ x hours', 'T_2m 99%', 20, 30, '', 'a ' + 'SAS')\n",
    "all('T_w', 0.999, r'$\\Delta$$\\mathregular{10^6 km^3}$ x hours', 'T_2m 99%', 30, 40, '', 'a ' + 'SAS')\n",
    "all('T_w', 0.999, r'$\\Delta$$\\mathregular{10^6 km^3}$ x hours', 'T_2m 99%', 40, 100, '', 'a ' + 'SAS')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef2b6ef5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}