test_GPS_OLS_WTHG-Copy1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime as dt\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "import glob\n",
    "from tidegauge_functions import read_GPS_SONEL\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "datum = 'dAB48_XXXXXXXXX_JPL14.neu'\n",
    "dir_in = 'data/GPS'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data/GPS\\dAB48_XXXXXXXXX_JPL14.neu\n"
     ]
    }
   ],
   "source": [
    "pattern = os.path.join(dir_in, datum)\n",
    "print(pattern)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "filenames = sorted(glob.glob(pattern))  # , key=alphanum_key)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "data/GPS\\dAB48_XXXXXXXXX_JPL14.neu\n",
      "            North   East  Vertical  NorthSTD  EastSTD  VerticalSTD\n",
      "Year                                                              \n",
      "2005-10-06  314.5  262.3     -67.3       1.2      0.7          3.7\n",
      "2005-10-07  313.4  266.7     -65.5       1.1      0.7          3.3\n",
      "2005-10-08  313.9  262.7     -63.4       1.0      0.6          3.2\n",
      "2005-10-09  315.5  265.1     -64.0       1.0      0.6          3.2\n",
      "2005-10-10  316.3  265.4     -64.7       1.1      0.7          3.2\n"
     ]
    }
   ],
   "source": [
    "for f, filepath in enumerate(filenames):\n",
    "    df = read_GPS_SONEL(filepath)\n",
    "    print(f'\\n\\n{filepath}')\n",
    "    print(df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.seasonal import seasonal_decompose\n",
    "def decompose(df, var, site, period):  \n",
    "    \n",
    "    decomposition = seasonal_decompose(df[var], freq = period)\n",
    "\n",
    "    trend = decomposition.trend\n",
    "    seasonal = decomposition.seasonal\n",
    "    residual = decomposition.resid\n",
    "    Amplitude = max(seasonal) - min(seasonal)\n",
    "    \n",
    "    decomposition.plot();\n",
    "    df['denoised'] = df[var] - residual\n",
    "    df['trend'] = trend\n",
    "    print(f'\\n\\n{filepath}')\n",
    "    print(df.head())\n",
    "    print(f'{site} Seasonal Amplitude =', Amplitude)\n",
    "    seasonal_amp = pd.DataFrame(list(zip(list({site}), list(Amplitude))),\n",
    "                            colums = ['Site_Name', 'Seasonal_Amplitude'])\n",
    "    seasonal_amp.to_csv('GPS_Seasonal_Amplitude_SONEL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Change Point Detection\n",
    "def change_pt_detection(df, var, points):\n",
    "    model = \"l1\"  \n",
    "    algo = rpt.Dynp(model=model, min_size=3, jump=5).fit(points)\n",
    "    my_bkps = algo.predict(n_bkps=10)\n",
    "\n",
    "    \n",
    "def change_pt_plot()\n",
    "    rpt.show.display(points, my_bkps, figsize=(10, 6))\n",
    "    plt.suptitle(f\"{site}Change Point Detection\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_calc_OLS(df, var, start, end):\n",
    "    df = df[start:end]\n",
    "    x, y = np.arange(len(df[var].dropna())), df[var].dropna()\n",
    "    x = sm.add_constant(x)\n",
    "    model = sm.OLS(y, x)\n",
    "    res = model.fit()\n",
    "    return res\n",
    "\n",
    "def test_plot_OLS(df, res, res2, site, var, trend_array, trend_array2, point1, point2):\n",
    "    df1 = df[(df.index <= point1)]\n",
    "    df2 = df[(df.index >= point2)]\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(12,6));\n",
    "    ax.plot(df1[var].dropna().index, df1[var].dropna().values, \n",
    "            label='trend', marker='.', linestyle='', color = 'darkgrey')\n",
    "    ax.plot(df1[var].dropna().index, [res1.params.x1*i + res1.params.const for i in np.arange(len(df1[var].dropna()))],\n",
    "           marker = '', linestyle = '-', color = 'blue')\n",
    "    # ax[1].plot(df['Vertical'].dropna().index, res.resid.values);\n",
    "    # ax[1].plot(df['Vertical'].dropna().index, np.abs(res.resid.values));\n",
    "    # ax[1].hlines(0, 0, len(res.resid), color='k');\n",
    "    # ax[1].set_title(\"Residuals\");\n",
    "    ax.plot(df2[var].dropna().index, df2[var].dropna().values, \n",
    "            label='trend', marker='.', linestyle='', color = 'darkgrey')\n",
    "    ax.plot(df2[var].dropna().index, [res2.params.x1*i + res2.params.const for i in np.arange(len(df2[var].dropna()))],\n",
    "           marker = '', linestyle = '-', color = 'red')\n",
    "    # ax[1].plot(df['Vertical'].dropna().index, res.resid.values);\n",
    "    # ax[1].plot(df['Vertical'].dropna().index, np.abs(res.resid.values));\n",
    "    # ax[1].hlines(0, 0, len(res.resid), color='k');\n",
    "    plt.axvline(dt.datetime(2013, 1, 5), color= 'red', linestyle='--')\n",
    "    ax.set_title(f\" Pre 2013 Earthquake Trend = {trend_array * 1:.2f} mm/yr, Post 2013 Earthquake Trend= {trend_array2 * 1:.2f} mm/yr\", fontsize=18);\n",
    "    # ax[1].set_title(\"Residuals\");\n",
    "    plt.suptitle(\"Port Alexander, Alaska, USA\", fontsize=20)\n",
    "    plt.xlabel('Date', fontsize=18)\n",
    "    plt.ylabel('Vertical Land Motion (mm)', fontsize=18)\n",
    "    plt.rc('xtick',labelsize=16)\n",
    "    plt.rc('ytick',labelsize=16)\n",
    "    plt.savefig(f'figs/test_GPS_OLS_split_{site}.png')\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for f, filepath in enumerate(filenames):\n",
    "    print(f'\\n\\n{filepath}')\n",
    "    print(f'f: {f}')\n",
    "    df = read_GPS_SONEL(filepath)\n",
    "    df = decompose(df, 'Vertical', filepath[-24:-20], 365)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def convert_trend_toyearly(df, res):\n",
    "    \n",
    "    period = df.index.year.value_counts().max()\n",
    "    yearlytrend = res.params.x1 * period\n",
    "    \n",
    "    return yearlytrend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = read_GPS_SONEL(filepath)\n",
    "df.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SONEL_trend_array = np.full(np.shape(filenames)[0], np.nan)\n",
    "SONEL_trend_array2 = np.full(np.shape(filenames)[0], np.nan)\n",
    "site_name_array = []\n",
    "\n",
    "\n",
    "for f, filepath in enumerate(filenames):\n",
    "    print(f'\\n\\n{filepath}')\n",
    "    print(f'f: {f}')\n",
    "    ## ToDo - extract and keep other important thing about each site from the file...Lat/Lon, name, etc.\n",
    "    \n",
    "    # Read in data\n",
    "    df = read_GPS_SONEL(filepath)\n",
    "    \n",
    "    # Get trend using linear regression\n",
    "    res1 = test_calc_OLS(df, 'Vertical', '2005-10-6', '2013-1-4')\n",
    "    SONEL_trend_array[f] = convert_trend_toyearly(df, res1)\n",
    "    \n",
    "    res2 = test_calc_OLS(df, 'Vertical', '2013-1-5','2020-1-1')\n",
    "    SONEL_trend_array2[f] = convert_trend_toyearly(df, res2)\n",
    "    \n",
    "    # Save site name\n",
    "    site_name_array.append(filepath[-24:-20])\n",
    "    \n",
    "    print(f\"Trend 1 = {SONEL_trend_array[f] * 1:.2f} mm/yr\")\n",
    "    print(res1.summary())\n",
    "    print(f\"Trend 2 = {SONEL_trend_array2[f] * 1:.2f} mm/yr\")\n",
    "    print(res2.summary())\n",
    "   \n",
    "    # Make Plot\n",
    "    test_plot_OLS(df, res1, res2, filepath[-18:-14], 'Vertical', SONEL_trend_array[f], SONEL_trend_array2[f], '2013-1-5','2013-1-5')"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}