notebooks/plot_sept_sip.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 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 calendar\n",
    "\n",
    "from esio import EsioData as ed\n",
    "from esio import ice_plot\n",
    "from esio import import_data\n",
    "\n",
    "\n",
    "import dask\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_workers=8)\n",
    "client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_plot_title(cmod, custom_names, E):\n",
    "    if cmod in custom_names:\n",
    "        new_title = custom_names[cmod]\n",
    "    elif cmod in E.model.keys():\n",
    "        new_title = E.model[cmod]['model_label']\n",
    "    else:\n",
    "        new_title = cmod\n",
    "    return new_title"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_subplot_title(cmod, custom_names, E, ax=None):\n",
    "\n",
    "    ax.set_title(get_plot_title(cmod, custom_names, E))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = ed.EsioData.load()\n",
    "# Directories\n",
    "all_models=['rasmesrl']\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':['SIP']}\n",
    "# Some models not to include\n",
    "MME_NO = ['hcmr']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Target time range\n",
    "valid_start = np.datetime64('2018-09-01')\n",
    "valid_end = np.datetime64('2018-09-30')\n",
    "\n",
    "# SIO Report Year_Month\n",
    "cyear = '2018'\n",
    "cmonth = 'June'\n",
    "# cmonth = 'July'\n",
    "# cmonth = 'August'\n",
    "\n",
    "# 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": [
    "# Get recent observations\n",
    "ds_81 = xr.open_mfdataset(E.obs['NSIDC_0081']['sipn_nc']+'_yearly/*.nc', concat_dim='time', autoclose=True, parallel=True)"
   ]
  },
  {
   "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','dampedAnomalyTrend'] # 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)\n",
    "# Load in regional data\n",
    "ds_region = xr.open_dataset(os.path.join(E.grid_dir, 'sio_2016_mask_Update.nc'))"
   ]
  },
  {
   "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 SIP data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "user_sip = [] # Dictionary to store user submitted (regridded) SIP data\n",
    "r_yr_mon = cyear+'_'+cmonth\n",
    "\n",
    "mon_2_int = {v: k for k,v in enumerate(calendar.month_name)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### npssipn (RASM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmod = 'RASM'\n",
    "if r_yr_mon=='2018_June':\n",
    "    sip_f = '/home/disk/sipn/upload/npssipn/forecast/Forecast/RASM_Sep2018_Junreportdata.nc'\n",
    "elif r_yr_mon=='2018_July':\n",
    "    sip_f = '/home/disk/sipn/upload/npssipn/forecast/Forecast/RASM_Sep2018_Julyreportdata.nc'\n",
    "elif r_yr_mon=='2018_August':\n",
    "    sip_f = '/home/disk/sipn/upload/npssipn/forecast/Forecast/RASM_Sep2018_Aug_reportdata.nc'\n",
    "else:\n",
    "    raise ValueError(r_yr_mon,\"not found, add sip_f path.\")\n",
    "\n",
    "# Open\n",
    "ds_user = xr.open_dataset(sip_f)\n",
    "sip_f = None\n",
    "\n",
    "# Regrid/format to sipn\n",
    "ds_user.rename({'lat':'nj', 'long':'ni','longitude':'lon','latitude':'lat'}, inplace=True);\n",
    "ds_user.set_coords(['lat','lon'], inplace=True)\n",
    "# Select SIP\n",
    "da_in = ds_user.SIP_ID\n",
    "mask = ds_user.Grid_Area.notnull()\n",
    "da_in = da_in.where(mask)\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",
    "# 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('2018-08-01')]), \n",
    "                                          dims='model', coords={'model':da_out.model})\n",
    "user_sip.append(da_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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/Modified_CanSIPS_Sep2018_JuneInit_SIP.nc'\n",
    "elif r_yr_mon=='2018_July':\n",
    "    sip_f='/home/disk/sipn/upload/ecsipn/July_2018/Modified_CanSIPS_Sep2018_JulyInit_SIP.nc'\n",
    "elif r_yr_mon=='2018_August':\n",
    "    sip_f='/home/disk/sipn/upload/ecsipn/August_2018/Modified_CanSIPS_Sep2018_AugustInit_SIP.nc'\n",
    "else:\n",
    "    raise ValueError(r_yr_mon,\"not found, add sip_f path.\")   \n",
    "    \n",
    "ds_user = xr.open_dataset(sip_f)\n",
    "sip_f = None\n",
    "ds_user.rename({'lat':'nj', 'lon':'ni','longitude':'lon','latitude':'lat'}, inplace=True);\n",
    "ds_user.set_coords(['lat','lon'], inplace=True)\n",
    "# Select SIP\n",
    "da_in = ds_user.sip.isel(time=0)\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",
    "da_in.plot()\n",
    "plt.figure()\n",
    "da_out.plot()\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('2018-08-01')]), \n",
    "                                          dims='model', coords={'model':da_out.model})\n",
    "user_sip.append(da_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### nrlsipn (NESM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cmod = 'NESM'\n",
    "\n",
    "# if r_yr_mon=='2018_June':\n",
    "#     sip_f = '/home/disk/sipn/nicway/data/model/usnavysipn/forecast/metrics_from_nrl/NESM_Sep2018_Maydata_SIP_common.nc'\n",
    "# elif r_yr_mon=='2018_July':\n",
    "#     sip_f=''\n",
    "# else:\n",
    "#     raise ValueError(r_yr_mon,\"not found, add sip_f path.\")  \n",
    "\n",
    "# ds_user = xr.open_dataset(sip_f)\n",
    "# sip_f = None\n",
    "# ds_user.rename({'Longitude':'lon','Latitude':'lat'}, inplace=True);\n",
    "# # Select SIP\n",
    "# da_in = ds_user.SIP.isel(MT=0)\n",
    "# da_in = da_in / 100 # scale to 0-1\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",
    "# # da_in.plot()\n",
    "# # plt.figure()\n",
    "# # da_out.plot()\n",
    "# # Store in dict of user submited SIP\n",
    "# da_out.coords['model'] = cmod\n",
    "# print(\"Not adding User submited NESM here, Joe wanted us to use our own calculated one\")\n",
    "# #user_sip.append(da_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Test plot\n",
    "# da_in.plot()\n",
    "# plt.figure()\n",
    "# da_out.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sean H."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmod = 'BBGLM\\n(Horvath)'\n",
    "\n",
    "if r_yr_mon=='2018_June':\n",
    "    sip_f = '/home/disk/sipn/upload/seansipn/forecast/2018_it_June_vt_Sept.nc'\n",
    "elif r_yr_mon=='2018_July':\n",
    "    sip_f='/home/disk/sipn/upload/seansipn/forecast/July/2018_it_July_vt_Sept.nc'\n",
    "elif r_yr_mon=='2018_August':\n",
    "    sip_f='/home/disk/sipn/upload/seansipn/forecast/August/2018_it_August_vt_Sept.nc'\n",
    "else:\n",
    "    raise ValueError(r_yr_mon,\"not found, add sip_f path.\")  \n",
    "\n",
    "ds_user = xr.open_dataset(sip_f)\n",
    "da_in = ds_user.mean(dim='valid_time').SIP\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",
    "da_in.plot()\n",
    "plt.figure()\n",
    "da_out.plot()\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('2018-08-13')]), \n",
    "                                          dims='model', coords={'model':da_out.model})\n",
    "user_sip.append(da_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Xioajun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmod = 'Lamont\\n(Yuan et al.)'\n",
    "\n",
    "if r_yr_mon=='2018_June':\n",
    "    sip_f = None # no forecast given\n",
    "elif r_yr_mon=='2018_July':\n",
    "    sip_f='/home/disk/sipn/nicway/data/model/Xioajun/forecast/native/Xiaojun_Yuan_Yuan_Sep2018_Julydata_SIC_Arctic_native.nc'\n",
    "elif r_yr_mon=='2018_August':\n",
    "    sip_f='/home/disk/sipn/nicway/data/model/Xioajun/forecast/native/Yuan_Sep2018_Augustdata_SIC_Arctic_native.nc'    \n",
    "else:\n",
    "    raise ValueError(r_yr_mon,\"not found, add sip_f path.\")  \n",
    "\n",
    "if sip_f:\n",
    "    ds_user = xr.open_dataset(sip_f)\n",
    "    \n",
    "    da_in = ds_user.rename({'X':'lon','Y':'lat'}).isel(T=0).ice\n",
    "    da_in = da_in/100\n",
    "    # SIC to SIP\n",
    "    da_in = (da_in>0.15).where(da_in.notnull())\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",
    "    da_in.plot()\n",
    "    plt.figure()\n",
    "    da_out.plot()\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('2018-07-15')]), \n",
    "                                              dims='model', coords={'model':da_out.model})\n",
    "    user_sip.append(da_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NASA GMAO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmod = 'GMAO'\n",
    "\n",
    "if r_yr_mon=='2018_June':\n",
    "    sip_f = None # no forecast given\n",
    "elif r_yr_mon=='2018_July':\n",
    "    sip_f='/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",
    "if sip_f:\n",
    "    ds_user = xr.open_dataset(sip_f)\n",
    "    \n",
    "    ds_user = ds_user.set_coords(['LAT','LON'])\n",
    "    \n",
    "    da_in = ds_user.rename({'LON':'lon','LAT':'lat','SIP':'sip'}).sip\n",
    "    \n",
    "    da_in = da_in/100 # 100% to 0-1\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",
    "    da_in.plot()\n",
    "    plt.figure()\n",
    "    da_out.plot()\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_sip.append(da_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Merge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "usr_SIP = xr.concat(user_sip, dim='model')\n",
    "usr_SIP = usr_SIP.rename({'nj':'x', 'ni':'y'})\n",
    "usr_SIP.name = 'SIP'\n",
    "usr_SIP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loop through SIPN SIC data and plot those models that have forecasts through target month\n",
    "- Use most recent init_time (if multiple)\n",
    "- Use monthly sept mean SIC for each ensemble member to calculate SIP"
   ]
  },
  {
   "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 13th of June)\n",
    "last_init_day = np.datetime64(datetime.datetime(int(cyear),mon_2_int[cmonth],13))\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 = 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": [
    "# 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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_inits_used = {}\n",
    "for cvar in variables:\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):\n",
    "            print(\"\")\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",
    "            start_time = timeit.default_timer()\n",
    "            ds_model = xr.open_mfdataset(model_forecast, \n",
    "                        chunks={'fore_time': 1, '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",
    "            print(\"    Loading took  \", (timeit.default_timer() - start_time)/60, \" minutes.\")\n",
    "\n",
    "            # Select var of interest (if available)\n",
    "            if cvar in ds_model.variables:\n",
    "#                     print('found ',cvar)\n",
    "                ds_model = ds_model[cvar]\n",
    "            else:\n",
    "                print('    cvar not found.')\n",
    "                continue\n",
    "                  \n",
    "            # Select init times\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",
    "            print(\"    Init time used is\",ds_model.init_time.values)\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",
    "            # Average over for_time\n",
    "            ds_model = ds_model.mean(dim='fore_time')\n",
    "            \n",
    "            start_time = timeit.default_timer()\n",
    "#                 print(\"Found data for model \", cmod, \". Plotting...\")    \n",
    "            if metric=='mean': # Calc ensemble mean\n",
    "                ds_model = ds_model.mean(dim='ensemble')\n",
    "            elif metric=='SIP': # Calc probability\n",
    "                ok_ens = ((ds_model.notnull().sum(dim='x').sum(dim='y'))>0) # select ensemble members with any data\n",
    "                ds_model = ((ds_model.where(ok_ens, drop=True)>=0.15) ).mean(dim='ensemble').where(ds_model.isel(ensemble=0).notnull())\n",
    "            elif metric=='anomaly': # Calc anomaly in reference to mean observed 1980-2010\n",
    "                ds_model = ds_model.mean(dim='ensemble') - da_obs_mean\n",
    "                # Add back lat/long (get dropped because of round off differences)\n",
    "                ds_model['lat'] = da_obs_mean.lat\n",
    "                ds_model['lon'] = da_obs_mean.lon\n",
    "            else:\n",
    "                raise ValueError('metric not implemented')\n",
    "            print(\"    Calc metric took  \", (timeit.default_timer() - start_time)/60, \" minutes.\")\n",
    "\n",
    "            # Build MME\n",
    "            if 'ensemble' in ds_model:\n",
    "                ds_model = ds_model.drop('ensemble')\n",
    "            if cmod not in MME_NO: # Exclude some models (bad) from MME\n",
    "                ds_model.coords['model'] = cmod\n",
    "                ds_model = ds_model.expand_dims('model')\n",
    "                ds_model.coords['init_date'] = xr.DataArray(ds_model.init_time.values, dims='model', coords={'model':ds_model.model})\n",
    "                MME_list.append(ds_model)\n",
    "                # Store init time\n",
    "                model_inits_used[cmod] = ds_model.init_time.values\n",
    "                print('    Added ',cmod,' to MME.')\n",
    "\n",
    "# Done with current it\n",
    "print(\"Took \", (timeit.default_timer() - start_time)/60, \" minutes.\")\n",
    "\n",
    "# Concat over all models\n",
    "ds_MME = xr.concat(MME_list, dim='model')\n",
    "ds_MME.name = 'SIP'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert to human readable\n",
    "inits_hr = {}\n",
    "for cmod, idate in model_inits_used.items():\n",
    "    inits_hr[get_plot_title(cmod, custom_names, E)] = pd.to_datetime(model_inits_used[cmod]).strftime('%Y-%m-%d')[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inits_hr"
   ]
  },
  {
   "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": [
    "# Take mean over init_time (different models init different times. mean doesn't change values)\n",
    "model_inits = ds_MME.init_time # Save init_times\n",
    "models_list = ds_MME.model\n",
    "ds_MME = ds_MME.mean(dim='init_time')\n",
    "ds_MME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Merge user SIP with SPIN SIP\n",
    "ds_SIP_All = xr.concat([ds_MME,usr_SIP], dim='model')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load into memory\n",
    "ds_SIP_All.load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save to disk\n",
    "if not os.path.exists('/home/disk/sipn/nicway/data/model/SIO/'+cyear+'/'+cmonth):\n",
    "    os.makedirs('/home/disk/sipn/nicway/data/model/SIO/'+cyear+'/'+cmonth)\n",
    "ds_SIP_All.to_netcdf('/home/disk/sipn/nicway/data/model/SIO/'+cyear+'/'+cmonth+'/SIP.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load in last results\n",
    "# ds_SIP_All = xr.open_mfdataset('/home/disk/sipn/nicway/data/model/SIO/'+cyear+'/'+cmonth+'/SIP.nc').SIP\n",
    "# cvar = 'sic'\n",
    "# metric = 'SIP'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_SIP_All"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Damped persistence trend has slightly different lat long (10^-6) so use lat lon from any other model here\n",
    "# print(\"hack for bad lat lon in peristence, fix upstream!\")\n",
    "# ds_SIP_All['lat'] = ds_SIP_All.lat.sel(model='gfdlsipn').drop('model')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_MME_avg = ds_SIP_All.mean(dim='model')\n",
    "ds_MME_std = ds_SIP_All.std(dim='model')"
   ]
  },
  {
   "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=10)\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",
    "\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",
    "\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=4, Nplots=15, extent=central_extent, central_longitude=0)\n",
    "\n",
    "for (i, cmod) in enumerate(ds_SIP_All.model.values):\n",
    "    print(cmod)\n",
    "    # Plot\n",
    "    add_subplot_title(cmod, custom_names, E, ax=axes[i])\n",
    "    p = ds_SIP_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",
    "    add_subplot_title(cmod, custom_names, E, ax=axes[i])\n",
    "\n",
    "\n",
    "# MME Mean\n",
    "i = i + 1\n",
    "\n",
    "\n",
    "pmme = ds_MME_avg.plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                                  transform=ccrs.PlateCarree(),\n",
    "                                  add_colorbar=False, \n",
    "                                  cmap=cmap_c,vmin=c_vmin, vmax=c_vmax)\n",
    "axes[i].set_title('Mean')\n",
    "\n",
    "\n",
    "# MME Standard deviation\n",
    "i = i + 1\n",
    "\n",
    "\n",
    "pmme = ds_MME_std.plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                                  transform=ccrs.PlateCarree(),\n",
    "                                  add_colorbar=False, \n",
    "                                  cmap=cmap_c,vmin=c_vmin, vmax=c_vmax)\n",
    "axes[i].set_title('σ')\n",
    "    \n",
    "\n",
    "# Make pretty\n",
    "f.subplots_adjust(bottom=0.05)\n",
    "# cbar_ax = f.add_axes([0.85, 0.15, 0.05, 0.7])\n",
    "cbar_ax = f.add_axes([0.25, 0.001, .5, 0.04]) #  [left, bottom, width, height] w\n",
    "if p:\n",
    "    cbar = f.colorbar(p, cax=cbar_ax, label=c_label, orientation='horizontal')\n",
    "    if metric=='anomaly':\n",
    "        cbar.set_ticks(np.arange(-1,1.1,0.2))\n",
    "    else:\n",
    "        cbar.set_ticks(np.arange(0,1.1,0.2))\n",
    "\n",
    "# Save to file\n",
    "f_out = os.path.join(fig_dir,'panArctic_'+metric+'_'+runType+'_'+cyear+'_'+cmonth+'.png')\n",
    "f.savefig(f_out,bbox_inches='tight', dpi=300)\n",
    "f_out = os.path.join(fig_dir,'panArctic_'+metric+'_'+runType+'_'+cyear+'_'+cmonth+'_lowRES.png')\n",
    "f.savefig(f_out,bbox_inches='tight', dpi=90)\n",
    "print(\"saved \", f_out)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare past mean forecasts\n",
    "ds_SIP_June = xr.open_mfdataset('/home/disk/sipn/nicway/data/model/SIO/'+'2018'+'/'+'June'+'/SIP.nc')\n",
    "ds_SIP_July = xr.open_mfdataset('/home/disk/sipn/nicway/data/model/SIO/'+'2018'+'/'+'July'+'/SIP.nc')\n",
    "ds_SIP_August = xr.open_mfdataset('/home/disk/sipn/nicway/data/model/SIO/'+'2018'+'/'+'August'+'/SIP.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# July\n",
    "# common_models = [x for x in ds_SIP_June.model.values if x not in ['awispin', 'nicosipn','ecmwfsipn']] # ecmwf not july, awi and nico no June\n",
    "# August\n",
    "common_models = [x for x in ds_SIP_August.model.values if x not in ['ecmwfsipn','GMAO']] # ecmwf and GMAO didn't submit aug"
   ]
  },
  {
   "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": [
    "# 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=1, Nplots=20, extent=central_extent, central_longitude=0)\n",
    "\n",
    "i = 0\n",
    "pmme1 = ds_SIP_June.mean(dim='model').SIP.plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                                  transform=ccrs.PlateCarree(),\n",
    "                                  add_colorbar=False, \n",
    "                                  cmap=cmap_c,vmin=c_vmin, vmax=c_vmax)\n",
    "axes[i].set_title('June Multi\\nModel Mean')\n",
    "\n",
    "i = 1\n",
    "pmme2 = ds_SIP_July.mean(dim='model').SIP.plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                                  transform=ccrs.PlateCarree(),\n",
    "                                  add_colorbar=False, \n",
    "                                  cmap=cmap_c,vmin=c_vmin, vmax=c_vmax)\n",
    "axes[i].set_title('July Multi\\nModel Mean')\n",
    "\n",
    "i = 2\n",
    "pmme3 = ds_SIP_August.mean(dim='model').SIP.plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                                  transform=ccrs.PlateCarree(),\n",
    "                                  add_colorbar=False, \n",
    "                                  cmap=cmap_c,vmin=c_vmin, vmax=c_vmax)\n",
    "axes[i].set_title('August Multi\\nModel Mean')\n",
    "\n",
    "\n",
    "i = 3\n",
    "cmap_diff = matplotlib.colors.LinearSegmentedColormap.from_list(\"\", [\"blue\",\"white\",\"red\",])\n",
    "cmap_diff.set_bad(color = 'lightgrey')\n",
    "pmme3 = (ds_SIP_August.sel(model=common_models).mean(dim='model')-ds_SIP_July.sel(model=common_models).mean(dim='model')).SIP.plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                                  transform=ccrs.PlateCarree(),\n",
    "                                  add_colorbar=False, \n",
    "                                  cmap=cmap_diff,vmin=-.2, vmax=.2)\n",
    "axes[i].set_title('Difference\\n(August - July)')\n",
    "\n",
    "\n",
    "f.subplots_adjust(bottom=0.02)\n",
    "cbar_ax = f.add_axes([0.25, 0.001, .5, 0.08]) #  [left, bottom, width, height] w\n",
    "cbar = f.colorbar(pmme2, cax=cbar_ax, label=c_label, orientation='horizontal')\n",
    "cbar.set_ticks(np.arange(0,1.1,0.2))\n",
    "\n",
    "cbar_ax_2 = f.add_axes([.93, 0.08, .025, 0.75]) #  [left, bottom, width, height] w\n",
    "cbar2 = f.colorbar(pmme3, cax=cbar_ax_2, label='SIP Change (-)', orientation='vertical')\n",
    "\n",
    "# Save to file\n",
    "f_out = os.path.join(fig_dir,'panArctic_'+metric+'_'+runType+'_'+cyear+'_Aug_July_diff.png')\n",
    "f.savefig(f_out,bbox_inches='tight', dpi=300)\n",
    "print(\"saved \", f_out)\n",
    "f_out = os.path.join(fig_dir,'panArctic_'+metric+'_'+runType+'_'+cyear+'_Aug_July_diff_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": [
    "# 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=10)\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",
    "\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",
    "# 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=4, Nplots=16, extent=central_extent, central_longitude=0)\n",
    "\n",
    "for (i, cmod) in enumerate(ds_SIP_All.model.values):\n",
    "    print(cmod)\n",
    "    # Plot\n",
    "    add_subplot_title(cmod, custom_names, E, ax=axes[i])\n",
    "    p = ds_SIP_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",
    "    add_subplot_title(cmod, custom_names, E, ax=axes[i])\n",
    "    \n",
    "    # Add init date\n",
    "    axes[i].annotate(pd.to_datetime(ds_SIP_All.sel(model=cmod).init_date.values).strftime('%m/%d'), xy=(.6, -.13), xycoords='axes fraction', )\n",
    "    \n",
    "\n",
    "\n",
    "# MME Mean\n",
    "i = i + 1\n",
    "# f.delaxes(axes[i])\n",
    "# i = i + 1\n",
    "# f.delaxes(axes[i])\n",
    "# i = i + 1\n",
    "\n",
    "pmme = ds_MME_avg.plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                                  transform=ccrs.PlateCarree(),\n",
    "                                  add_colorbar=False, \n",
    "                                  cmap=cmap_c,vmin=c_vmin, vmax=c_vmax)\n",
    "axes[i].set_title('Mean')\n",
    "\n",
    "\n",
    "# MME Standard deviation\n",
    "i = i + 1\n",
    "\n",
    "\n",
    "pmme = ds_MME_std.plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                                  transform=ccrs.PlateCarree(),\n",
    "                                  add_colorbar=False, \n",
    "                                  cmap=cmap_c,vmin=c_vmin, vmax=c_vmax)\n",
    "axes[i].set_title('σ')\n",
    "\n",
    "# Difference from last SIO\n",
    "i = i + 1\n",
    "cmap_diff = matplotlib.colors.LinearSegmentedColormap.from_list(\"\", [\"blue\",\"white\",\"red\",])\n",
    "cmap_diff.set_bad(color = 'lightgrey')\n",
    "pmme3 = (ds_SIP_August.sel(model=common_models).mean(dim='model')-ds_SIP_July.sel(model=common_models).mean(dim='model')).SIP.plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                                  transform=ccrs.PlateCarree(),\n",
    "                                  add_colorbar=False, \n",
    "                                  cmap=cmap_diff,vmin=-.2, vmax=.2)\n",
    "axes[i].set_title('Model-mean difference\\n(August - July)')\n",
    "    \n",
    "\n",
    "# Make pretty\n",
    "f.subplots_adjust(bottom=0.05)\n",
    "# cbar_ax = f.add_axes([0.85, 0.15, 0.05, 0.7])\n",
    "cbar_ax = f.add_axes([0.25, 0.001, .5, 0.04]) #  [left, bottom, width, height] w\n",
    "if p:\n",
    "    cbar = f.colorbar(p, cax=cbar_ax, label=c_label, orientation='horizontal')\n",
    "    if metric=='anomaly':\n",
    "        cbar.set_ticks(np.arange(-1,1.1,0.2))\n",
    "    else:\n",
    "        cbar.set_ticks(np.arange(0,1.1,0.2))\n",
    "        \n",
    "\n",
    "# cbar_ax_2 = f.add_axes([0.93, 0.07, .025, .15]) #  [left, bottom, width, height] w\n",
    "cbar_ax_2 = f.add_axes([0.715, 0.07, .025, .15]) #  [left, bottom, width, height] w\n",
    "\n",
    "cbar2 = f.colorbar(pmme3, cax=cbar_ax_2, label='SIP Change (-)', orientation='vertical')\n",
    "\n",
    "\n",
    "# Save to file\n",
    "f_out = os.path.join(fig_dir,'panArctic_'+metric+'_'+runType+'_'+cyear+'_'+cmonth+'_withDiff.png')\n",
    "f.savefig(f_out,bbox_inches='tight', dpi=300)\n",
    "f_out = os.path.join(fig_dir,'panArctic_'+metric+'_'+runType+'_'+cyear+'_'+cmonth+'_withDiff_lowRES.png')\n",
    "f.savefig(f_out,bbox_inches='tight', dpi=90)\n",
    "print(\"saved \", f_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## For Post-season, this figure overlays the observed ice iedge (>15%) on each model's SIP + std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get Observed mean sea ice edge location mean(SIC in sept) > 0.15\n",
    "cvar = 'sic'\n",
    "metric = 'SIP'\n",
    "obs_sept_SIP = ds_81[cvar].sel(time=slice(cyear+'-09-01',cyear+'-09-30')).mean(dim='time')>=0.15\n",
    "# Fill in pole hole\n",
    "obs_sept_SIP = obs_sept_SIP.where(obs_sept_SIP.hole_mask==0, other=1)\n",
    "obs_sept_SIP.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_ds = ds_SIP_June # pick month\n",
    "c_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TESTING CONTOUR PLOT\n",
    "# (f, axes) = ice_plot.polar_axis(extent=central_extent, central_longitude=0)\n",
    "\n",
    "\n",
    "# pmme = c_ds[metric].mean(dim='model').plot.pcolormesh(ax=axes, x='lon', y='lat', \n",
    "#                                   transform=ccrs.PlateCarree(),\n",
    "#                                   add_colorbar=False, \n",
    "#                                   cmap=cmap_c,vmin=c_vmin, vmax=c_vmax)\n",
    "\n",
    "\n",
    "# # pmme = obs_sept_SIP.plot.pcolormesh(ax=axes, x='lon', y='lat', \n",
    "# #                                   transform=ccrs.PlateCarree(),\n",
    "# #                                   add_colorbar=False, \n",
    "# #                                   cmap='binary',vmin=c_vmin, vmax=c_vmax,\n",
    "# #                                   alpha=0.1)\n",
    "\n",
    "\n",
    "\n",
    "# po = obs_sept_SIP.plot.contour(ax=axes, x='xm', y='ym',\n",
    "#                               transform=ccrs.NorthPolarStereo(central_longitude=-45),\n",
    "#                               colors=('k'),\n",
    "#                               linewidths=[3],\n",
    "#                               levels=[0.5]) #, label='Median ice edge 1981-2010')\n",
    "# ice_plot.remove_small_contours(po, thres=10**6)"
   ]
  },
  {
   "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=10)\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",
    "\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",
    "# 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=4, Nplots=14, extent=central_extent, central_longitude=0)\n",
    "\n",
    "for (i, cmod) in enumerate(c_ds.model.values):\n",
    "    print(cmod)\n",
    "    # Plot\n",
    "    add_subplot_title(cmod, custom_names, E, ax=axes[i])\n",
    "    p = c_ds[metric].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",
    "    \n",
    "    \n",
    "    po = obs_sept_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, custom_names, E, ax=axes[i])\n",
    "    \n",
    "    # Add init date\n",
    "    #axes[i].annotate(pd.to_datetime(c_ds.sel(model=cmod).init_date.values).strftime('%m/%d'), xy=(.6, -.13), xycoords='axes fraction', )\n",
    "    \n",
    "\n",
    "\n",
    "# MME Mean\n",
    "i = i + 1\n",
    "# f.delaxes(axes[i])\n",
    "# i = i + 1\n",
    "# f.delaxes(axes[i])\n",
    "# i = i + 1\n",
    "\n",
    "pmme = c_ds[metric].mean(dim='model').plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                                  transform=ccrs.PlateCarree(),\n",
    "                                  add_colorbar=False, \n",
    "                                  cmap=cmap_c,vmin=c_vmin, vmax=c_vmax)\n",
    "\n",
    "po = obs_sept_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",
    "axes[i].set_title('Mean')\n",
    "\n",
    "\n",
    "\n",
    "# MME Standard deviation\n",
    "i = i + 1\n",
    "\n",
    "\n",
    "pmme = c_ds[metric].std(dim='model').plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                                  transform=ccrs.PlateCarree(),\n",
    "                                  add_colorbar=False, \n",
    "                                  cmap=cmap_c,vmin=c_vmin, vmax=c_vmax)\n",
    "\n",
    "po = obs_sept_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",
    "axes[i].set_title('σ')\n",
    "\n",
    "                                                    \n",
    "# Make pretty\n",
    "# f.subplots_adjust(bottom=0.05)\n",
    "# cbar_ax = f.add_axes([0.85, 0.15, 0.05, 0.7])\n",
    "cbar_ax = f.add_axes([0.35, 0.20, .5, 0.04]) #  [left, bottom, width, height] w\n",
    "if p:\n",
    "    cbar = f.colorbar(p, cax=cbar_ax, label=c_label, orientation='horizontal')\n",
    "    if metric=='anomaly':\n",
    "        cbar.set_ticks(np.arange(-1,1.1,0.2))\n",
    "    else:\n",
    "        cbar.set_ticks(np.arange(0,1.1,0.2))\n",
    "        \n",
    "\n",
    "# cbar_ax_2 = f.add_axes([0.93, 0.07, .025, .15]) #  [left, bottom, width, height] w\n",
    "# cbar_ax_2 = f.add_axes([0.715, 0.07, .025, .15]) #  [left, bottom, width, height] w\n",
    "\n",
    "# cbar2 = f.colorbar(pmme3, cax=cbar_ax_2, label='SIP Change (-)', orientation='vertical')\n",
    "\n",
    "\n",
    "# Save to file\n",
    "f_out = os.path.join(fig_dir,'panArctic_'+metric+'_'+runType+'_'+cyear+'_'+cmonth+'_POSTSEASON.png')\n",
    "f.savefig(f_out,bbox_inches='tight', dpi=300)\n",
    "f_out = os.path.join(fig_dir,'panArctic_'+metric+'_'+runType+'_'+cyear+'_'+cmonth+'_POSTSEASON_lowRES.png')\n",
    "f.savefig(f_out,bbox_inches='tight', dpi=90)\n",
    "print(\"saved \", f_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now plot BSS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Force round off lat/lon differences to be the same\n",
    "obs_sept_SIP['lat'] = c_ds.lat\n",
    "obs_sept_SIP['lon'] = c_ds.lon\n",
    "ds_region['lat'] = c_ds.lat\n",
    "ds_region['lon'] = c_ds.lon\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mask only to Arctic seas (importand to do this before spatially averaging!)\n",
    "c_ds = c_ds.where(ds_region.mask.isin(ds_region.ocean_regions))\n",
    "obs_sept_SIP = obs_sept_SIP.where(ds_region.mask.isin(ds_region.ocean_regions))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('ticks')\n",
    "sns.set_context(\"talk\", font_scale=.8, rc={\"lines.linewidth\": 2.5})\n",
    "# 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=10)\n",
    "        #cmap_c = plt.get_cmap('jet') \n",
    "        cmap_c.set_bad(color = 'lightgrey')\n",
    "        c_label = 'Brier Skill Score (-)'\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",
    "\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 = 'Brier Skill Score (-)'\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",
    "# 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=14, extent=central_extent, central_longitude=0)\n",
    "\n",
    "for (i, cmod) in enumerate(c_ds.model.values):\n",
    "    print(cmod)\n",
    "    # Plot\n",
    "    add_subplot_title(cmod, custom_names, E, ax=axes[i])\n",
    "    c_BSS = ((c_ds[metric].sel(model=cmod)-obs_sept_SIP)**2)\n",
    "    p = c_BSS.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",
    "    \n",
    "    \n",
    "    po = obs_sept_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, custom_names, E, ax=axes[i])\n",
    "    \n",
    "    # Add init date\n",
    "    \n",
    "\n",
    "    axes[i].annotate('{0:.3f}'.format((c_BSS*ds_region.area).sum(dim=['x','y']).values/(9.4*10**6)),\n",
    "                     xy=(.25, -.16), \n",
    "                     xycoords='axes fraction', )\n",
    "    \n",
    "\n",
    "\n",
    "# MME Mean\n",
    "i = i + 1\n",
    "# f.delaxes(axes[i])\n",
    "# i = i + 1\n",
    "# f.delaxes(axes[i])\n",
    "# i = i + 1\n",
    "c_BSS = ((c_ds[metric].mean(dim='model')-obs_sept_SIP)**2)\n",
    "# Mask only to Arctic seas\n",
    "pmme = c_BSS.plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "                                  transform=ccrs.PlateCarree(),\n",
    "                                  add_colorbar=False, \n",
    "                                  cmap=cmap_c,vmin=c_vmin, vmax=c_vmax)\n",
    "\n",
    "po = obs_sept_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",
    "axes[i].set_title('Mean')\n",
    "\n",
    "axes[i].annotate('{0:.3f}'.format((c_BSS*ds_region.area).sum(dim=['x','y']).values/(9.4*10**6)),\n",
    "                     xy=(.25, -.16), \n",
    "                     xycoords='axes fraction', )\n",
    "\n",
    "\n",
    "\n",
    "# # MME Standard deviation\n",
    "# i = i + 1\n",
    "\n",
    "\n",
    "# pmme = c_ds[metric].std(dim='model').plot.pcolormesh(ax=axes[i], x='lon', y='lat', \n",
    "#                                   transform=ccrs.PlateCarree(),\n",
    "#                                   add_colorbar=False, \n",
    "#                                   cmap=cmap_c,vmin=c_vmin, vmax=c_vmax)\n",
    "\n",
    "# po = obs_sept_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",
    "# axes[i].set_title('σ')\n",
    "\n",
    "                                                    \n",
    "# Make pretty\n",
    "f.subplots_adjust(bottom=0.05)\n",
    "# cbar_ax = f.add_axes([0.85, 0.15, 0.05, 0.7])\n",
    "cbar_ax = f.add_axes([0.25, 0.001, .5, 0.04]) #  [left, bottom, width, height] w\n",
    "if p:\n",
    "    cbar = f.colorbar(p, cax=cbar_ax, label=c_label, orientation='horizontal')\n",
    "    if metric=='anomaly':\n",
    "        cbar.set_ticks(np.arange(-1,1.1,0.2))\n",
    "    else:\n",
    "        cbar.set_ticks(np.arange(0,1.1,0.2))\n",
    "        \n",
    "\n",
    "# cbar_ax_2 = f.add_axes([0.93, 0.07, .025, .15]) #  [left, bottom, width, height] w\n",
    "# cbar_ax_2 = f.add_axes([0.715, 0.07, .025, .15]) #  [left, bottom, width, height] w\n",
    "\n",
    "# cbar2 = f.colorbar(pmme3, cax=cbar_ax_2, label='SIP Change (-)', orientation='vertical')\n",
    "\n",
    "\n",
    "# Save to file\n",
    "f_out = os.path.join(fig_dir,'panArctic_BSS_'+runType+'_'+cyear+'_'+cmonth+'_POSTSEASON.png')\n",
    "f.savefig(f_out,bbox_inches='tight', dpi=300)\n",
    "f_out = os.path.join(fig_dir,'panArctic_BSS_'+runType+'_'+cyear+'_'+cmonth+'_POSTSEASON_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": [
    "# Combine each month datset\n",
    "ds_SIP_June['init'] = 'June'\n",
    "ds_SIP_July['init'] = 'July'\n",
    "ds_SIP_August['init'] = 'August'\n",
    "ds_SIP_ALL = xr.concat([ds_SIP_June,ds_SIP_July,ds_SIP_August.drop('init_date')],dim='init')\n",
    "ds_SIP_ALL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(((c_ds[metric]-obs_sept_SIP)**2)*ds_region.area).sum(dim=['x','y']).values/(9.4*10**6)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_SIP_ALL.model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define colors/lines etce for each model to use in plots\n",
    "cmap_c = itertools.cycle(sns.color_palette(\"Paired\", len(E.model.keys()) ))\n",
    "linecycler = itertools.cycle([\"-\",\"--\",\"-.\",\":\",\"--\"])\n",
    "markercycler = itertools.cycle([\"*\",\"o\",\"s\",\"v\",\"X\"])\n",
    "mod_cc = {}\n",
    "mod_ls = {}\n",
    "mod_mk = {}\n",
    "for (i, cmod) in enumerate(ds_SIP_ALL.model.values):\n",
    "    mod_cc[cmod] = next(cmap_c)\n",
    "    mod_ls[cmod] = next(linecycler)\n",
    "    mod_mk[cmod] = next(markercycler)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Inits_with_data = ds_SIP_ALL[metric].sel(model='Lamont\\n(Yuan et al.)').notnull().sum(dim=['x','y'])>0\n",
    "\n",
    "X = ((((ds_SIP_ALL[metric].sel(model='Lamont\\n(Yuan et al.)') - obs_sept_SIP)**2)*ds_region.area).sum(dim=['x','y'])/(9.4*10**6)).where(Inits_with_data)\n",
    "X.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot for June, July, August, the mean BSS scores\n",
    "sns.set_context(\"paper\", font_scale=2, rc={\"lines.linewidth\": 2.5})\n",
    "sns.set_style('whitegrid')\n",
    "(f,ax) = plt.subplots(1,1)\n",
    "f.set_size_inches((5,5))\n",
    "for (i, cmod) in enumerate(ds_SIP_ALL.model.values):\n",
    "    # Some models have msising inits, so remove those here (otherwise they sum to 0)\n",
    "    Inits_with_data = ds_SIP_ALL[metric].sel(model=cmod).notnull().sum(dim=['x','y'])>0\n",
    "\n",
    "    \n",
    "    c_BSS = ((((ds_SIP_ALL[metric].sel(model=cmod) - obs_sept_SIP)**2)*ds_region.area).sum(dim=['x','y'])/(9.4*10**6)).where(Inits_with_data)\n",
    "    ax.plot(c_BSS.init, c_BSS.values, \n",
    "            color=mod_cc[cmod],\n",
    "            marker=mod_mk[cmod],\n",
    "            linestyle=mod_ls[cmod],\n",
    "            markersize=10,\n",
    "            label=get_plot_title(cmod, custom_names, E))\n",
    "# MME\n",
    "MME_BSS = (((ds_SIP_ALL[metric].mean(dim='model') - obs_sept_SIP)**2)*ds_region.area).sum(dim=['x','y'])/(9.4*10**6)\n",
    "ax.plot(MME_BSS.init, MME_BSS.values, label='Mean', linewidth=8)\n",
    "f.legend(loc='lower right', bbox_to_anchor=(1.5, -0.04))\n",
    "ax.set_ylim([0,0.18])\n",
    "ax.set_ylabel('Brier Skill Score (-)')\n",
    "ax.set_xlabel('Initialization month')\n",
    "\n",
    "# Save to file\n",
    "f_out = os.path.join(fig_dir,'panArctic_Mean_BSS_'+cyear+'_POSTSEASON.png')\n",
    "f.savefig(f_out,bbox_inches='tight', dpi=300)\n",
    "f_out = os.path.join(fig_dir,'panArctic_Mean_BSS_'+cyear+'_POSTSEASON_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": [
    "# Need to add check -- Compare June to July models (Make sure they are different! (unless they are supposed to be unchanged))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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
}