bp_wetbulb.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "abdba210",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<p>To authorize access needed by Earth Engine, open the following\n",
       "        URL in a web browser and follow the instructions:</p>\n",
       "        <p><a href=https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=eNZ9jp3CZx8DOEc5GgFRg6Z80X62UnKvdA0GMQzzU3I&tc=0mwRBg9YXO3hlZ8yAegy2wjVdYtL38jeytHN1Ulquqg&cc=IGbtqgsh5IsgiWVyMVZdPhRqRV5Sm3qWfstQ1wRwOhM>https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=eNZ9jp3CZx8DOEc5GgFRg6Z80X62UnKvdA0GMQzzU3I&tc=0mwRBg9YXO3hlZ8yAegy2wjVdYtL38jeytHN1Ulquqg&cc=IGbtqgsh5IsgiWVyMVZdPhRqRV5Sm3qWfstQ1wRwOhM</a></p>\n",
       "        <p>The authorization workflow will generate a code, which you\n",
       "        should paste in the box below</p>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enter verification code: 4/1AfJohXlJADDa3treYN1Vg10GHhJLTDrssM288FmqD713roG8HdrkWKOoFbQ\n",
      "\n",
      "Successfully saved authorization token.\n"
     ]
    }
   ],
   "source": [
    "import ee\n",
    "%matplotlib inline\n",
    "import math\n",
    "import warnings\n",
    "import intake\n",
    "warnings.filterwarnings('ignore')\n",
    "from collections import defaultdict\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from scipy import stats\n",
    "import cftime\n",
    "import xarray as xr\n",
    "xr.set_options(display_style='html')\n",
    "\n",
    "import datetime, calendar\n",
    "\n",
    "plt.style.use(\"seaborn-darkgrid\")\n",
    "\n",
    "ee.Authenticate()\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "28f44297",
   "metadata": {},
   "outputs": [],
   "source": [
    "HIST_START = 1980\n",
    "HIST_END = 2014\n",
    "\n",
    "PERCENTILE_STARTYEAR = 1980\n",
    "PERCENTILE_ENDYEAR = 2019\n",
    "\n",
    "NUM_BEST_MODELS = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9951f379",
   "metadata": {},
   "outputs": [],
   "source": [
    "FUTURE_YEARS = [(2080, 2099)]\n",
    "FUTURE_SCENARIOS = ['ssp119', 'ssp126']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "24b29c18",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('citiesonly.csv', 'r', encoding='utf-8') as ifile:\n",
    "    CITYLATLON = {i[0]: (float(i[1]), float(i[2]), i[3]) for i in [j.split(',') + [idx] for idx, j in enumerate(ifile.readlines()[1:])]}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5d04f9a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "MODELS = {\n",
    "    'tasmax': ('GFDL-ESM4', 'CanESM5', 'MRI-ESM2-0', 'IPSL-CM6A-LR', 'EC-Earth3-Veg-LR'),\n",
    "    'tasmin': ('GFDL-ESM4', 'IPSL-CM6A-LR', 'CanESM5', 'MRI-ESM2-0', 'EC-Earth3-Veg-LR'),\n",
    "    'pr': ('GFDL-ESM4', 'IPSL-CM6A-LR', 'CanESM5', 'MRI-ESM2-0', 'EC-Earth3-Veg-LR'),\n",
    "    'hurs': ('GFDL-ESM4', 'CanESM5', 'MRI-ESM2-0', 'IPSL-CM6A-LR', 'EC-Earth3-Veg-LR')\n",
    "}\n",
    "\n",
    "MODELRUN = 'r1i1p1f1'\n",
    "\n",
    "MODELGRID = {\n",
    "    'GFDL-ESM4': 'gr1',\n",
    "    'CanESM5': 'gn',\n",
    "    'MRI-ESM2-0': 'gn',\n",
    "    'IPSL-CM6A-LR': 'gr',\n",
    "    'EC-Earth3-Veg-LR': 'gr'\n",
    "}\n",
    "\n",
    "YEARLENGTH = {\n",
    "    'GFDL-ESM4': 365,\n",
    "    'CanESM5': 365,\n",
    "    'MRI-ESM2-0': 366,\n",
    "    'IPSL-CM6A-LR': 366,\n",
    "    'EC-Earth3-Veg-LR': 366\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f32606e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "MODEL_FAMILY = {'UKESM1-0-LL': 'HadAM',\n",
    " 'NorESM2-MM': 'CCM',\n",
    " 'NorESM2-LM': 'CCM',\n",
    " 'MRI-ESM2-0': 'UCLA GCM',\n",
    " 'MPI-ESM1-2-LR': 'ECMWF',\n",
    " 'MPI-ESM1-2-HR': 'ECMWF',\n",
    " 'MIROC6': 'MIROC',\n",
    " 'MIROC-ES2L': 'MIROC',\n",
    " 'KIOST-ESM': 'GFDL',\n",
    " 'KACE-1-0-G': 'HadAM',\n",
    " 'IPSL-CM6A-LR': 'IPSL',\n",
    " 'INM-CM5-0': 'INM',\n",
    " 'INM-CM4-8': 'INM',\n",
    " 'HadGEM3-GC31-MM': 'HadAM',\n",
    " 'HadGEM3-GC31-LL': 'HadAM',\n",
    " 'GFDL-ESM4': 'GFDL',\n",
    " 'GFDL-CM4_gr2': 'GFDL',\n",
    " 'GFDL-CM4': 'GFDL',\n",
    " 'FGOALS-g3': 'CCM',\n",
    " 'EC-Earth3-Veg-LR': 'ECMWF',\n",
    " 'EC-Earth3': 'ECMWF',\n",
    " 'CanESM5': 'CanAM',\n",
    " 'CNRM-ESM2-1': 'ECMWF',\n",
    " 'CNRM-CM6-1': 'ECMWF',\n",
    " 'CMCC-ESM2': 'CCM',\n",
    " 'CMCC-CM2-SR5': 'CCM',\n",
    " 'BCC-CSM2-MR': 'CCM',\n",
    " 'ACCESS-ESM1-5': 'HadAM',\n",
    " 'ACCESS-CM2': 'HadAM',\n",
    " 'TaiESM1': 'CCM',\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "66e61532",
   "metadata": {},
   "outputs": [],
   "source": [
    "VARIABLES = {\n",
    "    'tasmax': {\n",
    "        'era_varname': 'maximum_2m_air_temperature',\n",
    "        'nex_transform': lambda x: x - 273.5,\n",
    "        'era_transform': lambda x: x - 273.5\n",
    "    },\n",
    "#    'tasmin': {\n",
    "#        'era_varname': 'minimum_2m_air_temperature',\n",
    "#        'nex_transform': lambda x: x - 273.5,\n",
    "#        'era_transform': lambda x: x - 273.5\n",
    "#    },\n",
    "#    'pr': {\n",
    "#        'era_varname': 'total_precipitation',\n",
    "#        'nex_transform': lambda x: x * 86400,\n",
    "#        'era_transform': lambda x: x * 1000\n",
    "#    },\n",
    "    'hurs': {\n",
    "        'era_varname': 'total_precipitation',\n",
    "        'nex_transform': lambda x: x,\n",
    "        'era_transform': lambda x: x\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8b53602c",
   "metadata": {},
   "outputs": [],
   "source": [
    "CAT_URL = \"https://storage.googleapis.com/cmip6/pangeo-cmip6.json\"\n",
    "COLL = intake.open_esm_datastore(CAT_URL)\n",
    "\n",
    "class Dataset:\n",
    "    def __init__(self, varname, model, scenario):\n",
    "        self.varname =varname\n",
    "        self.model = model\n",
    "        self.scenario = scenario\n",
    "        \n",
    "        print('Extracting {0} {1} {2}'.format(varname, scenario, model))\n",
    "        \n",
    "        run = MODELRUN\n",
    "        grid = MODELGRID[model]\n",
    "        \n",
    "        cat = COLL.search(experiment_id=scenario, table_id='day', variable_id=varname, source_id=model, member_id=run, grid_label=grid)\n",
    "        ds_dict = cat.to_dataset_dict(zarr_kwargs={'consolidated': True}, progressbar=False)\n",
    "        ds = list(ds_dict.values())[0][varname][0][0].sel(time=slice(['{0}-07-01'.format(min([i[0]-1 for i in FUTURE_YEARS])), '{0}-01-01'.format(HIST_START)][int(scenario=='historical')], ['{0}-12-31'.format(max([i[1] for i in FUTURE_YEARS])), '{0}-12-31'.format(HIST_END)][int(scenario=='historical')])).sel(lat=[i[0] for i in CITYLATLON.values()], lon=[[i[1], i[1]+360][int(i[1]<0)] for i in CITYLATLON.values()], method='nearest')\n",
    "        #ds = ds.select(lat=[i[0] for i in CITYLATLON.values()], lon=[[i[1], i[1]+360][int(i[1]<0)] for i in citylatlon.values()], method='nearest')\n",
    "        self.data = ds.to_numpy()\n",
    "\n",
    "        self.times = [str(d)[:10] for d in ds.time.data]\n",
    "        \n",
    "        del(ds)\n",
    "        \n",
    "    def get_timeseries(self, dates, loc_id):\n",
    "        timestart_idx = self.times.index(dates[0])\n",
    "        timeend_idx = self.times.index(dates[1]) + 1\n",
    "        return self.data[timestart_idx:timeend_idx, loc_id, loc_id]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b9af88e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calendardate_percentiles(nex_varname, q, latlon, sh_hem=False):\n",
    "    era_varname = VARIABLES[nex_varname]['era_varname']\n",
    "    hist_start = PERCENTILE_STARTYEAR\n",
    "    hist_end = PERCENTILE_ENDYEAR\n",
    "    allyears = []\n",
    "    for year in range(PERCENTILE_STARTYEAR, PERCENTILE_ENDYEAR):\n",
    "        allyears.append(VARIABLES[nex_varname]['era_transform'](get_eravar(era_varname, latlon, start_year=year, end_year=year, southern_hem=False)))\n",
    "    if not sh_hem:\n",
    "        return np.percentile(np.vstack(allyears), q, axis=0)\n",
    "    else:\n",
    "        res = np.percentile(np.vstack(allyears), q, axis=0)\n",
    "        return np.concatenate([res[152:], res[:152]])\n",
    "\n",
    "def wholeyear_percentile(nex_varname, q, latlon):\n",
    "    era_varname = VARIABLES[nex_varname]['era_varname']\n",
    "    hist_start = PERCENTILE_STARTYEAR\n",
    "    hist_end = PERCENTILE_ENDYEAR\n",
    "    allyears = []\n",
    "    for year in range(hist_start, hist_end):\n",
    "        allyears.append(VARIABLES[nex_varname]['era_transform'](get_eravar(era_varname, latlon, start_year=year, end_year=year, southern_hem=False)))\n",
    "    return np.percentile(np.concatenate(allyears).flatten(), q)\n",
    "\n",
    "def yearextreme_percentile(nex_varname, q, latlon, wantmax):\n",
    "    era_varname = VARIABLES[nex_varname]['era_varname']\n",
    "    hist_start = PERCENTILE_STARTYEAR\n",
    "    hist_end = PERCENTILE_ENDYEAR\n",
    "    allyears = []\n",
    "    for year in range(hist_start, hist_end):\n",
    "        allyears.append([np.min, np.max][int(wantmax)](VARIABLES[nex_varname]['era_transform'](get_eravar(era_varname, latlon, start_year=year, end_year=year, southern_hem=False))))\n",
    "    return np.percentile(np.array(allyears), q)\n",
    "\n",
    "def thresholdexceedance_mediancount(nex_varname, threshold, latlon, want_gte):\n",
    "    era_varname = VARIABLES[nex_varname]['era_varname']\n",
    "    data = VARIABLES[nex_varname]['era_transform'](get_eravar(era_varname, latlon, start_year=PERCENTILE_STARTYEAR, end_year=PERCENTILE_ENDYEAR, southern_hem=False))\n",
    "    if data.size % 365 != 0:\n",
    "        raise Exception('Data array length is not an integer multiple of 365')\n",
    "    byyear = data.reshape(data.size//365, 365)\n",
    "    if want_gte:\n",
    "        return np.median(np.sum(byyear >= threshold, axis=1))\n",
    "    else:\n",
    "        return np.median(np.sum(byyear <= threshold, axis=1))\n",
    "\n",
    "\n",
    "def get_rmsd(d1, d2):\n",
    "    c1 = seasonal_means(d1)\n",
    "    c2 = seasonal_means(d2)\n",
    "    return np.sqrt(np.mean(np.sum((c1 - c2)**2)))\n",
    "\n",
    "def count_runs(tf_array, min_runsize):\n",
    "    falses = np.zeros(tf_array.shape[0]).reshape((tf_array.shape[0],1))\n",
    "    extended_a = np.concatenate([[0], tf_array, [0]])\n",
    "    df = np.diff(extended_a)\n",
    "    starts = np.nonzero(df == 1)[0]\n",
    "    ends = np.nonzero(df == -1)[0]\n",
    "    count = 0\n",
    "    for idx in range(starts.size):\n",
    "        if ends[idx] - starts[idx] >= min_runsize:\n",
    "            count += 1\n",
    "    return count\n",
    "\n",
    "def removeLeapDays(arr, start_year, end_year, extralong=False, southern_hem=False):\n",
    "    if extralong:\n",
    "        indices = list(range(184))\n",
    "        jan1_idx = 184\n",
    "\n",
    "        for year in range(start_year, end_year+1):\n",
    "            indices += [jan1_idx + i for i in range(365)]\n",
    "            jan1_idx += 365\n",
    "            if calendar.isleap(year):\n",
    "                jan1_idx += 1\n",
    "        return arr[indices]\n",
    "    elif not southern_hem:\n",
    "        indices = []\n",
    "        jan1_idx = 0\n",
    "        for year in range(start_year, end_year+1):\n",
    "            indices += [jan1_idx + i for i in range(365)]\n",
    "            jan1_idx += 365\n",
    "            if calendar.isleap(year):\n",
    "                jan1_idx += 1\n",
    "        return arr[indices]\n",
    "    else:\n",
    "        indices = []\n",
    "        jul1_idx = 0\n",
    "        for year in range(start_year-1, end_year):\n",
    "            indices += [jul1_idx + i for i in range(365)]\n",
    "            jul1_idx += 365\n",
    "            if calendar.isleap(year):\n",
    "                jul1_idx += 1\n",
    "        return arr[indices]\n",
    "\n",
    "def get_eravar(varname, latlon, start_year, end_year, southern_hem=False, extralong=False):\n",
    "    model = 'ERA5'\n",
    "    dataset = ee.ImageCollection(\"ECMWF/ERA5/DAILY\")\n",
    "    gee_geom = ee.Geometry.Point((latlon[1], latlon[0]))\n",
    "    if extralong:\n",
    "        data_vars = dataset.select(varname).filter(ee.Filter.date('{0}-07-01'.format(start_year-1), '{0}-01-01'.format(end_year+1)))\n",
    "    elif not southern_hem:\n",
    "        data_vars = dataset.select(varname).filter(ee.Filter.date('{0}-01-01'.format(start_year), '{0}-01-01'.format(end_year+1)))\n",
    "    else:\n",
    "        data_vars = dataset.select(varname).filter(ee.Filter.date('{0}-07-01'.format(start_year-1), '{0}-07-01'.format(end_year)))\n",
    "    result = [i[4] for i in data_vars.getRegion(gee_geom, 2500, 'epsg:4326').getInfo()[1:]]\n",
    "    return removeLeapDays(np.array(result), start_year, end_year, extralong=extralong, southern_hem=southern_hem)\n",
    "\n",
    "def get_var(varname, model, loc_id, start_year, end_year, southern_hem=False, extralong=False, scenario='ssp585'):\n",
    "    scenario = [scenario, 'historical'][int(start_year < 2015)]\n",
    "    dataset = datasets[(varname, model, scenario)]\n",
    "    if extralong:\n",
    "        dates = ('{0}-07-01'.format(start_year-1), '{0}-12-31'.format(end_year))\n",
    "    elif not southern_hem:\n",
    "        dates = ('{0}-01-01'.format(start_year), '{0}-12-31'.format(end_year))\n",
    "    else:\n",
    "        dates = ('{0}-07-01'.format(start_year-1), '{0}-06-30'.format(end_year))\n",
    "    \n",
    "    ds = dataset.get_timeseries(dates, loc_id)\n",
    "    if YEARLENGTH[model] == 366:\n",
    "        return removeLeapDays(np.array(ds), start_year, end_year, extralong=extralong, southern_hem=southern_hem)\n",
    "    else:\n",
    "        return np.array(ds)\n",
    "\n",
    "    \n",
    "def quarters(d, start_year, end_year, southern_hem=False):\n",
    "    q2 = []  # 60-151\n",
    "    q3 = []  # 152-243\n",
    "    q4 = []  # 244-334\n",
    "    q1 = []  # 335-59\n",
    "    if not southern_hem:\n",
    "        jan1_idx = 365\n",
    "        for year in range(start_year, end_year):\n",
    "            tmp = np.concatenate((d[jan1_idx - 365 : jan1_idx - 365 + 60], d[jan1_idx + 335 : jan1_idx + 365]), axis=0)\n",
    "            q1.append(tmp)\n",
    "            q2.append(d[jan1_idx + 60 : jan1_idx + 152])\n",
    "            q3.append(d[jan1_idx + 152 : jan1_idx + 244])\n",
    "            q4.append(d[jan1_idx + 244 : jan1_idx + 335])\n",
    "\n",
    "            jan1_idx += 365 + [0, 0][int(False and calendar.isleap(year))]\n",
    "        mam_res = np.vstack(q2)\n",
    "        jja_res = np.vstack(q3)\n",
    "        son_res = np.vstack(q4)\n",
    "        djf_res = np.vstack(q1)\n",
    "    else:\n",
    "        jul1_idx = 365\n",
    "        for year in range(start_year, end_year+1):\n",
    "            tmp = np.concatenate((d[jul1_idx - 365 : jul1_idx - 365 + 60], d[jul1_idx + 335 : jul1_idx + 365]), axis=0)\n",
    "            q3.append(tmp)\n",
    "            q4.append(d[jul1_idx + 60 : jul1_idx + 152])\n",
    "            q1.append(d[jul1_idx + 152 : jul1_idx + 244])\n",
    "            q2.append(d[jul1_idx + 244 : jul1_idx + 335])\n",
    "\n",
    "            jul1_idx += 365 + [0, 0][int(False and calendar.isleap(year))]\n",
    "        mam_res = np.vstack(q4)\n",
    "        jja_res = np.vstack(q1)\n",
    "        son_res = np.vstack(q2)\n",
    "        djf_res = np.vstack(q3)\n",
    "    return mam_res, jja_res, son_res, djf_res\n",
    "    \n",
    "def seasonal_means(d):\n",
    "    q = quarters(d, HIST_START, HIST_END)\n",
    "    return np.array([np.mean(q[0], axis=1), np.mean(q[1], axis=1), np.mean(q[2], axis=1), np.mean(q[3], axis=1)])\n",
    "\n",
    "def calibration_function(hist_obs, hist_mod):\n",
    "# Calibration functions are P-P plots of historical and modeled values\n",
    "\n",
    "    source = np.sort(hist_obs.flatten())\n",
    "    target= np.sort(hist_mod.flatten())\n",
    "   \n",
    "    if (np.max(source) == 0 and np.min(source) == 0):\n",
    "        return np.arange(0, target.size) / target.size\n",
    "    if (np.max(target) == 0 and np.min(target) == 0):\n",
    "        return np.arange(0, source.size) / source.size\n",
    "    new_indices = []\n",
    "\n",
    "    for target_idx, target_value in enumerate(target):\n",
    "        if target_idx < len(source):\n",
    "            source_value = source[target_idx]\n",
    "            if source_value > target[-1]:\n",
    "                new_indices.append(target.size - 1)\n",
    "            else:\n",
    "                new_indices.append(np.argmax(target >= source_value))\n",
    "    return np.array(new_indices) / source.size\n",
    "\n",
    "def calibrate_component(uncalibrated_data, calibration_fxn):\n",
    "    N = len(uncalibrated_data)\n",
    "    unsorted_uncalib = [(i, idx) for idx, i in enumerate(uncalibrated_data)]\n",
    "    sorted_uncalib = sorted(unsorted_uncalib)\n",
    "    result = [0] * N\n",
    "    for j in range(N):\n",
    "        X_j = j / (N + 1)\n",
    "        Y_jprime = calibration_fxn[math.floor(X_j * len(calibration_fxn))]\n",
    "        jprime = math.floor(Y_jprime * (N + 1))\n",
    "        result[sorted_uncalib[j][1]] = sorted_uncalib[min(len(sorted_uncalib)-1, jprime)][0]\n",
    "    \n",
    "    return result\n",
    "\n",
    "def calibrate(uncalibrated_data, calibration_fxn):\n",
    "    mam = []\n",
    "    jja = []\n",
    "    son = []\n",
    "    djf = []\n",
    "    mam_idx = []\n",
    "    jja_idx = []\n",
    "    son_idx = []\n",
    "    djf_idx = []\n",
    "    for idx, i in enumerate(uncalibrated_data):\n",
    "        if idx % 365 >= 60 and idx % 365 < 152:\n",
    "            mam.append(uncalibrated_data[idx])\n",
    "            mam_idx.append(idx)\n",
    "        elif idx % 365 >= 152 and idx % 365 < 244:\n",
    "            jja.append(uncalibrated_data[idx])\n",
    "            jja_idx.append(idx)\n",
    "        elif idx % 365 >= 244 and idx % 365 < 335:\n",
    "            son.append(uncalibrated_data[idx])\n",
    "            son_idx.append(idx)\n",
    "        else:\n",
    "            djf.append(uncalibrated_data[idx])\n",
    "            djf_idx.append(idx)\n",
    "    \n",
    "    mam_calib = calibrate_component(np.array(mam), calibration_fxn[0])\n",
    "    jja_calib = calibrate_component(np.array(jja), calibration_fxn[1])\n",
    "    son_calib = calibrate_component(np.array(son), calibration_fxn[2])\n",
    "    djf_calib = calibrate_component(np.array(djf), calibration_fxn[3])\n",
    "    \n",
    "    result = [0] * len(uncalibrated_data)\n",
    "    for i in range(len(mam_idx)):\n",
    "        result[mam_idx[i]] = mam_calib[i]\n",
    "    for i in range(len(jja_idx)):\n",
    "        result[jja_idx[i]] = jja_calib[i]\n",
    "    for i in range(len(son_idx)):\n",
    "        result[son_idx[i]] = son_calib[i]\n",
    "    for i in range(len(djf_idx)):\n",
    "        result[djf_idx[i]] = djf_calib[i]\n",
    "\n",
    "    return np.array(result)\n",
    "\n",
    "def get_gamma(count, size):\n",
    "    return np.random.gamma(shape = count + 0.5, size=size)\n",
    "def get_beta(count, num, size):\n",
    "    return np.random.beta(a = count + 0.5, b = num - count + 0.5, size=size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "05aa86a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Hazard:\n",
    "    pass\n",
    "\n",
    "class Tempwave_simple(Hazard):\n",
    "    def __init__(self, varname, min_duration, threshold, want_gte=True):\n",
    "        if type(threshold) == np.ndarray and threshold.size % 365 != 0:\n",
    "            raise Exception('Comparison array length is not an integer multiple of 365')\n",
    "        self.varname = varname\n",
    "        self.want_gte = want_gte\n",
    "        self.min_duration = min_duration\n",
    "        self.threshold = threshold  # May be scalar or 365-long array\n",
    "        self.probmodel = 'Poisson'\n",
    "    def count(self, datalist):\n",
    "        data = datalist[0]\n",
    "        if type(self.threshold) in (float, int, np.float64, np.int32):\n",
    "            threshold = self.threshold\n",
    "        else:   # type is np array\n",
    "            threshold = np.array([])\n",
    "            while threshold.size < data.size:\n",
    "                threshold = np.concatenate([threshold, self.threshold])\n",
    "        if self.want_gte:\n",
    "            tf_array = data >= threshold\n",
    "        else:\n",
    "            tf_array = data <= threshold\n",
    "        return count_runs(tf_array, self.min_duration)\n",
    "\n",
    "def wetbulbtemp(T, RH):\n",
    "# JA Knox et al. 2017. Two simple and accurate approximations for wet-bulb\n",
    "# temperature in moist conditions, with forecasting applications. Bull. Am.\n",
    "# Meteorol. Soc. 98(9): 1897-1906. doi:10.1175/BAMS-D-16-0246.1\n",
    "    T = T.astype(np.float64)\n",
    "    rh_percent = RH.astype(np.float64)\n",
    "    return T * np.arctan(0.151977 * np.sqrt(rh_percent + 8.313659)) + np.arctan(T + rh_percent) - np.arctan(rh_percent - 1.676331) + ((0.00391838 * ((rh_percent)**(3/2))) * np.arctan(0.023101 * rh_percent)) - 4.686035\n",
    "\n",
    "\n",
    "def wetbulbtemp_chen(T, RH):\n",
    "    return -4391976 + (0.0198197 * RH) + (0.526359 * T) + (0.00730271 * RH * T) + (0.00024315 * RH * RH) - (0.0000258101 * T * RH * RH)\n",
    "\n",
    "class WetbulbHeatwave(Hazard):\n",
    "    def __init__(self, wbgt_threshold, min_duration):\n",
    "        self.varname = 'tasmax+hurs'\n",
    "        self.min_duration = min_duration\n",
    "        self.wbgt_threshold = wbgt_threshold\n",
    "        self.probmodel = 'Poisson'\n",
    "    def count(self, datalist):\n",
    "        data_t = datalist[0]\n",
    "        data_h = datalist[1]\n",
    "        data = wetbulbtemp(data_t, data_h)\n",
    "        tf_array = data >= self.wbgt_threshold\n",
    "        return count_runs(tf_array, self.min_duration)\n",
    "\n",
    "class WetbulbDays(Hazard):\n",
    "    def __init__(self, wbgt_threshold):\n",
    "        self.varname = 'tasmax+hurs'\n",
    "        self.wbgt_threshold = wbgt_threshold\n",
    "        self.probmodel = 'binomial'\n",
    "    def count(self, datalist):\n",
    "        data_t = datalist[0]\n",
    "        data_h = datalist[1]\n",
    "        data = wetbulbtemp(data_t, data_h)\n",
    "        byyear = data.reshape(data.size//365, 365)\n",
    "        return np.sum((np.max(byyear, axis=1) >= self.wbgt_threshold) * 1)\n",
    "    \n",
    "    def count_nc(self, datalist, targetcount):\n",
    "        data_t = datalist[0]\n",
    "        data_h = datalist[1]\n",
    "        data = wetbulbtemp(data_t, data_h)\n",
    "        if data.size % 365 != 0:\n",
    "            raise Exception('Data array length is not an integer multiple of 365')\n",
    "        byyear = data.reshape(data.size//365, 365)\n",
    "        return np.sum(np.abs(np.sum(byyear >= self.wbgt_threshold, axis=1) - targetcount) < 0.5)\n",
    "    \n",
    "class Heatwave_highlow(Hazard):\n",
    "    def __init__(self, hightemp, lowtemp, min_duration):\n",
    "        self.varname = 'tasmax+tasmin'\n",
    "        self.min_duration = min_duration\n",
    "        self.hightemp = hightemp\n",
    "        self.lowtemp = lowtemp\n",
    "        self.probmodel = 'Poisson'\n",
    "    def count(self, datalist):\n",
    "        data_tx = datalist[0]\n",
    "        data_tn = datalist[1]\n",
    "        if type(self.hightemp) in (float, int, np.float64, np.int32):\n",
    "            high_threshold = self.hightemp\n",
    "        else:   # type is np array\n",
    "            high_threshold = np.array([])\n",
    "            while high_threshold.size < data_tx.size:\n",
    "                high_threshold = np.concatenate([high_threshold, self.hightemp])\n",
    "        if type(self.lowtemp) in (float, int, np.float64, np.int32):\n",
    "            low_threshold = self.lowtemp\n",
    "        else:   # type is np array\n",
    "            low_threshold = np.array([])\n",
    "            while low_threshold.size < data_tn.size:\n",
    "                low_threshold = np.concatenate([low_threshold, self.lowtemp])\n",
    "        tf_array_tx = data_tx >= high_threshold\n",
    "        tf_array_tn = data_tn >= low_threshold\n",
    "        return count_runs(tf_array_tx * tf_array_tn, self.min_duration)\n",
    "\n",
    "class Threshold_simple(Hazard):\n",
    "    def __init__(self, varname, var_threshold, want_gte):\n",
    "        self.varname = varname\n",
    "        self.var_threshold = var_threshold\n",
    "        self.want_gte = want_gte\n",
    "        self.probmodel = 'binomial'\n",
    "    def count(self, datalist):\n",
    "        data = datalist[0]\n",
    "        if data.size % 365 != 0:\n",
    "            raise Exception('Data array length is not an integer multiple of 365')\n",
    "        byyear = data.reshape(data.size//365, 365)\n",
    "        if self.want_gte:\n",
    "            return np.sum(np.sum(byyear >= self.var_threshold, axis=1) >= self.count_threshold)\n",
    "        else:\n",
    "            return np.sum(np.sum(byyear <= self.var_threshold, axis=1) >= self.count_threshold)\n",
    "    def count_nc(self, datalist, targetcount):\n",
    "        data = datalist[0]\n",
    "        if data.size % 365 != 0:\n",
    "            raise Exception('Data array length is not an integer multiple of 365')\n",
    "        byyear = data.reshape(data.size//365, 365)\n",
    "        return np.sum(np.abs(np.sum(byyear >= self.var_threshold, axis=1) - targetcount) < 0.5)\n",
    "        \n",
    "\n",
    "class Hotdays_inrange(Hazard):\n",
    "    def __init__(self, hightemp, lowtemp):\n",
    "        self.varname = 'tasmax'\n",
    "        self.hightemp = hightemp\n",
    "        self.lowtemp = lowtemp\n",
    "        self.probmodel = 'binomial'\n",
    "    def count(self, datalist):\n",
    "        data = datalist[0]\n",
    "        if data.size % 365 != 0:\n",
    "            raise Exception('Data array length is not an integer multiple of 365')\n",
    "        tf_array_high = data <= self.hightemp\n",
    "        tf_array_low = data >= self.lowtemp\n",
    "        return runs(tf_array_high * tf_array_low, self.min_duration, 'count')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ae03622d",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Location:\n",
    "    def __init__(self, name, loc_id, latlon):\n",
    "        self.name = name\n",
    "        self.loc_id = loc_id\n",
    "        self.latlon = latlon\n",
    "        self.hist_observed = {}\n",
    "        self.hist_modeled = {scenario: {} for scenario in FUTURE_SCENARIOS}\n",
    "        self.best_models = {scenario: {} for scenario in FUTURE_SCENARIOS}\n",
    "        self.fut_modeled = {scenario: {} for scenario in FUTURE_SCENARIOS}\n",
    "        self.calib_fxns = {scenario: {} for scenario in FUTURE_SCENARIOS}\n",
    "        \n",
    "    def get_data(self, varname, future_start, future_end):\n",
    "    # Gets historical observations from ERA5 Daily Aggregate (from GEE)\n",
    "    # Goes through all NEX-GDDP-CMIP6 models in GEE and gets historical model outputs\n",
    "    # Chooses three best models based on quarterly RMSD\n",
    "\n",
    "        def relhum(T, Tdp):\n",
    "            T = T.astype(np.float64)\n",
    "            Tdp = Tdp.astype(np.float64)\n",
    "            numerator = np.exp(17.625 * Tdp / (243.04 + Tdp))\n",
    "            denominator = np.exp(17.625 * T / (243.04 + T))\n",
    "            return 100 * numerator / denominator\n",
    "        \n",
    "        #print('  Getting historical {0}'.format(varname))\n",
    "        if varname=='hurs':\n",
    "            era_dewpoint = VARIABLES['tasmax']['era_transform'](get_eravar('dewpoint_2m_temperature', self.latlon, HIST_START, HIST_END, southern_hem=False))\n",
    "            era_maxtemp = VARIABLES['tasmax']['era_transform'](get_eravar(VARIABLES['tasmax']['era_varname'], self.latlon, HIST_START, HIST_END, southern_hem=False))\n",
    "            hist_obs = relhum(era_maxtemp, era_dewpoint)\n",
    "        else:\n",
    "            hist_obs = VARIABLES[varname]['era_transform'](get_eravar(VARIABLES[varname]['era_varname'], self.latlon, HIST_START, HIST_END, southern_hem=False))\n",
    "        self.hist_observed[varname] = hist_obs\n",
    "        for scenario in FUTURE_SCENARIOS:\n",
    "            hist_mods = {}\n",
    "            rmsds = []\n",
    "            for model in MODELS[varname]:\n",
    "                #print('    Getting {0}'.format(model))\n",
    "                hist_mod = VARIABLES[varname]['nex_transform'](get_var(varname, model, self.loc_id, HIST_START, HIST_END, southern_hem=False, scenario=scenario))\n",
    "                hist_mods[model] = hist_mod\n",
    "                rmsds.append((get_rmsd(hist_obs, hist_mod), model))\n",
    "            rmsds.sort()\n",
    "\n",
    "            best_models = []\n",
    "            families = []\n",
    "            idx = 0\n",
    "            while len(best_models) < NUM_BEST_MODELS:\n",
    "                if not MODEL_FAMILY[rmsds[idx][1]] in families:\n",
    "                    best_models.append(rmsds[idx][1])\n",
    "                    families.append(MODEL_FAMILY[rmsds[idx][1]])\n",
    "                idx += 1\n",
    "\n",
    "            #for m in best_models:\n",
    "            #    print(m, [i[0] for i in rmsds if i[1]==m][0])\n",
    "            #best_models = []\n",
    "            #for idx in range(min(NUM_BEST_MODELS, len(MODELS[varname]))):\n",
    "            #    best_models.append(rmsds[idx][1])\n",
    "\n",
    "            self.hist_modeled[scenario][varname] = hist_mods\n",
    "            self.best_models[scenario][varname] = best_models\n",
    "\n",
    "        # Get calibration functions\n",
    "            #print('  Getting calibration functions')\n",
    "            self.calib_fxns[scenario][varname] = {}\n",
    "            hist_obs = self.hist_observed[varname]\n",
    "            hist_mod = self.hist_modeled[scenario][varname]\n",
    "            for model in self.best_models[scenario][varname]:\n",
    "                o_quarters = quarters(hist_obs, HIST_START, HIST_END)\n",
    "                m_quarters = quarters(hist_mod[model], HIST_START, HIST_END)\n",
    "                self.calib_fxns[scenario][varname][model] = [calibration_function(o_quarters[i].flatten(), m_quarters[i].flatten()) for i in range(4)]\n",
    "\n",
    "        # Get future model outputs\n",
    "            #print('  Getting future {0}'.format(varname))\n",
    "            self.fut_modeled[scenario][varname] = {\n",
    "                model: VARIABLES[varname]['nex_transform'](get_var(varname, model, self.loc_id, future_start, future_end, extralong=True, scenario=scenario)) for model in best_models\n",
    "            }\n",
    "        \n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "41ad76b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Estimate:\n",
    "    def __init__(self, location, hazard, future_start, future_end, sh_year):\n",
    "        self.location = location\n",
    "        self.hazard = hazard\n",
    "        self.future_start = future_start\n",
    "        self.future_end = future_end\n",
    "        self.sh_year = sh_year \n",
    "        self.estimate = {scenario: {} for scenario in FUTURE_SCENARIOS}\n",
    "        self.expected_estimate = {scenario: {} for scenario in FUTURE_SCENARIOS}\n",
    "        \n",
    "        for varname in self.hazard.varname.split('+'):\n",
    "            if not varname in list(self.location.hist_observed.keys()):\n",
    "                self.location.get_data(varname, future_start, future_end)\n",
    "        \n",
    "    def future_count(self):\n",
    "    # Gets future model outputs for three best models\n",
    "    # Calibrates based on stored model-specific quarterly calibration functions\n",
    "    # Calculates event count within future time series\n",
    "    # Draws 10,000 posterior rate parameters from appropriate Jeffrys prior distribution (parameterized with count)\n",
    "    # For each rate parameter, draw one event count\n",
    "        for scenario in FUTURE_SCENARIOS:\n",
    "            fut_mod = {}\n",
    "            varnames = self.hazard.varname.split('+')\n",
    "            for varname in varnames:\n",
    "                fut_mod[varname] = {}\n",
    "                for model in self.location.best_models[scenario][varname]:\n",
    "                    if self.sh_year:\n",
    "                        fut_mod[varname][model] = self.location.fut_modeled[scenario][varname][model][:-184]\n",
    "                    else:\n",
    "                        fut_mod[varname][model] = self.location.fut_modeled[scenario][varname][model][184:]\n",
    "            best_models = []\n",
    "            for idx in range(min(len(MODELS[varname]), NUM_BEST_MODELS)):\n",
    "                best_models.append('+'.join([self.location.best_models[scenario][varname][idx] for varname in varnames]))\n",
    "            posterior_rateparams = {}\n",
    "            posterior_draws = {}\n",
    "            estimate = {}\n",
    "            for modelplus in best_models:\n",
    "                calib_data = []\n",
    "                for idx, varname in enumerate(varnames):\n",
    "                    model = modelplus.split('+')[idx]\n",
    "                    calib_data.append(np.array(calibrate(fut_mod[varname][model], self.location.calib_fxns[scenario][varname][model])))\n",
    "                if self.sh_year:\n",
    "                    count = self.hazard.count(calib_data)\n",
    "                else:\n",
    "                    count = self.hazard.count([cd[152:-213] for cd in calib_data])\n",
    "                if self.hazard.probmodel == 'Poisson':\n",
    "                    posterior_rateparams[modelplus] = get_gamma(count, 10000)\n",
    "                    posterior_draws[modelplus] = np.random.poisson(posterior_rateparams[modelplus], 10000)\n",
    "                else:  # self.hazard.probmodel == 'binomial'\n",
    "                    posterior_rateparams[modelplus] = get_beta(count, self.future_end - self.future_start + 1, 10000)\n",
    "                    posterior_draws[modelplus] = np.random.binomial(self.future_end - self.future_start + 1, posterior_rateparams[modelplus], 10000)\n",
    "                self.estimate[scenario][modelplus] = [np.percentile(posterior_draws[modelplus], q)/(self.future_end - self.future_start + 1) for q in (25, 50, 75)] #np.sum(posterior_draws[modelplus]) / 10000\n",
    "\n",
    "            print()\n",
    "            print(scenario)\n",
    "            for modelplus in best_models:\n",
    "                if self.hazard.probmodel == 'Poisson':\n",
    "                    print('{0}: {1}'.format(modelplus, str([i for i in self.estimate[scenario][modelplus]])[1:-1]))\n",
    "                else:\n",
    "                    print('{0}: {1}'.format(modelplus, str(['{0:.1f}%'.format(i / (self.future_end - self.future_start + 1) * 100) for i in self.estimate[scenario][modelplus]])[1:-1]))\n",
    "    \n",
    "    def expected_count(self):\n",
    "        for scenario in FUTURE_SCENARIOS:\n",
    "            fut_mod = {}\n",
    "            varnames = self.hazard.varname.split('+')\n",
    "            for varname in varnames:\n",
    "                fut_mod[varname] = {}\n",
    "                for model in self.location.best_models[scenario][varname]:\n",
    "                    if self.sh_year:\n",
    "                        fut_mod[varname][model] = self.location.fut_modeled[scenario][varname][model][:-184]\n",
    "                    else:\n",
    "                        fut_mod[varname][model] = self.location.fut_modeled[scenario][varname][model][184:]\n",
    "            best_models = []\n",
    "            for idx in range(min(len(MODELS[varname]), NUM_BEST_MODELS)):\n",
    "                best_models.append('+'.join([self.location.best_models[scenario][varname][idx] for varname in varnames]))\n",
    "            estimates = {}\n",
    "            for modelplus in best_models:\n",
    "                res_sum = [0, 0, 0]\n",
    "                for targetcount in range(0, 365):\n",
    "                    \n",
    "                    calib_data = []\n",
    "                    for idx, varname in enumerate(varnames):\n",
    "                        model = modelplus.split('+')[idx]\n",
    "                        calib_data.append(np.array(calibrate(fut_mod[varname][model], self.location.calib_fxns[scenario][varname][model])))\n",
    "                    if self.sh_year:\n",
    "                        count = self.hazard.count_nc(calib_data, targetcount)\n",
    "                    else:\n",
    "                        count = self.hazard.count_nc([cd[152:-213] for cd in calib_data], targetcount)\n",
    "\n",
    "                    posterior_rateparams = get_beta(count, self.future_end - self.future_start + 1, 10000)\n",
    "                    posterior_draws = np.random.binomial(self.future_end - self.future_start + 1, posterior_rateparams, 10000)\n",
    "                    \n",
    "                    probs = [np.percentile(posterior_draws, q) / (self.future_end - self.future_start + 1) for q in (25, 50, 75)]\n",
    "                    res_sum = [res_sum[idx] + (targetcount * probs[idx]) for idx in (0, 1, 2)]\n",
    "                self.expected_estimate[scenario][modelplus] = [max(min(res_sum[idx], 365), 0) for idx in (0, 1, 2)]\n",
    "\n",
    "                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "b78f2d07",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Extracting tasmax historical GFDL-ESM4\n",
      "Extracting tasmax ssp119 GFDL-ESM4\n",
      "Extracting tasmax ssp126 GFDL-ESM4\n",
      "Extracting tasmax historical CanESM5\n",
      "Extracting tasmax ssp119 CanESM5\n",
      "Extracting tasmax ssp126 CanESM5\n",
      "Extracting tasmax historical MRI-ESM2-0\n",
      "Extracting tasmax ssp119 MRI-ESM2-0\n",
      "Extracting tasmax ssp126 MRI-ESM2-0\n",
      "Extracting tasmax historical IPSL-CM6A-LR\n",
      "Extracting tasmax ssp119 IPSL-CM6A-LR\n",
      "Extracting tasmax ssp126 IPSL-CM6A-LR\n",
      "Extracting tasmax historical EC-Earth3-Veg-LR\n",
      "Extracting tasmax ssp119 EC-Earth3-Veg-LR\n",
      "Extracting tasmax ssp126 EC-Earth3-Veg-LR\n",
      "Extracting hurs historical GFDL-ESM4\n",
      "Extracting hurs ssp119 GFDL-ESM4\n",
      "Extracting hurs ssp126 GFDL-ESM4\n",
      "Extracting hurs historical CanESM5\n",
      "Extracting hurs ssp119 CanESM5\n",
      "Extracting hurs ssp126 CanESM5\n",
      "Extracting hurs historical MRI-ESM2-0\n",
      "Extracting hurs ssp119 MRI-ESM2-0\n",
      "Extracting hurs ssp126 MRI-ESM2-0\n",
      "Extracting hurs historical IPSL-CM6A-LR\n",
      "Extracting hurs ssp119 IPSL-CM6A-LR\n",
      "Extracting hurs ssp126 IPSL-CM6A-LR\n",
      "Extracting hurs historical EC-Earth3-Veg-LR\n",
      "Extracting hurs ssp119 EC-Earth3-Veg-LR\n",
      "Extracting hurs ssp126 EC-Earth3-Veg-LR\n",
      "CPU times: total: 19min 7s\n",
      "Wall time: 1h 41min 41s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "datasets = {}\n",
    "for varname in VARIABLES:\n",
    "    for model in MODELS[varname]:\n",
    "        for scenario in ['historical'] + FUTURE_SCENARIOS:\n",
    "            datasets[(varname, model, scenario)] = Dataset(varname, model, scenario)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "44bf2b40",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "120\n",
      "Houston US\n",
      "\n",
      "Days wet-bulb > 30\n",
      "2080 2099\n",
      "ssp119 MRI-ESM2-0+IPSL-CM6A-LR: 2.5500000000000003,27.650000000000006,365\n",
      "ssp119 CanESM5+MRI-ESM2-0: 1.6,13.900000000000002,365\n",
      "ssp119 EC-Earth3-Veg-LR+EC-Earth3-Veg-LR: 0.0,0.05,365\n",
      "ssp126 MRI-ESM2-0+IPSL-CM6A-LR: 0.7000000000000001,26.35,365\n",
      "ssp126 CanESM5+MRI-ESM2-0: 1.6500000000000001,23.100000000000005,365\n",
      "ssp126 EC-Earth3-Veg-LR+EC-Earth3-Veg-LR: 0.0,0.0,365\n",
      "\n",
      "123\n",
      "Phoenix US\n",
      "\n",
      "Days wet-bulb > 30\n",
      "2080 2099\n",
      "ssp119 CanESM5+CanESM5: 0.6,1.0499999999999998,365\n",
      "ssp119 MRI-ESM2-0+EC-Earth3-Veg-LR: 1.6,10.000000000000002,365\n",
      "ssp119 EC-Earth3-Veg-LR+MRI-ESM2-0: 0.75,4.6,365\n",
      "ssp126 CanESM5+CanESM5: 0.8500000000000001,2.35,365\n",
      "ssp126 MRI-ESM2-0+EC-Earth3-Veg-LR: 1.35,5.55,365\n",
      "ssp126 EC-Earth3-Veg-LR+MRI-ESM2-0: 0.7500000000000001,2.6,365\n",
      "CPU times: total: 3min 28s\n",
      "Wall time: 4min 11s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "with open('wetbulb_chen.csv', 'w') as ofile:\n",
    "    ofile.write('Location,Hazard,Scenario,Model,Year_range,Q1,Q2,Q3\\n')\n",
    "    for name in CITYLATLON:\n",
    "        if True:\n",
    "            lat, lon, loc_id = CITYLATLON[name]\n",
    "            print()\n",
    "            print(loc_id)\n",
    "            print(name)\n",
    "\n",
    "            #high_95c = calendardate_percentiles('tasmax', 95, (lat, lon), sh_hem=lat < 0)\n",
    "            #Tx_95y = wholeyear_percentile('tasmax', 95, (lat, lon))\n",
    "            #median_95p_freq = thresholdexceedance_mediancount('tasmax', Tx_95y, (lat, lon), True)\n",
    "            test_hazards = [\n",
    "                {'name': 'Days wet-bulb > 30', 'obj': WetbulbDays(30), 'sh_year': lat < 0},\n",
    "            ]\n",
    "           # print('  5th pctl mintemp = {0:.1f}'.format(cold_5n))\n",
    "            loc = Location(name, loc_id, (lat, lon))\n",
    "\n",
    "            for haz in test_hazards:\n",
    "\n",
    "                for fut_start, fut_end in FUTURE_YEARS:\n",
    "                    print()\n",
    "                    print(haz['name'])\n",
    "                    print(fut_start, fut_end)\n",
    "                    est = Estimate(loc, haz['obj'], fut_start, fut_end, haz['sh_year'])\n",
    "                    if haz['obj'].probmodel == 'Poisson':\n",
    "                        est.future_count()\n",
    "                        for scenario in FUTURE_SCENARIOS:\n",
    "                            for model in est.estimate[scenario]:\n",
    "                                res = list(est.estimate[scenario][model])\n",
    "                                ofile.write('{0},{1},{2},{3},{4},{5},{6},{7}\\n'.format(loc.name, haz['name'], scenario, model, '{0}-{1}'.format(fut_start, fut_end), res[0], res[1], res[2]))\n",
    "\n",
    "                    else:\n",
    "                        est.expected_count()\n",
    "                        for scenario in FUTURE_SCENARIOS:\n",
    "                            for model in est.expected_estimate[scenario]:\n",
    "                                res = list(est.expected_estimate[scenario][model])\n",
    "                                print('{0} {1}: {2}'.format(scenario, model, ','.join([str(i) for i in res])))\n",
    "                                ofile.write('{0},{1},{2},{3},{4},{5},{6},{7}\\n'.format(loc.name, haz['name'], scenario, model, '{0}-{1}'.format(fut_start, fut_end), res[0], res[1], res[2]))\n",
    "            del(loc)\n",
    "            del(est)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35ff563c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19c815ed",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c42df46",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "c7dd43ed",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GFDL-CM4: min modeled value does not exceed observed 10th percentile  True\n",
      "GFDL-CM4: max modeled value does not exceed observed 90th percentile  True\n",
      "CanESM5: min modeled value does not exceed observed 10th percentile  True\n",
      "CanESM5: max modeled value does not exceed observed 90th percentile  True\n",
      "ACCESS-CM2: min modeled value does not exceed observed 10th percentile  True\n",
      "ACCESS-CM2: max modeled value does not exceed observed 90th percentile  True\n",
      "GFDL-CM4: min modeled value does not exceed observed 10th percentile  True\n",
      "GFDL-CM4: max modeled value does not exceed observed 90th percentile  True\n",
      "CanESM5: min modeled value does not exceed observed 10th percentile  True\n",
      "CanESM5: max modeled value does not exceed observed 90th percentile  True\n",
      "ACCESS-CM2: min modeled value does not exceed observed 10th percentile  True\n",
      "ACCESS-CM2: max modeled value does not exceed observed 90th percentile  True\n",
      "GFDL-CM4: min modeled value does not exceed observed 10th percentile  True\n",
      "GFDL-CM4: max modeled value does not exceed observed 90th percentile  True\n",
      "CanESM5: min modeled value does not exceed observed 10th percentile  True\n",
      "CanESM5: max modeled value does not exceed observed 90th percentile  True\n",
      "ACCESS-CM2: min modeled value does not exceed observed 10th percentile  True\n",
      "ACCESS-CM2: max modeled value does not exceed observed 90th percentile  True\n",
      "GFDL-CM4: min modeled value does not exceed observed 10th percentile  True\n",
      "GFDL-CM4: max modeled value does not exceed observed 90th percentile  True\n",
      "CanESM5: min modeled value does not exceed observed 10th percentile  False\n",
      "CanESM5: max modeled value does not exceed observed 90th percentile  True\n",
      "ACCESS-CM2: min modeled value does not exceed observed 10th percentile  False\n",
      "ACCESS-CM2: max modeled value does not exceed observed 90th percentile  True\n"
     ]
    }
   ],
   "source": [
    "for quarter in range(4):\n",
    "    obs_10 = np.percentile(quarters(hist_obs_tx, HIST_START, HIST_END)[quarter], 10)\n",
    "    obs_90 = np.percentile(quarters(hist_obs_tx, HIST_START, HIST_END)[quarter], 90)\n",
    "    for model in best_models_tx:\n",
    "        mod = quarters(hist_mods_tx[model] - 273.15, HIST_START, HIST_END)[quarter].flatten()\n",
    "        print('{0}: min modeled value does not exceed observed 10th percentile  {1}'.format(model, min(mod) <= obs_10))\n",
    "        print('{0}: max modeled value does not exceed observed 90th percentile  {1}'.format(model, max(mod) >= obs_90))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2847140a",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}