12_pr.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T02:28:54.587324Z",
     "start_time": "2020-05-16T02:28:45.731762Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "variable = 'pr'\n",
    "var_name = variable + \"/\"\n",
    "model_name_file = 'NorESM2-LM'\n",
    "import numpy as np\n",
    "import time\n",
    "import xarray as xr\n",
    "import os\n",
    "os.environ['KMP_DUPLICATE_LIB_OK']='True'\n",
    "import pickle as pkl\n",
    "import pickle\n",
    "import numpy as np\n",
    "from numpy import *\n",
    "import xarray as xr\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from datetime import datetime\n",
    "import sys\n",
    "from rpy2.robjects.packages import importr\n",
    "import rpy2.robjects.numpy2ri\n",
    "rpy2.robjects.numpy2ri.activate()\n",
    "qmap = importr('qmap')\n",
    "from rpy2.robjects import pandas2ri\n",
    "from rpy2.robjects import r\n",
    "utils = importr(\"utils\")\n",
    "#extremes = importr(\"extRemes\")\n",
    "#ismev = importr(\"ismev\")\n",
    "#eva = importr(\"eva\")\n",
    "gnFit = importr(\"gnFit\")\n",
    "from datetime import datetime\n",
    "statis = importr('stats')\n",
    "import time\n",
    "from datetime import datetime\n",
    "from matplotlib.offsetbox import AnchoredText\n",
    "from netCDF4 import num2date, date2num\n",
    "statis = importr('stats')\n",
    "import rpy2.robjects.numpy2ri\n",
    "rpy2.robjects.numpy2ri.activate()\n",
    "import glob\n",
    "import time\n",
    "from tqdm import tqdm\n",
    "import sys\n",
    "sys.setrecursionlimit(10000)\n",
    "import sys\n",
    "import time\n",
    "import traceback\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import glob\n",
    "import sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T02:28:54.592882Z",
     "start_time": "2020-05-16T02:28:54.589468Z"
    }
   },
   "outputs": [],
   "source": [
    "if sys.platform ==\"darwin\":\n",
    "\n",
    "    gcm_path = '/media/udit/cmip6/cmip6_raw_amar_regrided/'\n",
    "    obs_path = '/Volumes/cmip6/'\n",
    "    target_folder = '/Volumes/cmip6/'\n",
    "else:\n",
    "    gcm_path = '/media/udit/cmip6/cmip6_raw_amar_regrided/'\n",
    "    obs_path = '/media/udit/cmip6/'\n",
    "    target_folder = '/media/udit/cmip6/'\n",
    "    \n",
    "model_name = model_name_file+'/'\n",
    "\n",
    "scenarios = ['hist','ssp126','ssp245','ssp370','ssp585']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T02:28:54.605264Z",
     "start_time": "2020-05-16T02:28:54.594478Z"
    }
   },
   "outputs": [],
   "source": [
    "def bias_correction(x,y):\n",
    "    q_map = qmap.fitQmap(x, y,method=\"RQUANT\", qstep=0.01, wett_day=False)\n",
    "    qm1 = qmap.doQmap(y, q_map)\n",
    "    bias_corrected_output = {}\n",
    "    bias_corrected_output['params'] = q_map\n",
    "    bias_corrected_output['outputs'] = qm1\n",
    "    return bias_corrected_output\n",
    "\n",
    "def bias_correction_model(y,q_map):\n",
    "    qm1 = qmap.doQmap(y, q_map)\n",
    "    bias_corrected_output = {}\n",
    "    bias_corrected_output['outputs'] = qm1\n",
    "    return bias_corrected_output\n",
    "\n",
    "\n",
    "\n",
    "def save_as_pickled_object(obj, filepath):\n",
    "    \"\"\"\n",
    "    This is a defensive way to write pickle.write, allowing for very large files on all platforms\n",
    "    \"\"\"\n",
    "    max_bytes = 2**31 - 1\n",
    "    bytes_out = pickle.dumps(obj)\n",
    "    n_bytes = sys.getsizeof(bytes_out)\n",
    "    with open(filepath, 'wb') as f_out:\n",
    "        for idx in range(0, n_bytes, max_bytes):\n",
    "            f_out.write(bytes_out[idx:idx + max_bytes])\n",
    "\n",
    "\n",
    "def try_to_load_as_pickled_object_or_None(filepath):\n",
    "\n",
    "    max_bytes = 2**31 - 1\n",
    "    try:\n",
    "        input_size = os.path.getsize(filepath)\n",
    "        bytes_in = bytearray(0)\n",
    "        with open(filepath, 'rb') as f_in:\n",
    "            for _ in range(0, input_size, max_bytes):\n",
    "                bytes_in += f_in.read(max_bytes)\n",
    "        obj = pickle.loads(bytes_in)\n",
    "    except:\n",
    "        return None\n",
    "    return obj\n",
    "\n",
    "def read_netcdfs(files, dim, transform_func=None):\n",
    "    def process_one_path(path):\n",
    "        # use a context manager, to ensure the file gets closed after use\n",
    "        with xr.open_dataset(path) as ds:\n",
    "            # transform_func should do some sort of selection or\n",
    "            # aggregation\n",
    "            if transform_func is not None:\n",
    "                ds = transform_func(ds)\n",
    "            # load all data from the transformed dataset, to ensure we can\n",
    "            # use it after closing each original file\n",
    "            ds.load()\n",
    "            return ds\n",
    "\n",
    "    paths = sorted(glob.glob(files))\n",
    "    datasets = [process_one_path(p) for p in paths]\n",
    "    combined = xr.concat(datasets, dim)\n",
    "    return combined\n",
    "\n",
    "def preproc_hist(ds):\n",
    "    ds = ds.sel(time = slice('01-01-1951','2014-12-31'),lat = slice(5,38),lon = slice(60,100))\n",
    "    return ds\n",
    "def preproc_model(ds):\n",
    "    ds = ds.sel(lat = slice(5,38),lon = slice(60,100))\n",
    "    return ds\n",
    "\n",
    "\n",
    "\n",
    "def preproc_hist2(ds):\n",
    "    start = pd.to_datetime('1951-01-01', format='%Y%m%d', errors='ignore')\n",
    "    end = pd.to_datetime('2014-12-31', format='%Y%m%d', errors='ignore')\n",
    "    ds = ds.sel(time = slice(start,end),lat = slice(5,38),lon = slice(60,100))\n",
    "    return ds\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# here we suppose we only care about the combined mean of each file;\n",
    "# you might also use indexing operations like .sel to subset datasets\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T02:28:54.615583Z",
     "start_time": "2020-05-16T02:28:54.606640Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/media/udit/cmip6/cmip6_raw_amar_regrided/NorESM2-LM/\n"
     ]
    }
   ],
   "source": [
    "model_path = gcm_path + model_name\n",
    "model_path_historical = model_path\n",
    "model_path_RCP126 = model_path \n",
    "model_path_RCP245 = model_path\n",
    "model_path_RCP370 = model_path\n",
    "model_path_RCP585 = model_path\n",
    "\n",
    "print(model_path_historical)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T02:28:54.621905Z",
     "start_time": "2020-05-16T02:28:54.617607Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['/media/udit/cmip6/cmip6_raw_amar_regrided/NorESM2-LM/pr/pr_day_NorESM2-LM_historical_r1i1p1f1_gn_20150101-21001231.nc']\n"
     ]
    }
   ],
   "source": [
    "prcp_hist = model_path +var_name+ \"*historical*.nc\"\n",
    "file_lists = glob.glob(prcp_hist)\n",
    "print(file_lists)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T02:29:14.305442Z",
     "start_time": "2020-05-16T02:28:54.623112Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/udit/.local/lib/python3.6/site-packages/ipykernel_launcher.py:8: FutureWarning: In xarray version 0.15 the default behaviour of `open_mfdataset`\n",
      "will change. To retain the existing behavior, pass\n",
      "combine='nested'. To use future default behavior, pass\n",
      "combine='by_coords'. See\n",
      "http://xarray.pydata.org/en/stable/combining.html#combining-multi\n",
      "\n",
      "  \n",
      "/home/udit/anaconda3/envs/master/lib/python3.6/site-packages/xarray/backends/api.py:931: FutureWarning: Also `open_mfdataset` will no longer accept a `concat_dim` argument.\n",
      "To get equivalent behaviour from now on please use the new\n",
      "`combine_nested` function instead (or the `combine='nested'` option to\n",
      "`open_mfdataset`).The datasets supplied have global dimension coordinates. You may want\n",
      "to use the new `combine_by_coords` function (or the\n",
      "`combine='by_coords'` option to `open_mfdataset`) to order the datasets\n",
      "before concatenation. Alternatively, to continue concatenating based\n",
      "on the order the datasets are supplied in future, please use the new\n",
      "`combine_nested` function (or the `combine='nested'` option to\n",
      "open_mfdataset).\n",
      "  from_openmfds=True,\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/media/udit/cmip6/cmip6_raw_amar_regrided/NorESM2-LM/pr/*245*.nc\n",
      "/media/udit/cmip6/cmip6_raw_amar_regrided/NorESM2-LM/pr/*126*.nc\n",
      "['/media/udit/cmip6/cmip6_raw_amar_regrided/NorESM2-LM/pr/pr_day_NorESM2-LM_ssp245_r1i1p1f1_gn_20150101-21001231.nc']\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/udit/.local/lib/python3.6/site-packages/ipykernel_launcher.py:30: FutureWarning: In xarray version 0.15 the default behaviour of `open_mfdataset`\n",
      "will change. To retain the existing behavior, pass\n",
      "combine='nested'. To use future default behavior, pass\n",
      "combine='by_coords'. See\n",
      "http://xarray.pydata.org/en/stable/combining.html#combining-multi\n",
      "\n",
      "/home/udit/anaconda3/envs/master/lib/python3.6/site-packages/xarray/backends/api.py:931: FutureWarning: Also `open_mfdataset` will no longer accept a `concat_dim` argument.\n",
      "To get equivalent behaviour from now on please use the new\n",
      "`combine_nested` function instead (or the `combine='nested'` option to\n",
      "`open_mfdataset`).The datasets supplied have global dimension coordinates. You may want\n",
      "to use the new `combine_by_coords` function (or the\n",
      "`combine='by_coords'` option to `open_mfdataset`) to order the datasets\n",
      "before concatenation. Alternatively, to continue concatenating based\n",
      "on the order the datasets are supplied in future, please use the new\n",
      "`combine_nested` function (or the `combine='nested'` option to\n",
      "open_mfdataset).\n",
      "  from_openmfds=True,\n",
      "/home/udit/.local/lib/python3.6/site-packages/ipykernel_launcher.py:33: FutureWarning: In xarray version 0.15 the default behaviour of `open_mfdataset`\n",
      "will change. To retain the existing behavior, pass\n",
      "combine='nested'. To use future default behavior, pass\n",
      "combine='by_coords'. See\n",
      "http://xarray.pydata.org/en/stable/combining.html#combining-multi\n",
      "\n",
      "/home/udit/.local/lib/python3.6/site-packages/ipykernel_launcher.py:36: FutureWarning: In xarray version 0.15 the default behaviour of `open_mfdataset`\n",
      "will change. To retain the existing behavior, pass\n",
      "combine='nested'. To use future default behavior, pass\n",
      "combine='by_coords'. See\n",
      "http://xarray.pydata.org/en/stable/combining.html#combining-multi\n",
      "\n",
      "/home/udit/.local/lib/python3.6/site-packages/ipykernel_launcher.py:39: FutureWarning: In xarray version 0.15 the default behaviour of `open_mfdataset`\n",
      "will change. To retain the existing behavior, pass\n",
      "combine='nested'. To use future default behavior, pass\n",
      "combine='by_coords'. See\n",
      "http://xarray.pydata.org/en/stable/combining.html#combining-multi\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEWCAYAAACDoeeyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df5QcZZ3v8fdnEn4HFjABIxIDLKLIleBGQPEoguxFj15gj6DsysYrGtTlePHHCnLPWfGiK7Agi0cXb5QfwWX5KV5YVnQjgog/wAQhgHFFICISkyAkJAghmfneP+oZ6DQ90z3dVd1VNZ/XOXVmuqq66js9Pd95+vs89ZQiAjMzq6ehQQdgZmbFcZI3M6sxJ3kzsxpzkjczqzEneTOzGnOSNzOrMSf5SUTSLEnrJU0ZdCxm1h9O8jUmabmkt40+johHImJaRAwPMq7xSDpZ0mJJGyRd2mL7ByX9Jv2z+q6klzVsO0PSxrRtdNmzYfsCSf8laUTS+zuIZY6kJZL+lL7Oadj2Vkm3SForaXkHx5KksyX9MS3nSFLaNl3Sj9P6NZJ+KumQti+WWQec5K1sHgM+D1zcvEHSW4B/BI4CdgYeBq5o2u2q9I9sdHmoYds9wEeBu9oFIWlL4HrgX4GdgIXA9Wk9wNMpxr/v8OeaDxwN7A+8FngncFLath74ADAjnets4N8lTe3w2GZjcpKvKUnfBGaRJYv1kj4tabakGE0ekm6V9HlJP0n7/Lukl0i6XNJTkn4uaXbDMV8laZGkJ1KL+Li8446I6yLi/wF/bLH5XcA1EXF/RDwHnAm8WdJeHR77qxFxM/BsB7sfCkwF/jkiNkTElwEBh6Vj3RkR3wQeGvsQm5kHnBcRj0bE74HzgPenYz0bEf8VESPpHMNkyX7nDo9tNiYn+ZqKiBOAR4B3pRbtOWPs+l7gBGA3YC/gp8AlZAlmGfBZAEnbAYuAfwN2AY4H/kXSa1odVNK/pNJDq2Vplz+W0tL4GGC/hnXvSv+E7pf0kS7PA/AaYGlsPu/H0rS+2+Pd0/D4nuZjpdflWeAG4BsRsarLc5k9z0neLomIByNiLXAT8GBEfD8iNgHXAAek/d4JLI+ISyJiU0TcBXwLeHerg0bERyNixzGW13YZ63eA4yS9VtI2wD8AAWybtl8NvJqs7PEh4B8kHd/luaYBa5vWrQW2z+l4a4Fpo3V5gPS67AD8NXB7l+cx24yTvK1s+P6ZFo+npe9fARzU2CIH/gZ4aX/ChFRq+SzZP5ffAsuBdcCjafsvI+KxiBiOiJ8AFzDGP6FmTZ21s8jq5Ds07bZDOl+7Y53ecKyvpdXNx9sBWN/0SWG0dHMFcJqk/TuJ3Ww8TvL1lucUo78DftjUIp8WES1LIpK+1pQ4G5f7uw0i1dX3johdyJL9VOC+sXZn8/LOeMdt7Kx9BLgfeG1jS5usw7Rt7BHxjw3H+nBafT9Zp+uo/dscawtgz3G2m3XESb7eVpJforgReKWkEyRtkZbXS3p1q50j4sNNibNxGbOuLWmqpK2BKcAUSVs3dBRvLWm/NBxxFrAAuCAinkzbj5K0U9p+IPAxshEyo8feMh1bwBbpeGP9DdxK1gH6MUlbSTo5rf9BOtZQOtYW2UNt3TDyppXLgE9I2k3ZsM9PApemYx0s6U0pvm0knQrsCtwxzvHMOhMRXmq6kA01fARYA3wKmE3Wup2att8KfLBh/88DlzY8fhvwm4bH+wD/AawmG/3yA2BOzjGfkWJsXM5I23Yk6/x8GvgD8EVgSsNzr0hxrQd+BXys6di3tjj2oePEcgCwhKxsdRdwQMO2Q1sc69ZxjiXgHOCJtJwDKG17C1lH7Lq07YfAmwf9/vFSj2X0TWZmZjXkco2ZWY05yZuZ1ZiTvJlZjTnJm5nVWCUmQJo+fXrMnj170GHU0q+XjD/1yiv/Ys+2+zTuazZoS5YseTwiZvRyjP/+1u3ij0+0n6x1ydIN34uII3s5V9EqkeRnz57N4sWLBx1GLR0xdOy42xctvqbtPo37mg2apN/2eow/PjHMnd+b1Xa/KTMfmN7ruYpWiSRvvek0Sff63OZ9F4046Vs1BTDCyKDDyIWTfA31ktTzNBqHk71VTRBsLO+9dSbESd7MrAW35K10ytKCb9YYl1v1VgVBMFyT2QCc5M3MWhjJdRLXwXGSr4GytuCL4k8GVrQAhnNI8mmm0tuArcjy7bUR8VlJewBXkt2B7S7ghMhuaZk7XwxlZtbCCNF26cAG4LCI2B+YAxwp6WCym7WfHxF7A08CJxb1c7glX3LjtdIXjXQ+hr0sxoq3uUXe6c/V6fHMJiKAjTnU5COb5nd9erhFWoLshvB/ndYvJJti+8KeT9iCk3yJTDRhVy3Bj6dOP4tVXxC5lGsAJE0huy/BnwNfBR4E1kR2H2XIbl+5Wy4na6Gwck26U86dku6RdL+kz6X1l0p6WNLdaZlTVAxmZl0JGO5gAaZLWtywzH/RobJ7Ds8BXg4cSHaz+RZnLEaRLfnRWtR6SVsAt0u6KW37+4i4tsBzV4Jbr2bllF3x2pHHI2JuR8eMWCPpVuBgYEdJU1Nr/uXAY91F2l5hLfnItKpFmZmVnBjuYGl7FGmGpB3T99uQ3VJzGXAL8O602zwa7kWct0JH10iaIuluYBWwKCJGb0z8BUlLJZ0vaasxnjt/9CPQ6tWriwyz744YOtateLMSyzpe1XbpwEzgFklLgZ+T5cEbgVPJbuz+G+AlwEVF/SyFdrxGxDAwJ/0n+7ak/YDPkN2EeUtgAdkP+39aPHdB2s7cuXP9CcA64lE1lodsnHxHSXz840QsJbshfPP6h8jq84Xryzj5iFgD3AocGRErUilnA3AJffpBzcwmYiTUdqmCwlrykmYAG1Nnw2gt6mxJMyNihSQBRwP3FRVD2bhEU5zGFryviLVe5dWSL4MiyzUzgYVpjOgQcHVE3CjpB+kfgIC7gQ8XGIOZ2YQFYrgmEwIUluTHqUUdVtQ5y8yt+GL59bW8VaUc046veDUzaxKI52LKoMPIRT0+j5iNwS1860Z2MdRQ26UK3JI3M2vBHa9mFeEbjNtERYjhqEZLvR0n+T4ZL7E0D/lzicFs8Ebckre8uGXZO7+Glqes47Ue6bEeP4WZWY5GO17rwEneKs0teCvKsMfJW95ci8+fpziwbviKVzOzmhvx6BrLi1vw3fNrZ0XIJihzkrceOUH11+jr7bKNtROIjTWZ1sBJ3sysSQS+GMq65xb8YLlFb+3JF0OZmdVV4Ja8dcEt+HJxi97G445XM7OaCqpzD9d2nOT7xK14s+oIYKPnrjEzqyt5Pnkzs7oKfMWrNXE5xiajOt+QxS15M7OaipBb8ubWe1108nusUwu1U928v8d6TvMdz8r+emYdr71PayBpd+Ay4KXACLAgIi6QdAbwIWB12vX0iPhOzydswUnezOxFcrvH6ybgkxFxl6TtgSWSFqVt50fEuXmcZDxO8mYdaL5wqkqt0k4V+cl0rGO3Wl+G1zPreO29Jh8RK4AV6ft1kpYBu/V84AlwkjebgFZJqcgrZ9sl3l7OOchy43jnLktnbodXvE6XtLjh8YKIWNBqR0mzgQOAO4BDgJMl/S2wmKy1/2RPAY+hsJ4FSVtLulPSPZLul/S5tH4PSXdIekDSVZK2LCoGM7NujF7x2m4BHo+IuQ3LWAl+GvAt4JSIeAq4ENgLmEPW0j+vqJ+lyJb8BuCwiFgvaQvgdkk3AZ8gq0VdKelrwIlkP7BZrfSjpDOR1nhjDFUaNHDE0LEDac3ndSPvlP++BVweEdcBRMTKhu1fB27M5WQtFNaSj8z69HCLtARwGHBtWr8QOLqoGMzMuhEBG0eG2i7tSBJwEbAsIr7UsH5mw27HAPfl/kMkhdbkJU0BlgB/DnwVeBBYExGb0i6PMkYnhKT5wHyAWbNmFRmmWS46qTN30ppuHm6Ylyq13pv1u06flWtyaQMfApwA3Cvp7rTudOB4SXPIGr7LgZPyOFkrhSb5iBgG5kjaEfg28OpWu43x3AXAAoC5c+e23MfMrCh5XPEaEbdDywMVMia+lb6MromINZJuBQ4GdpQ0NbXmXw481u75v17y0ECHWlW5BWTl0sl7ye+39oq+F0BeQyjLoLAkL2kGsDEl+G2AtwFnA7cA7wauBOYB13d7Dv8xmE1uxeUAT2vQiZnAwlSXHwKujogbJf0SuFLS54FfkHVKjOuVf7En3FVgpGZmTXyP1zYiYinZwP/m9Q8BBxZ1XjOzXmWja3qfu6YMKnPFa6vLyc3MiuDb/5mZ1ZzLNWZmNeXRNWZmNefRNQPi2ryZFS1CbHKSH6xWF0E48ZtNXqM5IZsupncu15iZ1ZRr8iXlUo6Z5cVJ3syspjxO3sysRIqYqMzj5M3MaioCNnVwU5AqcJI3s8oqcrpxl2vMzGrKNXkzsxIo8uYh4SRfXs2/cA+pNKuusRJ40X/X7ng1M6upCNfkzcwK1a4EU+w9nsWwR9eUl8szZtYr1+TNzHJUbMt8Yjx3TQm59W5muYmsLl8HlU3yTupmViSPrjEzy0GZyjSjIqeOV0m7A5cBLwVGgAURcYGknYGrgNnAcuC4iHiy5xO2UI/uYzOznEW0XzqwCfhkRLwaOBj4O0n7AqcBN0fE3sDN6XEhnOTNzFqIUNul/TFiRUTclb5fBywDdgOOAham3RYCRxf0Y7hcY2bWLGup51uTlzQbOAC4A9g1IlZk54oVknbJ9WQNnOTNzFrocAjldEmLGx4viIgFzTtJmgZ8CzglIp7K6z60nSgsyY/T4XAG8CFgddr19Ij4TlFxmJl1o8Oa++MRMXe8HSRtQZbgL4+I69LqlZJmplb8TGBVT8GOo8iW/GiHw12StgeWSFqUtp0fEecWeO7cLBq5ZrPhmq1GAng4p9nElXFUzahAjOQzukbARcCyiPhSw6YbgHnAWenr9T2fbAyFJflUbxqtOa2TNNrhkLuJvFnaJexujt/8j8DMqi+na6EOAU4A7pV0d1p3Ollyv1rSicAjQGEJpC81+aYOh0OAkyX9LbCYrLVfyPhQM7Ou5NTxGhG3w5hXVR3e8wk6UHiSb9HhcCFwJtk/yjOB84APtHjefGA+wKxZs1503G4/6pX5I6KZlUhNpjUodJx8qw6HiFgZEcMRMQJ8HTiw1XMjYkFEzI2IuTNmzCgyTDOzF8ljnHwZFDm6pmWHw2iPcnp4DHBfUTH006DuXmNWRWX/RB3AyEg1kng7RZZrxupwOF7SHLLXcTlwUoExDIyTu1mFBVCRlno7RY6uGavDwWPizaz0PNWwmU0KU7bfnnhuIwAjzz2XrYyRro5V9jLNZpzkzczqqjodq+04yRdktMXi2rxVzZSdd8q+GUlN2eFhGMoSntLXGB7/GI0t9sr+Dbglb2ZWUwHh0TXWCc91Y1UjZZfPxFBWdx9es+75bUNbbpltG23ld1Cbr1QdfjNO8mZm9eVyjXXL9XortdH6e4tZGJ8fXaNJcFM5J3nrlZO9ldGm1Y+336nLIZSV4YuhzMzqzRdDWW6aO6bGatl73noru+p2srbg0TVmZvUlt+StKLVqDZlVUeCOV8tfq1sTujxjNghyx6sVy8nd+mm8G9ZP2veiW/JWNe64tWaNybzXMmHtyow1GSXa0RUNks7uZJ2ZWS2MjpNvt1RAp5etHdFi3dvzDMSyllDtWkNWiE7eJ837tHvcSyx1fN8q2i99iUMaknRct88fN8lL+oike4F9JC1tWB4GlnZ7UjOz0osOln6EETECnNzt89vV5P8NuAn4InBaw/p1EfFEtye18RVZO/eonerqpiWeR+t9rPdjHVvvJbZI0qeAq4CnR1d2kofHTfIRsRZYCxwPIGkXYGtgmqRpEfFIL1Fb/zm5W15G30t1TfYluxjqA2SfHT7atH7Pdk/saHSNpHcBXwJeBqwCXgEsA14zoTCtY25xW1mT56QYpRWUbVqDfckS/JvIovsR8LVOntjpEMrPAwcD34+IAyS9ldS6t2LllezLmjDsxarwu6pCjD3LqSUv6WLgncCqiNgvrTsD+BCwOu12ekR8Z5zDLASeAr6cHh+f1rXtkO00yW+MiD+mXt6hiLjFQyjNilH3MkhV5FiuuRT4CnBZ0/rzI+LcDo+xT0Ts3/D4Fkn3dPLETpP8GknTgNuAyyWtAjZ1+FzLQaczVY61f/P62n/ctlxNyn84OSX5iLhN0uweD/MLSQdHxM8AJB0E/LiTJ3Y6Tv4o4Bng48B3gQeBd3URqJlZNXQ2hHK6pMUNy/wJnOHkNCT9Ykk7tdn3IOAnkpZLWg78FHiLpHsljTucvaOWfEQ83fBwYSfPsWLldQm6W/RmLzaBi50ej4i5XZziQuBMsn8VZwLnkY2gGcuRXZwDaJPkJa2j9YcWARERO4zz3N3JalAvJZsFYkFEXCBpZ7KxnrOB5cBxEfFkV9Fb15zczdoocHRNRKwc/V7S14Eb2+z/227P1W6c/PbdHpisZv/JiLhL0vbAEkmLgPcDN0fEWZJOI7vI6tQezmNWK5Oy/l1CRY6TlzQzIlakh8cA9xV1rsJmoUw/wIr0/TpJy4DdyOr7h6bdFgK34iRv5uReNvkNobyCLOdNl/Qo8FngUElz0lmWAyflc7YX68tUw6ln+QDgDmDX0f9gEbEiXUVrZlYeOU5AFhGtrim6KJ+jt1d4kk9DL78FnBIRT0md1blSL/V8gFmzZhUXoNmAuQVfUuWa1qBrnQ6h7IqkLcgS/OURcV1avVLSzLR9Jtk0CS8SEQsiYm5EzJ0xY0aRYZqZvYhG2i9VUFhLXlmT/SJgWUR8qWHTDcA84Kz09fqiYjArG7fard+KLNccApwA3Cvp7rTudLLkfrWkE4FHAI/lM7PyqUm5psjRNbeTjadv5fCizmtWJm65V1Qf7/xUNN/I26wgTvAV5yRvZlZjTvJWZZ67pjhuwVefqM7omXac5M3Mmrkmb2ZWc07yZmY15iRvZo1ci68Xl2vMJjkn9Zpzkjczq6nw6BqzScst+EnCLXkzs/pyTd5qwRdFTYxb8ZOIk7yZWU0FTvJmk4Fb7pOTcLnGrNac3M1J3syszpzkzcxqzEnezKymPAulWT25Fm/Pc5I3M6svT2tgVhNuvVsrLtdYrUymK1+d1K2tHC+GknQx8E5gVUTsl9btDFwFzAaWA8dFxJP5nHFzQ0Uc1Mys8qKDpTOXAkc2rTsNuDki9gZuTo8L4Za8TRpuwVun8rziNSJukzS7afVRwKHp+4XArcCp+Zxxc07yZmYtaKSjLD9d0uKGxwsiYkEHz9s1IlYARMQKSbt0E2MnnORtM4tGrqlMXV5TpgAwtP32AAyvWfP8NrfarSedl2Mej4i5xQbTG9fkzcxaULRferBS0kyA9HVVHjG3UliSl3SxpFWS7mtYd4ak30u6Oy3vKOr8NokMKVvM8pRfx2srNwDz0vfzgOt7Oto4iizXXAp8Bbisaf35EXFugee1HlVlOGUMDwMw/OTabIWGWDR81QAjsjrJq+NV0hVknazTJT0KfBY4C7ha0onAI0Bhf2yFJfkxepTNzKohv9E1x4+x6fB8zjC+QdTkT5a0NJVzdhprJ0nzJS2WtHj16tX9jM+qJkZeWMzyENm0Bu2WKuh3kr8Q2AuYA6wAzhtrx4hYEBFzI2LujBkz+hWfVdwRQ8eWvsxk5Tc6Tr7Ajte+6esQyohYOfq9pK8DN/bz/DYxVanNmxUiKpLF2+hrS350yFByDHDfWPuamQ2SW/JtjNGjfKikOWRdGsuBk4o6v+XHLXqbdHKcoGzQihxd06pH+aKizmdmlqeqdKy242kNzMxacJI3M6uroDYdr07yZmYtVKVjtR0neeuYO2BtUnGSN7PJrvEffp2md87zpiGD5iRvZtYsotObhpSek7yZWSv1yPFO8mZmrbhcY2ZWVwG4XGNmVmP1yPG+x6tNXJ1GUZiNxROUmZnVmEfXmJWQP2UMzhFDx9bn9fcslDbZ+epXg/q+D7KLoeqR5Z3kzWzC6nql62Y8C6VZedQ20VREHV9/t+TNzOrKNXmzTF1rsja+OrbcN+e5a8xsEpsUNXmXa8xe4Ba91Urkd/s/ScuBdcAwsCki5uZz5M44yZvZhNW29d4o35b8WyPi8TwP2CkneTOzVupRrfHcNZavRSPXTI5WntWeRkbaLsB0SYsblvktDhXAf0paMsb2Qrklb2bWLOj0YqjHO6ixHxIRj0naBVgk6VcRcVuvIXbKLXkzsyYiULRfOhERj6Wvq4BvAwcWGPqLOMmbmbUS0X5pQ9J2krYf/R74S+C+giPfTGFJXtLFklZJuq9h3c6SFkl6IH3dqajzm5n1JIckD+wK3C7pHuBO4D8i4ruFxt2kyJr8pcBXgMsa1p0G3BwRZ0k6LT0+tcAYbECKHjfvzl0rVOc1+fEPE/EQsH/vR+peYS351LHwRNPqo4CF6fuFwNFFnd/MrBcdjq4pvX7X5HeNiBUA6esuY+0oaf7osKTVq1f3LUArP7firXgdlGoqMu1BaYdQRsQCYAHA3Llzq/Fq2kC1Kg35H4J1JahMEm+n3y35lZJmAqSvq/p8fjOzzox0sFRAv5P8DcC89P084Po+n98mGU+YZt3Ka5z8oBU5hPIK4KfAPpIelXQicBZwhKQHgCPSYzOz8nFNfnwRcfwYmw4v6pxmrcz+v+cCsPykTw04EquMCBiuSD2mjdJ2vFo9LBq5pueSSS+dpw9ceFBP57ZJrCIt9Xac5K208hgZo+025hCJTUpO8mbVsse/fhGAh9/3mQFHYqUXgO/xataZQd4a8M92+hMAa5/ctu/ntioLCNfkzUphov88zvnl2wH49L43FRGO1UHgjlezieq0RZ/nVapv3e0BALbd/TkA/mzqM7kd22rONXkzsxpzkjfrTqsWfVFzzBw87UEA/ttWjwGwMXyfHOtEdS52asdJ3gYm78T+0D+9AYCRXTcA8Ma9H2LnKesBeMXU7K1+49MvTHw60Em+rdwCqMhUwu04ydukNfpJwjNVWktuyZuVy9BzAmB4U1aSWfXMNHYYehaAbbQVwPMt+3P22m8AEVp1eFoDM7P6CgiPkzcrhz984o0ATM2ue+K51JLfdupGth3KpjV4Jv3Bzt4ie3zMssf7HKVVjq94NTMb2xFv+sLz3y+6/X8PMJIuuSZvVg5Kn6pHR0dqQ/bNlkPD/OhPewMwZdtfA7BmZGsAthva0N8grVoiPLrGrCxGk7xSwyu2HgZg9TPbseK5HQH4mfYA4OENMwD41bqXPv/8eXv3KdBJJqZo0CH0xi15M7O6CmJ4eNBB5MJJ3irvnn/+OAD7nHk+ANqYlWs2DE/lR6v2AuB27QnAqnXbA/D02m1eOMAb+xXp5PL9H54+6BC656mGzcxqzkMozcpJG7Na8ONrp7FpwxQA4k/ZW33KU9njKcCDn/rEQOKz8gsgcmrJSzoSuIDsbfeNiDgrlwN3yLM1mZk1i3TTkHZLG5KmAF8F3g7sCxwvad+Co9+MW/JWGzE6mCN9jYCRZ7O3+JZPvPBWf+C0j/c5MquinDpeDwR+ExEPAUi6EjgK+GUeB++EogLDhCStBn7b59NOB8p2WWTZYnI87ZUtprLFA/nH9IqImNHLASR9lyyudrYGnm14vCAiFjQc593AkRHxwfT4BOCgiDi5l/gmohIt+V5/Yd2QtDgi5vb7vOMpW0yOp72yxVS2eKCcMUXEkTkdqtXFAn1tWbsmb2ZWnEeB3Rsevxx4rJ8BOMmbmRXn58DekvaQtCXwXuCGfgZQiXLNgCxov0vflS0mx9Ne2WIqWzxQzphyERGbJJ0MfI9sCOXFEXF/P2OoRMermZl1x+UaM7Mac5I3M6sxJ3lA0o6SrpX0K0nLJL1B0hmSfi/p7rS8o4/x7NNw3rslPSXpFEk7S1ok6YH0dacBxzOw1yjF9XFJ90u6T9IVkrZOHVx3pNfoqtTZNch4LpX0cMNrNKdf8aSY/leK535Jp6R1A3kfjRPPQN9HdeeaPCBpIfCjiPhGSgrbAqcA6yPi3AHHNgX4PXAQ8HfAExFxlqTTgJ0i4tQBxvM/GdBrJGk34HZg34h4RtLVwHeAdwDXRcSVkr4G3BMRFw4wnkOBGyPi2qJjaBHTfsCVZFddPgd8F/gI8CEG8D4aJ56/oQR/a3U16VvyknYA3gxcBBARz0XEmsFGtZnDgQcj4rdkl0MvTOsXAkcPOJ5BmwpsI2kq2T/mFcBhwGhC7fdr1BxPX8dDt/Bq4GcR8aeI2AT8EDiGwb2PxorHCjTpkzywJ7AauETSLyR9Q9J2advJkpZKurifH2mbvBe4In2/a0SsAEhfdxlwPDCg1ygifg+cCzxCltzXAkuANSmBQHYhym6Diici/jNt/kJ6jc6XtFU/4knuA94s6SWStiX7lLM7g3sfjRUPlONvrZac5LPW1+uACyPiAOBp4DTgQmAvYA7ZH+15/Q4slY7+B3BNv8/dSot4BvYapURwFLAH8DJgO7KZ/pr1pR7ZKh5J7wM+A7wKeD2wM9C38lpELAPOBhaRlUbuATaN+6TBxDPwv7U6c5LPWnuPRsQd6fG1wOsiYmVEDEfECPB1sjpiv70duCsiVqbHKyXNBEhfVw0yngG/Rm8DHo6I1RGxEbiO7B5PO6ZyCfT3EvKW8UTEishsAC6hz++jiLgoIl4XEW8GngAeYIDvo1bxlORvrbYmfZKPiD8Av5O0T1p1OPDL0T+C5Biyj5r9djybl0ZuAOal7+cB1w8yngG/Ro8AB0vaVpJIvzfgFuDdaZ9+vkat4lnWkExFVvvu6/tI0i7p6yzgr8h+fwN7H7WKpyR/a7Xl0TVAGtb2DWBL4CGyUSNfJvv4GMBy4KTROmafYtoW+B2wZ0SsTeteAlwNzCJLKsdGxBMDjOebDPY1+hzwHrKP/L8APkhWg7+SrDTyC+B9qRU9qHhuAmaQzUZ4N/DhiFjfj3hSTD8CXgJsBD4RETcP+H3UKp6Bvo/qzknezKzGJn25xsyszpzkzcxqzEnezKzGnOTNzGrMSd7MrMac5K3UJPVtuKFZHTnJm5nVmJO8VYIy/5TmIr9X0nvS+kMl3aoX7nATtXAAAAC9SURBVAdwebq61MzwjbytOv6K7KrI/YHpwM8l3Za2HQC8hmyemh8Dh5DN7W426bklb1XxJuCKNJHVSrK5yF+ftt0ZEY+mCa7uBmYPKEaz0nGSt6oYrwTTODfNMP6EavY8J3mrituA90iaImkG2d287hxwTGal5xaPVcW3gTeQ3WgigE9HxB8kvWqwYZmVm2ehNDOrMZdrzMxqzEnezKzGnOTNzGrMSd7MrMac5M3MasxJ3sysxpzkzcxq7P8DwEDSLiTjLi0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "if variable == \"pr\":\n",
    "    observed = (obs_path +\"final_kelvin.nc\")\n",
    "    prcp_hist = model_path +var_name+ \"*historical*.nc\"\n",
    "    file_lists = glob.glob(prcp_hist)\n",
    "    try:\n",
    "        model_hist = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_hist)\n",
    "    except:\n",
    "        model_hist = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_hist2)\n",
    "    \n",
    "    observed = xr.open_dataset(observed)\n",
    "    try:\n",
    "        observed = observed.rename({'prcp':'pr'})\n",
    "    except:\n",
    "        pass\n",
    "    observed = observed['pr']\n",
    "    observed.isel(time = 1).plot()\n",
    "    model_hist.pr.values = model_hist.pr.values*86400\n",
    "    model_hist = model_hist['pr']\n",
    "    lats = observed.lat.values\n",
    "    lons = observed.lon.values\n",
    "    prcp_245 = model_path +var_name+ \"*245*.nc\"\n",
    "    prcp_126 = model_path +var_name+ \"*126*.nc\"\n",
    "    prcp_585 = model_path +var_name+ \"*585*.nc\"\n",
    "    prcp_370 = model_path +var_name+ \"*370*.nc\"\n",
    "#     print(\"reading 45\")\n",
    "    print(prcp_245)\n",
    "    print(prcp_126)\n",
    "    file_lists = glob.glob(prcp_245)\n",
    "    print(file_lists)\n",
    "    model_245 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "\n",
    "    file_lists = glob.glob(prcp_126)\n",
    "    model_126 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "    \n",
    "    file_lists = glob.glob(prcp_370)\n",
    "    model_370 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "    \n",
    "    file_lists = glob.glob(prcp_585)\n",
    "    model_585 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "    \n",
    "    \n",
    "    \n",
    "    model_245.pr.values = model_245.pr.values*86400\n",
    "    model_585.pr.values = model_585.pr.values*86400\n",
    "    model_126.pr.values = model_126.pr.values*86400\n",
    "    model_370.pr.values = model_370.pr.values*86400\n",
    "    \n",
    "    \n",
    "    model_245 = model_245['pr']\n",
    "    model_585 = model_585['pr']\n",
    "    model_126 = model_126['pr']\n",
    "    model_370 = model_370['pr']\n",
    "    model_245.close()\n",
    "    model_585.close()\n",
    "    model_126.close()\n",
    "    model_370.close()\n",
    "    model_hist.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T02:29:14.314617Z",
     "start_time": "2020-05-16T02:29:14.306894Z"
    }
   },
   "outputs": [],
   "source": [
    "if variable == \"tasmax\" or variable == \"tasmin\":\n",
    "    observed = (obs_path +\"final_kelvin.nc\")\n",
    "    observed = xr.open_dataset(observed)\n",
    "    observed.close()\n",
    "    observed = observed[variable]\n",
    "    temp_hist = model_path +var_name+ \"*historical*.nc\"\n",
    "    print(temp_hist)\n",
    "    file_lists = glob.glob(temp_hist)\n",
    "    print(file_lists)\n",
    "    print(\"hello\")\n",
    "    try:\n",
    "        model_hist = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_hist)\n",
    "    except:\n",
    "        model_hist = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_hist2)\n",
    "    \n",
    "    model_hist[variable].values = model_hist[variable].values\n",
    "    model_hist = model_hist[variable]\n",
    "    lats = observed.lat.values\n",
    "    lons = observed.lon.values\n",
    "    temp_245 = model_path +var_name+ \"*245*.nc\"\n",
    "    temp_126 = model_path +var_name+ \"*126*.nc\"\n",
    "    temp_585 = model_path +var_name+ \"*585*.nc\"\n",
    "    temp_370 = model_path +var_name+ \"*370*.nc\"\n",
    "#     print(\"reading 45\")\n",
    "    file_lists = glob.glob(temp_245)\n",
    "    print(file_lists)\n",
    "    model_245 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "\n",
    "    file_lists = glob.glob(temp_126)\n",
    "    model_126 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "    \n",
    "    file_lists = glob.glob(temp_370)\n",
    "    model_370 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "    \n",
    "    file_lists = glob.glob(temp_585)\n",
    "    model_585 = xr.open_mfdataset( file_lists, concat_dim='time', preprocess=preproc_model)\n",
    "    \n",
    "    \n",
    "    \n",
    "    model_245[variable].values = model_245[variable].values\n",
    "    model_585[variable].values = model_585[variable].values\n",
    "    model_126[variable].values = model_126[variable].values\n",
    "    model_370[variable].values = model_370[variable].values\n",
    "    model_245 = model_245[variable]\n",
    "    model_585 = model_585[variable]\n",
    "    model_126 = model_126[variable]\n",
    "    model_370 = model_370[variable]\n",
    "    model_245.close()\n",
    "    model_585.close()\n",
    "    model_126.close()\n",
    "    model_370.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T02:30:51.949342Z",
     "start_time": "2020-05-16T02:29:14.316505Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reading hist\n",
      "reading 45\n",
      "reading 26\n",
      "reading 85\n",
      "reading 70\n"
     ]
    }
   ],
   "source": [
    "print(\"reading hist\")\n",
    "model_hist = model_hist.interp(lat = lats, lon = lons)\n",
    "print(\"reading 45\")\n",
    "\n",
    "model_245 = model_245.interp(lat = lats, lon = lons)\n",
    "print(\"reading 26\")\n",
    "\n",
    "model_126 = model_126.interp(lat = lats, lon = lons)\n",
    "print(\"reading 85\")\n",
    "\n",
    "model_585 = model_585.interp(lat = lats, lon = lons)\n",
    "print(\"reading 70\")\n",
    "model_370 = model_370.interp(lat = lats, lon = lons)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T02:30:51.954664Z",
     "start_time": "2020-05-16T02:30:51.951970Z"
    }
   },
   "outputs": [],
   "source": [
    "# obs_test = observed.isel(lat = 56, lon  = 56)\n",
    "# hist_test = model_hist.isel(lat = 56,lon = 56)\n",
    "# test_245 = model_245.isel(lat = 56,lon = 56)\n",
    "\n",
    "# #test_245 = bias_correction_model(y_245,q_map)\n",
    "# temp = bias_correction(obs_test.values,hist_test.values)\n",
    "# q_map = temp['params']\n",
    "# #test_245_corr = bias_correction_model(test_245.values,q_map)\n",
    "\n",
    "# hist_corr = temp['outputs']\n",
    "\n",
    "# # plt.plot(obs_test[0:365])\n",
    "# # plt.plot(hist_test[0:365])\n",
    "# # plt.plot(hist_corr[0:365])\n",
    "# plt.plot(test_245[0:365])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T02:31:01.842931Z",
     "start_time": "2020-05-16T02:30:51.956648Z"
    }
   },
   "outputs": [],
   "source": [
    "bias_corrected_results_hist = np.zeros([len(model_hist.time.values),\n",
    "                                           len(model_hist.lat.values),\n",
    "                                           len(model_hist.lon.values)])\n",
    "bias_corrected_results_hist[:] = np.nan\n",
    "\n",
    "\n",
    "bias_corrected_results_585 = np.zeros([len(model_585.time.values),\n",
    "                                           len(model_585.lat.values),\n",
    "                                           len(model_585.lon.values)])\n",
    "bias_corrected_results_585[:] = np.nan\n",
    "\n",
    "\n",
    "\n",
    "bias_corrected_results_126 = np.zeros([len(model_126.time.values),\n",
    "                                           len(model_126.lat.values),\n",
    "                                           len(model_126.lon.values)])\n",
    "bias_corrected_results_126[:] = np.nan\n",
    "\n",
    "bias_corrected_results_245 = np.zeros([len(model_245.time.values),\n",
    "                                           len(model_245.lat.values),\n",
    "                                           len(model_245.lon.values)])\n",
    "bias_corrected_results_245[:] = np.nan\n",
    "\n",
    "bias_corrected_results_370 = np.zeros([len(model_370.time.values),\n",
    "                                           len(model_370.lat.values),\n",
    "                                           len(model_370.lon.values)])\n",
    "bias_corrected_results_370[:] = np.nan\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T02:31:09.456463Z",
     "start_time": "2020-05-16T02:31:01.844674Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(124, 147, 24106)\n"
     ]
    }
   ],
   "source": [
    "model_hist_values = model_hist.values\n",
    "hist_dict = {}\n",
    "hist_dict['time'] = model_hist.time.values\n",
    "hist_dict['lon'] =  model_hist.lon.values\n",
    "hist_dict['lat'] =  model_hist.lat.values\n",
    "\n",
    "model245_values = model_245.values\n",
    "rcp245_dict = {}\n",
    "rcp245_dict['time'] = model_245.time.values\n",
    "rcp245_dict['lon'] =  model_245.lon.values\n",
    "rcp245_dict['lat'] =  model_245.lat.values\n",
    "\n",
    "model585_values = model_585.values\n",
    "\n",
    "rcp585_dict = {}\n",
    "rcp585_dict['time'] = model_585.time.values\n",
    "rcp585_dict['lon'] =  model_585.lon.values\n",
    "rcp585_dict['lat'] =  model_585.lat.values\n",
    "\n",
    "\n",
    "\n",
    "model126_values = model_126.values\n",
    "\n",
    "rcp126_dict = {}\n",
    "rcp126_dict['time'] = model_126.time.values\n",
    "rcp126_dict['lon'] =  model_126.lon.values\n",
    "rcp126_dict['lat'] =  model_126.lat.values\n",
    "\n",
    "model370_values = model_370.values\n",
    "\n",
    "rcp370_dict = {}\n",
    "rcp370_dict['time'] = model_370.time.values\n",
    "rcp370_dict['lon'] =  model_370.lon.values\n",
    "rcp370_dict['lat'] =  model_370.lat.values\n",
    "\n",
    "observation_attr_values = observed.values\n",
    "print(observation_attr_values.shape)\n",
    "#observation_attr_values = observation_attr_values.reshape((-1,len(lats),len(lons)))\n",
    "#print(observation_attr_values.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T02:31:09.461587Z",
     "start_time": "2020-05-16T02:31:09.458178Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(23360, 124, 147)\n",
      "(124, 147, 24106)\n",
      "(31390, 124, 147)\n"
     ]
    }
   ],
   "source": [
    "print(model_hist_values.shape)\n",
    "print(observation_attr_values.shape)\n",
    "print(model126_values.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T03:47:32.532077Z",
     "start_time": "2020-05-16T02:31:09.462880Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.125 81.125\n",
      "7.375 79.875\n",
      "7.375 81.125\n",
      "8.625 77.375\n",
      "8.625 81.125\n",
      "9.875 77.375\n",
      "9.875 78.625\n",
      "11.125 76.125\n",
      "11.125 77.375\n",
      "11.125 78.625\n",
      "12.375 76.125\n",
      "12.375 77.375\n",
      "12.375 78.625\n",
      "12.375 79.875\n",
      "13.625 74.875\n",
      "13.625 76.125\n",
      "13.625 77.375\n",
      "13.625 78.625\n",
      "13.625 79.875\n",
      "14.875 74.875\n",
      "14.875 76.125\n",
      "14.875 77.375\n",
      "14.875 78.625\n",
      "14.875 79.875\n",
      "16.125 73.625\n",
      "16.125 74.875\n",
      "16.125 76.125\n",
      "16.125 77.375\n",
      "16.125 78.625\n",
      "16.125 79.875\n",
      "16.125 81.125\n",
      "17.375 73.625\n",
      "17.375 74.875\n",
      "17.375 76.125\n",
      "17.375 77.375\n",
      "17.375 78.625\n",
      "17.375 79.875\n",
      "17.375 81.125\n",
      "17.375 82.375\n",
      "18.625 73.625\n",
      "18.625 74.875\n",
      "18.625 76.125\n",
      "18.625 77.375\n",
      "18.625 78.625\n",
      "18.625 79.875\n",
      "18.625 81.125\n",
      "18.625 82.375\n",
      "18.625 83.625\n",
      "19.875 73.625\n",
      "19.875 74.875\n",
      "19.875 76.125\n",
      "19.875 77.375\n",
      "19.875 78.625\n",
      "19.875 79.875\n",
      "19.875 81.125\n",
      "19.875 82.375\n",
      "19.875 83.625\n",
      "19.875 84.875\n",
      "19.875 86.125\n",
      "21.125 71.125\n",
      "21.125 73.625\n",
      "21.125 74.875\n",
      "21.125 76.125\n",
      "21.125 77.375\n",
      "21.125 78.625\n",
      "21.125 79.875\n",
      "21.125 81.125\n",
      "21.125 82.375\n",
      "21.125 83.625\n",
      "21.125 84.875\n",
      "21.125 86.125\n",
      "22.375 69.875\n",
      "22.375 71.125\n",
      "22.375 72.375\n",
      "22.375 73.625\n",
      "22.375 74.875\n",
      "22.375 76.125\n",
      "22.375 77.375\n",
      "22.375 78.625\n",
      "22.375 79.875\n",
      "22.375 81.125\n",
      "22.375 82.375\n",
      "22.375 83.625\n",
      "22.375 84.875\n",
      "22.375 86.125\n",
      "22.375 87.375\n",
      "22.375 88.625\n",
      "22.375 89.875\n",
      "22.375 91.125\n",
      "22.375 92.375\n",
      "23.625 68.625\n",
      "23.625 69.875\n",
      "23.625 71.125\n",
      "23.625 72.375\n",
      "23.625 73.625\n",
      "23.625 74.875\n",
      "23.625 76.125\n",
      "23.625 77.375\n",
      "23.625 78.625\n",
      "23.625 79.875\n",
      "23.625 81.125\n",
      "23.625 82.375\n",
      "23.625 83.625\n",
      "23.625 84.875\n",
      "23.625 86.125\n",
      "23.625 87.375\n",
      "23.625 88.625\n",
      "23.625 89.875\n",
      "23.625 91.125\n",
      "23.625 92.375\n",
      "24.875 67.375\n",
      "24.875 68.625\n",
      "24.875 69.875\n",
      "24.875 71.125\n",
      "24.875 72.375\n",
      "24.875 73.625\n",
      "24.875 74.875\n",
      "24.875 76.125\n",
      "24.875 77.375\n",
      "24.875 78.625\n",
      "24.875 79.875\n",
      "24.875 81.125\n",
      "24.875 82.375\n",
      "24.875 83.625\n",
      "24.875 84.875\n",
      "24.875 86.125\n",
      "24.875 87.375\n",
      "24.875 88.625\n",
      "24.875 89.875\n",
      "24.875 91.125\n",
      "24.875 92.375\n",
      "24.875 93.625\n",
      "26.125 62.375\n",
      "26.125 63.625\n",
      "26.125 64.875\n",
      "26.125 66.125\n",
      "26.125 67.375\n",
      "26.125 68.625\n",
      "26.125 69.875\n",
      "26.125 71.125\n",
      "26.125 72.375\n",
      "26.125 73.625\n",
      "26.125 74.875\n",
      "26.125 76.125\n",
      "26.125 77.375\n",
      "26.125 78.625\n",
      "26.125 79.875\n",
      "26.125 81.125\n",
      "26.125 82.375\n",
      "26.125 83.625\n",
      "26.125 84.875\n",
      "26.125 86.125\n",
      "26.125 87.375\n",
      "26.125 88.625\n",
      "26.125 89.875\n",
      "26.125 91.125\n",
      "26.125 92.375\n",
      "26.125 93.625\n",
      "26.125 94.875\n",
      "27.375 63.625\n",
      "27.375 64.875\n",
      "27.375 66.125\n",
      "27.375 67.375\n",
      "27.375 68.625\n",
      "27.375 69.875\n",
      "27.375 71.125\n",
      "27.375 72.375\n",
      "27.375 73.625\n",
      "27.375 74.875\n",
      "27.375 76.125\n",
      "27.375 77.375\n",
      "27.375 78.625\n",
      "27.375 79.875\n",
      "27.375 81.125\n",
      "27.375 82.375\n",
      "27.375 83.625\n",
      "27.375 84.875\n",
      "27.375 86.125\n",
      "27.375 87.375\n",
      "27.375 88.625\n",
      "27.375 89.875\n",
      "27.375 91.125\n",
      "27.375 92.375\n",
      "27.375 93.625\n",
      "27.375 94.875\n",
      "27.375 96.125\n",
      "28.625 62.375\n",
      "28.625 63.625\n",
      "28.625 64.875\n",
      "28.625 66.125\n",
      "28.625 67.375\n",
      "28.625 68.625\n",
      "28.625 69.875\n",
      "28.625 71.125\n",
      "28.625 72.375\n",
      "28.625 73.625\n",
      "28.625 74.875\n",
      "28.625 76.125\n",
      "28.625 77.375\n",
      "28.625 78.625\n",
      "28.625 79.875\n",
      "28.625 81.125\n",
      "28.625 82.375\n",
      "28.625 83.625\n",
      "28.625 84.875\n",
      "28.625 86.125\n",
      "28.625 87.375\n",
      "28.625 88.625\n",
      "28.625 89.875\n",
      "28.625 91.125\n",
      "28.625 92.375\n",
      "28.625 93.625\n",
      "28.625 94.875\n",
      "28.625 96.125\n",
      "28.625 97.375\n",
      "29.875 67.375\n",
      "29.875 68.625\n",
      "29.875 69.875\n",
      "29.875 71.125\n",
      "29.875 72.375\n",
      "29.875 73.625\n",
      "29.875 74.875\n",
      "29.875 76.125\n",
      "29.875 77.375\n",
      "29.875 78.625\n",
      "29.875 79.875\n",
      "29.875 81.125\n",
      "29.875 82.375\n",
      "29.875 83.625\n",
      "29.875 84.875\n",
      "29.875 86.125\n",
      "29.875 87.375\n",
      "29.875 88.625\n",
      "29.875 89.875\n",
      "29.875 91.125\n",
      "29.875 92.375\n",
      "29.875 93.625\n",
      "29.875 94.875\n",
      "29.875 96.125\n",
      "31.125 67.375\n",
      "31.125 68.625\n",
      "31.125 69.875\n",
      "31.125 71.125\n",
      "31.125 72.375\n",
      "31.125 73.625\n",
      "31.125 74.875\n",
      "31.125 76.125\n",
      "31.125 77.375\n",
      "31.125 78.625\n",
      "31.125 79.875\n",
      "31.125 81.125\n",
      "32.375 68.625\n",
      "32.375 69.875\n",
      "32.375 71.125\n",
      "32.375 72.375\n",
      "32.375 73.625\n",
      "32.375 74.875\n",
      "32.375 76.125\n",
      "32.375 77.375\n",
      "32.375 78.625\n",
      "32.375 79.875\n",
      "32.375 81.125\n",
      "33.625 69.875\n",
      "33.625 71.125\n",
      "33.625 72.375\n",
      "33.625 73.625\n",
      "33.625 74.875\n",
      "33.625 76.125\n",
      "33.625 77.375\n",
      "33.625 78.625\n",
      "33.625 79.875\n",
      "33.625 81.125\n",
      "34.875 68.625\n",
      "34.875 69.875\n",
      "34.875 71.125\n",
      "34.875 72.375\n",
      "34.875 73.625\n",
      "34.875 74.875\n",
      "34.875 76.125\n",
      "34.875 77.375\n",
      "34.875 78.625\n",
      "34.875 79.875\n",
      "36.125 72.375\n",
      "36.125 73.625\n",
      "36.125 74.875\n",
      "36.125 76.125\n"
     ]
    }
   ],
   "source": [
    "correct_params = []\n",
    "for i,lat in enumerate(lats):\n",
    "    for j,lon in enumerate(lons):\n",
    "        params_dict = {}\n",
    "\n",
    "        if np.isnan(model_hist_values[0,i,j]) or np.isnan(observation_attr_values[i,j,0]):\n",
    "            bias_corrected_results_hist[:,i,j] = np.nan\n",
    "            params_dict['lat'] = lat\n",
    "            params_dict['lon'] = lon\n",
    "            params_dict['params'] = np.nan\n",
    "        else:\n",
    "            try:\n",
    "\n",
    "                y = model_hist_values[:,i,j]\n",
    "                x = observation_attr_values[i,j,:]\n",
    "                \n",
    "                \n",
    "                y_245 = model245_values[:,i,j]\n",
    "                y_585 = model585_values[:,i,j]\n",
    "                y_126 = model126_values[:,i,j]\n",
    "                y_370 = model370_values[:,i,j]\n",
    "                \n",
    "                temp = bias_correction(x,y)\n",
    "                q_map = temp['params']\n",
    "                temp_245 = bias_correction_model(y_245,q_map)\n",
    "                temp_585 = bias_correction_model(y_585,q_map)\n",
    "                temp_126 = bias_correction_model(y_126,q_map)\n",
    "                temp_370 = bias_correction_model(y_370,q_map)\n",
    "                bias_corrected_results_245[:,i,j] = temp_245['outputs']\n",
    "                bias_corrected_results_585[:,i,j] = temp_585['outputs']\n",
    "                bias_corrected_results_126[:,i,j] = temp_126['outputs']\n",
    "                bias_corrected_results_370[:,i,j] = temp_370['outputs']\n",
    "                bias_corrected_results_hist[:,i,j] = temp['outputs']\n",
    "\n",
    "                if i%5==0 and j%5==0:\n",
    "                    print(lat,lon)\n",
    "                params_dict['lat'] = lat\n",
    "                params_dict['lon'] = lon\n",
    "                params_dict['params'] = temp['params']\n",
    "\n",
    "            except:\n",
    "                bias_corrected_results_hist[:,i,j] = np.nan\n",
    "                bias_corrected_results_245[:,i,j] = np.nan\n",
    "                bias_corrected_results_126[:,i,j] = np.nan\n",
    "                bias_corrected_results_585[:,i,j] = np.nan\n",
    "                bias_corrected_results_370[:,i,j] = np.nan\n",
    "                params_dict['lat'] = lat\n",
    "                params_dict['lon'] = lon\n",
    "                params_dict['params'] = np.nan\n",
    "\n",
    "        correct_params.append(params_dict)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T03:47:32.711795Z",
     "start_time": "2020-05-16T03:47:32.534226Z"
    }
   },
   "outputs": [],
   "source": [
    "ds_hist = xr.Dataset({variable: (('time', 'lat','lon'), bias_corrected_results_hist)},\n",
    "                        coords={'lat': lats,\n",
    "                                'lon': lons,\n",
    "                                'time':hist_dict['time'] })\n",
    "\n",
    "\n",
    "ds_ssp585 = xr.Dataset({variable: (('time', 'lat','lon'), bias_corrected_results_585)},\n",
    "                        coords={'lat': lats,\n",
    "                                'lon': lons,\n",
    "                                'time':rcp585_dict['time'] })\n",
    "\n",
    "ds_ssp126 = xr.Dataset({variable: (('time', 'lat','lon'), bias_corrected_results_126)},\n",
    "                        coords={'lat': lats,\n",
    "                                'lon': lons,\n",
    "                                'time':rcp126_dict['time'] })\n",
    "\n",
    "ds_ssp245 = xr.Dataset({variable: (('time', 'lat','lon'), bias_corrected_results_245)},\n",
    "                        coords={'lat': lats,\n",
    "                                'lon': lons,\n",
    "                                'time':rcp245_dict['time'] })\n",
    "\n",
    "ds_ssp370 = xr.Dataset({variable: (('time', 'lat','lon'), bias_corrected_results_370)},\n",
    "                        coords={'lat': lats,\n",
    "                                'lon': lons,\n",
    "                                'time':rcp370_dict['time'] })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T03:48:30.097953Z",
     "start_time": "2020-05-16T03:47:32.713050Z"
    }
   },
   "outputs": [],
   "source": [
    "for scenario in scenarios:\n",
    "    name = \"ds_\"+ str(scenario)\n",
    "    ds = globals()[name]\n",
    "    name_netcdf =  model_name_file + \"_\" + variable + \"_\" + str(scenario) + '.nc'\n",
    "    name_pickle =  model_name_file + \"_\" + variable + \"_\" + str(scenario) + '.pkl'\n",
    "        \n",
    "    try:\n",
    "        os.makedirs(target_folder + \"/Bias_corrected_kelvin/\")\n",
    "        netcdf_corr = target_folder + \"/Bias_corrected_kelvin/\"+name_netcdf\n",
    "        pickle_corr = target_folder + \"/Bias_corrected_kelvin/\"+name_pickle\n",
    "        ds.to_netcdf(netcdf_corr)\n",
    "        if scenario == \"hist\":\n",
    "            save_as_pickled_object(correct_params, pickle_corr)\n",
    "\n",
    "\n",
    "    except FileExistsError:\n",
    "        netcdf_corr = target_folder + \"/Bias_corrected_kelvin/\"+name_netcdf\n",
    "        pickle_corr = target_folder + \"/Bias_corrected_kelvin/\"+name_pickle\n",
    "        ds.to_netcdf(netcdf_corr)\n",
    "        if scenario == \"hist\":\n",
    "            save_as_pickled_object(correct_params, pickle_corr)    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-16T03:48:32.079585Z",
     "start_time": "2020-05-16T03:48:30.101271Z"
    }
   },
   "outputs": [],
   "source": [
    "%reset -f"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (master)",
   "language": "python",
   "name": "master"
  },
  "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}