bias-vs-trend_ensemble_get-data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pre-proj: check bias versus trends in multimodel ensemble\n",
    "## Get and store data to common grid 1°x1° from topographic file\n",
    "conda env: `phd_v3` (in `envs/phd`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n",
      "3.8.5 | packaged by conda-forge | (default, Jul 24 2020, 01:25:15) \n",
      "[GCC 7.5.0]\n"
     ]
    }
   ],
   "source": [
    "# To reload external files automatically (ex: utils)\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import xarray as xr\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import proplot as plot # New plot library (https://proplot.readthedocs.io/en/latest/)\n",
    "plot.rc['savefig.dpi'] = 300 # 1200 is too big! #https://proplot.readthedocs.io/en/latest/basics.html#Creating-figures\n",
    "from scipy import stats\n",
    "import xesmf as xe # For regridding (https://xesmf.readthedocs.io/en/latest/)\n",
    "\n",
    "import sys\n",
    "sys.path.insert(1, '/home/mlalande/notebooks/utils') # to include my util file in previous directory\n",
    "import utils as u # my personal functions\n",
    "u.check_python_version()\n",
    "# u.check_virtual_memory()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "period_past = slice('1979','2014')\n",
    "seasons = ['Annual', 'DJFMA', 'JJAS']\n",
    "# seasons = ['Annual']\n",
    "\n",
    "# Make a extended version for regridding properly on the edges\n",
    "latlim, lonlim = u.get_domain_HMA()\n",
    "latlim_ext, lonlim_ext = slice(latlim.start-5, latlim.stop+5), slice(lonlim.start-5, lonlim.stop+5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load topography"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = xr.open_dataset('GMTED2010_15n240_1000deg.nc').drop_dims('nbounds').swap_dims(\n",
    "    {'nlat': 'latitude', 'nlon': 'longitude'}).drop({'nlat', 'nlon'}).rename(\n",
    "    {'latitude': 'lat', 'longitude': 'lon'}).sel(lat=latlim_ext, lon=lonlim_ext)\n",
    "elevation = ds.elevation\n",
    "elevation_std = ds.elevation_stddev"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['BCC-CSM2-MR',\n",
       " 'BCC-ESM1',\n",
       " 'CAS-ESM2-0',\n",
       " 'CESM2',\n",
       " 'CESM2-FV2',\n",
       " 'CESM2-WACCM',\n",
       " 'CESM2-WACCM-FV2',\n",
       " 'CNRM-CM6-1',\n",
       " 'CNRM-CM6-1-HR',\n",
       " 'CNRM-ESM2-1',\n",
       " 'CanESM5',\n",
       " 'CanESM5-CanOE',\n",
       " 'GFDL-CM4',\n",
       " 'GISS-E2-1-G',\n",
       " 'GISS-E2-1-H',\n",
       " 'HadGEM3-GC31-LL',\n",
       " 'HadGEM3-GC31-MM',\n",
       " 'IPSL-CM6A-LR',\n",
       " 'MIROC-ES2L',\n",
       " 'MIROC6',\n",
       " 'MPI-ESM1-2-HR',\n",
       " 'MPI-ESM1-2-LR',\n",
       " 'MRI-ESM2-0',\n",
       " 'NorESM2-LM',\n",
       " 'SAM0-UNICON',\n",
       " 'TaiESM1',\n",
       " 'UKESM1-0-LL']"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list_models = u.get_model_names()\n",
    "list_models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## tas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Near-Surface Air Temperature'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var = 'tas'\n",
    "table = u.get_table(var)\n",
    "label, units, levels, cmap, extend, \\\n",
    "    levels_diff, cmap_diff, extend_diff, \\\n",
    "    levels_bias, cmap_bias, extend_bias = u.get_var_infos(var)\n",
    "label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IPSL-CM6A-LR (1/10)\n",
      "['r10i1p1f1', 'r11i1p1f1', 'r12i1p1f1', 'r13i1p1f1', 'r14i1p1f1', 'r15i1p1f1', 'r16i1p1f1', 'r17i1p1f1', 'r18i1p1f1', 'r19i1p1f1', 'r1i1p1f1', 'r20i1p1f1', 'r21i1p1f1', 'r22i1p1f1', 'r23i1p1f1', 'r24i1p1f1', 'r25i1p1f1', 'r26i1p1f1', 'r27i1p1f1', 'r28i1p1f1', 'r29i1p1f1', 'r2i1p1f1', 'r30i1p1f1', 'r31i1p1f1', 'r32i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "32\n",
      "Reuse existing file: bilinear_28x25_35x60.nc\n",
      "\n",
      "\n",
      "MIROC-ES2L (2/10)\n",
      "['r10i1p1f2', 'r1i1p1f2', 'r2i1p1f2', 'r3i1p1f2', 'r4i1p1f2', 'r5i1p1f2', 'r6i1p1f2', 'r7i1p1f2', 'r8i1p1f2', 'r9i1p1f2']\n",
      "10\n",
      "Reuse existing file: bilinear_13x21_35x60.nc\n",
      "\n",
      "\n",
      "MIROC6 (3/10)\n",
      "['r10i1p1f1', 'r11i1p1f1', 'r12i1p1f1', 'r13i1p1f1', 'r14i1p1f1', 'r15i1p1f1', 'r16i1p1f1', 'r17i1p1f1', 'r18i1p1f1', 'r19i1p1f1', 'r1i1p1f1', 'r20i1p1f1', 'r21i1p1f1', 'r22i1p1f1', 'r23i1p1f1', 'r24i1p1f1', 'r25i1p1f1', 'r26i1p1f1', 'r27i1p1f1', 'r28i1p1f1', 'r29i1p1f1', 'r2i1p1f1', 'r30i1p1f1', 'r31i1p1f1', 'r32i1p1f1', 'r33i1p1f1', 'r34i1p1f1', 'r35i1p1f1', 'r36i1p1f1', 'r37i1p1f1', 'r38i1p1f1', 'r39i1p1f1', 'r3i1p1f1', 'r40i1p1f1', 'r41i1p1f1', 'r42i1p1f1', 'r43i1p1f1', 'r44i1p1f1', 'r45i1p1f1', 'r46i1p1f1', 'r47i1p1f1', 'r48i1p1f1', 'r49i1p1f1', 'r4i1p1f1', 'r50i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "50\n",
      "Reuse existing file: bilinear_25x42_35x60.nc\n",
      "\n",
      "\n",
      "MPI-ESM1-2-HR (4/10)\n",
      "['r10i1p1f1', 'r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "10\n",
      "Reuse existing file: bilinear_37x64_35x60.nc\n",
      "\n",
      "\n",
      "MPI-ESM1-2-LR (5/10)\n",
      "['r10i1p1f1', 'r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "10\n",
      "Reuse existing file: bilinear_19x32_35x60.nc\n",
      "\n",
      "\n",
      "MRI-ESM2-0 (6/10)\n",
      "['r1i1000p1f1', 'r1i1p1f1', 'r1i2p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1']\n",
      "7\n",
      "Reuse existing file: bilinear_32x54_35x60.nc\n",
      "\n",
      "\n",
      "NorESM2-LM (7/10)\n",
      "['r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1']\n",
      "3\n",
      "Reuse existing file: bilinear_18x25_35x60.nc\n",
      "\n",
      "\n",
      "SAM0-UNICON (8/10)\n",
      "['r1i1p1f1']\n",
      "1\n",
      "Reuse existing file: bilinear_37x49_35x60.nc\n",
      "\n",
      "\n",
      "TaiESM1 (9/10)\n",
      "['r1i1p1f1']\n",
      "1\n",
      "Reuse existing file: bilinear_37x49_35x60.nc\n",
      "\n",
      "\n",
      "UKESM1-0-LL (10/10)\n",
      "['r10i1p1f2', 'r11i1p1f2', 'r12i1p1f2', 'r16i1p1f2', 'r17i1p1f2', 'r18i1p1f2', 'r19i1p1f2', 'r1i1p1f2', 'r2i1p1f2', 'r3i1p1f2', 'r4i1p1f2', 'r5i1p1f3', 'r6i1p1f3', 'r7i1p1f3', 'r8i1p1f2', 'r9i1p1f2']\n",
      "16\n",
      "Reuse existing file: bilinear_28x32_35x60.nc\n",
      "\n",
      "\n",
      "CPU times: user 1min 37s, sys: 12.8 s, total: 1min 50s\n",
      "Wall time: 6min 40s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "list_models = u.get_model_names()\n",
    "list_models = list_models[17:]\n",
    "# list_models = ['EC-Earth3']\n",
    "\n",
    "for i, model in enumerate(list_models):\n",
    "    print(model + ' ('+str(i+1)+'/'+str(len(list_models))+')')\n",
    "    institute, grid, realization, calendar = u.get_model_infos(model, var)\n",
    "    \n",
    "    #################\n",
    "    ### Load past ###\n",
    "    #################\n",
    "    \n",
    "    realizations = !ls {'/bdd/CMIP6/CMIP/'+institute+'/'+model+'/historical/'}\n",
    "    temp = []\n",
    "    no_file = []\n",
    "    \n",
    "    for i, realization in enumerate(realizations):\n",
    "        try:\n",
    "            da = xr.open_mfdataset(\n",
    "                '/bdd/CMIP6/CMIP/'+institute+'/'+model+'/historical/'+realization+'/'+table+'/'+var+'/'+grid+'/latest/*.nc', \n",
    "                combine='by_coords'\n",
    "            ).sel(time=period_past, lat=latlim_ext, lon=lonlim_ext)[var]\n",
    "            \n",
    "            temp.append(da)\n",
    "            \n",
    "        except OSError:\n",
    "            no_file.append(realization)\n",
    "            \n",
    "    for nf in no_file:\n",
    "        realizations.remove(nf)\n",
    "        \n",
    "    print(realizations)\n",
    "    \n",
    "    da_past = xr.concat(temp, pd.Index(realizations, name='realization'))\n",
    "    \n",
    "    print(da_past.realization.size)\n",
    "    \n",
    "    np.testing.assert_equal((int(period_past.stop) - int(period_past.start) + 1)*12, da_past.time.size)\n",
    "    \n",
    "    # tas (K -> °C)\n",
    "    #################################################################################\n",
    "\n",
    "    da_past -= 273.15\n",
    "        \n",
    "    #################################################################################\n",
    "    \n",
    "    # Regrid\n",
    "    regridder = xe.Regridder(da_past, elevation, 'bilinear', periodic=False, reuse_weights=True)\n",
    "    da_past_regrid = regridder(da_past)\n",
    "    \n",
    "    if 'height' in da_past_regrid.coords: da_past_regrid = da_past_regrid.drop('height')\n",
    "    \n",
    "    # Store data on CICLAD\n",
    "    path = '/data/mlalande/CMIP6_HMA_paper'\n",
    "    da_past_regrid.sel(lat=latlim, lon=lonlim).to_netcdf(path+'/'+var+'/'+var+'_'+table+'_'+model+'_historical_ens_'+grid+'_197901-201412_HMA.nc')\n",
    "        \n",
    "    print('\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## snc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Snow Cover Extent'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var = 'snc'\n",
    "table = u.get_table(var)\n",
    "label, units, levels, cmap, extend, \\\n",
    "    levels_diff, cmap_diff, extend_diff, \\\n",
    "    levels_bias, cmap_bias, extend_bias = u.get_var_infos(var)\n",
    "label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BCC-CSM2-MR (1/27)\n",
      "['r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1']\n",
      "3\n",
      "Reuse existing file: bilinear_32x54_35x60.nc\n",
      "\n",
      "\n",
      "BCC-ESM1 (2/27)\n",
      "['r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1']\n",
      "3\n",
      "Reuse existing file: bilinear_13x21_35x60.nc\n",
      "\n",
      "\n",
      "CAS-ESM2-0 (3/27)\n",
      "['r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1']\n",
      "4\n",
      "Reuse existing file: bilinear_24x42_35x60.nc\n",
      "\n",
      "\n",
      "CESM2 (4/27)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'snc' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'snc' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'snc' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'snc' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'snc' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['r10i1p1f1', 'r11i1p1f1', 'r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "11\n",
      "Reuse existing file: bilinear_37x49_35x60.nc\n",
      "\n",
      "\n",
      "CESM2-FV2 (5/27)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'snc' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'snc' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1']\n",
      "3\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/dask/utils.py:31: RuntimeWarning: All-NaN slice encountered\n",
      "  return func(*args, **kwargs)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reuse existing file: bilinear_18x25_35x60.nc\n",
      "\n",
      "\n",
      "CESM2-WACCM (6/27)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'snc' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1']\n",
      "3\n",
      "Reuse existing file: bilinear_37x49_35x60.nc\n",
      "\n",
      "\n",
      "CESM2-WACCM-FV2 (7/27)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'snc' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'snc' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'snc' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1']\n",
      "3\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/dask/utils.py:31: RuntimeWarning: All-NaN slice encountered\n",
      "  return func(*args, **kwargs)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reuse existing file: bilinear_18x25_35x60.nc\n",
      "\n",
      "\n",
      "CNRM-CM6-1 (8/27)\n",
      "['r10i1p1f2', 'r11i1p1f2', 'r12i1p1f2', 'r13i1p1f2', 'r14i1p1f2', 'r15i1p1f2', 'r16i1p1f2', 'r17i1p1f2', 'r18i1p1f2', 'r19i1p1f2', 'r1i1p1f2', 'r20i1p1f2', 'r21i1p1f2', 'r22i1p1f2', 'r23i1p1f2', 'r24i1p1f2', 'r25i1p1f2', 'r26i1p1f2', 'r27i1p1f2', 'r28i1p1f2', 'r29i1p1f2', 'r2i1p1f2', 'r30i1p1f2', 'r3i1p1f2', 'r4i1p1f2', 'r5i1p1f2', 'r6i1p1f2', 'r7i1p1f2', 'r8i1p1f2', 'r9i1p1f2']\n",
      "30\n",
      "Reuse existing file: bilinear_25x42_35x60.nc\n",
      "\n",
      "\n",
      "CNRM-CM6-1-HR (9/27)\n",
      "['r1i1p1f2']\n",
      "1\n",
      "Reuse existing file: bilinear_70x121_35x60.nc\n",
      "\n",
      "\n",
      "CNRM-ESM2-1 (10/27)\n",
      "['r10i1p1f2', 'r11i1p1f2', 'r1i1p1f2', 'r2i1p1f2', 'r3i1p1f2', 'r4i1p1f2', 'r5i1p1f2', 'r6i1p1f2', 'r8i1p1f2', 'r9i1p1f2']\n",
      "10\n",
      "Reuse existing file: bilinear_25x42_35x60.nc\n",
      "\n",
      "\n",
      "CanESM5 (11/27)\n",
      "['r10i1p1f1', 'r10i1p2f1', 'r11i1p1f1', 'r11i1p2f1', 'r12i1p1f1', 'r12i1p2f1', 'r13i1p1f1', 'r13i1p2f1', 'r14i1p1f1', 'r14i1p2f1', 'r15i1p1f1', 'r15i1p2f1', 'r16i1p1f1', 'r16i1p2f1', 'r17i1p1f1', 'r17i1p2f1', 'r18i1p1f1', 'r18i1p2f1', 'r19i1p1f1', 'r19i1p2f1', 'r1i1p1f1', 'r1i1p2f1', 'r20i1p1f1', 'r20i1p2f1', 'r21i1p1f1', 'r21i1p2f1', 'r22i1p1f1', 'r22i1p2f1', 'r23i1p1f1', 'r23i1p2f1', 'r24i1p1f1', 'r24i1p2f1', 'r25i1p1f1', 'r25i1p2f1', 'r26i1p2f1', 'r27i1p2f1', 'r28i1p2f1', 'r29i1p2f1', 'r2i1p1f1', 'r2i1p2f1', 'r30i1p2f1', 'r31i1p2f1', 'r32i1p2f1', 'r33i1p2f1', 'r34i1p2f1', 'r35i1p2f1', 'r36i1p2f1', 'r37i1p2f1', 'r38i1p2f1', 'r39i1p2f1', 'r3i1p1f1', 'r3i1p2f1', 'r40i1p2f1', 'r4i1p1f1', 'r4i1p2f1', 'r5i1p1f1', 'r5i1p2f1', 'r6i1p1f1', 'r6i1p2f1', 'r7i1p1f1', 'r7i1p2f1', 'r8i1p1f1', 'r8i1p2f1', 'r9i1p1f1', 'r9i1p2f1']\n",
      "65\n",
      "Reuse existing file: bilinear_13x21_35x60.nc\n",
      "\n",
      "\n",
      "CanESM5-CanOE (12/27)\n",
      "['r1i1p2f1', 'r2i1p2f1', 'r3i1p2f1']\n",
      "3\n",
      "Reuse existing file: bilinear_13x21_35x60.nc\n",
      "\n",
      "\n",
      "GFDL-CM4 (13/27)\n",
      "['r1i1p1f1']\n",
      "1\n",
      "Reuse existing file: bilinear_35x48_35x60.nc\n",
      "\n",
      "\n",
      "GISS-E2-1-G (14/27)\n",
      "['r102i1p1f1', 'r10i1p1f1', 'r10i1p1f2', 'r10i1p3f1', 'r10i1p5f1', 'r11i1p1f2', 'r1i1p1f1', 'r1i1p1f2', 'r1i1p1f3', 'r1i1p3f1', 'r1i1p5f1', 'r2i1p1f1', 'r2i1p1f2', 'r2i1p1f3', 'r2i1p3f1', 'r2i1p5f1', 'r3i1p1f1', 'r3i1p1f2', 'r3i1p1f3', 'r3i1p3f1', 'r3i1p5f1', 'r4i1p1f1', 'r4i1p1f2', 'r4i1p1f3', 'r4i1p3f1', 'r4i1p5f1', 'r5i1p1f1', 'r5i1p1f2', 'r5i1p1f3', 'r5i1p3f1', 'r6i1p1f1', 'r6i1p1f2', 'r6i1p3f1', 'r6i1p5f1', 'r7i1p1f1', 'r7i1p1f2', 'r7i1p3f1', 'r7i1p5f1', 'r8i1p1f1', 'r8i1p1f2', 'r8i1p3f1', 'r8i1p5f1', 'r9i1p1f1', 'r9i1p1f2', 'r9i1p3f1', 'r9i1p5f1']\n",
      "46\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/dask/utils.py:31: RuntimeWarning: All-NaN slice encountered\n",
      "  return func(*args, **kwargs)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reuse existing file: bilinear_18x24_35x60.nc\n",
      "\n",
      "\n",
      "GISS-E2-1-H (15/27)\n",
      "['r10i1p1f1', 'r1i1p1f1', 'r1i1p1f2', 'r1i1p3f1', 'r1i1p5f1', 'r2i1p1f1', 'r2i1p1f2', 'r2i1p3f1', 'r2i1p5f1', 'r3i1p1f1', 'r3i1p1f2', 'r3i1p3f1', 'r3i1p5f1', 'r4i1p1f1', 'r4i1p1f2', 'r4i1p3f1', 'r5i1p1f1', 'r5i1p1f2', 'r5i1p3f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "23\n",
      "Reuse existing file: bilinear_18x24_35x60.nc\n",
      "\n",
      "\n",
      "HadGEM3-GC31-LL (16/27)\n",
      "['r1i1p1f3', 'r2i1p1f3', 'r3i1p1f3', 'r4i1p1f3']\n",
      "4\n",
      "Reuse existing file: bilinear_28x32_35x60.nc\n",
      "\n",
      "\n",
      "HadGEM3-GC31-MM (17/27)\n",
      "['r1i1p1f3', 'r2i1p1f3', 'r3i1p1f3', 'r4i1p1f3']\n",
      "4\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/dask/utils.py:31: RuntimeWarning: All-NaN slice encountered\n",
      "  return func(*args, **kwargs)\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/dask/core.py:121: RuntimeWarning: All-NaN slice encountered\n",
      "  return func(*(_execute_task(a, cache) for a in args))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reuse existing file: bilinear_63x72_35x60.nc\n",
      "\n",
      "\n",
      "IPSL-CM6A-LR (18/27)\n",
      "['r10i1p1f1', 'r11i1p1f1', 'r12i1p1f1', 'r13i1p1f1', 'r14i1p1f1', 'r15i1p1f1', 'r16i1p1f1', 'r17i1p1f1', 'r18i1p1f1', 'r19i1p1f1', 'r1i1p1f1', 'r20i1p1f1', 'r21i1p1f1', 'r22i1p1f1', 'r23i1p1f1', 'r24i1p1f1', 'r25i1p1f1', 'r26i1p1f1', 'r27i1p1f1', 'r28i1p1f1', 'r29i1p1f1', 'r2i1p1f1', 'r30i1p1f1', 'r31i1p1f1', 'r32i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "32\n",
      "Reuse existing file: bilinear_28x25_35x60.nc\n",
      "\n",
      "\n",
      "MIROC-ES2L (19/27)\n",
      "['r10i1p1f2', 'r1i1p1f2', 'r2i1p1f2', 'r3i1p1f2', 'r4i1p1f2', 'r5i1p1f2', 'r6i1p1f2', 'r7i1p1f2', 'r8i1p1f2', 'r9i1p1f2']\n",
      "10\n",
      "Reuse existing file: bilinear_13x21_35x60.nc\n",
      "\n",
      "\n",
      "MIROC6 (20/27)\n",
      "['r10i1p1f1', 'r11i1p1f1', 'r12i1p1f1', 'r13i1p1f1', 'r14i1p1f1', 'r15i1p1f1', 'r16i1p1f1', 'r17i1p1f1', 'r18i1p1f1', 'r19i1p1f1', 'r1i1p1f1', 'r20i1p1f1', 'r21i1p1f1', 'r22i1p1f1', 'r23i1p1f1', 'r24i1p1f1', 'r25i1p1f1', 'r26i1p1f1', 'r27i1p1f1', 'r28i1p1f1', 'r29i1p1f1', 'r2i1p1f1', 'r30i1p1f1', 'r31i1p1f1', 'r32i1p1f1', 'r33i1p1f1', 'r34i1p1f1', 'r35i1p1f1', 'r36i1p1f1', 'r37i1p1f1', 'r38i1p1f1', 'r39i1p1f1', 'r3i1p1f1', 'r40i1p1f1', 'r41i1p1f1', 'r42i1p1f1', 'r43i1p1f1', 'r44i1p1f1', 'r45i1p1f1', 'r46i1p1f1', 'r47i1p1f1', 'r48i1p1f1', 'r49i1p1f1', 'r4i1p1f1', 'r50i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "50\n",
      "Reuse existing file: bilinear_25x42_35x60.nc\n",
      "\n",
      "\n",
      "MPI-ESM1-2-HR (21/27)\n",
      "['r10i1p1f1', 'r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "10\n",
      "Reuse existing file: bilinear_37x64_35x60.nc\n",
      "\n",
      "\n",
      "MPI-ESM1-2-LR (22/27)\n",
      "['r10i1p1f1', 'r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "10\n",
      "Reuse existing file: bilinear_19x32_35x60.nc\n",
      "\n",
      "\n",
      "MRI-ESM2-0 (23/27)\n",
      "['r1i1p1f1', 'r1i2p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1']\n",
      "6\n",
      "Reuse existing file: bilinear_32x54_35x60.nc\n",
      "\n",
      "\n",
      "NorESM2-LM (24/27)\n",
      "['r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1']\n",
      "3\n",
      "Reuse existing file: bilinear_18x25_35x60.nc\n",
      "\n",
      "\n",
      "SAM0-UNICON (25/27)\n",
      "['r1i1p1f1']\n",
      "1\n",
      "Reuse existing file: bilinear_37x49_35x60.nc\n",
      "\n",
      "\n",
      "TaiESM1 (26/27)\n",
      "['r1i1p1f1']\n",
      "1\n",
      "Reuse existing file: bilinear_37x49_35x60.nc\n",
      "\n",
      "\n",
      "UKESM1-0-LL (27/27)\n",
      "['r10i1p1f2', 'r16i1p1f2', 'r17i1p1f2', 'r18i1p1f2', 'r19i1p1f2', 'r1i1p1f2', 'r2i1p1f2', 'r3i1p1f2', 'r4i1p1f2', 'r5i1p1f3', 'r6i1p1f3', 'r7i1p1f3', 'r8i1p1f2', 'r9i1p1f2']\n",
      "14\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/dask/utils.py:31: RuntimeWarning: All-NaN slice encountered\n",
      "  return func(*args, **kwargs)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reuse existing file: bilinear_28x32_35x60.nc\n",
      "\n",
      "\n",
      "CPU times: user 3min 17s, sys: 42.6 s, total: 3min 59s\n",
      "Wall time: 20min 46s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "list_models = u.get_model_names()\n",
    "# list_models = list_models[17:]\n",
    "# list_models = ['EC-Earth3']\n",
    "\n",
    "for i, model in enumerate(list_models):\n",
    "    print(model + ' ('+str(i+1)+'/'+str(len(list_models))+')')\n",
    "    institute, grid, realization, calendar = u.get_model_infos(model, var)\n",
    "    \n",
    "    #################\n",
    "    ### Load past ###\n",
    "    #################\n",
    "    \n",
    "    realizations = !ls {'/bdd/CMIP6/CMIP/'+institute+'/'+model+'/historical/'}\n",
    "    temp = []\n",
    "    no_file = []\n",
    "    \n",
    "    for i, realization in enumerate(realizations):\n",
    "        try:\n",
    "            da = xr.open_mfdataset(\n",
    "                '/bdd/CMIP6/CMIP/'+institute+'/'+model+'/historical/'+realization+'/'+table+'/'+var+'/'+grid+'/latest/*.nc', \n",
    "                combine='by_coords'\n",
    "            ).sel(time=period_past, lat=latlim_ext, lon=lonlim_ext)[var]\n",
    "            \n",
    "            temp.append(da)\n",
    "            \n",
    "        except OSError:\n",
    "            no_file.append(realization)\n",
    "            \n",
    "    for nf in no_file:\n",
    "        realizations.remove(nf)\n",
    "        \n",
    "    print(realizations)\n",
    "    \n",
    "    da_past = xr.concat(temp, pd.Index(realizations, name='realization'))\n",
    "    \n",
    "    print(da_past.realization.size)\n",
    "    \n",
    "    np.testing.assert_equal((int(period_past.stop) - int(period_past.start) + 1)*12, da_past.time.size)\n",
    "    \n",
    "    # snc (0-1 to % when needed and take out wrong values)\n",
    "    #################################################################################\n",
    "    \n",
    "    # Some models have values between 0 and 1 even if units in %    \n",
    "    if da_past.max().values <= 1.1:\n",
    "        da_past = da_past*100\n",
    "    \n",
    "    # Some models have wrong values (+-0.1 for numeric precision)\n",
    "    da_past = da_past.where( (da_past >= -0.1) & (da_past <= 100.1) )\n",
    "    \n",
    "    #################################################################################\n",
    "    \n",
    "    # Regrid\n",
    "    regridder = xe.Regridder(da_past, elevation, 'bilinear', periodic=False, reuse_weights=True)\n",
    "    da_past_regrid = regridder(da_past)\n",
    "    \n",
    "    if 'height' in da_past_regrid.coords: da_past_regrid = da_past_regrid.drop('height')\n",
    "    \n",
    "    # Store data on CICLAD\n",
    "    path = '/data/mlalande/CMIP6_HMA_paper'\n",
    "    da_past_regrid.sel(lat=latlim, lon=lonlim).to_netcdf(path+'/'+var+'/'+var+'_'+table+'_'+model+'_historical_ens_'+grid+'_197901-201412_HMA.nc')\n",
    "        \n",
    "    print('\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## pr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Total Precipitation'"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var = 'pr'\n",
    "table = u.get_table(var)\n",
    "label, units, levels, cmap, extend, \\\n",
    "    levels_diff, cmap_diff, extend_diff, \\\n",
    "    levels_bias, cmap_bias, extend_bias = u.get_var_infos(var)\n",
    "label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CESM2 (1/24)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'pr' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'pr' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'pr' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'pr' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'pr' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['r10i1p1f1', 'r11i1p1f1', 'r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "11\n",
      "Reuse existing file: bilinear_37x49_35x60.nc\n",
      "\n",
      "\n",
      "CESM2-FV2 (2/24)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'pr' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'pr' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'pr' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1']\n",
      "3\n",
      "Reuse existing file: bilinear_18x25_35x60.nc\n",
      "\n",
      "\n",
      "CESM2-WACCM (3/24)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'pr' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1']\n",
      "3\n",
      "Reuse existing file: bilinear_37x49_35x60.nc\n",
      "\n",
      "\n",
      "CESM2-WACCM-FV2 (4/24)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'pr' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'pr' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n",
      "/data/mlalande/miniconda3/envs/phd_v3/lib/python3.8/site-packages/xarray/conventions.py:490: SerializationWarning: variable 'pr' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1']\n",
      "3\n",
      "Reuse existing file: bilinear_18x25_35x60.nc\n",
      "\n",
      "\n",
      "CNRM-CM6-1 (5/24)\n",
      "['r10i1p1f2', 'r11i1p1f2', 'r12i1p1f2', 'r13i1p1f2', 'r14i1p1f2', 'r15i1p1f2', 'r16i1p1f2', 'r17i1p1f2', 'r18i1p1f2', 'r19i1p1f2', 'r1i1p1f2', 'r20i1p1f2', 'r21i1p1f2', 'r22i1p1f2', 'r24i1p1f2', 'r25i1p1f2', 'r26i1p1f2', 'r27i1p1f2', 'r28i1p1f2', 'r29i1p1f2', 'r2i1p1f2', 'r30i1p1f2', 'r3i1p1f2', 'r4i1p1f2', 'r5i1p1f2', 'r6i1p1f2', 'r7i1p1f2', 'r8i1p1f2', 'r9i1p1f2']\n",
      "29\n",
      "Reuse existing file: bilinear_25x42_35x60.nc\n",
      "\n",
      "\n",
      "CNRM-CM6-1-HR (6/24)\n",
      "['r1i1p1f2']\n",
      "1\n",
      "Reuse existing file: bilinear_70x121_35x60.nc\n",
      "\n",
      "\n",
      "CNRM-ESM2-1 (7/24)\n",
      "['r10i1p1f2', 'r11i1p1f2', 'r1i1p1f2', 'r2i1p1f2', 'r3i1p1f2', 'r4i1p1f2', 'r5i1p1f2', 'r7i1p1f2', 'r8i1p1f2', 'r9i1p1f2']\n",
      "10\n",
      "Reuse existing file: bilinear_25x42_35x60.nc\n",
      "\n",
      "\n",
      "CanESM5 (8/24)\n",
      "['r10i1p1f1', 'r10i1p2f1', 'r11i1p1f1', 'r11i1p2f1', 'r12i1p1f1', 'r12i1p2f1', 'r13i1p1f1', 'r13i1p2f1', 'r14i1p1f1', 'r14i1p2f1', 'r15i1p1f1', 'r15i1p2f1', 'r16i1p1f1', 'r16i1p2f1', 'r17i1p1f1', 'r17i1p2f1', 'r18i1p1f1', 'r18i1p2f1', 'r19i1p1f1', 'r19i1p2f1', 'r1i1p1f1', 'r1i1p2f1', 'r20i1p1f1', 'r20i1p2f1', 'r21i1p1f1', 'r21i1p2f1', 'r22i1p1f1', 'r22i1p2f1', 'r23i1p1f1', 'r23i1p2f1', 'r24i1p1f1', 'r24i1p2f1', 'r25i1p1f1', 'r25i1p2f1', 'r26i1p2f1', 'r27i1p2f1', 'r28i1p2f1', 'r29i1p2f1', 'r2i1p1f1', 'r2i1p2f1', 'r30i1p2f1', 'r31i1p2f1', 'r32i1p2f1', 'r33i1p2f1', 'r34i1p2f1', 'r35i1p2f1', 'r36i1p2f1', 'r37i1p2f1', 'r38i1p2f1', 'r39i1p2f1', 'r3i1p1f1', 'r3i1p2f1', 'r40i1p2f1', 'r4i1p1f1', 'r4i1p2f1', 'r5i1p1f1', 'r5i1p2f1', 'r6i1p1f1', 'r6i1p2f1', 'r7i1p1f1', 'r7i1p2f1', 'r8i1p1f1', 'r8i1p2f1', 'r9i1p1f1', 'r9i1p2f1']\n",
      "65\n",
      "Reuse existing file: bilinear_13x21_35x60.nc\n",
      "\n",
      "\n",
      "CanESM5-CanOE (9/24)\n",
      "['r1i1p2f1', 'r2i1p2f1', 'r3i1p2f1']\n",
      "3\n",
      "Reuse existing file: bilinear_13x21_35x60.nc\n",
      "\n",
      "\n",
      "GFDL-CM4 (10/24)\n",
      "['r1i1p1f1']\n",
      "1\n",
      "Reuse existing file: bilinear_35x48_35x60.nc\n",
      "\n",
      "\n",
      "GISS-E2-1-G (11/24)\n",
      "['r101i1p1f1', 'r102i1p1f1', 'r10i1p1f1', 'r10i1p1f2', 'r10i1p3f1', 'r10i1p5f1', 'r11i1p1f2', 'r1i1p1f1', 'r1i1p1f2', 'r1i1p1f3', 'r1i1p3f1', 'r1i1p5f1', 'r2i1p1f1', 'r2i1p1f2', 'r2i1p1f3', 'r2i1p3f1', 'r2i1p5f1', 'r3i1p1f1', 'r3i1p1f2', 'r3i1p1f3', 'r3i1p3f1', 'r3i1p5f1', 'r4i1p1f1', 'r4i1p1f2', 'r4i1p1f3', 'r4i1p3f1', 'r4i1p5f1', 'r5i1p1f1', 'r5i1p1f2', 'r5i1p1f3', 'r5i1p3f1', 'r6i1p1f1', 'r6i1p1f2', 'r6i1p3f1', 'r6i1p5f1', 'r7i1p1f1', 'r7i1p1f2', 'r7i1p3f1', 'r7i1p5f1', 'r8i1p1f1', 'r8i1p1f2', 'r8i1p3f1', 'r8i1p5f1', 'r9i1p1f1', 'r9i1p1f2', 'r9i1p3f1', 'r9i1p5f1']\n",
      "47\n",
      "Reuse existing file: bilinear_18x24_35x60.nc\n",
      "\n",
      "\n",
      "GISS-E2-1-H (12/24)\n",
      "['r10i1p1f1', 'r1i1p1f1', 'r1i1p1f2', 'r1i1p3f1', 'r1i1p5f1', 'r2i1p1f1', 'r2i1p1f2', 'r2i1p3f1', 'r2i1p5f1', 'r3i1p1f1', 'r3i1p1f2', 'r3i1p3f1', 'r3i1p5f1', 'r4i1p1f1', 'r4i1p1f2', 'r4i1p3f1', 'r4i1p5f1', 'r5i1p1f1', 'r5i1p1f2', 'r5i1p3f1', 'r5i1p5f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "25\n",
      "Reuse existing file: bilinear_18x24_35x60.nc\n",
      "\n",
      "\n",
      "HadGEM3-GC31-LL (13/24)\n",
      "['r1i1p1f3', 'r2i1p1f3', 'r3i1p1f3', 'r4i1p1f3']\n",
      "4\n",
      "Reuse existing file: bilinear_28x32_35x60.nc\n",
      "\n",
      "\n",
      "HadGEM3-GC31-MM (14/24)\n",
      "['r1i1p1f3', 'r2i1p1f3', 'r3i1p1f3', 'r4i1p1f3']\n",
      "4\n",
      "Reuse existing file: bilinear_63x72_35x60.nc\n",
      "\n",
      "\n",
      "IPSL-CM6A-LR (15/24)\n",
      "['r10i1p1f1', 'r11i1p1f1', 'r12i1p1f1', 'r13i1p1f1', 'r14i1p1f1', 'r15i1p1f1', 'r16i1p1f1', 'r17i1p1f1', 'r18i1p1f1', 'r19i1p1f1', 'r1i1p1f1', 'r20i1p1f1', 'r21i1p1f1', 'r22i1p1f1', 'r23i1p1f1', 'r24i1p1f1', 'r25i1p1f1', 'r26i1p1f1', 'r27i1p1f1', 'r28i1p1f1', 'r29i1p1f1', 'r2i1p1f1', 'r30i1p1f1', 'r31i1p1f1', 'r32i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "32\n",
      "Reuse existing file: bilinear_28x25_35x60.nc\n",
      "\n",
      "\n",
      "MIROC-ES2L (16/24)\n",
      "['r10i1p1f2', 'r1i1p1f2', 'r2i1p1f2', 'r3i1p1f2', 'r4i1p1f2', 'r5i1p1f2', 'r6i1p1f2', 'r7i1p1f2', 'r8i1p1f2', 'r9i1p1f2']\n",
      "10\n",
      "Reuse existing file: bilinear_13x21_35x60.nc\n",
      "\n",
      "\n",
      "MIROC6 (17/24)\n",
      "['r10i1p1f1', 'r11i1p1f1', 'r12i1p1f1', 'r13i1p1f1', 'r14i1p1f1', 'r15i1p1f1', 'r16i1p1f1', 'r17i1p1f1', 'r18i1p1f1', 'r19i1p1f1', 'r1i1p1f1', 'r20i1p1f1', 'r21i1p1f1', 'r22i1p1f1', 'r23i1p1f1', 'r24i1p1f1', 'r25i1p1f1', 'r26i1p1f1', 'r27i1p1f1', 'r28i1p1f1', 'r29i1p1f1', 'r2i1p1f1', 'r30i1p1f1', 'r31i1p1f1', 'r32i1p1f1', 'r33i1p1f1', 'r34i1p1f1', 'r35i1p1f1', 'r36i1p1f1', 'r37i1p1f1', 'r38i1p1f1', 'r39i1p1f1', 'r3i1p1f1', 'r40i1p1f1', 'r41i1p1f1', 'r42i1p1f1', 'r43i1p1f1', 'r44i1p1f1', 'r45i1p1f1', 'r46i1p1f1', 'r47i1p1f1', 'r48i1p1f1', 'r49i1p1f1', 'r4i1p1f1', 'r50i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "50\n",
      "Reuse existing file: bilinear_25x42_35x60.nc\n",
      "\n",
      "\n",
      "MPI-ESM1-2-HR (18/24)\n",
      "['r10i1p1f1', 'r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "10\n",
      "Reuse existing file: bilinear_37x64_35x60.nc\n",
      "\n",
      "\n",
      "MPI-ESM1-2-LR (19/24)\n",
      "['r10i1p1f1', 'r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "10\n",
      "Reuse existing file: bilinear_19x32_35x60.nc\n",
      "\n",
      "\n",
      "MRI-ESM2-0 (20/24)\n",
      "['r1i1000p1f1', 'r1i1p1f1', 'r1i2p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1']\n",
      "7\n",
      "Reuse existing file: bilinear_32x54_35x60.nc\n",
      "\n",
      "\n",
      "NorESM2-LM (21/24)\n",
      "['r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1']\n",
      "3\n",
      "Reuse existing file: bilinear_18x25_35x60.nc\n",
      "\n",
      "\n",
      "SAM0-UNICON (22/24)\n",
      "['r1i1p1f1']\n",
      "1\n",
      "Reuse existing file: bilinear_37x49_35x60.nc\n",
      "\n",
      "\n",
      "TaiESM1 (23/24)\n",
      "['r1i1p1f1']\n",
      "1\n",
      "Reuse existing file: bilinear_37x49_35x60.nc\n",
      "\n",
      "\n",
      "UKESM1-0-LL (24/24)\n",
      "['r10i1p1f2', 'r11i1p1f2', 'r12i1p1f2', 'r16i1p1f2', 'r17i1p1f2', 'r18i1p1f2', 'r19i1p1f2', 'r1i1p1f2', 'r2i1p1f2', 'r3i1p1f2', 'r4i1p1f2', 'r5i1p1f3', 'r6i1p1f3', 'r7i1p1f3', 'r8i1p1f2', 'r9i1p1f2']\n",
      "16\n",
      "Reuse existing file: bilinear_28x32_35x60.nc\n",
      "\n",
      "\n",
      "CPU times: user 2min 34s, sys: 1min, total: 3min 35s\n",
      "Wall time: 26min 21s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "list_models = u.get_model_names()\n",
    "list_models = list_models[3:]\n",
    "# list_models = ['EC-Earth3']\n",
    "\n",
    "for i, model in enumerate(list_models):\n",
    "    print(model + ' ('+str(i+1)+'/'+str(len(list_models))+')')\n",
    "    institute, grid, realization, calendar = u.get_model_infos(model, var)\n",
    "    \n",
    "    #################\n",
    "    ### Load past ###\n",
    "    #################\n",
    "    \n",
    "    realizations = !ls {'/bdd/CMIP6/CMIP/'+institute+'/'+model+'/historical/'}\n",
    "    temp = []\n",
    "    no_file = []\n",
    "    \n",
    "    for i, realization in enumerate(realizations):\n",
    "        try:\n",
    "            if model == 'CESM2' and realization == 'r9i1p1f1':\n",
    "                da = xr.open_mfdataset(\n",
    "                    '/bdd/CMIP6/CMIP/'+institute+'/'+model+'/historical/'+realization+'/'+table+'/'+var+'/'+grid+'/latest/*185001-201412.nc', \n",
    "                    combine='by_coords'\n",
    "                ).sel(time=period_past, lat=latlim_ext, lon=lonlim_ext)[var]\n",
    "            else:\n",
    "                da = xr.open_mfdataset(\n",
    "                    '/bdd/CMIP6/CMIP/'+institute+'/'+model+'/historical/'+realization+'/'+table+'/'+var+'/'+grid+'/latest/*.nc', \n",
    "                    combine='by_coords'\n",
    "                ).sel(time=period_past, lat=latlim_ext, lon=lonlim_ext)[var]\n",
    "            \n",
    "            temp.append(da)\n",
    "            \n",
    "        except OSError:\n",
    "            no_file.append(realization)\n",
    "            \n",
    "    for nf in no_file:\n",
    "        realizations.remove(nf)\n",
    "        \n",
    "    print(realizations)\n",
    "    \n",
    "    da_past = xr.concat(temp, pd.Index(realizations, name='realization'))\n",
    "    \n",
    "    print(da_past.realization.size)\n",
    "    \n",
    "    np.testing.assert_equal((int(period_past.stop) - int(period_past.start) + 1)*12, da_past.time.size)\n",
    "    \n",
    "    # pr (mm/s -> mm/day)\n",
    "    #################################################################################\n",
    "    \n",
    "    da_past = da_past * 86400\n",
    "    \n",
    "    #################################################################################\n",
    "    \n",
    "    # Regrid\n",
    "    regridder = xe.Regridder(da_past, elevation, 'bilinear', periodic=False, reuse_weights=True)\n",
    "    da_past_regrid = regridder(da_past)\n",
    "    \n",
    "    if 'height' in da_past_regrid.coords: da_past_regrid = da_past_regrid.drop('height')\n",
    "    \n",
    "    # Store data on CICLAD\n",
    "    path = '/data/mlalande/CMIP6_HMA_paper'\n",
    "    da_past_regrid.sel(lat=latlim, lon=lonlim).to_netcdf(path+'/'+var+'/'+var+'_'+table+'_'+model+'_historical_ens_'+grid+'_197901-201412_HMA.nc')\n",
    "        \n",
    "    print('\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/bdd/CMIP6/CMIP/NCAR/CESM2/historical/r9i1p1f1/Amon/pr/gn/'"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'/bdd/CMIP6/CMIP/'+institute+'/'+model+'/historical/'+realization+'/'+table+'/'+var+'/'+grid+'/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ta 500 hPa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Air Temperature'"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var = 'ta'\n",
    "level = 500\n",
    "table = u.get_table(var)\n",
    "label, units, levels, cmap, extend, \\\n",
    "    levels_diff, cmap_diff, extend_diff, \\\n",
    "    levels_bias, cmap_bias, extend_bias = u.get_var_infos(var)\n",
    "label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IPSL-CM6A-LR (1/10)\n",
      "['r10i1p1f1', 'r11i1p1f1', 'r12i1p1f1', 'r13i1p1f1', 'r14i1p1f1', 'r15i1p1f1', 'r16i1p1f1', 'r17i1p1f1', 'r18i1p1f1', 'r19i1p1f1', 'r1i1p1f1', 'r20i1p1f1', 'r21i1p1f1', 'r22i1p1f1', 'r23i1p1f1', 'r24i1p1f1', 'r25i1p1f1', 'r26i1p1f1', 'r27i1p1f1', 'r28i1p1f1', 'r29i1p1f1', 'r2i1p1f1', 'r30i1p1f1', 'r31i1p1f1', 'r32i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "32\n",
      "Reuse existing file: bilinear_28x25_35x60.nc\n",
      "\n",
      "\n",
      "MIROC-ES2L (2/10)\n",
      "['r10i1p1f2', 'r1i1p1f2', 'r2i1p1f2', 'r3i1p1f2', 'r4i1p1f2', 'r5i1p1f2', 'r6i1p1f2', 'r8i1p1f2', 'r9i1p1f2']\n",
      "9\n",
      "Reuse existing file: bilinear_13x21_35x60.nc\n",
      "\n",
      "\n",
      "MIROC6 (3/10)\n",
      "['r10i1p1f1', 'r11i1p1f1', 'r12i1p1f1', 'r13i1p1f1', 'r14i1p1f1', 'r15i1p1f1', 'r16i1p1f1', 'r17i1p1f1', 'r18i1p1f1', 'r19i1p1f1', 'r1i1p1f1', 'r20i1p1f1', 'r21i1p1f1', 'r22i1p1f1', 'r23i1p1f1', 'r24i1p1f1', 'r25i1p1f1', 'r26i1p1f1', 'r27i1p1f1', 'r28i1p1f1', 'r29i1p1f1', 'r2i1p1f1', 'r30i1p1f1', 'r31i1p1f1', 'r32i1p1f1', 'r33i1p1f1', 'r34i1p1f1', 'r35i1p1f1', 'r36i1p1f1', 'r37i1p1f1', 'r38i1p1f1', 'r39i1p1f1', 'r3i1p1f1', 'r40i1p1f1', 'r41i1p1f1', 'r42i1p1f1', 'r43i1p1f1', 'r44i1p1f1', 'r45i1p1f1', 'r46i1p1f1', 'r47i1p1f1', 'r48i1p1f1', 'r49i1p1f1', 'r4i1p1f1', 'r50i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "50\n",
      "Reuse existing file: bilinear_25x42_35x60.nc\n",
      "\n",
      "\n",
      "MPI-ESM1-2-HR (4/10)\n",
      "['r10i1p1f1', 'r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "10\n",
      "Reuse existing file: bilinear_37x64_35x60.nc\n",
      "\n",
      "\n",
      "MPI-ESM1-2-LR (5/10)\n",
      "['r10i1p1f1', 'r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1', 'r6i1p1f1', 'r7i1p1f1', 'r8i1p1f1', 'r9i1p1f1']\n",
      "10\n",
      "Reuse existing file: bilinear_19x32_35x60.nc\n",
      "\n",
      "\n",
      "MRI-ESM2-0 (6/10)\n",
      "['r1i1000p1f1', 'r1i1p1f1', 'r1i2p1f1', 'r2i1p1f1', 'r3i1p1f1', 'r4i1p1f1', 'r5i1p1f1']\n",
      "7\n",
      "Reuse existing file: bilinear_32x54_35x60.nc\n",
      "\n",
      "\n",
      "NorESM2-LM (7/10)\n",
      "['r1i1p1f1', 'r2i1p1f1', 'r3i1p1f1']\n",
      "3\n",
      "Reuse existing file: bilinear_18x25_35x60.nc\n",
      "\n",
      "\n",
      "SAM0-UNICON (8/10)\n",
      "['r1i1p1f1']\n",
      "1\n",
      "Reuse existing file: bilinear_37x49_35x60.nc\n",
      "\n",
      "\n",
      "TaiESM1 (9/10)\n",
      "['r1i1p1f1']\n",
      "1\n",
      "Reuse existing file: bilinear_37x49_35x60.nc\n",
      "\n",
      "\n",
      "UKESM1-0-LL (10/10)\n",
      "['r10i1p1f2', 'r12i1p1f2', 'r16i1p1f2', 'r17i1p1f2', 'r18i1p1f2', 'r19i1p1f2', 'r1i1p1f2', 'r2i1p1f2', 'r3i1p1f2', 'r4i1p1f2', 'r5i1p1f3', 'r6i1p1f3', 'r7i1p1f3', 'r8i1p1f2', 'r9i1p1f2']\n",
      "15\n",
      "Reuse existing file: bilinear_28x32_35x60.nc\n",
      "\n",
      "\n",
      "CPU times: user 14min 30s, sys: 3min 50s, total: 18min 21s\n",
      "Wall time: 47min 57s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "list_models = u.get_model_names()\n",
    "list_models = list_models[17:]\n",
    "# list_models = ['EC-Earth3']\n",
    "\n",
    "for i, model in enumerate(list_models):\n",
    "    print(model + ' ('+str(i+1)+'/'+str(len(list_models))+')')\n",
    "    institute, grid, realization, calendar = u.get_model_infos(model, var)\n",
    "    \n",
    "    #################\n",
    "    ### Load past ###\n",
    "    #################\n",
    "    \n",
    "    realizations = !ls {'/bdd/CMIP6/CMIP/'+institute+'/'+model+'/historical/'}\n",
    "    temp = []\n",
    "    no_file = []\n",
    "    \n",
    "    for i, realization in enumerate(realizations):\n",
    "        try:\n",
    "            da = xr.open_mfdataset(\n",
    "                '/bdd/CMIP6/CMIP/'+institute+'/'+model+'/historical/'+realization+'/'+table+'/'+var+'/'+grid+'/latest/*.nc', \n",
    "                combine='by_coords'\n",
    "            ).sel(time=period_past, lat=latlim_ext, lon=lonlim_ext, plev=level*100)[var]\n",
    "            \n",
    "            temp.append(da)\n",
    "            \n",
    "        except OSError:\n",
    "            no_file.append(realization)\n",
    "            \n",
    "    for nf in no_file:\n",
    "        realizations.remove(nf)\n",
    "        \n",
    "    print(realizations)\n",
    "    \n",
    "    da_past = xr.concat(temp, pd.Index(realizations, name='realization'))\n",
    "    \n",
    "    print(da_past.realization.size)\n",
    "    \n",
    "    np.testing.assert_equal((int(period_past.stop) - int(period_past.start) + 1)*12, da_past.time.size)\n",
    "    \n",
    "    # ta (K -> °C)\n",
    "    #################################################################################\n",
    "\n",
    "    da_past -= 273.15\n",
    "        \n",
    "    #################################################################################\n",
    "    \n",
    "    # Regrid\n",
    "    regridder = xe.Regridder(da_past, elevation, 'bilinear', periodic=False, reuse_weights=True)\n",
    "    da_past_regrid = regridder(da_past)\n",
    "    \n",
    "    if 'height' in da_past_regrid.coords: da_past_regrid = da_past_regrid.drop('height')\n",
    "    \n",
    "    # Store data on CICLAD\n",
    "    path = '/data/mlalande/CMIP6_HMA_paper'\n",
    "    da_past_regrid.sel(lat=latlim, lon=lonlim).to_netcdf(path+'/'+var+str(level)+'/'+var+str(level)+'_'+table+'_'+model+'_historical_ens_'+grid+'_197901-201412_HMA.nc')\n",
    "        \n",
    "    print('\\n')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:phd_v3] *",
   "language": "python",
   "name": "conda-env-phd_v3-py"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}