PROMICE_processing_tools.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import Counter\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_data(file, year):\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    Loading PROMICE data for a given station and all or given year(s)\n",
    "    into a DataFrame. \n",
    "    \n",
    "    \n",
    "    INTPUTS:\n",
    "        file: Path to the desired file containing PROMICE data [string]\n",
    "        year: Year to import. If 'all', all the years are imported [int,string]\n",
    "    \n",
    "    OUTPUTS:\n",
    "        promice_data: Dataframe containing PROMICE data for the desired settings [DataFrame]\n",
    "        \n",
    "    '''\n",
    "        \n",
    "    #extract site name\n",
    "    global site\n",
    "    site=file.split('/')[-1].split('_')[0]+'_'+file.split('/')[-1].split('_')[1]\n",
    "    \n",
    "    if site[:3]=='MIT' or site[:3]=='EGP' or site[:3]=='CEN':\n",
    "        site=site.split('_')[0]\n",
    "    \n",
    "    #load data\n",
    "    promice_data=pd.read_csv(file, delim_whitespace=True)\n",
    "    \n",
    "    #set invalid values (-999) to nan \n",
    "    promice_data[promice_data==-999.0]=np.nan\n",
    "    \n",
    "    #only keep selected year(s) if needed\n",
    "    global yr; yr=year\n",
    "    \n",
    "    if year!='all':\n",
    "        promice_data=promice_data[promice_data.Year==year]\n",
    "    elif isinstance(year, list):\n",
    "        promice_data[promice_data.Year.isin(year)]\n",
    "    \n",
    "    if promice_data.empty:\n",
    "        print('ERROR: Selected year not available')\n",
    "        return\n",
    "    \n",
    "    \n",
    "    return promice_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def BIC_processing(promice_data, fig_save=False, fig_path=None, visualisation=False):\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    Processing, filtering and exclusion of ice ablation and albedo timeseries \n",
    "    spanning the onset of bare ice conditions.\n",
    "    \n",
    "    \n",
    "    INPUTS:\n",
    "        promice_data: Dataset imported using load_data() [DataFrame]\n",
    "        visualisation: If True, displays raw and processed DPT as well as air temperature, \n",
    "                       boom height and albedo time series (default: False) [boolean]\n",
    "        fig_save: If True, saves figure generated with visualisation option\n",
    "                  (default: False) [boolean]\n",
    "        fig_path: Path where to save figure generated with visualisation option\n",
    "                  (default: None) [boolean]\n",
    "    \n",
    "    OUTPUTS:\n",
    "        promice_data_proc: Processed ice ablation and albedo time \n",
    "                           series [DataFrame]\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    \n",
    "    def DPT_processing(df,year):\n",
    "        \n",
    "        '''\n",
    "        \n",
    "        Processing of Depth Pressure Sensor (DPT) time series using PROMICE 2019-08-02\n",
    "        datasets.        \n",
    "        \n",
    "        INPUTS:\n",
    "            promice_data: Dataframe imported using load_data() [DataFrame]\n",
    "            year: Year of promice_data to process [int]\n",
    "            \n",
    "        OUTPUTS:\n",
    "            DPT_proc: Processed ice ablation time series [pandas.series]\n",
    "            albedo: Albedo time series associated with ice ablation [pandas.series]\n",
    "            DPT_flag: Assessement of the ease to determine bare ice appearance\n",
    "                      from ice ablation (0=no data, 1=, 2=low confidence, \n",
    "                      3= high confidence) [int]\n",
    "            albedo_flag: Determines if an albedo time series will be excluded (0)\n",
    "                         or not (1) [int]\n",
    "            BID: Bare Ice Day, day of bare ice appearance based on ice ablation [float]\n",
    "        \n",
    "        '''\n",
    "        \n",
    "        #load ice ablation and associated time\n",
    "        global z, doy\n",
    "        z=df[\"DepthPressureTransducer_Cor(m)\"].copy()\n",
    "        doy=df[\"DayOfYear\"]\n",
    "        \n",
    "        \n",
    "        def assign_nans(days):\n",
    "            if type(days)==int:\n",
    "                z[days==doy]=np.nan\n",
    "            else:\n",
    "                for d in days:\n",
    "                    z[d==doy]=np.nan\n",
    "        \n",
    "        #initialize albedo_flag\n",
    "        albedo_flag=1\n",
    "        \n",
    "        #manually process DPT measurements for a given site and year and\n",
    "        #identify theoretical ice ablation onset (IAO) \n",
    "        global no_proc\n",
    "        no_proc=0\n",
    "        \n",
    "        if site=='NUK_U':\n",
    "            if year in [2007,2014,2015]: DPT_flag=1\n",
    "            elif year in [2008,2009]: DPT_flag=0\n",
    "            elif year==2010:\n",
    "                z[doy>203]-=2.6; z[doy>204]-=8.8; z[doy>=256]=np.nan\n",
    "                DPT_flag=3; IAO=121\n",
    "            elif year==2011:\n",
    "                assign_nans(235)\n",
    "                z[doy>=236]-=2; z[doy>=276]=np.nan\n",
    "                DPT_flag=3; IAO=156; albedo_flag=0\n",
    "            elif year==2012:\n",
    "                z[(doy>=260)&(z>=-3)]=np.nan\n",
    "                DPT_flag=3; IAO=148\n",
    "            elif year==2013:\n",
    "                z[doy>202]-=2; z[doy>203]-=5.8\n",
    "                DPT_flag=3; IAO=152\n",
    "            elif year==2016:\n",
    "                z[doy>=263]=np.nan; z[(doy<=70)&(z<=-0.5)]=np.nan\n",
    "                z[(doy>38)&(doy<43)]=np.nan\n",
    "                assign_nans(np.arange(71,76)); assign_nans([32,92,93])\n",
    "                DPT_flag=3; IAO=160\n",
    "            elif year==2017:\n",
    "                z[(doy>=250)&(z>=-0.3)]=np.nan\n",
    "                DPT_flag=3; IAO=204\n",
    "            elif year==2018:\n",
    "                z[doy<=103]=np.nan; z[doy>=212]-=8\n",
    "                assign_nans([123,124,210,211])\n",
    "                z[doy>=210]-=1.25; DPT_flag=3; IAO=161\n",
    "            elif year==2019:\n",
    "                DPT_flag=3; IAO=128; z[doy<97]=np.nan\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='NUK_L':\n",
    "            if year==2007: DPT_flag=1\n",
    "            elif year==2008:\n",
    "                z[doy>212]-=7.1; z[doy>211]-=5; assign_nans(156)\n",
    "                DPT_flag=3; IAO=121\n",
    "            elif year==2009:\n",
    "                assign_nans([90,92,95,102,105,142,154,159,169,174]); z[doy>236]-=0.5\n",
    "                DPT_flag=3; IAO=136\n",
    "            elif year==2010:\n",
    "                assign_nans(207); z[doy>206]-=11.1; z[doy>264]=np.nan\n",
    "                DPT_flag=3; IAO=107\n",
    "            elif year==2011:\n",
    "                assign_nans([210,238,239,240]); z[doy>209]-=1.15\n",
    "                z[(doy>=237)&(doy<=239)]-=0.4; z[(doy>=241)&(doy<=242)]-=0.2\n",
    "                z[doy>280]=np.nan; DPT_flag=3; IAO=156; albedo_flag=0\n",
    "            elif year==2012:\n",
    "                assign_nans(241); z[doy>241]-=12; DPT_flag=3; IAO=150\n",
    "            elif year==2013:\n",
    "                assign_nans([23,24,204]); z[doy>24]+=2.6; z[doy>204]-=5.25\n",
    "                DPT_flag=3; IAO=155; albedo_flag=0\n",
    "            elif year==2014:\n",
    "                z[doy>206]-=2.5; z[doy>207]-=3.35; DPT_flag=3; IAO=121\n",
    "                assign_nans([54,55])\n",
    "            elif year in [2015,2016]:\n",
    "                DPT_flag=3\n",
    "                if year==2015:\n",
    "                    IAO=159\n",
    "                else:\n",
    "                    IAO=100\n",
    "            elif year==2017:\n",
    "                assign_nans([117,118,219,220,233,234,235]); z[doy>117]-=14.92; z[doy>=142]-=0.5\n",
    "                DPT_flag=3; IAO=142\n",
    "            elif year==2018:\n",
    "                assign_nans([194,196,211]); z[doy>59]-=14.75\n",
    "                z[(doy>=197)&(doy<=210)]+=7.35; z[(doy>=194)&(doy<=211)]+=0.27\n",
    "                z[(doy>=50)&(doy<=70)]=np.nan\n",
    "                DPT_flag=2\n",
    "            elif year==2019:\n",
    "                DPT_flag=3; IAO=97; albedo_flag=0\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='NUK_N':\n",
    "            if year in [2010,2013]: DPT_flag=1\n",
    "            elif year==2011:\n",
    "                z[doy>=280]=np.nan; DPT_flag=3; IAO=164\n",
    "            elif year==2012:   \n",
    "                z[doy>=240]=np.nan; DPT_flag=3; IAO=158\n",
    "                z[(doy>=90)&(doy<=156)&(z<0)]=np.nan; assign_nans(96); \n",
    "                z[doy==168]-=0.19; assign_nans(np.arange(145,152))\n",
    "            elif year==2014:\n",
    "                DPT_flag=3; IAO=165\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='NUK_K':\n",
    "            if year==2014:\n",
    "                z[(doy>=225)&(doy<=238)]+=1.19; assign_nans(239); DPT_flag=2\n",
    "            elif year==2015:\n",
    "                DPT_flag=3; IAO=200\n",
    "            elif year==2016:\n",
    "                assign_nans([23,24,124])\n",
    "                DPT_flag=3; IAO=166\n",
    "            elif year==2017:\n",
    "                DPT_flag=3; IAO=204\n",
    "            elif year in [2018,2019]:\n",
    "                DPT_flag=3; IAO=176\n",
    "                if year==2019:\n",
    "                    assign_nans([94,127,128,129])\n",
    "                elif year==2018:\n",
    "                    assign_nans(177); z[doy>177]+=0.07\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='UPE_L':\n",
    "            if year==2009:\n",
    "                z[(doy>=120)]=np.nan; DPT_flag=1\n",
    "            elif year==2010:\n",
    "                assign_nans(range(271,275)); DPT_flag=3; IAO=129\n",
    "            elif year==2011:\n",
    "                z[(doy>=262)&(doy<=311)]=np.nan; DPT_flag=3; IAO=157\n",
    "            elif year==2012:\n",
    "                z[doy>150]-=0.3; z[doy>151]+=0.2; z[doy>152]+=0.1; z[doy>=225]-=0.1\n",
    "                z[doy>=226]-=0.15; z[(doy>=258)&(doy<=308)]=np.nan; DPT_flag=3; IAO=148\n",
    "                assign_nans([149,150,152,153])\n",
    "            elif year==2013:\n",
    "                z[doy<=70]=np.nan; assign_nans([74,214,215]); z[(doy>=73)&(doy<=75)]=np.nan\n",
    "                z[doy>=214]-=15.3\n",
    "                DPT_flag=3; IAO=160\n",
    "            elif year==2014:\n",
    "                z[(doy>=266)&(doy<=270)]=np.nan; DPT_flag=3; IAO=158\n",
    "            elif year==2015:\n",
    "                assign_nans(213); z[doy>=214]-=0.32; z[(doy>=253)&(doy<=261)]=np.nan\n",
    "                DPT_flag=3; IAO=151\n",
    "            elif year==2016: DPT_flag=0\n",
    "            elif year==2017: DPT_flag=1\n",
    "            elif year==2018:\n",
    "                DPT_flag=3; IAO=167\n",
    "            elif year==2019:\n",
    "                DPT_flag=3; IAO=124; z[doy<99]=np.nan\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='UPE_U':\n",
    "            if year==2009: DPT_flag=1\n",
    "            elif (year>=2010)&(year<=2012):\n",
    "                DPT_flag=3\n",
    "                if year==2010:\n",
    "                    IAO=143; assign_nans(114)\n",
    "                elif year==2011:\n",
    "                    IAO=159\n",
    "                elif year==2012:\n",
    "                    IAO=153\n",
    "            elif year in [2013,2018]: DPT_flag=2\n",
    "            elif year==2014:\n",
    "                DPT_flag=3; IAO=171\n",
    "            elif year==2015:\n",
    "                z[doy>=216]-= 5.48; assign_nans(215); DPT_flag=2\n",
    "            elif year==2016:\n",
    "                DPT_flag=3; IAO=183\n",
    "            elif year==2017:\n",
    "                if year==2017:\n",
    "                    IAO=190\n",
    "                    DPT_flag=3\n",
    "            elif year==2019:\n",
    "                DPT_flag=3; IAO=161\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='KPC_L':\n",
    "            if year in [2008,2010,2012]: DPT_flag=1\n",
    "            elif year==2009:\n",
    "                DPT_flag=0\n",
    "            elif year==2011: DPT_flag=0\n",
    "            elif year==2013:\n",
    "                DPT_flag=3\n",
    "                IAO=163\n",
    "            elif year==2014:\n",
    "                DPT_flag=3\n",
    "                IAO=165\n",
    "            elif year==2015:\n",
    "                DPT_flag=3\n",
    "                IAO=172\n",
    "            elif year==2016:\n",
    "                z[doy>209]-=6.297; assign_nans([210,211]); DPT_flag=3; IAO=166\n",
    "                z[doy>211]-=0.1\n",
    "            elif year==2017:\n",
    "                DPT_flag=3; IAO=160\n",
    "            elif year==2018:\n",
    "                DPT_flag=2;\n",
    "            elif year==2019:\n",
    "                z[doy>192]-=4.478; assign_nans(193); DPT_flag=3; IAO=163\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='KPC_U': \n",
    "            if (year>=2008)&(year<=2019): DPT_flag=1\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='KAN_L':\n",
    "            if year==2008: DPT_flag=1\n",
    "            elif year==2009:\n",
    "                z[(doy>=250)&(z>=-3.2)]=np.nan; z[(doy>=260)&(np.isnan(z))]=-3.5 \n",
    "                assign_nans(276); DPT_flag=3; IAO=136\n",
    "            elif year==2010:\n",
    "                z[(doy>=260)&(z>=-5.1)]=np.nan; z[doy>=280]-=0.1; z[(doy>=260)&(np.isnan(z))]=-5.45 \n",
    "                assign_nans([136,284]); DPT_flag=3; IAO=121\n",
    "            elif year==2011:\n",
    "                z[doy>155]-=4; assign_nans(155); DPT_flag=3; IAO=156\n",
    "            elif year==2012:\n",
    "                z[doy>=234]-=11.56; z[(doy>=282)&(doy<=298)]=np.nan; assign_nans([234,235])\n",
    "                DPT_flag=3; IAO=149\n",
    "            elif year==2013:\n",
    "                DPT_flag=3; IAO=141\n",
    "            elif year==2014:\n",
    "                z[(doy>=290)&(doy<=319)]=np.nan; DPT_flag=3; IAO=146; z[doy>133]+=0.07\n",
    "            elif year==2015:\n",
    "                z[doy>=119]-=7.08; assign_nans([118,278]); z[doy>=188]-=0.29\n",
    "                DPT_flag=3; IAO=162\n",
    "            elif year==2016:\n",
    "                z[(doy>198)]-=1.6; z[(doy>=260)&(z>=-4.5)]=np.nan; z[(doy>=260)&(z<=-5.5)]=np.nan\n",
    "                z[(doy>=273)&(doy<=280)]=np.nan; assign_nans(198); DPT_flag=3; IAO=100\n",
    "            elif year==2017:\n",
    "                z[doy>244]-=3.6; z[doy>58]-=z[doy==59]; assign_nans(244); DPT_flag=3; IAO=123\n",
    "            elif year==2018:\n",
    "                assign_nans(240); z[doy>=241]-=10.65; DPT_flag=3; IAO=137\n",
    "            elif year==2019:\n",
    "                z[doy>=248]-=1.3; DPT_flag=3; IAO=115\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='KAN_M':\n",
    "            if (year>=2008)&(year<=2011): DPT_flag=0\n",
    "            elif year==2012:\n",
    "                DPT_flag=3; IAO=169\n",
    "            elif year==2013:\n",
    "                DPT_flag=3; IAO=173; assign_nans(183)\n",
    "            elif year in [2014,2015]: DPT_flag=1\n",
    "            elif year==2016:\n",
    "                DPT_flag=3; IAO=156\n",
    "            elif year==2017:\n",
    "                z[(doy>100)&(doy<135)&(z>0.15)]=np.nan; DPT_flag=3; IAO=202\n",
    "            elif year==2018:\n",
    "                z[doy>=238]=np.nan; DPT_flag=2\n",
    "            elif year==2019:\n",
    "                z[doy>=247]=np.nan; DPT_flag=3; IAO=172\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='QAS_M':\n",
    "            if year==2016: DPT_flag=2\n",
    "            elif year==2017:\n",
    "                assign_nans(range(141,146)); assign_nans([151,152,187,188])\n",
    "                z[(doy>=196)&(doy<=236)]=np.nan; DPT_flag=1\n",
    "            elif year==2018:\n",
    "                z[doy>242]-=5.0 ; assign_nans(range(241,244)); DPT_flag=3; IAO=192\n",
    "                z[doy>=244]-=0.2\n",
    "            elif year==2019:\n",
    "                z[doy>=242]-=4.92; DPT_flag=3; IAO=178\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='QAS_A':\n",
    "            if (year>=2012)&(year<=2015): DPT_flag=0\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='QAS_L':\n",
    "            if year in [2007,2009,2013]: DPT_flag=1\n",
    "            elif year==2008:\n",
    "                z[doy>=219]-=7; z[doy>=277]+=5.1; assign_nans([77,78,92,103,110,220,277])\n",
    "                z[doy>=220]-=1.9; DPT_flag=3; IAO=121\n",
    "            elif year==2010:\n",
    "                z[(doy>=127)&(doy<=128)]=np.nan; assign_nans(147)\n",
    "                z[(doy>=129)&(doy<=146)&(doy<=148)]-=3.5; z[(doy>=144)&(doy<=150)]=np.nan\n",
    "                z[doy>=129]-=1.3; DPT_flag=2; IAO=129\n",
    "            elif year==2011:\n",
    "                z[doy>=222]=np.nan; DPT_flag=3; IAO=155\n",
    "            elif year==2012:\n",
    "                assign_nans([133,134]); z[doy>132]-=0.55; z[doy>230]-=0.8\n",
    "                DPT_flag=3; IAO=151\n",
    "            elif year==2014:\n",
    "                assign_nans([125,235]); z[doy>124]-=1.48; z[doy>=236]-=7.9; z[doy>159]-=0.08\n",
    "                DPT_flag=3; IAO=149\n",
    "            elif year==2015:\n",
    "                DPT_flag=3; IAO=171\n",
    "            elif year==2016:\n",
    "                assign_nans(224); z[doy>=225]-=13.3; DPT_flag=3; IAO=133\n",
    "            elif year==2017:\n",
    "                DPT_flag=3; IAO=152; z[doy>=152]-=0.05\n",
    "            elif year==2018:\n",
    "                assign_nans([123,238]); z[doy>122]-=1; z[doy>132]-=0.9; z[doy>130]+=0.2\n",
    "                z[doy>126]-=0.2; z[doy>=239]-=8.2; DPT_flag=3; IAO=155\n",
    "                assign_nans(np.arange(123,133)); assign_nans(np.arange(133,137))\n",
    "            elif year==2019:\n",
    "                z[doy>142]=np.nan; z[doy>=177]-=1; z[doy>=140]-=0.4; DPT_flag=1 \n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='QAS_U':\n",
    "            if year in [2008,2009,2013]: DPT_flag=1 \n",
    "            elif year==2010:\n",
    "                DPT_flag=3; IAO=175\n",
    "            elif year in [2011,2015,2018]: DPT_flag=2\n",
    "            elif year==2012:\n",
    "                DPT_flag=3; IAO=196; albedo_flag=0\n",
    "            elif year==2014:\n",
    "                DPT_flag=3; IAO=206\n",
    "            elif year==2016:\n",
    "                z[doy>224]=np.nan; DPT_flag=3; IAO=202\n",
    "            elif year==2017:\n",
    "                z[doy<=143]=np.nan\n",
    "                DPT_flag=3; IAO=220\n",
    "            elif year==2019:\n",
    "                z[doy>=243]-=2.24; z[(doy>=241)&(doy<=242)]=np.nan; z[doy<=100]=np.nan\n",
    "                DPT_flag=3; IAO=200\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='THU_L':\n",
    "            if year in [2010,2011,2012,2013,2018]: DPT_flag=1\n",
    "            elif year==2014:\n",
    "                z[doy>255]=np.nan; DPT_flag=3; IAO=177\n",
    "            elif year==2015:\n",
    "                z[doy>=20]+=0.1; z[doy>=182]-=0.3; z[doy>=260]-=0.2; z[doy>250]=np.nan\n",
    "                DPT_flag=3; IAO=163\n",
    "            elif year==2016:\n",
    "                z[doy>=204]-=0.2; DPT_flag=3; IAO=180\n",
    "            elif year==2017:\n",
    "                DPT_flag=3; IAO=194\n",
    "            elif year==2019:\n",
    "                z[(doy>43)&(doy<59)]=np.nan\n",
    "                z[(doy<100)&(z<-0.2)]=np.nan; DPT_flag=3; IAO=161; albedo_flag=0\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='THU_U':\n",
    "            if year in [2010,2018]: DPT_flag=1\n",
    "            elif year in [2011,2013,2015,2019]: DPT_flag=2\n",
    "            elif year==2012:\n",
    "                DPT_flag=3; IAO=191\n",
    "            elif year==2014:\n",
    "                z[doy>=258]+=0.85; assign_nans(258); DPT_flag=2\n",
    "            elif year==2016:\n",
    "                z[doy>=204]+=2.5; assign_nans(203); z[z<-0.5]=np.nan; DPT_flag=2\n",
    "            elif year==2017:\n",
    "                z[z<-0.06]=np.nan; DPT_flag=1\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='THU_U2':\n",
    "            if year==2018:\n",
    "                assign_nans(216); DPT_flag=2\n",
    "            elif year==2019: DPT_flag=1 #instrument broken\n",
    "        \n",
    "        elif site=='SCO_L':\n",
    "            if year==2008: DPT_flag=2\n",
    "            elif (year>=2009)&(year<=2013):\n",
    "                DPT_flag=3\n",
    "                if year==2009:\n",
    "                    IAO=143\n",
    "                elif year==2010:\n",
    "                    IAO=144; albedo_flag=0\n",
    "                elif year==2011:\n",
    "                    IAO=158\n",
    "                elif year==2012: \n",
    "                    IAO=150\n",
    "                elif year==2013: \n",
    "                    IAO=151\n",
    "            elif year==2014:\n",
    "                z[doy>=222]-=14.79; assign_nans(221); DPT_flag=2\n",
    "            elif year in [2015,2016]: \n",
    "                if year==2015:\n",
    "                    DPT_flag=1\n",
    "                else:\n",
    "                    assign_nans(np.arange(115,120)); DPT_flag=3; IAO=158\n",
    "            elif year==2017:\n",
    "                assign_nans(np.arange(110,116))\n",
    "                z[doy>=219]-=9.88; assign_nans(218); DPT_flag=3; IAO=146\n",
    "            elif year==2018:\n",
    "                DPT_flag=3; IAO=148\n",
    "            elif year==2019: DPT_flag=1\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='SCO_U':\n",
    "            if year in [2008,2018]: DPT_flag=2\n",
    "            elif year==2009:\n",
    "                z[doy>=316]=np.nan; DPT_flag=3; IAO=184; albedo_flag=0\n",
    "            elif year==2010: DPT_flag=1\n",
    "            elif year==2011:\n",
    "                DPT_flag=3; IAO=161\n",
    "            elif year==2012:\n",
    "                z[doy>=242]-=1.47; assign_nans(241); DPT_flag=3; IAO=150\n",
    "            elif year in [2013,2014]: \n",
    "                DPT_flag=3\n",
    "                if year==2013:\n",
    "                    z[(doy>22)&(doy<26)]=np.nan; IAO=155\n",
    "                else:\n",
    "                    IAO=162\n",
    "            elif year==2015:\n",
    "                assign_nans([101,102,103,113])\n",
    "                DPT_flag=3; IAO=170\n",
    "            elif year==2016:\n",
    "                assign_nans([47,60,61])\n",
    "                DPT_flag=3; IAO=160\n",
    "            elif year==2017:\n",
    "                z[doy>=217]-=12.08; assign_nans(216); DPT_flag=3; IAO=153\n",
    "                z[(doy>49)&(doy<54)]=np.nan; assign_nans(np.arange(82,93))\n",
    "            elif year==2019:\n",
    "                DPT_flag=3; IAO=160\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='TAS_A':\n",
    "            if year in [2013,2015,2018]: DPT_flag=1\n",
    "            elif year==2014:\n",
    "                z[(doy>41)&(doy<75)]=np.nan;\n",
    "                DPT_flag=3; IAO=188\n",
    "            elif year==2016: \n",
    "                DPT_flag=3; IAO=176\n",
    "            elif year==2017:\n",
    "                DPT_flag=3; IAO=210\n",
    "            elif year==2019:\n",
    "                DPT_flag=3; IAO=195; albedo_flag=0\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='TAS_L':\n",
    "            if year in range(2007,2010): DPT_flag=0\n",
    "            elif year in [2010,2011,2014,2015,2016]: DPT_flag=1\n",
    "            elif year==2012: \n",
    "                z[doy>=252]-=2.86; assign_nans(251); z[doy<=83]=np.nan\n",
    "                DPT_flag=3; IAO=153\n",
    "            elif year==2013:\n",
    "                z[doy>=235]-=np.nan; DPT_flag=2; IAO=140\n",
    "            elif year==2017: \n",
    "                z[doy>=209]-=10.36; assign_nans([207,208]); DPT_flag=2\n",
    "            elif year==2018:\n",
    "                z[doy>=272]+=1.29; assign_nans(271); DPT_flag=3; IAO=156\n",
    "            elif year==2019:\n",
    "                z[doy>=121]-=1.135; DPT_flag=3; IAO=143; albedo_flag=0\n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='TAS_U':\n",
    "            if year==2008:\n",
    "                DPT_flag=3; IAO=167; assign_nans(143)\n",
    "            elif year in [2009,2012]: DPT_flag=1\n",
    "            elif year==2010:\n",
    "                z[doy>=219]-=8.07; assign_nans(218); DPT_flag=2\n",
    "            elif year==2011:\n",
    "                z[doy>=203]-=1.558; assign_nans([202,222]); DPT_flag=3; IAO=161; albedo_flag=0\n",
    "            elif year==2013:\n",
    "                DPT_flag=3; IAO=167\n",
    "            elif year==2014: \n",
    "                z[doy>=216]-=1.29; assign_nans([138,215]); DPT_flag=3; IAO=136\n",
    "                z[doy>=138]+=1.27; assign_nans(137);\n",
    "            elif year==2015: DPT_flag=2 \n",
    "            else: no_proc=1\n",
    "        \n",
    "        elif site=='MIT':\n",
    "            if year==2009:\n",
    "                z[doy>=224]-=6.6; assign_nans([124,223]); DPT_flag=3; IAO=222; albedo_flag=0\n",
    "            elif year==2010:\n",
    "                DPT_flag=3; IAO=190\n",
    "            elif year==2011:\n",
    "                assign_nans(224); DPT_flag=3; IAO=203; z[doy>223]-=0.08\n",
    "            elif year==2012:\n",
    "                z[doy>=250]-=13.70; assign_nans(249); DPT_flag=2\n",
    "            elif year in [2013,2015,2016]: DPT_flag=1\n",
    "            elif year==2014:\n",
    "                DPT_flag=3; IAO=188\n",
    "            elif year==2017:\n",
    "                assign_nans(207); DPT_flag=3; IAO=208; albedo_flag=0\n",
    "            elif year==2018: DPT_flag=2\n",
    "            elif year==2019:\n",
    "                z[doy<12]+=4.4\n",
    "                DPT_flag=3; IAO=186; albedo_flag=0\n",
    "            else: no_proc=1\n",
    "        \n",
    "        \n",
    "        #stations in the accumulation area are not used for this application\n",
    "        elif site=='EGP' or site=='CEN' or site=='KAN_U':\n",
    "            print('WARNING: No processing available for %s' %site)\n",
    "            DPT_flag=0\n",
    "        \n",
    "        if no_proc==1:\n",
    "            print('WARNING: No processing available for %s %s' %(site,year))\n",
    "            DPT_flag=0\n",
    "        \n",
    "        if DPT_flag!=3:\n",
    "            IAO=np.nan\n",
    "                \n",
    "            \n",
    "        albedo=df['Albedo_theta<70d']\n",
    "        \n",
    "        #correction for measurement platform obstruction of the radiometer field \n",
    "        #of view after Aoki et al (2011) that increases average PROMICE BBA by 0.034\n",
    "        albedo+=0.034\n",
    "        \n",
    "        if albedo_flag==0:\n",
    "            albedo[:]=np.nan\n",
    "        if np.sum(np.isnan(albedo))==len(albedo):\n",
    "            albedo_flag=0\n",
    "        \n",
    "        #adjust DPT to have a null pre-melt season ice ablation\n",
    "        z-=np.nanmean(z[(doy<IAO)&(doy>IAO-45)]) \n",
    "        \n",
    "        #determine the start of significant ice ablation (>6cm after IAO)\n",
    "        if DPT_flag==3 and albedo_flag==1:\n",
    "            try:\n",
    "                indx=np.where(z[doy>IAO]<-0.06)[0][0]\n",
    "                BID=doy[doy>IAO].iloc[indx]\n",
    "            except IndexError:\n",
    "                BID=np.nan\n",
    "                albedo_flag=0\n",
    "        else:\n",
    "            BID=np.nan\n",
    "                \n",
    "        \n",
    "        return z, albedo, DPT_flag, BID, albedo_flag\n",
    "    \n",
    "    \n",
    "    def plot_variables():\n",
    "        \n",
    "        fs=13\n",
    "        mpl.rc('xtick', labelsize=fs); mpl.rc('ytick', labelsize=fs); mpl.rc('lines', markersize=3)\n",
    "        \n",
    "        plt.figure(figsize=(10,15))\n",
    "        ax1=plt.subplot(411)\n",
    "        ln1=ax1.plot(df_y_init[\"DayOfYear\"],df_y_init[\"DepthPressureTransducer_Cor(m)\"],'ro-', \n",
    "                     label='Raw', zorder=3)\n",
    "        ax1.set_ylabel('Ice ablation (meters)',fontsize=fs,color='b')\n",
    "        ax2=ax1.twinx()\n",
    "        ln2=ax2.plot(df_y[\"DayOfYear\"],DPT_proc,'go-', label='Processed', zorder=2)\n",
    "        ax2.set_ylabel('Ice ablation (meters)',fontsize=fs,color='g')\n",
    "        ax2.axhline(0,color='gray',LineStyle='--', zorder=1)\n",
    "        ax2.legend(loc='upper left')\n",
    "        lns = ln1+ln2\n",
    "        labs = [l.get_label() for l in lns]\n",
    "        plt.legend(lns, labs)\n",
    "        ax1.get_xaxis().set_visible(False)\n",
    "        \n",
    "        \n",
    "        ax2=plt.subplot(412, sharex=ax1)\n",
    "        ax2.plot(df_y_init[\"DayOfYear\"],df_y_init['Albedo_theta<70d'],'o-',color='purple')\n",
    "        ax2.set_ylabel('Albedo (unitless)',fontsize=fs,color='purple')\n",
    "        ax2.get_xaxis().set_visible(False)\n",
    "        \n",
    "        ax3=plt.subplot(413, sharex=ax1)\n",
    "        ax3.plot(df_y_init[\"DayOfYear\"], HeightSensorBoom_m,'o-',color='orange')\n",
    "        ax3.set_ylabel('Boom height (meters)',fontsize=fs,color='orange')\n",
    "        ax3.get_xaxis().set_visible(False)\n",
    "        \n",
    "        ax4=plt.subplot(414, sharex=ax1)\n",
    "        ax4.plot(df_y_init[\"DayOfYear\"],df_y_init['AirTemperature(C)'],'ro-')\n",
    "        ax4.set_ylabel('Air temperature (°C)',fontsize=fs,color='red')\n",
    "        ax4.set_xlabel('Day of year (DOY)',fontsize=fs)\n",
    "        ax4.axhline(0,color='gray',LineStyle='--')\n",
    "        plt.suptitle('%s %s' %(site,y),fontsize=fs+6)\n",
    "        \n",
    "\n",
    "        if fig_save:\n",
    "            plt.savefig(fig_path+site+'_'+str(y)+'_DPT_proc.png',dpi=300, bbox_inches='tight')\n",
    "            \n",
    "            \n",
    "    \n",
    "    if yr=='all' or isinstance(yr, list):\n",
    "        \n",
    "        #search for available years\n",
    "        if yr=='all':\n",
    "            years=list(Counter(promice_data.Year))\n",
    "        else:\n",
    "            years=yr\n",
    "        \n",
    "        #run BIC_processing() for each year\n",
    "        for i,y in enumerate(years):\n",
    "            \n",
    "            df_y=promice_data[promice_data.Year==int(y)]\n",
    "            df_y_init=df_y.copy()\n",
    "            \n",
    "            DPT_proc, albedo, DPT_flag, BID, albedo_flag=DPT_processing(df=df_y,year=int(y))\n",
    "            DPT_flags=np.zeros(len(DPT_proc))\n",
    "            DPT_flags[:]=DPT_flag\n",
    "            BIDs=np.zeros(len(DPT_proc))\n",
    "            BIDs[:]=BID\n",
    "            albedo_flags=np.zeros(len(DPT_proc))\n",
    "            albedo_flags[:]=albedo_flag\n",
    "            \n",
    "            HeightSensorBoom_m=df_y['HeightSensorBoom(m)']\n",
    "            \n",
    "            #reallocate processed variables to outputs\n",
    "            if i==0:\n",
    "                promice_data_proc=pd.DataFrame({'DPT_proc':DPT_proc, 'DPT_flag':DPT_flags,\n",
    "                                                'BID':BIDs,\n",
    "                                                'DOY': df_y[\"DayOfYear\"], 'Year':df_y[\"Year\"],\n",
    "                                                'Air_temperature_C':df_y['AirTemperature(C)'],\n",
    "                                                'Albedo_theta_inf_70d':albedo,\n",
    "                                                'Albedo_flag':albedo_flags,\n",
    "                                                'HeightSensorBoom_m': HeightSensorBoom_m})\n",
    "            \n",
    "            else:\n",
    "                promice_data_proc_sc=pd.DataFrame({'DPT_proc':DPT_proc, 'DPT_flag':DPT_flags,\n",
    "                                                   'BID':BIDs,\n",
    "                                                   'DOY': df_y[\"DayOfYear\"], 'Year':df_y[\"Year\"],\n",
    "                                                   'Air_temperature_C':df_y['AirTemperature(C)'],\n",
    "                                                   'Albedo_theta_inf_70d':albedo,\n",
    "                                                   'Albedo_flag':albedo_flags,\n",
    "                                                   'Height_sensor_boom_m': HeightSensorBoom_m})\n",
    "                \n",
    "                promice_data_proc=promice_data_proc.append(promice_data_proc_sc)\n",
    "                \n",
    "            \n",
    "            if visualisation: plot_variables()\n",
    "            \n",
    "            \n",
    "    \n",
    "    else:\n",
    "        \n",
    "        df_y=promice_data.copy()\n",
    "        \n",
    "        #run BIC_processing() for the selected year\n",
    "        DPT_proc, albedo, DPT_flag, BID, albedo_flag=DPT_processing(df=df_y,year=int(y))\n",
    "        DPT_flags=np.zeros(len(DPT_proc))\n",
    "        DPT_flags[:]=DPT_flag\n",
    "        BIDs=np.zeros(len(DPT_proc))\n",
    "        BIDs[:]=BID\n",
    "        albedo_flags=np.zeros(len(DPT_proc))\n",
    "        albedo_flags[:]=albedo_flag\n",
    "    \n",
    "        \n",
    "        HeightSensorBoom_m=df_y['HeightSensorBoom(m)']\n",
    "        \n",
    "        if y==2015 and site=='THU_L':\n",
    "            HeightSensorBoom_m[doy==184]=np.nan\n",
    "        elif y==2012 and site=='THU_U':\n",
    "            HeightSensorBoom_m[doy==172]=np.nan; HeightSensorBoom_m[doy==211]=np.nan\n",
    "        \n",
    "        \n",
    "        #reallocate processed variables to outputs\n",
    "        promice_data_proc=pd.DataFrame({'DPT_proc':DPT_proc, 'DPT_flag':DPT_flags,\n",
    "                                        'BID':BIDs,\n",
    "                                        'DOY': df_y[\"DayOfYear\"], 'Year':df_y[\"Year\"],\n",
    "                                        'Air_temperature_C':df_y['AirTemperature(C)'],\n",
    "                                        'Albedo_theta_inf_70d':albedo,\n",
    "                                        'Albedo_flag':albedo_flags,\n",
    "                                        'HeightSensorBoom_m': HeightSensorBoom_m})\n",
    "        \n",
    "        if visualisation: plot_variables()\n",
    "\n",
    "    \n",
    "    return promice_data_proc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def BIC_composite(inpath, dt=45,verbose=True, visualisation=True, save_pkl=False):\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    Composites of air temperature, snow height, ice ablation and albedo time series \n",
    "    centered on bare ice appearance and spanning ± dt days. \n",
    "    \n",
    "    \n",
    "    INPUTS:\n",
    "        inpath: Path to the folder containing PROMICE .txt files [string]\n",
    "        dt: Number of days to consider around bare ice appearance \n",
    "            (default: 45) [int]\n",
    "        verbose: If True, prints and counts station years (default: True) [boolean]\n",
    "        visualisation: If True, displays composite (default: True) [boolean]\n",
    "        save_pkl: If True, saves composite dictionnary into a .pkl file \n",
    "                  (default: False) [boolean]\n",
    "            \n",
    "    OUTPUTS:\n",
    "        composite: Dictionnary containing air_temperature_degrees, ice_ablation_meters,\n",
    "                   snow_height_meters, albedo_unitless keys and associated dataframes\n",
    "                   of station years (per row) spanning BID±dt (along columns). [dictionnary]\n",
    "                     \n",
    "    '''\n",
    "    \n",
    "    files=glob.glob(inpath+'*day_v03_upd.txt')\n",
    "    sites=[]\n",
    "    \n",
    "    for f in files:\n",
    "        site=f.split(os.sep)[1].split('_')[0]+'_'+f.split(os.sep)[1].split('_')[1]\n",
    "        if site[:3]=='MIT' or site[:3]=='CEN' or site[:3]=='EGP':\n",
    "            site=site.split('_')[0]\n",
    "        sites.append(site)\n",
    "    \n",
    "    sites=list(Counter(sites))\n",
    "    \n",
    "    variables={'air_temperature_degrees':'Air_temperature_C',\n",
    "               'ice_ablation_meters':'DPT_proc',\n",
    "               'snow_height_meters':'Height_sensor_boom_m',\n",
    "               'albedo_unitless':'Albedo_theta_inf_70d'} \n",
    "    \n",
    "    #compute timespan associated with ±dt\n",
    "    time_span=2*dt+1\n",
    "    \n",
    "    \n",
    "    composite={var:np.zeros(time_span) for var in variables}\n",
    "    \n",
    "    station_years=[]\n",
    "    stations=[]\n",
    "\n",
    "    for s in sites:\n",
    "        \n",
    "        #stations in the accumulation area are not used for this application\n",
    "        if s=='CEN' or s=='EGP' or s=='KAN_U':\n",
    "            continue\n",
    "        else:\n",
    "            df=load_data(file=inpath+s+'_day_v03_upd.txt', year='all')\n",
    "            df_corr=BIC_processing(df)\n",
    "            \n",
    "            years=list(Counter(df_corr.Year))\n",
    "            \n",
    "            for year in years:\n",
    "                df_yr=df_corr[df_corr.Year==int(year)]\n",
    "                \n",
    "                if df_yr.DPT_flag.iloc[0]==3 and df_yr.Albedo_flag.iloc[0]==1:\n",
    "                    station_years.append(s)\n",
    "                    stations.append(s+' '+str(year))\n",
    "                    BID=int(df_yr.BID.iloc[0])\n",
    "                    df_yr_comp=df_yr[(df_yr.DOY>=BID-dt)&(df_yr.DOY<=BID+dt)]\n",
    "                    \n",
    "                    if verbose: print(s,year)\n",
    "                    \n",
    "                    #append nan values if time_span higher than data set coverage \n",
    "                    if len(df_yr_comp)!=time_span:\n",
    "                        to_append=np.zeros((time_span-len(df_yr_comp),np.shape(df_yr_comp)[1]))\n",
    "                        to_append[:]=np.nan\n",
    "                        df_yr_comp=df_yr_comp.append(pd.DataFrame(to_append))\n",
    "                        \n",
    "                    mean_boom_height=np.nanmean(df_yr_comp.Height_sensor_boom_m[dt+10:dt+30])\n",
    "                    df_yr_comp.Height_sensor_boom_m=mean_boom_height-df_yr_comp.Height_sensor_boom_m\n",
    "                    \n",
    "                    composite={key:np.vstack((composite[key], df_yr_comp[value])) \n",
    "                               for key, value in variables.items()}\n",
    "                    \n",
    "    composite={var:composite[var][1:,:] for var in variables}\n",
    "    \n",
    "    \n",
    "    if verbose: \n",
    "        nb_stations=len(list(Counter(station_years)))\n",
    "        print('%s station years from %s stations' %(len(station_years), nb_stations))\n",
    "    \n",
    "    if visualisation:\n",
    "        \n",
    "        var_names=[]\n",
    "        [var_names.append(key) for key in composite.keys()]\n",
    "        plt.style.use('default') \n",
    "        plt.figure(figsize=(15,15))\n",
    "        op=0\n",
    "        \n",
    "        for i,var in enumerate(composite.values()):\n",
    "            \n",
    "            if var_names[i] in ['ice_ablation_meters','albedo_unitless',\n",
    "                                'air_temperature_degrees','snow_height_meters']:  \n",
    "                op+=1\n",
    "                fs=20\n",
    "                ax=plt.subplot(2,2,op)\n",
    "                plt.plot(np.nanmean(var,axis=0),'k-',zorder=1)\n",
    "                plt.plot(np.nanmean(var,axis=0),'ko',alpha=0.3,markersize=5)\n",
    "                plt.fill_between(np.arange(0,len(var.T)),np.nanmean(var,axis=0)-np.nanstd(var,axis=0),\n",
    "                                 np.nanmean(var,axis=0)+np.nanstd(var,axis=0),color='gray',alpha=0.2)\n",
    "                plt.xticks(np.arange(0,time_span,10),np.arange(-dt,time_span-dt,10),fontsize=fs-2)\n",
    "                plt.yticks(fontsize=fs-2)\n",
    "                ylims=ax.get_ylim()\n",
    "                xlims=ax.get_xlim()\n",
    "                plt.xlabel('bare ice day',fontsize=fs)\n",
    "                \n",
    "                if var_names[i]=='albedo_unitless':\n",
    "                    \n",
    "                    plt.hlines(np.nanmean(var,axis=0)[dt],LineStyle='--',color='red',xmin=xlims[0],xmax=dt,\n",
    "                           label='%.3f±%.3f' %(np.nanmean(var,axis=0)[dt],np.nanstd(var,axis=0)[dt]),zorder=3)\n",
    "                    plt.ylabel('%s, %s' %(var_names[i].split('_')[0],var_names[i].split('_')[1]),fontsize=fs)\n",
    "                    plt.legend(fontsize=fs-5, loc='upper right')\n",
    "                    plt.vlines(dt,ymin=ylims[0],ymax=np.nanmean(var,axis=0)[dt],LineStyle='--',color='black',zorder=3)\n",
    "                \n",
    "                else:\n",
    "                    \n",
    "                    plt.hlines(np.nanmean(var,axis=0)[dt],LineStyle='--',color='black',xmin=xlims[0],xmax=dt, zorder=3)\n",
    "                    plt.axhline(0,color='gray',zorder=2, LineStyle='--')\n",
    "                    plt.ylabel('%s, %s' %(var_names[i].split('_')[0]+' '+var_names[i].split('_')[1],var_names[i].split('_')[2]),fontsize=fs)\n",
    "                    plt.vlines(dt,ymin=ylims[0],ymax=np.nanmean(var,axis=0)[dt],LineStyle='--',color='black',zorder=3)\n",
    "                \n",
    "                plt.ylim([ylims[0],ylims[1]]) \n",
    "                plt.xlim([0,2*dt])\n",
    "                plt.subplots_adjust(hspace=0.25,wspace=0.25)\n",
    "\n",
    "\n",
    "    #convert station years to index-friendly format\n",
    "    stations = [st.replace(' ','_') for st in stations]\n",
    "    \n",
    "    #load Dataframes into output dictionnary \n",
    "    composite={var:pd.DataFrame(composite[var],index=stations,columns=np.arange(-dt,time_span-dt))\n",
    "               for var in composite}\n",
    "    \n",
    "    \n",
    "    if save_pkl:\n",
    "        \n",
    "        file_name='PROMICE_composite_%dBID'%dt\n",
    "        \n",
    "        f = open(inpath+file_name+'.pkl','wb')\n",
    "        pickle.dump(composite,f)\n",
    "        f.close()\n",
    "    \n",
    "    \n",
    "    return composite"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}