notebooks/Eval_weekly_forecasts.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paper or Website plots\n",
    "PaperPlots = False # Paper and website have different formating. True for running for paper figures. False for running website plots."
   ]
  },
  {
   "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",
    "if not PaperPlots:\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 numpy.ma as ma\n",
    "import pandas as pd\n",
    "import struct\n",
    "import os\n",
    "import xarray as xr\n",
    "import glob\n",
    "import datetime\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",
    "import json\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",
    "import dask\n",
    "import timeit\n",
    "from dask.distributed import Client\n",
    "\n",
    "# General plotting settings\n",
    "sns.set_style('whitegrid')\n",
    "sns.set_context(\"talk\", font_scale=.8, rc={\"lines.linewidth\": 2.5})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Update_Evaluation_Maps():\n",
    "\n",
    "    client = Client(n_workers=8)\n",
    "    client\n",
    "\n",
    "    ### Paper Figure Parameters\n",
    "\n",
    "    # Time period options\n",
    "    Last_valid_time = np.datetime64('2018-11-30') # last day in valid time space to include in analysis. \n",
    "    # Use  np.datetime64('2018-11-30') to reproduce Paper.\n",
    "\n",
    "    # Domain options to use for calculating metrics over\n",
    "    customDomainName = 'ArcticSeas' \n",
    "    # Options:\n",
    "    # 'ArcticSeas' will use the Arctic Seas domain shown in Supplemental information\n",
    "    # 'X' where X is a model name to use as the domain (i.e. X='rasmesrl' because it is the smallest domain)\n",
    "\n",
    "    ### Code\n",
    "\n",
    "    metrics_all = ['anomaly','mean','SIP']\n",
    "    runType = 'forecast'\n",
    "    variables = ['sic']\n",
    "\n",
    "    # Only use 'sic' here\n",
    "    cvar = variables[0]\n",
    "\n",
    "    # Get path data\n",
    "    E = ed.EsioData.load()\n",
    "    grid_dir = E.grid_dir\n",
    "\n",
    "    # Load in regional data\n",
    "    ds_region = xr.open_dataset(os.path.join(grid_dir, 'sio_2016_mask_Update.nc'))\n",
    "\n",
    "    # Define fig dir and make if doesn't exist (Hard coded here)\n",
    "    if PaperPlots:\n",
    "        fig_dir = os.path.join('/home/disk/sipn/nicway/Nic/figures', 'model', 'MME', cvar, 'BSS')\n",
    "    else:\n",
    "        fig_dir = os.path.join('/home/disk/sipn/nicway/public_html/sipn/figures/model/Eval')\n",
    "\n",
    "    if not os.path.exists(fig_dir):\n",
    "        os.makedirs(fig_dir)\n",
    "\n",
    "    # Load from Zarr\n",
    "    ds_m = xr.open_zarr(os.path.join(E.data_dir,'model/zarr',cvar+'.zarr'))\n",
    "\n",
    "    # Select only initlizations from January 2018 to November 2018 (for GRL paper)\n",
    "    if PaperPlots:\n",
    "        ds_m = ds_m.where(ds_m.valid_end<=Last_valid_time, drop=True)\n",
    "    # print(\"Warning: using ALL data\")\n",
    "\n",
    "    # Drop models we don't want\n",
    "    ds_ALL = ds_m.where(ds_m.model!='dampedAnomaly', drop=True)\n",
    "\n",
    "    def remove_list(X, L):\n",
    "        for cl in L:\n",
    "            X.remove(cl)\n",
    "\n",
    "    # Get recent observations\n",
    "    ds_81 = xr.open_mfdataset(E.obs['NSIDC_0081']['sipn_nc']+'_yearly/*.nc', concat_dim='time', autoclose=True, parallel=True)\n",
    "\n",
    "    # Hack to decode strings\n",
    "    ds_ALL['model'] = [s.decode(\"utf-8\") for s in ds_ALL.model.values]\n",
    "\n",
    "    # # Drop models that we don't evaluate (i.e. monthly means)\n",
    "    exl_mods = ['awispin','nicosipn','szapirosipn']\n",
    "\n",
    "    # Order models to plot in\n",
    "    models_first = ['climatology','dampedAnomalyTrend','MME']\n",
    "    X = list(ds_ALL.model.values)\n",
    "    [X.remove(a) for a in models_first]\n",
    "    model_plot_order = models_first + X\n",
    "    model_plot_order.remove('Observed')\n",
    "    # model_plot_order\n",
    "\n",
    "\n",
    "\n",
    "    # Get custom domain mask if we have defined it (i.e. not None)\n",
    "    if customDomainName=='ArcticSeas':\n",
    "        customDomain = None\n",
    "        fig_ext = '' # Figure extention to separte domain names\n",
    "    else:\n",
    "        customDomain = ds_ALL.sel(model=customDomainName).isel(fore_time=0,init_end=40).SIP.notnull()\n",
    "        # Drop coods\n",
    "        customDomain = customDomain.drop(['fore_time','init_end','init_start','valid_end','valid_start','model'])\n",
    "        if customDomain.sum().values==0:\n",
    "            raise ValueError(\"Warning: you need to choose a fore_time and init_end that has model data represenative of its extent\")\n",
    "        fig_ext = customDomainName+'_' # Figure extention to separte domain names\n",
    "\n",
    "    # Calculate the IIEE (Integrated Ice Edge Error)\n",
    "    l = []\n",
    "    for cmod in ds_ALL.model.values: \n",
    "        c_SIP_IIEE = metrics.IIEE(da_mod=ds_ALL.sel(model=cmod)['mean'], \n",
    "                                          da_obs=ds_ALL.sel(model='Observed')['mean'], \n",
    "                                          region=ds_region, \n",
    "                                          testplots=False,\n",
    "                                          customDomain=customDomain)\n",
    "        c_SIP_IIEE.coords['model'] = cmod\n",
    "        l.append(c_SIP_IIEE)\n",
    "    SIP_IIEE = xr.concat(l, dim='model')\n",
    "    SIP_IIEE\n",
    "\n",
    "    # For SIP, calculate the Brier Skill Score for panArctic \n",
    "    # BS = (SIP_model - SIP_observed)^2\n",
    "    #\n",
    "    l = []\n",
    "    for cmod in ds_ALL.model.values: \n",
    "        c_SIP_BSS = metrics.BrierSkillScore(da_mod_sip=ds_ALL.sel(model=cmod).SIP, \n",
    "                                          da_obs_ip=ds_ALL.sel(model='Observed').SIP, \n",
    "                                          region=ds_region, \n",
    "                                          testplots=False,\n",
    "                                          customDomain=customDomain)\n",
    "        c_SIP_BSS.coords['model'] = cmod\n",
    "        l.append(c_SIP_BSS)\n",
    "    SIP_BSS = xr.concat(l, dim='model')\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    def add_subplot_title(cmod, E, ax=None, BSS_val=''):\n",
    "        if cmod in E.model.keys():\n",
    "            ax.set_title((E.model[cmod]['model_label']).rstrip('*')+'\\n('+BSS_val+')')\n",
    "        else:\n",
    "            ax.set_title(cmod)\n",
    "\n",
    "    # add missing info for climatology\n",
    "    E.model_color['climatology'] = (0,0,0)\n",
    "    E.model_linestyle['climatology'] = '--'\n",
    "    E.model_marker['climatology'] = '*'\n",
    "    E.model['climatology'] = {'model_label':'Climatology\\nTrend'}\n",
    "\n",
    "\n",
    "\n",
    "    # Aggregate over space (x,y), including all pixels in valid Arctic seas (masked above with BrierSkillScore())\n",
    "    BSS_agg = SIP_BSS.mean(dim=['x','y'])\n",
    "    BSS_agg.load() # Compute and load result into memory\n",
    "\n",
    "    ### At what lead time is the MME significantly (95%) better than the Damped Anomaly?\n",
    "\n",
    "    from scipy import stats\n",
    "\n",
    "    model_1 = 'MME'\n",
    "    model_2 = 'dampedAnomalyTrend'\n",
    "    alphaval = 0.05 # 95%\n",
    "\n",
    "    t_all=[]\n",
    "    p_all=[]\n",
    "    cv_all=[]\n",
    "    for ft in np.arange(0,BSS_agg.fore_time.size):\n",
    "        x = BSS_agg.sel(model=model_1).isel(fore_time=ft)\n",
    "        y = BSS_agg.sel(model=model_2).isel(fore_time=ft)\n",
    "\n",
    "        x = x.where(x.notnull() & y.notnull(), drop=True)\n",
    "        y = y.where(x.notnull() & y.notnull(), drop=True)\n",
    "\n",
    "        df = x.size + y.size -2\n",
    "\n",
    "        t, p = stats.ttest_ind(x, y, equal_var=False)\n",
    "        cv = stats.t.ppf(1.0 - alphaval, df)\n",
    "\n",
    "        t_all.append(t)\n",
    "        p_all.append(p)\n",
    "        cv_all.append(cv)\n",
    "\n",
    "\n",
    "    plt.figure()\n",
    "    plt.plot(BSS_agg.fore_time.values.astype('timedelta64[D]').astype(int)/7,\n",
    "            abs(np.array(t_all)),'-k*', label='t-value')\n",
    "    plt.plot(BSS_agg.fore_time.values.astype('timedelta64[D]').astype(int)/7,\n",
    "            cv_all,'-r*', label='critical t-value')\n",
    "    plt.ylabel('t-value')\n",
    "    plt.legend()\n",
    "\n",
    "\n",
    "    plt.figure()\n",
    "    plt.plot(BSS_agg.fore_time.values.astype('timedelta64[D]').astype(int)/7,\n",
    "            p_all,'-k*',label='p value')\n",
    "    plt.plot(BSS_agg.fore_time.values.astype('timedelta64[D]').astype(int)/7,\n",
    "            np.ones(len(p_all))*alphaval,'-r*',label='Alpha p value')\n",
    "    plt.ylabel('p-value')\n",
    "    plt.legend()\n",
    "\n",
    "    # Read off what lead time in weeks where the black line crosses the red line\n",
    "\n",
    "    ### Plot BS spatial plots valid for Summer months (Figure 3)\n",
    "\n",
    "    # Remove some select models\n",
    "    enough_init = model_plot_order.copy()\n",
    "    [enough_init.remove(a) for a in exl_mods]\n",
    "    print(enough_init)\n",
    "\n",
    "    # Define lead time to plot\n",
    "    week_lead_time = 4\n",
    "\n",
    "    period_dict = {'June':[np.datetime64('2018-06-01'),np.datetime64('2018-07-01')],\n",
    "                  'July':[np.datetime64('2018-07-01'),np.datetime64('2018-08-01')],\n",
    "                  'August':[np.datetime64('2018-08-01'),np.datetime64('2018-09-01')],\n",
    "                  'September':[np.datetime64('2018-09-01'),np.datetime64('2018-10-01')]}\n",
    "\n",
    "    for period in period_dict:\n",
    "        print(period)\n",
    "        start_date_map = period_dict[period][0]\n",
    "        end_date_map = period_dict[period][1]    \n",
    "\n",
    "        # Select time slice of valid\n",
    "        BS_Sept = SIP_BSS.where( (SIP_BSS.valid_start>=start_date_map) & \n",
    "                      (SIP_BSS.valid_start<=end_date_map), drop=True)\n",
    "\n",
    "        # Average over valid time \n",
    "        SIP_BSS_init_avg = BS_Sept.sel(model=enough_init).mean(dim='init_end')\n",
    "\n",
    "        sns.set_context(\"talk\", font_scale=.8, rc={\"lines.linewidth\": 2.5})\n",
    "\n",
    "        # Set up color maps\n",
    "        cmap_c = matplotlib.colors.LinearSegmentedColormap.from_list(\"\", [\"white\",\"orange\",\"red\",\"#990000\"], N=10)\n",
    "        cmap_c.set_bad(color = 'lightgrey')\n",
    "        c_label = 'BS (0=best, 1=worst)'\n",
    "        c_vmin = 0\n",
    "        c_vmax = 1\n",
    "\n",
    "        for ft in [SIP_BSS_init_avg.fore_time.values[week_lead_time]]:  \n",
    "\n",
    "            # Grab current lead time\n",
    "            c_ft_ds = SIP_BSS_init_avg.sel(fore_time=ft)\n",
    "\n",
    "            # Drop models without any data\n",
    "            c_ft_ds = c_ft_ds.where(c_ft_ds.notnull().sum(dim=['x','y'])>0, drop=True)\n",
    "\n",
    "            nrows = np.int(np.ceil(np.sqrt(c_ft_ds.model.size))) - 1\n",
    "            ncols = int(np.ceil(c_ft_ds.model.size/nrows))\n",
    "            Nplots = c_ft_ds.model.size + 1\n",
    "            #print(Nplots)\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=ncols, nrows=nrows, Nplots=Nplots, \n",
    "                                                  extent=central_extent, central_longitude=0)\n",
    "\n",
    "\n",
    "            # Get Observed mean sea ice edge location mean(SIC in sept) > 0.15\n",
    "            obs_SIP = ds_81['sic'].sel(time=slice(start_date_map,end_date_map)).mean(dim='time')>=0.15\n",
    "            # Fill in pole hole\n",
    "            obs_SIP = obs_SIP.where(obs_SIP.hole_mask==0, other=1)\n",
    "\n",
    "            for (i, cmod) in enumerate(c_ft_ds.model.values):\n",
    "                if cmod in c_ft_ds.model.values:\n",
    "                    # Plot\n",
    "                    add_subplot_title(cmod, E, ax=axes[i])\n",
    "                    p = c_ft_ds.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",
    "\n",
    "                    po = obs_SIP.plot.contour(ax=axes[i], x='xm', y='ym',\n",
    "                                          transform=ccrs.NorthPolarStereo(central_longitude=-45),\n",
    "                                          colors=('k'),\n",
    "                                          linewidths=[1],\n",
    "                                          levels=[0.5]) #, label='Median ice edge 1981-2010')\n",
    "                    ice_plot.remove_small_contours(po, thres=10**6)\n",
    "\n",
    "                    add_subplot_title(cmod, E, ax=axes[i], BSS_val='{0:.3f}'.format(c_ft_ds.sel(model=cmod).mean(dim=['x','y']).load().item()))\n",
    "\n",
    "            # Make pretty\n",
    "            cbar_ax = f.add_axes([0.2, 0.05, .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(np.arange(-1,1.1,0.2))\n",
    "\n",
    "            # Set title of all plots\n",
    "            lead_time_days = str(ft.astype('timedelta64[D]').astype(int))\n",
    "            #print(lead_time_days)\n",
    "\n",
    "            if not PaperPlots: # only add for website plots\n",
    "                cbar_ax.text(0.35, 1.1, 'Wayand et al. (in review)', fontsize=12)\n",
    "\n",
    "            valid_start_str = pd.to_datetime(start_date_map).strftime('%Y-%m-%d')\n",
    "            valid_end_str = pd.to_datetime(end_date_map).strftime('%Y-%m-%d')\n",
    "\n",
    "            # Save to file\n",
    "            f_out = os.path.join(fig_dir,fig_ext+'BSS_Avg_Valid_'+valid_start_str+'_to_'+valid_end_str+'_'+lead_time_days.zfill(3)+'_day_lead_time.png')\n",
    "            f.savefig(f_out,bbox_inches='tight', dpi=300)\n",
    "\n",
    "\n",
    "\n",
    "    ### Plot Brier Score vs lead time\n",
    "\n",
    "\n",
    "    min_N_samples = 10 # Min number of samples to allow for mean\n",
    "    BSS_agg_init = BSS_agg.mean(dim='init_end')\n",
    "    sns.set_style(\"whitegrid\")\n",
    "\n",
    "    sns.set_context(\"talk\", font_scale=1.5, rc={\"lines.linewidth\": 2.5})\n",
    "\n",
    "    # Get sample size of for each lead time\n",
    "    for_sample = BSS_agg.sel(model='MME').notnull().sum(dim='init_end')\n",
    "    for_sample\n",
    "\n",
    "    # Use threshold of sample size to cut off lead times\n",
    "    max_lead = for_sample.where(for_sample>=min_N_samples,drop=True).fore_time.max().values.astype('timedelta64[D]').astype(int)/7\n",
    "\n",
    "    f = plt.figure(figsize=(10,10))\n",
    "    NM = 10\n",
    "    ax1 = plt.subplot2grid((NM, NM), (0, 0), colspan=NM, rowspan=NM-1)\n",
    "    ax2 = plt.subplot2grid((NM, NM), (NM-1, 0), colspan=NM, rowspan=1)\n",
    "\n",
    "    for cmod in model_plot_order:\n",
    "        if cmod in exl_mods:\n",
    "            continue\n",
    "        # Get model plotting specs\n",
    "        cc = E.model_color[cmod]\n",
    "        cl = E.model_linestyle[cmod]\n",
    "        cm = E.model_marker[cmod]\n",
    "        if cmod=='rasmesrl':\n",
    "            cflag = '*'\n",
    "        else:\n",
    "            cflag = ''\n",
    "\n",
    "        if cmod in ['MME','dampedAnomalyTrend','climatology']:\n",
    "            lw=5\n",
    "        else:\n",
    "            lw = 2\n",
    "\n",
    "        ax1.plot(BSS_agg_init.fore_time.values.astype('timedelta64[D]').astype(int)/7,\n",
    "                BSS_agg_init.sel(model=cmod).values, label=E.model[cmod]['model_label'].rstrip('*')+cflag,\n",
    "                color=cc,\n",
    "                linestyle=cl,\n",
    "                linewidth=lw,\n",
    "                marker=cm)\n",
    "    ax1.legend(loc='lower right', bbox_to_anchor=(1.4, 0))\n",
    "    ax1.set_ylabel('Pan-Arctic BS (-)')\n",
    "    ax1.set_xlim([-0.5,max_lead])\n",
    "    ax1.set_xticklabels([''])\n",
    "\n",
    "    # second axis\n",
    "    ax2.plot(for_sample.fore_time.values.astype('timedelta64[D]').astype(int)/7,\n",
    "             for_sample.values, '-ko')\n",
    "    ax2.set_ylabel('#\\nweeks')\n",
    "    ax2.set_xlabel('Lead time (Weeks)')\n",
    "    ax2.set_xlim(ax1.get_xlim());\n",
    "\n",
    "    ax2.set_ylim([0,for_sample.max()+5]);\n",
    "    ax2.set_yticks(np.arange(0,for_sample.max()+5,15));\n",
    "\n",
    "    if not PaperPlots: # only add for website plots\n",
    "        ax1.text(17, 0.025, 'Wayand et al. (in review)', fontsize=12)\n",
    "\n",
    "    # Save to file\n",
    "    f_out = os.path.join(fig_dir,fig_ext+'BSS_by_lead_time_PanArctic.png')\n",
    "    f.savefig(f_out,bbox_inches='tight', dpi=300)\n",
    "\n",
    "    ### Plot Brier Score vs lead time (Figure 1)\n",
    "\n",
    "    min_N_samples = 10 # Min number of samples to allow for mean\n",
    "    BSS_agg_init = BSS_agg.mean(dim='init_end')\n",
    "    sns.set_style(\"whitegrid\")\n",
    "\n",
    "    sns.set_context(\"talk\", font_scale=1.5, rc={\"lines.linewidth\": 2.5})\n",
    "\n",
    "    # Get sample size of for each lead time\n",
    "    for_sample = BSS_agg.sel(model='MME').notnull().sum(dim='init_end')\n",
    "    for_sample\n",
    "\n",
    "    # Use threshold of sample size to cut off lead times\n",
    "    max_lead = for_sample.where(for_sample>=min_N_samples,drop=True).fore_time.max().values.astype('timedelta64[D]').astype(int)/7\n",
    "\n",
    "    f = plt.figure(figsize=(15,10))\n",
    "    NM = 10\n",
    "    ax1 = plt.subplot2grid((NM, NM), (0, 0), colspan=NM, rowspan=NM-1)\n",
    "    ax2 = plt.subplot2grid((NM, NM), (NM-1, 0), colspan=NM, rowspan=1)\n",
    "\n",
    "    for cmod in model_plot_order:\n",
    "        if cmod in exl_mods:\n",
    "            continue\n",
    "        # Get model plotting specs\n",
    "        cc = E.model_color[cmod]\n",
    "        cl = E.model_linestyle[cmod]\n",
    "        cm = E.model_marker[cmod]\n",
    "        if cmod=='rasmesrl':\n",
    "            cflag = '*'\n",
    "        else:\n",
    "            cflag = ''\n",
    "\n",
    "        if cmod in ['MME','dampedAnomalyTrend','climatology']:\n",
    "            lw=5\n",
    "        else:\n",
    "            lw = 2\n",
    "\n",
    "        ax1.plot(BSS_agg_init.fore_time.values.astype('timedelta64[D]').astype(int)/7,\n",
    "                BSS_agg_init.sel(model=cmod).values, label=E.model[cmod]['model_label'].rstrip('*')+cflag,\n",
    "                color=cc,\n",
    "                linestyle=cl,\n",
    "                linewidth=lw,\n",
    "                marker=cm)\n",
    "    ax1.legend(loc='lower right', bbox_to_anchor=(1.285, -0.1))\n",
    "    ax1.set_ylabel('Pan-Arctic BS (-)')\n",
    "    ax1.set_xlim([-0.5,max_lead])\n",
    "    ax1.set_xticklabels([''])\n",
    "    ax1.set_ylim([0.02,0.11])\n",
    "\n",
    "    # second axis\n",
    "    ax2.plot(for_sample.fore_time.values.astype('timedelta64[D]').astype(int)/7,\n",
    "             for_sample.values, '-ko')\n",
    "    ax2.set_ylabel('#\\nweeks   ')\n",
    "    ax2.set_xlabel('Lead time (Weeks)')\n",
    "    ax2.set_xlim(ax1.get_xlim());\n",
    "\n",
    "    ax2.set_ylim([0,for_sample.max()+5]);\n",
    "    ax2.set_yticks(np.arange(0,for_sample.max()+5,15));\n",
    "\n",
    "    if not PaperPlots: # only add for website plots\n",
    "        ax1.text(17, 0.025, 'Wayand et al. (in review)', fontsize=12)\n",
    "\n",
    "    # Save to file\n",
    "    f_out = os.path.join(fig_dir,fig_ext+'BSS_by_lead_time_PanArctic_New.png')\n",
    "    f.savefig(f_out,bbox_inches='tight', dpi=300)\n",
    "\n",
    "    ### Plot the IIEE with lead time (SI)\n",
    "\n",
    "    SIP_IIEE.load()\n",
    "\n",
    "    min_N_samples = 10 # Min number of samples to allow for mean\n",
    "    SIP_IIEE_init = SIP_IIEE.mean(dim='init_end')\n",
    "    sns.set_style(\"whitegrid\")\n",
    "\n",
    "    sns.set_context(\"talk\", font_scale=1.5, rc={\"lines.linewidth\": 2.5})\n",
    "\n",
    "    # Get sample size of for each lead time\n",
    "    for_sample = SIP_IIEE.sel(model='MME').notnull().sum(dim='init_end')\n",
    "\n",
    "    # Use threshold of sample size to cut off lead times\n",
    "    max_lead = for_sample.where(for_sample>=min_N_samples,drop=True).fore_time.max().values.astype('timedelta64[D]').astype(int)/7\n",
    "\n",
    "    f = plt.figure(figsize=(10,10))\n",
    "    NM = 10\n",
    "    ax1 = plt.subplot2grid((NM, NM), (0, 0), colspan=NM, rowspan=NM-1)\n",
    "    ax2 = plt.subplot2grid((NM, NM), (NM-1, 0), colspan=NM, rowspan=1)\n",
    "\n",
    "    for cmod in model_plot_order:\n",
    "        if cmod in exl_mods:\n",
    "            continue\n",
    "        # Get model plotting specs\n",
    "        cc = E.model_color[cmod]\n",
    "        cl = E.model_linestyle[cmod]\n",
    "        cm = E.model_marker[cmod]\n",
    "        if cmod=='rasmesrl':\n",
    "            cflag = '*'\n",
    "        else:\n",
    "            cflag = ''\n",
    "        if cmod in ['MME','dampedAnomalyTrend','climatology']:\n",
    "            lw=5\n",
    "        else:\n",
    "            lw = 2\n",
    "\n",
    "        ax1.plot(SIP_IIEE_init.fore_time.values.astype('timedelta64[D]').astype(int)/7,\n",
    "                SIP_IIEE_init.sel(model=cmod).values, label=E.model[cmod]['model_label'].rstrip('*')+cflag,\n",
    "                color=cc,\n",
    "                linestyle=cl,\n",
    "                linewidth=lw,\n",
    "                marker=cm)\n",
    "    ax1.legend(loc='lower right', bbox_to_anchor=(1.4, 0))\n",
    "    ax1.set_ylabel('IIEE (Millions of km$^2$)')\n",
    "    ax1.set_xlim([-0.5,max_lead])\n",
    "    ax1.set_xticklabels([''])\n",
    "\n",
    "    # second axis\n",
    "    ax2.plot(for_sample.fore_time.values.astype('timedelta64[D]').astype(int)/7,\n",
    "             for_sample.values, '-ko')\n",
    "    ax2.set_ylabel('#\\nweeks')\n",
    "    ax2.set_xlabel('Lead time (Weeks)')\n",
    "    ax2.set_xlim(ax1.get_xlim());\n",
    "\n",
    "    ax2.set_ylim([0,for_sample.max()+5]);\n",
    "    ax2.set_yticks(np.arange(0,for_sample.max()+5,15));\n",
    "\n",
    "    if not PaperPlots: # only add for website plots\n",
    "        ax1.text(np.datetime64('2018-01-01'), 0.025, 'Wayand et al. (in review)', fontsize=12)\n",
    "\n",
    "    # Save to file\n",
    "    f_out = os.path.join(fig_dir,fig_ext+'IIEE_by_lead_time_PanArctic.png')\n",
    "    f.savefig(f_out,bbox_inches='tight', dpi=300)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    ### Define DA methods for each model\n",
    "\n",
    "    # copy past info from Table 1\n",
    "    DA_dict = {\n",
    "    'modcansipns_3':'SIC (NG)', \n",
    "    'modcansipns_4':'SIC (NG)',\n",
    "    'ecmwfsipn':'SIC (4DVAR)', \n",
    "    'ecmwf':'SIC (4DVAR)',\n",
    "    'yopp':'SIC (4DVAR)',\n",
    "    'gfdlsipn':'No Sea Ice DA',\n",
    "    'metreofr':'SIC (EnKF)',\n",
    "    'szapirosipn':'No Sea Ice DA',\n",
    "    'ncep-exp-bias-corr':'SIC (NG)',\n",
    "    'noaasipn':'SIC (NG)',\n",
    "    'usnavysipn':'SIC (3DVAR)',\n",
    "    'usnavyncep':'SIC (3DVAR)',\n",
    "    'usnavygofs':'SIC (3DVAR)',\n",
    "    'rasmesrl':'SIC (DI), SIT* (DI)',\n",
    "    'uclsipn':'No Sea Ice DA',\n",
    "    'ukmetofficesipn':'SIC (3DVAR)',\n",
    "    'ukmo':'SIC (3DVAR)',\n",
    "    'ncep':'SIC (NG)',\n",
    "    'kma':'SIC (3DVAR)'\n",
    "    }\n",
    "\n",
    "    # Add models\n",
    "    DA_dict['climatology'] = 'No Sea Ice DA'\n",
    "    DA_dict['dampedAnomalyTrend'] = 'SIC (DI)'\n",
    "    DA_dict['MME'] = 'MME'\n",
    "    DA_dict\n",
    "\n",
    "    DA_options = sorted(list(set(DA_dict.values())))\n",
    "    dict(zip(DA_options,np.arange(len(DA_options))))\n",
    "\n",
    "    DA_options = [DA_options[1],  DA_options[4], DA_options[5], DA_options[7], DA_options[2], DA_options[3], DA_options[6],DA_options[0],] # Reorder from simple to complex\n",
    "    DA_options_i = np.arange(len(DA_options))\n",
    "    DA_options_dict = dict(zip(DA_options,DA_options_i))\n",
    "    DA_options_dict\n",
    "\n",
    "    ### Plot BS by DA method (Figure 4)\n",
    "\n",
    "    # In place a multi lead times\n",
    "\n",
    "    # Lead times to plot\n",
    "    leads2plot = [0,1,2,3,4] # indices\n",
    "    sns.set_style(\"whitegrid\")\n",
    "\n",
    "    sns.set_context(\"talk\", font_scale=1, rc={\"lines.linewidth\": 2.5})\n",
    "    f, axes = plt.subplots(1, 1, figsize=(8, 5))\n",
    "\n",
    "    for cmod in BSS_agg.model.values:\n",
    "        if cmod in DA_dict.keys():\n",
    "            # Get model plotting specs\n",
    "            cc = E.model_color[cmod]\n",
    "            cl = E.model_linestyle[cmod]\n",
    "            cm = E.model_marker[cmod]\n",
    "            if cmod=='MME':\n",
    "                lw=4\n",
    "            else:\n",
    "                lw=2\n",
    "\n",
    "            BSS_init = BSS_agg.sel(model=cmod).isel(fore_time=leads2plot).mean(dim='init_end').load()\n",
    "\n",
    "            #rand_jit = np.random.randint(-100,100)/1000*2\n",
    "            c_x = np.linspace(DA_options_dict[DA_dict[cmod]],\n",
    "                              DA_options_dict[DA_dict[cmod]]+0.9,\n",
    "                              len(leads2plot))\n",
    "            #print(c_x)\n",
    "\n",
    "            axes.plot(c_x, \n",
    "                    BSS_init.values,\n",
    "                    color=cc,\n",
    "                    linestyle='-',\n",
    "                    linewidth=lw,\n",
    "                    marker=cm,\n",
    "                    label=E.model[cmod]['model_label'])\n",
    "        else:\n",
    "            print(cmod,\"not in dict\")\n",
    "    axes.set_xticks(DA_options_i)\n",
    "    axes.set_xticklabels(DA_options, rotation='45', ha='right')\n",
    "\n",
    "    plt.legend(loc='lower right', bbox_to_anchor=(1.36, -.25))\n",
    "    plt.ylabel('Pan-Arctic BS (-)')\n",
    "\n",
    "    if not PaperPlots: # only add for website plots\n",
    "        axes.text(5, 0.085, 'Wayand et al. (in review)', fontsize=12)\n",
    "\n",
    "    # Save to file\n",
    "    f_out = os.path.join(fig_dir,fig_ext+'BSS_week_Multi_by_DA_Type.png')\n",
    "    f.savefig(f_out,bbox_inches='tight', dpi=300)\n",
    "\n",
    "\n",
    "\n",
    "    ### Plot BS by initialization time (Figure 2)\n",
    "\n",
    "    sns.set_style(\"whitegrid\")\n",
    "\n",
    "    for ft_i in [4]:\n",
    "        BSS_agg_fore = BSS_agg.isel(fore_time=ft_i)\n",
    "\n",
    "        sns.set_context(\"talk\", font_scale=1.5, rc={\"lines.linewidth\": 2.5})\n",
    "        f = plt.figure(figsize=(10,10))\n",
    "        for cmod in model_plot_order:\n",
    "            # Get model plotting specs\n",
    "            cc = E.model_color[cmod]\n",
    "            cl = E.model_linestyle[cmod]\n",
    "            cm = E.model_marker[cmod]\n",
    "            if cmod in ['MME','dampedAnomalyTrend','climatology']:\n",
    "                lw=5\n",
    "            else:\n",
    "                lw = 2\n",
    "            if BSS_agg_fore.sel(model=cmod).notnull().sum().values==0:\n",
    "                continue # Don't plot\n",
    "            plt.plot(BSS_agg_fore.init_end.values,\n",
    "                    BSS_agg_fore.sel(model=cmod).values, label=E.model[cmod]['model_label'],\n",
    "                    color=cc,\n",
    "                    linestyle=cl,\n",
    "                    linewidth=lw,\n",
    "                    marker=cm)\n",
    "        plt.legend(loc='lower right', bbox_to_anchor=(1.4, -0.1))\n",
    "        plt.ylabel('Pan-Arctic BS (-)')\n",
    "        plt.xlabel('Initialization date')\n",
    "        #plt.title(BSS_agg_fore.fore_time.values.astype('timedelta64[D]').astype(int))\n",
    "        f.autofmt_xdate()\n",
    "\n",
    "\n",
    "        if not PaperPlots: # only add for website plots\n",
    "            plt.text(np.datetime64('2018-01-01'), 0.16, 'Wayand et al. (in review)', fontsize=12)\n",
    "        # Save to file\n",
    "        f_out = os.path.join(fig_dir,fig_ext+'BSS_by_init_time_'+str(BSS_agg_fore.fore_time.values.astype('timedelta64[D]').astype(int))+'_days.png')\n",
    "        f.savefig(f_out,bbox_inches='tight', dpi=300)\n",
    "        \n",
    "        print(\"Finished Eval_Weekly\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if __name__ == '__main__':\n",
    "    # Start up Client\n",
    "    client = Client(n_workers=8)\n",
    "    \n",
    "    # Call function\n",
    "    Update_Evaluation_Maps()"
   ]
  }
 ],
 "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
}