cube_construct_stats.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import required toolboxes and create dropdown list for dataset selection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "Created on Fri Jun 15 13:28:31 2018\n",
    "\n",
    "@author: rthomas\n",
    "\"\"\"\n",
    "import datetime\n",
    "from erddapy import ERDDAP\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "\n",
    "dtype = widgets.Dropdown(\n",
    "    options=['wave_spectral', 'wave_zero', 'weather'],\n",
    "    disabled=False,\n",
    ")\n",
    "print(\"Select dataset to use:\")\n",
    "display(dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the following cell to download the selected dataset from the Marine Institute ERDDAP instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "typed = dtype.value\n",
    "print(typed)\n",
    "\n",
    "# Set ERDDAP server details\n",
    "s = 'https://erddap.marine.ie/erddap'\n",
    "p = 'tabledap'\n",
    "r = 'csv'\n",
    "\n",
    "# Set global variables for the parameters without 'buoy_id', 'latitude' and 'longitude' in the list\n",
    "now = datetime.date.today()\n",
    "now_string = now.strftime('%Y-%m-%d')\n",
    "\n",
    "# Set variables based on data type (typed)\n",
    "if typed.lower() == 'wave_spectral':\n",
    "    metadata = ['station_id','time','buoy_id']\n",
    "    dataset_id = 'IWaveBNetwork_spectral'\n",
    "    syear = 2019    \n",
    "    master_params = ['PeakPeriod',\n",
    "                     'PeakDirection',\n",
    "                     'PeakSpread',\n",
    "                     'SignificantWaveHeight',\n",
    "                     'EnergyPeriod',\n",
    "                     'MeanWavePeriod_Tm01',\n",
    "                     'MeanWavePeriod_Tm02',\n",
    "                     'qcflag']\n",
    "\n",
    "elif typed.lower() == 'wave_zero':\n",
    "    metadata = ['station_id','time','buoy_id']\n",
    "    dataset_id = 'IWaveBNetwork_zerocrossing'\n",
    "    syear = 2019    \n",
    "    master_params = ['Hmax',\n",
    "                     'HmaxPeriod',\n",
    "                     'Havg',\n",
    "                     'Tavg',\n",
    "                     'qcflag']\n",
    "\n",
    "elif typed.lower() == 'weather':\n",
    "    metadata = ['station_id','time']\n",
    "    dataset_id = 'IWBNetwork'\n",
    "    syear = 2019\n",
    "    master_params = ['AtmosphericPressure',\n",
    "                     'WindDirection',\n",
    "                     'WindSpeed',\n",
    "                     'Gust',\n",
    "                     'WaveHeight',\n",
    "                     'WavePeriod',\n",
    "                     'MeanWaveDirection',\n",
    "                     'Hmax',\n",
    "                     'AirTemperature',\n",
    "                     'DewPoint',\n",
    "                     'SeaTemperature',\n",
    "                     'RelativeHumidity',\n",
    "                     'QC_Flag']\n",
    "\n",
    "# Get data for dataset\n",
    "e = ERDDAP(server='https://erddap.marine.ie/erddap',\n",
    "           protocol='tabledap',)\n",
    "e.response = 'csv'\n",
    "e.dataset_id = dataset_id\n",
    "e.constraints = {\n",
    "    'time>=': '%s-01-01T00:00:00Z' % (syear),\n",
    "    'time<=': '%sT00:00:00Z' % (now.strftime('%Y-%m-%d')),\n",
    "}\n",
    "e.variables = metadata + master_params        \n",
    "df = e.to_pandas(index_col=\"time (UTC)\", parse_dates=True)\n",
    "\n",
    "col_names = metadata + master_params\n",
    "col_names.remove('time')\n",
    "\n",
    "df.columns = col_names\n",
    "\n",
    "# Make a copy of the unaltered data download\n",
    "data_full = df.copy(deep=True)\n",
    "print(\"Full resolution data downloaded. Available as 'data_full'.\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the following cell to see full resolution data dataframe structure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_full.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate data availability as a percentage of expected data per day for each variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Utilise quality control flags to clean data set\n",
    "# Code to be added...\n",
    "\n",
    "# Take a working copy of the downloaded data\n",
    "data = df.copy(deep=True)\n",
    "\n",
    "# Add columns for date variable\n",
    "data['Date'] = data.index.date\n",
    "\n",
    "if typed!='weather':\n",
    "    df_summ = data.groupby(['station_id','Date','buoy_id']).count().reset_index(level=['station_id','Date','buoy_id'])\n",
    "    df_summ.loc[:,'res'] = 48\n",
    "    df_summ.loc[df_summ['buoy_id']==13,'res'] = 24\n",
    "else:\n",
    "    df_summ = data.groupby(['station_id','Date']).count().reset_index(level=['station_id','Date'])\n",
    "    df_summ.loc[:,'res'] = 24\n",
    "    \n",
    "# Create master availability dataframe to hold count converted to percentage of expected data points per day\n",
    "df_avail = pd.DataFrame()  \n",
    "\n",
    "for stn in df_summ.station_id.unique().tolist():\n",
    "    df_stn = df_summ[df_summ['station_id']==stn].copy(deep=True)\n",
    "    for param in master_params:\n",
    "        df_stn[param] = df_stn[param].astype(float)\n",
    "    df_stn['res'] = df_stn['res'].astype(float)\n",
    "    df_stn.loc[:,master_params] = df_stn.loc[:,master_params].div(df_stn['res'], axis='index').mul(100.0)\n",
    "    df_fulldates = pd.DataFrame(index = pd.date_range(df_stn.Date.min() - datetime.timedelta(days=df_stn.Date.min().day-1), df_stn.Date.max()))\n",
    "    df_fulldates['Date'] = df_fulldates.index.date\n",
    "    df_fulldates['Year'] = df_fulldates.index.year\n",
    "    df_fulldates['Month'] = df_fulldates.index.month\n",
    "    df_fulldates['DOY'] = df_fulldates.index.dayofyear\n",
    "    df_fulldates = df_fulldates.merge(df_stn, how='outer', left_on='Date', right_on='Date').fillna(0)\n",
    "    df_fulldates.station_id = stn\n",
    "\n",
    "    df_avail = pd.concat([df_avail,df_fulldates])\n",
    "\n",
    "# Set indices and tidy up dataframe\n",
    "df_avail = df_avail.set_index(['station_id', 'Date','Year','Month','DOY'])\n",
    "\n",
    "if typed != 'weather':\n",
    "    df_avail = df_avail.drop(['buoy_id','res','qcflag'], axis=1)\n",
    "else:\n",
    "    df_avail = df_avail.drop(['res','QC_Flag'], axis=1)\n",
    "df_avail.columns = pd.MultiIndex.from_product([df_avail.columns, ['avail']])\n",
    "\n",
    "print(\"Daily availability generated. Available as 'data_avail'.\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the following cell to see data availability dataframe structure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_avail.head()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate daily summary statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split out parameter types for different summary statistics\n",
    "\n",
    "params = []\n",
    "param_dir = []\n",
    "for item in master_params:\n",
    "    if 'qc' not in item.lower():\n",
    "        if 'Dir' in item:\n",
    "            param_dir.append(item)\n",
    "        else:\n",
    "            params.append(item)\n",
    "#%% Take a copy of the data\n",
    "data = df.copy()\n",
    "data['Date'] = data.index.date\n",
    "data['Year'] = data.index.year\n",
    "\n",
    "# Get north and east components for directional measurements\n",
    "param_comp = []\n",
    "for dirtn in param_dir:\n",
    "    data['%s_n' % (dirtn)] = np.cos(data[dirtn]*np.pi/180)\n",
    "    param_comp.append('%s_n' % (dirtn))\n",
    "    data['%s_e' % (dirtn)] = np.sin(data[dirtn]*np.pi/180)\n",
    "    param_comp.append('%s_e' % (dirtn))\n",
    "\n",
    "# Resample for summary statistics for non-directional measurements\n",
    "daily = data.groupby(['station_id','Date','Year'])[params].agg(['min','max','mean','std'])\n",
    "\n",
    "if len(param_dir)!=0:\n",
    "    # Resample for mean and std for directional measurement components (north and east)\n",
    "    data2 = data.groupby(['station_id','Date','Year'])[param_comp].agg(['mean','std'])\n",
    "\n",
    "    # Recalculate direction mean and std from averaged components (north and east)\n",
    "    # Add directly into daily dataframe\n",
    "    for dirtn in param_dir:\n",
    "        daily[(dirtn, 'mean')] = (360 + np.arctan2(data2[('%s_e' % (dirtn), 'mean')], data2[('%s_n' % (dirtn), 'mean')]) * 180/np.pi) % 360\n",
    "        daily[(dirtn, 'std')] = (360 + np.arctan2(data2[('%s_e' % (dirtn), 'std')], data2[('%s_n' % (dirtn), 'std')]) * 180/np.pi) % 360\n",
    "        daily[(dirtn, 'max')] = np.nan\n",
    "        daily[(dirtn, 'min')] = np.nan\n",
    "\n",
    "# Sort daily dataframe\n",
    "daily = daily[sorted(daily.columns.tolist())]\n",
    "data_daily = daily\n",
    "print(\"Daily statistics generated. Available as 'data_daily'.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the following cell to see the daily summary statistics dataframe structure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "data_daily.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot data with interactive widgets:\n",
    "#### This cell needs to be run twice to show plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from ipywidgets import interact\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "plot_daily = data_daily.loc[data_daily.index.levels[0].tolist()[0]]\n",
    "\n",
    "def plotting(station, year, xaxis, yaxis, stat):\n",
    "    idx = pd.IndexSlice\n",
    "\n",
    "    plot_daily = data_daily.loc[station]\n",
    "    \n",
    "    if xaxis in ('PeakDirection','WindDirection','MeanWaveDirection'):\n",
    "        statx = 'mean'\n",
    "    else:\n",
    "        statx = stat\n",
    "    if yaxis in ('PeakDirection','WindDirection','MeanWaveDirection'):\n",
    "        staty = 'mean'\n",
    "    else:\n",
    "        staty = stat\n",
    "        \n",
    "    x = plot_daily.loc[idx[:,year],idx[[xaxis],statx]].values\n",
    "    y = plot_daily.loc[idx[:,year],idx[[yaxis],staty]].values\n",
    "\n",
    "    if xaxis in ('PeakDirection','WindDirection','MeanWaveDirection'):\n",
    "        ax = plt.subplot(111, projection='polar')\n",
    "        ax.set_theta_direction(-1)      # Set degrees to match compass rose\n",
    "        ax.set_theta_zero_location(\"N\") # Set zero to north on the plot\n",
    "        ax.set_rlabel_position(180)     # Set labels to area on plot with minimal data\n",
    "        ax = plt.scatter(x*np.pi/180, y)\n",
    "        plt.title(\"Location: %s Year: %s\" % (station,year))\n",
    "        plt.xlabel(\"%s (%s) [unit]\" % (yaxis, staty))\n",
    "\n",
    "    else:\n",
    "        l = plt.plot(x, y,'o')\n",
    "        plt.setp(l, markersize=5)\n",
    "        plt.setp(l, markerfacecolor='C0')\n",
    "        plt.title(\"Location: %s Year: %s\" % (station,year))\n",
    "        plt.xlabel(\"%s (%s) [unit]\" % (xaxis, statx))\n",
    "        plt.ylabel(\"%s (%s) [unit]\" % (yaxis, staty))\n",
    "        \n",
    "    plt.show()\n",
    "\n",
    "interact(plotting,\n",
    "    station = widgets.Dropdown(\n",
    "        options=data_daily.index.levels[0].tolist(),\n",
    "        description='Location:',\n",
    "        disabled=False),\n",
    "\n",
    "    year = widgets.IntSlider(\n",
    "        min = int(plot_daily.index.levels[1].min()),\n",
    "        max = int(plot_daily.index.levels[1].max()),\n",
    "        description = 'Year:'),\n",
    "\n",
    "    xaxis = widgets.Dropdown(\n",
    "        options=data_daily.columns.levels[0].tolist(),\n",
    "        description='x-axis:',\n",
    "        value = data_daily.columns.levels[0].tolist()[0],\n",
    "        disabled=False),\n",
    "\n",
    "    yaxis = widgets.Dropdown(\n",
    "        options=data_daily.columns.levels[0].tolist(),\n",
    "        description='y-axis:',\n",
    "        value = data_daily.columns.levels[0].tolist()[1],\n",
    "        disabled=False),\n",
    "\n",
    "    stat = widgets.RadioButtons(\n",
    "        options=data_daily.columns.levels[1].tolist(),\n",
    "        description='Statistic:',\n",
    "        disabled=False),\n",
    "        )\n",
    "\n",
    "def avail_map(station1,parameter1,interval):\n",
    "    plot_avail = df_avail.loc[station1]\n",
    "    plot_avail.columns = plot_avail.columns.droplevel(level=1)\n",
    "    \n",
    "    if interval == 'Day':\n",
    "        plot_avail_set = plot_avail\n",
    "        ycol = 'DOY'\n",
    "        tag=False\n",
    "    elif interval == 'Month':\n",
    "        plot_avail_set = plot_avail.groupby(['Year','Month'])[df_avail.columns.levels[0].tolist()].agg(['mean'])\n",
    "        ycol='Month'\n",
    "        tag=True\n",
    "        plot_avail_set.columns = plot_avail_set.columns.droplevel(level=1)\n",
    "\n",
    "    plot_avail_set = plot_avail_set.reset_index()\n",
    "    pivot = plot_avail_set.pivot(index='Year', columns=ycol, values=parameter1)\n",
    "    cmap = sns.cm.rocket_r\n",
    "    ax = sns.heatmap(pivot\n",
    "                    ,annot=None\n",
    "                    ,cmap = cmap)\n",
    "    plt.title(\"%s data availability (as percentage per %s): %s\" % (station1, interval, parameter1))\n",
    "    plt.show()\n",
    "    \n",
    "interact(avail_map,    \n",
    "    station1 = widgets.Dropdown(\n",
    "        options=df_avail.index.levels[0].tolist(),\n",
    "        description='Location:',\n",
    "        disabled=False),\n",
    "         \n",
    "    parameter1 = widgets.Dropdown(\n",
    "        options=df_avail.columns.levels[0].tolist(),\n",
    "        description='Parameter:',\n",
    "        disabled=False),\n",
    "    \n",
    "    interval = widgets.RadioButtons(\n",
    "        options=['Month','Day'],\n",
    "        description='Interval:',\n",
    "        disabled=False),\n",
    "        )\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Commands to save data as csv files. \n",
    "\n",
    "Will not work from MyBinder container.\n",
    "\n",
    "Input details of the working directory to save the csv files and uncomment the commands (remove the # from the start of each line) to make the cells run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#working_directory = os.path.normpath('C:/Downloads')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#datafile = os.path.join(working,'%s_data_raw.csv' % (typed, typed))\n",
    "#availfile = os.path.join(working,'%s_data_availability.csv' % (typed, typed))\n",
    "#dailyfile = os.path.join(working,'%s_daily_summary.csv' % (typed, typed))\n",
    "\n",
    "#data_full.to_csv(datafile, sep=',')\n",
    "#data_avail.to_csv(availfile, sep=',')\n",
    "#data_daily.to_csv(dailyfile, sep=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_full.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### For the wave spectral datasert only: Run following cell to calculate Power percentiles for each site and overlay percentile lines on wave height period scatter plots. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "for stn in data_full['station_id'].unique().tolist():\n",
    "    stn_data = data_full[data_full['station_id']==stn].loc[:,['SignificantWaveHeight','EnergyPeriod']]\n",
    "\n",
    "    stn_data['P'] = 0.49 * stn_data.loc[:,'SignificantWaveHeight']**2 * stn_data.loc[:,'EnergyPeriod']\n",
    "\n",
    "    P_mean = stn_data.loc[:,'P'].mean()\n",
    "    P_max = stn_data.loc[:,'P'].max()\n",
    "    P_min = stn_data.loc[:,'P'].min()\n",
    "    P_90 = np.percentile(stn_data.loc[:,'P'].values,90)\n",
    "    P_99 = np.percentile(stn_data.loc[:,'P'].values,99)\n",
    "\n",
    "    print(\"%s: Power - mean = %s kW/m\" % (stn, round(P_mean,2)))\n",
    "    print(\"%s: Power - max = %s kW/m\" % (stn, round(P_max,2)))\n",
    "    print(\"%s: Power - min = %s kW/m\" % (stn, round(P_min,2)))\n",
    "    print(\"%s: Power - 90th percentile = %s kW/m\" % (stn, round(P_90,2)))\n",
    "    print(\"%s: Power - 99th percentile = %s kW/m\" % (stn, round(P_99,2)))\n",
    "\n",
    "    x = stn_data.loc[:,'EnergyPeriod'].values\n",
    "    y = stn_data.loc[:,'SignificantWaveHeight'].values\n",
    "    x99 = range(1,26)\n",
    "    y99 = []\n",
    "    y90 = []\n",
    "    for i in x99:\n",
    "        y99.append(P_99/(0.49*i))\n",
    "        y90.append(P_90/(0.49*i))\n",
    "\n",
    "    l = plt.plot(x, y,'o')\n",
    "    plt.setp(l, markersize=5)\n",
    "    plt.setp(l, markerfacecolor='C0')\n",
    "    m = plt.plot(x99, y99,linestyle='solid')\n",
    "    plt.setp(m, markersize=10)\n",
    "    plt.setp(m, color='red')\n",
    "    m = plt.plot(x99, y90,linestyle='solid')\n",
    "    plt.setp(m, markersize=10)\n",
    "    plt.setp(m, color='green')\n",
    "    #plt.axis([0, 25, 0, 25])\n",
    "    plt.title(\"Location: %s\" % (stn))\n",
    "    plt.xlabel(\"Energy Period [s]\" )\n",
    "    plt.ylabel(\"Significant Wave Hieght [m]\" )\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ctd_env312_notebooks",
   "language": "python",
   "name": "ctd_env312_notebooks"
  },
  "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.12.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}