notebooks/plot_IFD.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "\n",
    "This code is part of the SIPN2 project focused on improving sub-seasonal to seasonal predictions of Arctic Sea Ice. \n",
    "If you use this code for a publication or presentation, please cite the reference in the README.md on the\n",
    "main page (https://github.com/NicWayand/ESIO). \n",
    "\n",
    "Questions or comments should be addressed to nicway@uw.edu\n",
    "\n",
    "Copyright (c) 2018 Nic Wayand\n",
    "\n",
    "GNU General Public License v3.0\n",
    "\n",
    "\n",
    "'''\n",
    "\n",
    "'''\n",
    "Plot forecast maps with all available models.\n",
    "'''\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload\n",
    "import matplotlib\n",
    "#matplotlib.use('Agg')\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import OrderedDict\n",
    "import itertools\n",
    "import numpy as np\n",
    "import timeit\n",
    "import numpy.ma as ma\n",
    "import pandas as pd\n",
    "import struct\n",
    "import os\n",
    "import xarray as xr\n",
    "import xesmf as xe\n",
    "import glob\n",
    "import datetime\n",
    "import calendar\n",
    "import cartopy.crs as ccrs\n",
    "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n",
    "import seaborn as sns\n",
    "np.seterr(divide='ignore', invalid='ignore')\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "\n",
    "from esio import EsioData as ed\n",
    "from esio import ice_plot\n",
    "from esio import import_data\n",
    "from esio import metrics\n",
    "\n",
    "\n",
    "import dask\n",
    "#dask.config.set(scheduler='threads')  # overwrite default with threaded scheduler\n",
    "#from dask.distributed import Client\n",
    "\n",
    "\n",
    "# General plotting settings\n",
    "sns.set_style('ticks')\n",
    "sns.set_context(\"talk\", font_scale=.8, rc={\"lines.linewidth\": 2.5})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# client = Client()\n",
    "# client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = ed.EsioData.load()\n",
    "# Directories\n",
    "runType='forecast'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting Info\n",
    "runType = 'forecast'\n",
    "variables = ['sic']\n",
    "metrics_all = {'sic':['IFD']}\n",
    "# Some models are terrible/have serious issues, so don't include in MME\n",
    "MME_NO = ['hcmr']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Target time range\n",
    "valid_start = np.datetime64('2018-06-01')\n",
    "valid_end = np.datetime64('2018-09-30')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SIO Report Year_Month\n",
    "cyear = '2018'\n",
    "cmonth = 'August'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Models that create an Ensemble using different init days\n",
    "it_target_size = {'usnavysipn':10, 'ukmetofficesipn':21} # Define the number of past days to use TODO: assumes daily inits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define custom mod names\n",
    "custom_names = {'usnavysipn':'NESM','gfdlsipn':'GFDL/NOAA','noaasipn':'NCEP CPC','uclsipn':'UCL',\n",
    "                'ukmetofficesipn':'Met Office','ecmwfsipn':'ECMWF-c3s',\n",
    "               'nicosipn':'Nico Sun','awispin':'AWI\\n(Kauker et al.)','szapirosipn':'MPAS-CESM'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define models to plot\n",
    "No_No_models = ['MME_NEW','piomas','MME','modcansipns_3','modcansipns_4','ecmwfsipn','usnavygofs']\n",
    "non_SIPN_Report_models  = ['ecmwf','ukmo','kma','ncep','usnavyncep','metreofr'] # Models like s2s\n",
    "models_2_plot = list(E.model.keys())\n",
    "models_2_plot = [x for x in models_2_plot if np.any(x not in No_No_models) ] # remove some models\n",
    "models_2_plot = [x for x in models_2_plot if np.any(x not in non_SIPN_Report_models) ] # remove some more models\n",
    "models_2_plot = [x for x in models_2_plot if E.icePredicted[x]] # Only predictive models\n",
    "models_2_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stero_grid_file = E.obs['NSIDC_0051']['grid']\n",
    "obs_grid = import_data.load_grid_info(stero_grid_file, model='NSIDC')\n",
    "# Ensure latitude is within bounds (-90 to 90)\n",
    "# Have to do this because grid file has 90.000001\n",
    "obs_grid['lat_b'] = obs_grid.lat_b.where(obs_grid.lat_b < 90, other = 90)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Regridding Options\n",
    "method='nearest_s2d' # ['bilinear', 'conservative', 'nearest_s2d', 'nearest_d2s', 'patch']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load in User submited IFD data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "user_ifd = [] # Dictionary to store user submitted (regridded) SIP data\n",
    "r_yr_mon = cyear+'_'+cmonth\n",
    "mon_2_int = {v: k for k,v in enumerate(calendar.month_name)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### cansips (Mod CanSIPS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmod = 'Mod CanSIPS'\n",
    "\n",
    "if r_yr_mon=='2018_June':\n",
    "    sip_f = '/home/disk/sipn/upload/ecsipn/June_2018_fix/ifd_bc_i2018m06_Mod_CanSIPS.nc'\n",
    "elif r_yr_mon=='2018_July':\n",
    "    sip_f = '/home/disk/sipn/upload/ecsipn/July_2018/ifd_bc_i2018m07_Mod_CanSIPS.nc'\n",
    "elif r_yr_mon=='2018_August':\n",
    "    sip_f = '/home/disk/sipn/upload/ecsipn/August_2018/ifd_bc_i2018m08_Mod_CanSIPS.nc'\n",
    "else:\n",
    "    raise ValueError(r_yr_mon,\"not found, add sip_f path.\")\n",
    "\n",
    "ds_user = xr.open_dataset(sip_f, drop_variables='time')\n",
    "ds_user.rename({'longitude':'lon','latitude':'lat'}, inplace=True);\n",
    "ds_user.set_coords(['lat','lon'], inplace=True)\n",
    "da_in = ds_user.ifd\n",
    "print(da_in.max().values)\n",
    "print(da_in.min().values)\n",
    "\n",
    "# Set never melts to first valeus (to be consistent with how sipn does it)\n",
    "# da_in = da_in.where(da_in<da_in.max().values, other=da_in.min().values)\n",
    "# Set pernial to 275\n",
    "da_in = da_in.where(da_in<da_in.max().values, other=275)\n",
    "\n",
    "# Set pixels that are ice free at init, to -9999\n",
    "da_in = da_in.where(da_in!=da_in.min().values, other=152) #pd.to_datetime('2018-06-01').timetuple().tm_yday)\n",
    "\n",
    "# Calculate regridding matrix\n",
    "regridder = xe.Regridder(da_in, obs_grid, method, periodic=False)\n",
    "# Regrid \n",
    "da_out = regridder(da_in)\n",
    "# Remove weight file\n",
    "regridder.clean_weight_file() \n",
    "\n",
    "# Store in dict of user submited SIP\n",
    "da_out.coords['model'] = cmod\n",
    "da_out = da_out.expand_dims('model')\n",
    "# da_out.coords['init_date'] = xr.DataArray(np.array([np.datetime64(datetime.date(int(cyear), 1, 1) + datetime.timedelta(da_in.min().values - 1))]), \n",
    "#                                           dims='model', coords={'model':da_out.model})\n",
    "da_out.coords['init_date'] = xr.DataArray(np.array([np.datetime64('2018-08-01')]), \n",
    "                                          dims='model', coords={'model':da_out.model})\n",
    "user_ifd.append(da_out)\n",
    "\n",
    "plt.figure()\n",
    "da_in.plot()\n",
    "plt.figure()\n",
    "da_out.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.array([np.datetime64(datetime.date(int(cyear), 1, 1) + datetime.timedelta(da_in.min().values - 1))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da_out.model.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmod"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GMAO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmod = 'GMAO'\n",
    "\n",
    "if r_yr_mon=='2018_June':\n",
    "    sip_f = ''\n",
    "elif r_yr_mon=='2018_July':\n",
    "    sip_f = '/home/disk/sipn/nicway/data/model/gmao/gmao_Sep2018_Julydata_IFD_thatssmigrid.nc4'\n",
    "    #maskf = '/home/disk/sipn/nicway/data/model/gmao/gmao_Sep2018_Julydata_SIP_thatssmigrid.nc4'\n",
    "elif r_yr_mon=='2018_August':\n",
    "    sip_f = None \n",
    "else:\n",
    "    raise ValueError(r_yr_mon,\"not found, add sip_f path.\")\n",
    "\n",
    "\n",
    "if sip_f:\n",
    "    ds_user = xr.open_dataset(sip_f)\n",
    "    landmask = ds_user.IFD.notnull()\n",
    "    \n",
    "    # Mask\n",
    "    ds_user = ds_user.set_coords(['LAT','LON'])\n",
    "\n",
    "    da_in = ds_user.rename({'LON':'lon','LAT':'lat'}).IFD\n",
    "    \n",
    "    # Set never melts to first valeus (to be consistent with how sipn does it)\n",
    "    da_in = da_in.where(da_in<da_in.max().values, other=275)\n",
    "    \n",
    "    da_in = da_in.where(landmask)\n",
    "    \n",
    "\n",
    "    # Calculate regridding matrix\n",
    "    regridder = xe.Regridder(da_in, obs_grid, method, periodic=True)\n",
    "    # Regrid \n",
    "    da_out = regridder(da_in)\n",
    "    # Remove weight file\n",
    "    regridder.clean_weight_file() \n",
    "    \n",
    "    da_in.plot()\n",
    "    plt.figure()\n",
    "    da_out.plot()\n",
    "    \n",
    "    # Store in dict of user submited SIP\n",
    "    da_out.coords['model'] = cmod\n",
    "    da_out = da_out.expand_dims('model')\n",
    "    da_out.coords['init_date'] = xr.DataArray(np.array([np.datetime64(datetime.date(int(cyear), 1, 1) + datetime.timedelta(da_in.min().values - 1))]), \n",
    "                                              dims='model', coords={'model':da_out.model})\n",
    "    user_ifd.append(da_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Concat them together by model\n",
    "usr_IFD = xr.concat(user_ifd, dim='model')\n",
    "usr_IFD = usr_IFD.rename({'nj':'x', 'ni':'y'})\n",
    "usr_IFD.name = 'IFD'\n",
    "usr_IFD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the last init day as the 11th of the current month's report (i.e. June would be the 11th of June)\n",
    "last_init_day = np.datetime64(datetime.datetime(int(cyear),mon_2_int[cmonth],11))\n",
    "# Define the earliest init day as 1 year before last_init_day ( some models init on 1st of month (i.e. UCL))\n",
    "first_init_day = np.datetime64(datetime.datetime(int(cyear),1,1))-np.timedelta64(2,'D')    # last_init_day-np.timedelta64(365,'D')\n",
    "print(\"Looking for init times between\",first_init_day,\"and\",last_init_day)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for cvar in variables:\n",
    "    \n",
    "    # Define fig dir and make if doesn't exist\n",
    "    fig_dir = os.path.join(E.fig_dir, 'model', 'SIO_Reports', cyear, cmonth)\n",
    "    if not os.path.exists(fig_dir):\n",
    "        os.makedirs(fig_dir)\n",
    "       \n",
    "    # Loop through variable of interest + any metrics (i.e. SIP) based on that\n",
    "    for metric in metrics_all[cvar]:\n",
    "\n",
    "        MME_list = []      \n",
    "    \n",
    "        # For each model\n",
    "        for (i, cmod) in enumerate(models_2_plot): # ['nicosipn']\n",
    "            print(cmod)\n",
    "\n",
    "            # Load in Model\n",
    "            model_forecast = os.path.join(E.model[cmod][runType]['sipn_nc'], '*.nc') \n",
    "\n",
    "            # Check we have files \n",
    "            files = glob.glob(model_forecast)\n",
    "            if not files:\n",
    "                #print(\"Skipping model\", cmod, \"no forecast files found.\")\n",
    "                continue # Skip this model\n",
    "            \n",
    "            ds_model = xr.open_mfdataset(model_forecast, \n",
    "                        chunks={'init_time': 1, 'nj': 304, 'ni': 448}, \n",
    "                                         concat_dim='init_time', autoclose=True, parallel=True)\n",
    "            ds_model.rename({'nj':'x', 'ni':'y'}, inplace=True)\n",
    "\n",
    "            # Restrictions for calculating IFD\n",
    "            \n",
    "            # 1) Must be Daily or less\n",
    "            c_dt = (ds_model.fore_time[1]-ds_model.fore_time[0]).values.astype('timedelta64[D]').astype('int')\n",
    "            if c_dt > 1:\n",
    "                print(\"    fore_time timestep greater than 1 day, skipping.\")\n",
    "                continue\n",
    "            \n",
    "            # Select var of interest (if available)\n",
    "            if cvar in ds_model.variables:\n",
    "                ds_model = ds_model[cvar]\n",
    "            else:\n",
    "                print('    cvar not found.')\n",
    "                continue\n",
    "                  \n",
    "            # Use lagged ensemble\n",
    "            if cmod in list(it_target_size.keys()): # Uses a lagged ensemble... so grab all inits and make them ensembles\n",
    "                \n",
    "                # NESM (Joe) wants us to use the past months inits (June for July report)\n",
    "                if cmod == 'usnavysipn':\n",
    "                    # Replace last_init_day inline below (note the \"-1\" meaning last month)\n",
    "                    ds_model = ds_model.sel(init_time=slice(first_init_day,  np.datetime64(datetime.datetime(int(cyear),mon_2_int[cmonth] - 1, 11))  )) # First get all time between start and stop (up to and including the 10th)\n",
    "                else:\n",
    "                    ds_model = ds_model.sel(init_time=slice(first_init_day,last_init_day)) # First get all time between start and stop (up to and including the 10th)\n",
    "                \n",
    "                if ds_model.init_time.size==0:\n",
    "                    print(\"    No init_times found in requested start range\",first_init_day,\"to\",last_init_day,\". So skipping...\")\n",
    "                    continue\n",
    "                    \n",
    "                print(\"    Init times used are\",ds_model.init_time[-1*it_target_size[cmod]:].values)    \n",
    "                \n",
    "                ds_model = ds_model.sel(init_time=ds_model.init_time[-1*it_target_size[cmod]:]) # Get nearest date to 10th\n",
    "                avg_init_time = ds_model.init_time.astype('int').max().astype('datetime64[ns]')\n",
    "                ds_model = ds_model.stack(hybrid_ensemble=('ensemble', 'init_time')).reset_index('hybrid_ensemble')\n",
    "                ds_model.coords['hybrid_ensemble'] = np.arange(1,ds_model.hybrid_ensemble.size+1,1)\n",
    "                ds_model = ds_model.drop(['ensemble','init_time'])\n",
    "                ds_model.coords['init_time'] = avg_init_time\n",
    "                ds_model = ds_model.rename({'hybrid_ensemble':'ensemble'})\n",
    "                print(ds_model)\n",
    "\n",
    "            # Normal Ensemble\n",
    "            else:\n",
    "                # Find init time closest and earlier to the 10th of the month\n",
    "                ds_model = ds_model.sel(init_time=slice(first_init_day,last_init_day)) # First get all time between start and stop (up to and including the 10th)\n",
    "                if ds_model.init_time.size==0:\n",
    "                    print(\"    No init_times found in requested start range\",first_init_day,\"to\",last_init_day,\". So skipping...\")\n",
    "                    continue\n",
    "                ds_model = ds_model.sel(init_time=last_init_day, method='nearest') # Get nearest date to 10th\n",
    "                \n",
    "            print(\"    Init time used is\", ds_model.init_time.values)\n",
    "            \n",
    "            \n",
    "            # Get Valid time\n",
    "            ds_model = import_data.get_valid_time(ds_model.expand_dims('init_time'))\n",
    "                           \n",
    "            # Check if we have any valid times in range of target dates\n",
    "            ds_model = ds_model.where((ds_model.valid_time>=valid_start) & (ds_model.valid_time<=valid_end), drop=True) \n",
    "            if ds_model.fore_time.size == 0:\n",
    "                print(\"    No fore_time found for target period.\")\n",
    "                continue\n",
    "                \n",
    "            # 2) Must have a forecast that goes until end of sept.\n",
    "            last_f = ds_model.valid_time.max()\n",
    "            if last_f < valid_end - np.timedelta64(20,'D'): # allow 5 days wiggle room\n",
    "                print('    Last Valid time', last_f.values, ' was not within 20 days of target end ', valid_end)\n",
    "                continue\n",
    "                \n",
    "            # Averave SIC over ensembles\n",
    "            ds_model = ds_model.mean(dim='ensemble')\n",
    "            ds_model.load()\n",
    "                \n",
    "            # Calculate IFD\n",
    "            print(\"    Calculating IFD\")\n",
    "#             if cmod=='szapirosipn':\n",
    "#                 xr.exit()\n",
    "            DOY_first_vt = [x.timetuple().tm_yday for x in pd.to_datetime([ds_model.valid_time.min().values])]\n",
    "            ds_ifd = metrics.calc_IFD_10day(ds_model, sic_threshold=0.5, DOY_s=DOY_first_vt[0], time_dim='fore_time', \n",
    "                                            Nday=10, default_ice_free=pd.to_datetime('2018-06-01').timetuple().tm_yday)\n",
    "\n",
    "            # Build MME\n",
    "            ds_idf = ds_ifd.isel(init_time=0).drop('init_time')\n",
    "\n",
    "            if cmod not in MME_NO: # Exclude some models (bad) from MME\n",
    "                ds_ifd.coords['model'] = cmod\n",
    "                ds_ifd = ds_ifd.expand_dims('model')\n",
    "                ds_ifd.coords['init_date'] = xr.DataArray(ds_model.init_time.values, dims='model', coords={'model':ds_ifd.model})\n",
    "                MME_list.append(ds_ifd)\n",
    "                print('    Added ', cmod, ' to MME.')\n",
    "\n",
    "\n",
    "# Concat over all models\n",
    "ds_MME = xr.concat(MME_list, dim='model')\n",
    "ds_MME.name = 'IFD'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DOY_first_vt = [x.timetuple().tm_yday for x in pd.to_datetime([ds_model.valid_time.min().values])]\n",
    "# ds_ifd = metrics.calc_IFD_10day(ds_model, sic_threshold=0.5, DOY_s=DOY_first_vt[0], time_dim='fore_time', \n",
    "#                                 Nday=10, default_ice_free=pd.to_datetime('2018-06-01').timetuple().tm_yday)\n",
    "# print(ds_ifd.where(ds_ifd<275).max().values)\n",
    "# plt.figure()\n",
    "# ds_ifd.isel(init_time=0).plot(vmin=153,vmax=274)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Take mean\n",
    "model_inits = ds_MME.init_time.values\n",
    "ds_MME = ds_MME.mean(dim='init_time')\n",
    "ds_MME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Merge user SIP with SPIN SIP\n",
    "ds_IFD_All = xr.concat([ds_MME, usr_IFD], dim='model')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load into memory\n",
    "ds_IFD_All.load();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save to disk\n",
    "ds_IFD_All.to_netcdf('/home/disk/sipn/nicway/data/model/SIO/'+cyear+'/'+cmonth+'/IFD.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get mean and std of all models (set -9999 (ice free at init) and perenial ice (275) to missing))\n",
    "ds_MME_avg = ds_IFD_All.mean(dim='model') #.where((ds_IFD_All>1) & (ds_IFD_All<275))\n",
    "ds_MME_std = ds_IFD_All.std(dim='model')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_vmin = [x.timetuple().tm_yday for x in [pd.datetime(2018,6,1)]][0]\n",
    "c_vmax = [x.timetuple().tm_yday for x in [pd.datetime(2018,10,1)]][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get dates for plot\n",
    "last_IFD = pd.to_datetime((np.datetime64('2018-01-01') + np.timedelta64(int(c_vmax),'D')))\n",
    "date_str = pd.to_datetime(datetime.datetime(2018, 6, 1)).strftime('%Y-%m-%d')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get monthly ticks\n",
    "DOY_months = [x.timetuple().tm_yday for x in pd.date_range(date_str,'2018-10-30', freq='MS')]\n",
    "DOY_labels = [calendar.month_abbr[mn]+' 1' for mn in np.arange(6,last_IFD.month+1)]\n",
    "assert len(DOY_months)==len(DOY_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_IFD_All.loc[dict(model='Mod CanSIPS')] = ds_IFD_All.sel(model='Mod CanSIPS').where(ds_IFD_All.isel(model=0).notnull())\n",
    "ds_MME_std = ds_MME_std.where(ds_IFD_All.isel(model=0).notnull())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up plotting info\n",
    "if cvar=='sic':\n",
    "    if metric=='mean':\n",
    "        cmap_c = matplotlib.colors.ListedColormap(sns.color_palette(\"Blues_r\", 10))\n",
    "        cmap_c.set_bad(color = 'lightgrey')\n",
    "        c_label = 'Sea Ice Concentration (-)'\n",
    "        c_vmin = 0\n",
    "        c_vmax = 1\n",
    "    elif metric=='SIP':\n",
    "        cmap_c = matplotlib.colors.LinearSegmentedColormap.from_list(\"\", [\"white\",\"orange\",\"red\",\"#990000\"])\n",
    "        #cmap_c = plt.get_cmap('jet') \n",
    "        cmap_c.set_bad(color = 'lightgrey')\n",
    "        c_label = 'Sea Ice Probability (-)'\n",
    "        c_vmin = 0\n",
    "        c_vmax = 1\n",
    "    elif metric=='anomaly':\n",
    "#                         cmap_c = matplotlib.colors.ListedColormap(sns.color_palette(\"coolwarm\", 9))\n",
    "        cmap_c = matplotlib.colors.LinearSegmentedColormap.from_list(\"\", [\"red\",\"white\",\"blue\"])\n",
    "        cmap_c.set_bad(color = 'lightgrey')\n",
    "        c_label = 'SIC Anomaly to 1980-2010 Mean'\n",
    "        c_vmin = -1\n",
    "        c_vmax = 1\n",
    "    elif metric=='IFD':\n",
    "#         cmap_c = matplotlib.colors.LinearSegmentedColormap.from_list(\"\", [\"white\",\"orange\",\"red\",\"#990000\"], N=8)\n",
    "#         cmap_c.set_bad(color = 'lightgrey')\n",
    "        \n",
    "        cfile = '/home/disk/sipn/nicway/python/ESIO/notebooks/color_redhues.mat'\n",
    "        import scipy.io\n",
    "        matfile = scipy.io.loadmat(cfile)\n",
    "        redhues = matfile['redhues']\n",
    "        # add perenial ice color (black)\n",
    "        #redhues_b = np.concatenate((redhues,np.array([[0,0,0]])),axis=0)\n",
    "        cmap_c = matplotlib.colors.LinearSegmentedColormap.from_list('redhues', redhues, N=8)\n",
    "        cmap_c.set_bad(color = 'lightgrey')\n",
    "        \n",
    "        #cmap_std = matplotlib.colors.LinearSegmentedColormap.from_list('redhues', redhues, N=8)\n",
    "        cmap_std = matplotlib.colors.ListedColormap(sns.color_palette(\"Blues\", 10))\n",
    "        cmap_std.set_bad(color = 'lightgrey')\n",
    "        \n",
    "        cmap_b =  matplotlib.colors.ListedColormap(sns.color_palette(\"binary\", 2))\n",
    "        cmap_b.set_bad(color = 'None')\n",
    "        \n",
    "        cmap_obs =  matplotlib.colors.ListedColormap(sns.color_palette(\"Greys\", 2))\n",
    "        cmap_obs.set_bad(color = 'None')\n",
    "        \n",
    "        c_label = 'Ice Free Day'\n",
    "\n",
    "\n",
    "elif cvar=='hi':\n",
    "    if metric=='mean':\n",
    "        cmap_c = matplotlib.colors.ListedColormap(sns.color_palette(\"Reds_r\", 10))\n",
    "        cmap_c.set_bad(color = 'lightgrey')\n",
    "        c_label = 'Sea Ice Thickness (m)'\n",
    "        c_vmin = 0\n",
    "        c_vmax = None\n",
    "else:\n",
    "    raise ValueError(\"cvar not found.\") \n",
    "\n",
    "\n",
    "    \n",
    "def add_subplot_title(cmod, custom_names, E, ax=None):\n",
    "    if cmod in custom_names:\n",
    "        ax.set_title(custom_names[cmod])\n",
    "    elif cmod in E.model.keys():\n",
    "        ax.set_title(E.model[cmod]['model_label'])\n",
    "    else:\n",
    "        ax.set_title(cmod)\n",
    "\n",
    "        \n",
    "# New Plot\n",
    "central_extent = [-3850000*0.6, 3725000*0.6, -5325000*0.45, 5850000*0.45] # (x0, x1, y0, y1\n",
    "(f, axes) = ice_plot.multi_polar_axis(ncols=4, nrows=3, Nplots=ds_IFD_All.model.size + 1 + 1 + 1, extent=central_extent, central_longitude=0)\n",
    "\n",
    "for (i, cmod) in enumerate(ds_IFD_All.model.values):\n",
    "    print(cmod)\n",
    "    # Plot\n",
    "    add_subplot_title(cmod, custom_names, E, ax=axes[i])\n",
    "    p = ds_IFD_All.sel(model=cmod).plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                          transform=ccrs.PlateCarree(),\n",
    "                          add_colorbar=False,\n",
    "                          cmap=cmap_c,\n",
    "                          vmin=c_vmin, vmax=c_vmax)\n",
    "    # Plot perenial ice as black\n",
    "    cper_mask = ds_IFD_All.sel(model=cmod)\n",
    "    pb = cper_mask.where(cper_mask==275).plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                          transform=ccrs.PlateCarree(),\n",
    "                          add_colorbar=False,\n",
    "                          cmap=cmap_b,\n",
    "                          vmin=0, vmax=1)    \n",
    "    \n",
    "    add_subplot_title(cmod, custom_names, E, ax=axes[i])\n",
    "\n",
    "    # Add init date\n",
    "    axes[i].annotate(pd.to_datetime(ds_IFD_All.sel(model=cmod).init_date.values).strftime('%m/%d'), xy=(.6, -.15), xycoords='axes fraction', )\n",
    "    \n",
    "    \n",
    "# Add IFD from Observations\n",
    "i = i + 1\n",
    "\n",
    "ds_81 = xr.open_mfdataset(E.obs['NSIDC_0081']['sipn_nc']+'_yearly/*.nc', concat_dim='time', autoclose=True, parallel=True)\n",
    "ds_81 = ds_81.sel(time=slice(first_init_day+np.timedelta64(2,'D'),datetime.datetime.now()))\n",
    "DOY_first_vt = [x.timetuple().tm_yday for x in pd.to_datetime([ds_81.time.min().values])]\n",
    "obs_ifd = metrics.calc_IFD_10day(ds_81.sic, sic_threshold=0.5, DOY_s=DOY_first_vt[0], time_dim='time', \n",
    "                                Nday=10, default_ice_free=pd.to_datetime('2018-06-01').timetuple().tm_yday)\n",
    "\n",
    "obs_h = obs_ifd.plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                      transform=ccrs.PlateCarree(),\n",
    "                      add_colorbar=False,\n",
    "                      cmap=cmap_c,\n",
    "                      vmin=c_vmin, vmax=c_vmax)\n",
    "# Mask\n",
    "pb_obs = obs_ifd.where(obs_ifd==275).plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                      transform=ccrs.PlateCarree(),\n",
    "                      add_colorbar=False,\n",
    "                      cmap=cmap_b,\n",
    "                      vmin=0, vmax=0.5)  \n",
    "cd = datetime.datetime.now()\n",
    "cd = datetime.datetime(cd.year, cd.month, cd.day)\n",
    "axes[i].set_title('Observed\\n('+pd.to_datetime(cd).strftime('%Y-%m-%d')+')')    \n",
    "    \n",
    "    \n",
    "# MME Standard deviation\n",
    "i = i + 1\n",
    "\n",
    "\n",
    "pstd = ds_MME_std.plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                                  transform=ccrs.PlateCarree(),\n",
    "                                  add_colorbar=False, \n",
    "                                  cmap=cmap_std,vmin=ds_MME_std.min().values, vmax=ds_MME_std.max().values)\n",
    "axes[i].set_title('σ')\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Make pretty\n",
    "f.subplots_adjust(bottom=0.05)\n",
    "cbar_ax = f.add_axes([0.25, 0.001, .5, 0.04]) #  [left, bottom, width, height] w\n",
    "cbar = f.colorbar(p, cax=cbar_ax, label=c_label, orientation='horizontal')\n",
    "cbar.set_ticks(DOY_months)\n",
    "cbar.set_ticklabels(DOY_labels)\n",
    "\n",
    "#f.subplots_adjust(right=0.05)\n",
    "# cbar_ax_2 = f.add_axes([0.925, 0.1, .025, .28]) #  [left, bottom, width, height] w\n",
    "cbar_ax_2 = f.add_axes([0.55, 0.07, .025, .20]) #  [left, bottom, width, height] w\n",
    "\n",
    "cbar2 = f.colorbar(pstd, cax=cbar_ax_2, label='days', orientation='vertical')\n",
    "# cbar.set_ticks(DOY_months)\n",
    "# cbar.set_ticklabels(DOY_labels)\n",
    "\n",
    "# Save to file\n",
    "f_out = os.path.join(fig_dir,'panArctic_'+metric+'_'+runType+'_'+cyear+'_'+cmonth+'_plus_std.png')\n",
    "f.savefig(f_out,bbox_inches='tight', dpi=300)\n",
    "f_out = os.path.join(fig_dir,'panArctic_'+metric+'_'+runType+'_'+cyear+'_'+cmonth+'_plus_std_lowRes.png')\n",
    "f.savefig(f_out,bbox_inches='tight', dpi=90)\n",
    "print(\"saved \", f_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "3.6.4 esio",
   "language": "python",
   "name": "esio"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}