notebooks/Calc_Weekly_Model_Metrics.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",
    "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",
    "warnings.simplefilter(action='ignore', category=RuntimeWarning) # not good to supress but they divide by nan are annoying\n",
    "#warnings.simplefilter(action='ignore', category=UserWarning) # https://github.com/pydata/xarray/issues/2273\n",
    "import json\n",
    "from esio import EsioData as ed\n",
    "from esio import ice_plot\n",
    "from esio import import_data\n",
    "import subprocess\n",
    "import dask\n",
    "from dask.distributed import Client\n",
    "import timeit\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": [
    "#client = Client()\n",
    "#client\n",
    "dask.config.set(scheduler='threads')  # overwrite default with threaded scheduler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#def Update_PanArctic_Maps():\n",
    "# Plotting Info\n",
    "runType = 'forecast'\n",
    "variables = ['sic']\n",
    "metrics_all = {'sic':['anomaly','mean','SIP'], 'hi':['mean']}\n",
    "#metrics_all = {'sic':['SIP']}\n",
    "updateAll = False\n",
    "\n",
    "# Define Init Periods here, spaced by 7 days (aprox a week)\n",
    "# Now\n",
    "cd = datetime.datetime.now()\n",
    "cd = datetime.datetime(cd.year, cd.month, cd.day) # Set hour min sec to 0. \n",
    "# Hardcoded start date (makes incremental weeks always the same)\n",
    "start_t = datetime.datetime(1950, 1, 1) # datetime.datetime(1950, 1, 1)\n",
    "# Params for this plot\n",
    "Ndays = 7 # time period to aggregate maps to (default is 7)\n",
    "\n",
    "#Npers = 40 # number of periods agg (from current date)\n",
    "init_start_date = np.datetime64('2018-01-01')\n",
    "\n",
    "init_slice = np.arange(start_t, cd, datetime.timedelta(days=Ndays)).astype('datetime64[ns]')\n",
    "# init_slice = init_slice[-Npers:] # Select only the last Npers of periods (weeks) since current date\n",
    "init_slice = init_slice[init_slice>=init_start_date] # Select only the inits after init_start_date\n",
    "\n",
    "print(init_slice[0],init_slice[-1])\n",
    "print('')\n",
    "\n",
    "# Forecast times\n",
    "weeks = pd.to_timedelta(np.arange(0,52,1), unit='W')\n",
    "#months = pd.to_timedelta(np.arange(2,12,1), unit='M')\n",
    "#years = pd.to_timedelta(np.arange(1,2), unit='Y') - np.timedelta64(1, 'D') # need 364 not 365\n",
    "#slices = weeks.union(months).union(years).round('1d')\n",
    "da_slices = xr.DataArray(weeks, dims=('fore_time'))\n",
    "da_slices.fore_time.values.astype('timedelta64[D]')\n",
    "print(da_slices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#############################################################\n",
    "# Load in Observed and non-dynamic model Data\n",
    "#############################################################\n",
    "\n",
    "E = ed.EsioData.load()\n",
    "mod_dir = E.model_dir\n",
    "\n",
    "# Get median ice edge by DOY\n",
    "median_ice_fill = xr.open_mfdataset(os.path.join(E.obs_dir, 'NSIDC_0051', 'agg_nc', 'ice_edge.nc')).sic\n",
    "# Get mean sic by DOY\n",
    "mean_1980_2010_sic = xr.open_dataset(os.path.join(E.obs_dir, 'NSIDC_0051', 'agg_nc', 'mean_1980_2010_sic.nc')).sic\n",
    "# Get average sip by DOY\n",
    "mean_1980_2010_SIP = xr.open_dataset(os.path.join(E.obs_dir, 'NSIDC_0051', 'agg_nc', 'hist_SIP_1980_2010.nc')).sip    \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",
    "# Define models to plot\n",
    "models_2_plot = list(E.model.keys())\n",
    "models_2_plot = [x for x in models_2_plot if x not in ['piomas','MME','MME_NEW','uclsipn','hcmr']] # remove some models\n",
    "models_2_plot = [x for x in models_2_plot if E.icePredicted[x]] # Only predictive models\n",
    "# models_2_plot = ['dampedAnomalyTrend']\n",
    "models_2_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def is_in_time_range(x):\n",
    "    \n",
    "#     if x.sel(init_time=slice(time_bds[0],time_bds[1])).init_time.size>0: # We have some time in the time range\n",
    "#         return x\n",
    "#     else:\n",
    "#         return []\n",
    "# time_bds = [init_slice[0],init_slice[-1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###########################################################\n",
    "#          Loop through each dynamical model              #\n",
    "###########################################################\n",
    "\n",
    "# Plot all Models\n",
    "for cmod in models_2_plot:\n",
    "    print(cmod)\n",
    "\n",
    "    # Load in Model\n",
    "    # Find only files that have current year and month in filename (speeds up loading)\n",
    "    all_files = os.path.join(E.model[cmod][runType]['sipn_nc'], '*.nc') \n",
    "\n",
    "    # Check we have files \n",
    "    files = glob.glob(all_files)\n",
    "    if not files:\n",
    "        continue # Skip this model\n",
    "\n",
    "    # Get list of variablse we want to drop\n",
    "    drop_vars = [x for x in xr.open_dataset(sorted(files)[-1],autoclose=True).data_vars if x not in variables]\n",
    "    \n",
    "    # Load in model   \n",
    "    ds_model_ALL = xr.open_mfdataset(sorted(files), \n",
    "                                 chunks={ 'fore_time': 1,'init_time': 1,'nj': 304, 'ni': 448},  \n",
    "                                 concat_dim='init_time', autoclose=True, \n",
    "                                 parallel=True, drop_variables=drop_vars)\n",
    "                                 # preprocess=lambda x : is_in_time_range(x)) # 'fore_time': 1, ,\n",
    "    ds_model_ALL.rename({'nj':'x', 'ni':'y'}, inplace=True)\n",
    "    \n",
    "    # Sort by init_time\n",
    "    ds_model_ALL = ds_model_ALL.sortby('init_time')\n",
    "        \n",
    "    # Get Valid time\n",
    "    ds_model_ALL = import_data.get_valid_time(ds_model_ALL)\n",
    "    \n",
    "    # For each variable\n",
    "    for cvar in variables:\n",
    "\n",
    "        # For each init time period\n",
    "        for it in init_slice: \n",
    "            it_start = it-np.timedelta64(Ndays,'D') + np.timedelta64(1,'D') # Start period for init period (it is end of period). Add 1 day because when\n",
    "            # we select using slice(start,stop) it is inclusive of end points. So here we are defining the start of the init AND the start of the valid time.\n",
    "            # So we need to add one day, so we don't double count.\n",
    "            print(it_start,\"to\",it)\n",
    "\n",
    "            # For each forecast time we haven't plotted yet\n",
    "            #ft_to_plot = ds_status.sel(init_time=it)\n",
    "            #ft_to_plot = ft_to_plot.where(ft_to_plot.isnull(), drop=True).fore_time\n",
    "\n",
    "            for ft in da_slices.values: \n",
    "\n",
    "                cdoy_end = pd.to_datetime(it + ft).timetuple().tm_yday # Get current day of year end for valid time\n",
    "                cdoy_start = pd.to_datetime(it_start + ft).timetuple().tm_yday  # Get current day of year end for valid time\n",
    "\n",
    "                # Get datetime64 of valid time start and end\n",
    "                valid_start = it_start + ft\n",
    "                valid_end = it + ft\n",
    "\n",
    "                # Loop through variable of interest + any metrics (i.e. SIP) based on that\n",
    "                for metric in metrics_all[cvar]:\n",
    "\n",
    "                    # File paths and stuff\n",
    "                    out_metric_dir = os.path.join(E.model['MME_NEW'][runType]['sipn_nc'], cvar, metric)\n",
    "                    if not os.path.exists(out_metric_dir):\n",
    "                        os.makedirs(out_metric_dir) \n",
    "                        \n",
    "                    out_init_dir = os.path.join(out_metric_dir, pd.to_datetime(it).strftime('%Y-%m-%d'))\n",
    "                    if not os.path.exists(out_init_dir):\n",
    "                        os.makedirs(out_init_dir)\n",
    "                        \n",
    "                    out_mod_dir = os.path.join(out_init_dir, cmod)\n",
    "                    if not os.path.exists(out_mod_dir):\n",
    "                        os.makedirs(out_mod_dir)     \n",
    "                        \n",
    "                    out_nc_file = os.path.join(out_mod_dir, pd.to_datetime(it+ft).strftime('%Y-%m-%d')+'_'+cmod+'.nc')\n",
    "\n",
    "                    # Only update if either we are updating All or it doesn't yet exist\n",
    "                    # OR, its one of the last 3 init times \n",
    "                    if updateAll | (os.path.isfile(out_nc_file)==False) | np.any(it in init_slice[-3:]):\n",
    "                        #print(\"    Updating...\")\n",
    "\n",
    "                        # Select init period and fore_time of interest\n",
    "                        ds_model = ds_model_ALL.sel(init_time=slice(it_start, it))\n",
    "\n",
    "                        # Check we found any init_times in range\n",
    "                        if ds_model.init_time.size==0:\n",
    "                            #print('init_time not found.')\n",
    "                            continue\n",
    "\n",
    "                        # Select var of interest (if available)\n",
    "                        if cvar in ds_model.variables:\n",
    "                            ds_model = ds_model[cvar]\n",
    "                        else:\n",
    "                            #print('cvar not found.')\n",
    "                            continue\n",
    "\n",
    "                        # 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 and init_times\n",
    "                        ds_model = ds_model.mean(dim=['fore_time','init_time'])\n",
    "\n",
    "                        if metric=='mean': # Calc ensemble mean\n",
    "                            ds_model = ds_model.mean(dim='ensemble')\n",
    "\n",
    "                        elif metric=='SIP': # Calc probability\n",
    "                            # Remove ensemble members having missing data\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",
    "\n",
    "                        elif metric=='anomaly': # Calc anomaly in reference to mean observed 1980-2010\n",
    "                            # Get climatological mean\n",
    "                            da_obs_mean = mean_1980_2010_sic.isel(time=slice(cdoy_start,cdoy_end)).mean(dim='time')\n",
    "                            # Calc anomaly\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",
    "\n",
    "                        # drop ensemble if still present\n",
    "                        if 'ensemble' in ds_model:\n",
    "                            ds_model = ds_model.drop('ensemble')\n",
    "\n",
    "                        ds_model.coords['model'] = cmod\n",
    "                        if 'xm' in ds_model:\n",
    "                            ds_model = ds_model.drop(['xm','ym']) #Dump coords we don't use\n",
    "\n",
    "                        # Add Coords info\n",
    "                        ds_model.name = metric\n",
    "                        ds_model.coords['model'] = cmod\n",
    "                        ds_model.coords['init_start'] = it_start\n",
    "                        ds_model.coords['init_end'] = it\n",
    "                        ds_model.coords['valid_start'] = it_start+ft\n",
    "                        ds_model.coords['valid_end'] = it+ft\n",
    "                        ds_model.coords['fore_time'] = ft\n",
    "                        \n",
    "                        # Save to file\n",
    "                        ds_model.to_netcdf(out_nc_file)\n",
    "\n",
    "                        # Clean up for current model\n",
    "                        ds_model = None\n"
   ]
  },
  {
   "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": [
    "cmod = 'climatology'\n",
    "\n",
    "all_files = os.path.join(mod_dir,cmod,runType,'sipn_nc', str(cd.year)+'*.nc')\n",
    "files = glob.glob(all_files)\n",
    "\n",
    "obs_clim_model = xr.open_mfdataset(sorted(files), \n",
    "        chunks={'time': 30, 'x': 304, 'y': 448},  \n",
    "         concat_dim='time', autoclose=True, parallel=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from esio import metrics\n",
    "# ds_region = xr.open_mfdataset(os.path.join(E.grid_dir, 'sio_2016_mask_Update.nc'))\n",
    "\n",
    "# X = obs_clim_model.sic.sel(time=slice('2018-09-01','2018-09-30'))\n",
    "\n",
    "# metrics.calc_extent(da=X, region=ds_region).mean(dim='time').values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###########################################################\n",
    "#          climatology  trend                             #\n",
    "###########################################################\n",
    "\n",
    "cmod = 'climatology'\n",
    "\n",
    "all_files = os.path.join(mod_dir,cmod,runType,'sipn_nc', str(cd.year)+'*.nc')\n",
    "files = glob.glob(all_files)\n",
    "\n",
    "obs_clim_model = xr.open_mfdataset(sorted(files), \n",
    "        chunks={'time': 30, 'x': 304, 'y': 448},  \n",
    "         concat_dim='time', autoclose=True, parallel=True)\n",
    "\n",
    "# For each variable\n",
    "for cvar in variables:\n",
    "\n",
    "    # For each init time period\n",
    "    for it in init_slice: \n",
    "        it_start = it-np.timedelta64(Ndays,'D') + np.timedelta64(1,'D') # Start period for init period (it is end of period). Add 1 day because when\n",
    "        # we select using slice(start,stop) it is inclusive of end points. So here we are defining the start of the init AND the start of the valid time.\n",
    "        # So we need to add one day, so we don't double count.\n",
    "        print(it_start,\"to\",it)\n",
    "\n",
    "        for ft in da_slices.values: \n",
    "\n",
    "            cdoy_end = pd.to_datetime(it + ft).timetuple().tm_yday # Get current day of year end for valid time\n",
    "            cdoy_start = pd.to_datetime(it_start + ft).timetuple().tm_yday  # Get current day of year end for valid time\n",
    "\n",
    "            # Get datetime64 of valid time start and end\n",
    "            valid_start = it_start + ft\n",
    "            valid_end = it + ft\n",
    "\n",
    "            # Loop through variable of interest + any metrics (i.e. SIP) based on that\n",
    "            for metric in metrics_all[cvar]:\n",
    "\n",
    "                # File paths and stuff\n",
    "                out_metric_dir = os.path.join(E.model['MME_NEW'][runType]['sipn_nc'], cvar, metric)\n",
    "                if not os.path.exists(out_metric_dir):\n",
    "                    os.makedirs(out_metric_dir) \n",
    "\n",
    "                out_init_dir = os.path.join(out_metric_dir, pd.to_datetime(it).strftime('%Y-%m-%d'))\n",
    "                if not os.path.exists(out_init_dir):\n",
    "                    os.makedirs(out_init_dir)\n",
    "\n",
    "                out_mod_dir = os.path.join(out_init_dir, cmod)\n",
    "                if not os.path.exists(out_mod_dir):\n",
    "                    os.makedirs(out_mod_dir)     \n",
    "\n",
    "                out_nc_file = os.path.join(out_mod_dir, pd.to_datetime(it+ft).strftime('%Y-%m-%d')+'_'+cmod+'.nc')\n",
    "\n",
    "                # Only update if either we are updating All or it doesn't yet exist\n",
    "                # OR, its one of the last 3 init times \n",
    "                if updateAll | (os.path.isfile(out_nc_file)==False) | np.any(it in init_slice[-3:]):\n",
    "                    #print(\"    Updating...\")\n",
    "\n",
    "                    # Check if we have any valid times in range of target dates\n",
    "                    ds_model = obs_clim_model[cvar].where((obs_clim_model.time>=valid_start) & (obs_clim_model.time<=valid_end), drop=True) \n",
    "                    if 'time' in ds_model.lat.dims:\n",
    "                        ds_model.coords['lat'] = ds_model.lat.isel(time=0).drop('time') # Drop time from lat/lon dims (not sure why?)\n",
    "\n",
    "                    # If we have any time\n",
    "                    if ds_model.time.size > 0:\n",
    "\n",
    "                        # Average over time\n",
    "                        ds_model = ds_model.mean(dim='time')\n",
    "\n",
    "                        if metric=='mean': # Calc ensemble mean\n",
    "                            ds_model = ds_model\n",
    "                        elif metric=='SIP': # Calc probability\n",
    "                            # Issue of some ensemble members having missing data\n",
    "                            ocnmask = ds_model.notnull()\n",
    "                            ds_model = (ds_model>=0.15).where(ocnmask)\n",
    "                        elif metric=='anomaly': # Calc anomaly in reference to mean observed 1980-2010\n",
    "                            # Get climatological mean\n",
    "                            da_obs_mean = mean_1980_2010_sic.isel(time=slice(cdoy_start,cdoy_end)).mean(dim='time')\n",
    "                            # Get anomaly\n",
    "                            ds_model = ds_model - 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",
    "\n",
    "                        # Drop un-needed coords to match model format\n",
    "                        if 'doy' in ds_model.coords:\n",
    "                            ds_model = ds_model.drop(['doy'])\n",
    "                        if 'xm' in ds_model.coords:\n",
    "                            ds_model = ds_model.drop(['xm'])\n",
    "                        if 'ym' in ds_model.coords:\n",
    "                            ds_model = ds_model.drop(['ym'])\n",
    "                    \n",
    "                        # Add Coords info\n",
    "                        ds_model.name = metric\n",
    "                        ds_model.coords['model'] = cmod\n",
    "                        ds_model.coords['init_start'] = it_start\n",
    "                        ds_model.coords['init_end'] = it\n",
    "                        ds_model.coords['valid_start'] = it_start+ft\n",
    "                        ds_model.coords['valid_end'] = it+ft\n",
    "                        ds_model.coords['fore_time'] = ft\n",
    "                        \n",
    "                        # Save to file\n",
    "                        ds_model.to_netcdf(out_nc_file)\n",
    "\n",
    "                        # Clean up for current model\n",
    "                        ds_model = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "############################################################################\n",
    "#                               OBSERVATIONS                               #\n",
    "############################################################################\n",
    "\n",
    "cmod = 'Observed'\n",
    "\n",
    "updateAll = True # We ALWAYS want to update all observations, because each day we get new obs that can be used to evaluate forecasts from up to 12 months ago\n",
    "\n",
    "# For each variable\n",
    "for cvar in variables:\n",
    "\n",
    "    # For each init time period\n",
    "    for it in init_slice: \n",
    "        it_start = it-np.timedelta64(Ndays,'D') + np.timedelta64(1,'D') # Start period for init period (it is end of period). Add 1 day because when\n",
    "        # we select using slice(start,stop) it is inclusive of end points. So here we are defining the start of the init AND the start of the valid time.\n",
    "        # So we need to add one day, so we don't double count.\n",
    "        print(it_start,\"to\",it)\n",
    "\n",
    "        for ft in da_slices.values: \n",
    "\n",
    "            cdoy_end = pd.to_datetime(it + ft).timetuple().tm_yday # Get current day of year end for valid time\n",
    "            cdoy_start = pd.to_datetime(it_start + ft).timetuple().tm_yday  # Get current day of year end for valid time\n",
    "\n",
    "            # Get datetime64 of valid time start and end\n",
    "            valid_start = it_start + ft\n",
    "            valid_end = it + ft\n",
    "\n",
    "            # Loop through variable of interest + any metrics (i.e. SIP) based on that\n",
    "            for metric in metrics_all[cvar]:\n",
    "\n",
    "                # File paths and stuff\n",
    "                out_metric_dir = os.path.join(E.model['MME_NEW'][runType]['sipn_nc'], cvar, metric)\n",
    "                if not os.path.exists(out_metric_dir):\n",
    "                    os.makedirs(out_metric_dir) \n",
    "\n",
    "                out_init_dir = os.path.join(out_metric_dir, pd.to_datetime(it).strftime('%Y-%m-%d'))\n",
    "                if not os.path.exists(out_init_dir):\n",
    "                    os.makedirs(out_init_dir)\n",
    "\n",
    "                out_mod_dir = os.path.join(out_init_dir, cmod)\n",
    "                if not os.path.exists(out_mod_dir):\n",
    "                    os.makedirs(out_mod_dir)     \n",
    "\n",
    "                out_nc_file = os.path.join(out_mod_dir, pd.to_datetime(it+ft).strftime('%Y-%m-%d')+'_'+cmod+'.nc')\n",
    "\n",
    "                # Only update if either we are updating All or it doesn't yet exist\n",
    "                # OR, its one of the last 3 init times \n",
    "                if updateAll | (os.path.isfile(out_nc_file)==False) | np.any(it in init_slice[-3:]):\n",
    "                    #print(\"    Updating...\")\n",
    "\n",
    "                    # Check if we have any valid times in range of target dates\n",
    "                    ds_model = da_obs_c = ds_81[cvar].sel(time=slice(valid_start, valid_end))\n",
    "                    \n",
    "                    if 'time' in ds_model.lat.dims:\n",
    "                        ds_model.coords['lat'] = ds_model.lat.isel(time=0).drop('time') # Drop time from lat/lon dims (not sure why?)\n",
    "\n",
    "                    # Check we have all observations for this week (7)\n",
    "                    if ds_model.time.size == 7:\n",
    "\n",
    "                        if metric=='mean':\n",
    "                            ds_model = ds_model.mean(dim='time') #ds_81.sic.sel(time=(it + ft))\n",
    "                        elif metric=='SIP':\n",
    "                            ds_model = (ds_model >= 0.15).mean(dim='time').astype('int').where(ds_model.isel(time=0).notnull())\n",
    "                        elif metric=='anomaly':\n",
    "                            da_obs_VT = ds_model.mean(dim='time')\n",
    "                            da_obs_mean = mean_1980_2010_sic.isel(time=slice(cdoy_start,cdoy_end)).mean(dim='time')\n",
    "                            ds_model = da_obs_VT - da_obs_mean\n",
    "                        else:\n",
    "                            raise ValueError('Not implemented')\n",
    "\n",
    "                        # Drop coords we don't need\n",
    "                        ds_model = ds_model.drop(['hole_mask','xm','ym'])\n",
    "                        if 'time' in ds_model:\n",
    "                            ds_model = ds_model.drop('time')\n",
    "\n",
    "                        # Add Coords info\n",
    "                        ds_model.name = metric\n",
    "                        ds_model.coords['model'] = cmod\n",
    "                        ds_model.coords['init_start'] = it_start\n",
    "                        ds_model.coords['init_end'] = it\n",
    "                        ds_model.coords['valid_start'] = it_start+ft\n",
    "                        ds_model.coords['valid_end'] = it+ft\n",
    "                        ds_model.coords['fore_time'] = ft\n",
    "\n",
    "                        # Write to disk\n",
    "                        ds_model.to_netcdf(out_nc_file)\n",
    "\n",
    "                        # Clean up for current model\n",
    "                        ds_model = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Finished Calculating Weekly means\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from dask.distributed import Client\n",
    "# client = Client(n_workers=8)\n",
    "# client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cvar = 'sic' # hard coded for now\n",
    "# # Load in all data and write to Zarr\n",
    "# # Load in all metrics for given variable\n",
    "# print(\"Loading in weekly metrics...\")\n",
    "# ds_m = import_data.load_MME_by_init_end(E=E, \n",
    "#                                         runType=runType, \n",
    "#                                         variable=cvar, \n",
    "#                                         metrics=metrics_all[cvar])\n",
    "\n",
    "# # Drop models that we don't evaluate (i.e. monthly means)\n",
    "# models_keep = [x for x in ds_m.model.values if x not in ['noaasipn','modcansipns_3','modcansipns_4']]\n",
    "# ds_m = ds_m.sel(model=models_keep)\n",
    "# # Get list of dynamical models that are not observations\n",
    "# dynamical_Models = [x for x in ds_m.model.values if x not in ['Observed','climatology','dampedAnomaly','dampedAnomalyTrend']]\n",
    "# # # Get list of all models\n",
    "# # all_Models = [x for x in ds_m.model.values if x not in ['Observed']]\n",
    "# # Add MME\n",
    "# MME_avg = ds_m.sel(model=dynamical_Models).mean(dim='model') # only take mean over dynamical models\n",
    "# MME_avg.coords['model'] = 'MME'\n",
    "# ds_ALL = xr.concat([ds_m, MME_avg], dim='model')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# ###### ADD METADATA #################\n",
    "\n",
    "# ## Add coordinate system info\n",
    "# ds_ALL.coords['crs'] = xr.DataArray('crs')\n",
    "# ds_ALL['crs'].attrs = {\n",
    "#     'comment': '(https://nsidc.org/data/polar-stereo/ps_grids.html or https://nsidc.org/data/oib/epsg_3413.html) This is a container variable that describes the grid_mapping used by the data in this file. This variable does not contain any data; only information about the geographic coordinate system',\n",
    "#     'grid_mapping_name': 'polar_stereographic',\n",
    "#     'straight_vertical_longitude_from_pole':'-45',\n",
    "#     'latitude_of_projection_origin': '90.0',\n",
    "#     'standard_parallel':'70',\n",
    "#     'false_easting':'0',\n",
    "#     'false_northing':'0'\n",
    "#     }\n",
    "\n",
    "# # Add time coords\n",
    "# ds_ALL.coords['init_start'] = ds_ALL.init_end - np.timedelta64(Ndays,'D') + np.timedelta64(1,'D')\n",
    "# ds_ALL['init_start'].attrs = {\n",
    "#     'comment':        'Start date for weekly average period',\n",
    "#     'long_name':      'Start date for weekly average period',\n",
    "#     'standard_name':  \"start_init_date\"}\n",
    "\n",
    "# ds_ALL['init_end'].attrs = {\n",
    "#     'comment':        'End date for weekly average period',\n",
    "#     'long_name':      'End date for weekly average period',\n",
    "#     'standard_name':  \"end_init_date\"}\n",
    "\n",
    "# ds_ALL['fore_time'].attrs = {\n",
    "#     'comment':        'Forecast lead time',\n",
    "#     'long_name':      'Forecast lead time',\n",
    "#     'standard_name':  \"forecast_lead_time\"}\n",
    "\n",
    "# # Add Valid time (start and end period)\n",
    "# ds_ALL = import_data.get_valid_time(ds_ALL, init_dim='init_end', fore_dim='fore_time')\n",
    "# ds_ALL.rename({'valid_time':'valid_end'}, inplace=True);\n",
    "# ds_ALL.coords['valid_start'] = ds_ALL.valid_end - np.timedelta64(Ndays,'D') + np.timedelta64(1,'D')\n",
    "\n",
    "# # Add attributes\n",
    "# ds_ALL['valid_end'].attrs = {\n",
    "#     'comment':        'End Valid date for weekly average period',\n",
    "#     'long_name':      'End Valid date for weekly average period',\n",
    "#     'standard_name':  \"end_valid_date\"}\n",
    "\n",
    "# ds_ALL['valid_start'].attrs = {\n",
    "#     'comment':        'Start Valid date for weekly average period',\n",
    "#     'long_name':      'Start Valid date for weekly average period',\n",
    "#     'standard_name':  \"start_valid_date\"}\n",
    "\n",
    "# # Add Variable attributes\n",
    "# ds_ALL['SIP'].attrs = {\n",
    "#     'comment':        'Sea ice probability, calculated by averaging across ensemble members predictions of sea ice concentration >= 0.15',\n",
    "#     'grid_mapping':   'crs',\n",
    "#     'long_name':      'Sea ice probability',\n",
    "#     'standard_name':  \"sea_ice_probability\",\n",
    "#     'units':          'fraction'}\n",
    "\n",
    "# ds_ALL['anomaly'].attrs = {\n",
    "#     'comment':        'Anomaly of the forecasted sea ice concentration mean (ensemble average) compared to the 1980 to 2010 Observed Climatology',\n",
    "#     'grid_mapping':   'crs',\n",
    "#     'long_name':      'Anomaly',\n",
    "#     'standard_name':  \"anomaly\",\n",
    "#     'units':          'fraction'}\n",
    "\n",
    "# ds_ALL['mean'].attrs = {\n",
    "#     'comment':        'Mean of the forecasted sea ice concentration (ensemble average)',\n",
    "#     'grid_mapping':   'crs',\n",
    "#     'long_name':      'Sea ice concentration',\n",
    "#     'standard_name':  \"sea_ice_concentration\",\n",
    "#     'units':          'fraction'}\n",
    "\n",
    "# # Dataset Attributes\n",
    "# ds_ALL.attrs = {\n",
    "# 'comment':                         'Weekly mean sea ice concentration forecasted by multiple models as well as observed by remotly sensed passive microwave sensors.',\n",
    "# 'contact':                         'nicway@uw.edu',\n",
    "# 'creator_email':                   'nicway@uw.edu',\n",
    "# 'creator_name':                    'Nicholas Wayand, University of Washington',\n",
    "# 'creator_url':                     'https://atmos.uw.edu/sipn/',\n",
    "# 'date_created':                    '2018-12-03T00:00:00',\n",
    "# 'date_modified':                   datetime.datetime.now().strftime('%Y-%m-%dT%H:%M:%S'),\n",
    "# 'geospatial_lat_max':              str(float(ds_ALL.lat.max().values)),\n",
    "# 'geospatial_lat_min':              str(float(ds_ALL.lat.min().values)),\n",
    "# 'geospatial_lat_resolution':       '~25km',\n",
    "# 'geospatial_lat_units':            'degrees_north',\n",
    "# 'geospatial_lon_max':              str(float(ds_ALL.lon.max().values)),\n",
    "# 'geospatial_lon_min':              str(float(ds_ALL.lon.min().values)),\n",
    "# 'geospatial_lon_resolution':       '~25km',\n",
    "# 'geospatial_lon_units':            'degrees_east',\n",
    "# 'history':                         datetime.datetime.now().strftime('%Y-%m-%dT%H:%M:%S')+': updated by Nicholas Wayand',\n",
    "# 'institution':                     'UW, SIPN, ARCUS',\n",
    "# 'keywords':                        'Arctic > Sea ice concentration > Prediction',\n",
    "# 'product_version':                 '1.0',\n",
    "# 'project':                         'Sea Ice Prediction Network Phase II',\n",
    "# 'references':                      'Wayand, N.E., Bitz, C.M., and E. Blanchard-Wrigglesworth, (in review). A year-round sub-seasonal to seasonal sea ice prediction portal. Submited to Geophysical Research letters.',\n",
    "# 'source':                          'Numerical model predictions and Passive microwave measurments.',\n",
    "# 'summary':                         'Dataset is updated daily with weekly sea ice forecasts',\n",
    "# 'time_coverage_end':               pd.to_datetime(ds_ALL.valid_end.max().values).strftime('%Y-%m-%dT%H:%M:%S'),\n",
    "# 'time_coverage_start':             pd.to_datetime(ds_ALL.init_start.min().values).strftime('%Y-%m-%dT%H:%M:%S'),\n",
    "# 'title':                           'SIPN2 Sea ice Concentration Forecasts and Observations.'\n",
    "# }\n",
    "\n",
    "# ####################################\n",
    "# print(ds_ALL)\n",
    "\n",
    "\n",
    "# # Rechunk from ~1MB to 100MB chunks\n",
    "# # Chunk along fore_time and init_end\n",
    "# ds_ALL = ds_ALL.chunk({'fore_time': 10, 'init_end': 10, 'model': 1, 'x': 304, 'y': 448})\n",
    "\n",
    "# # Save to Zarr\n",
    "# print(\"Saving to Zarr...\")\n",
    "# ds_ALL.to_zarr(os.path.join(E.data_dir,'model/zarr', cvar+'.zarr'), mode='w')\n",
    "# print(\"Finished updating Weekly SIC metrics and saved to Zar\")\n",
    "# ds_ALL=None # Flush memory\n",
    "# MME_avg=None\n",
    "# ds_m=None"
   ]
  },
  {
   "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
}