rolling_mean_dascore.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rolling Mean Processing of DAS Data\n",
    "\n",
    "This Jupyter Notebook is created for performing rolling mean processing on a  [spool](https://dascore.org/tutorial/concepts.html#:~:text=read%20the%20docs!-,Data%20structures,-DASCore%20has%20two) of distributed acoustic sensing (DAS) data. It uses [DASCore](https://dascore.org/) package and the ```lf_das.py``` script.\n",
    "\n",
    "\n",
    "<svg width=\"100%\" height=\"1\">\n",
    "  <line x1=\"0\" y1=\"0\" x2=\"100%\" y2=\"0\" style=\"stroke:rgb(0,0,0);stroke-width:2\" />\n",
    "</svg>\n",
    "\n",
    "\n",
    "#### Notes: \n",
    "1. Before using this notebook, make sure you have included the ```lf_das.py``` script in the current directory with this notebook and successfully installed DASCore using ```pip``` (recommended) or ```conda```:\n",
    "    ```python\n",
    "    pip install dascore\n",
    "    ```\n",
    "    or\n",
    "    ```python\n",
    "    conda install dascore -c conda-forge\n",
    "    ```   \n",
    "2. Please find all supported I/O [here](https://dascore.org/#:~:text=specialized%20analysis/visualization.-,Supported%20file%20formats,-name).\n",
    " \n",
    "\n",
    "Current DASCore version: 0.0.13 (tested)\n",
    "\n",
    "Date: 09/07/2023\n",
    "\n",
    "\n",
    "Contact: [Ahmad Tourei](https://github.com/ahmadtourei/)\n",
    "\n",
    "ahmadtourei@gmail.com"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import libraries\n",
    "import warnings\n",
    "warnings.simplefilter('ignore')\n",
    "\n",
    "import dascore as dc\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from time import time\n",
    "from dascore.units import s # s as seconds\n",
    "from lf_das import _get_filename\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get a spool of data to work on"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define data path (spool of data) and output folder \n",
    "data_path = '/mnt/h/data'\n",
    "output_data_folder =  '/mnt/h/results'\n",
    "output_figure_folder = '/mnt/h/figures'\n",
    "\n",
    "# get the sorted spool of data form the defined data path (on first run, it will index the patches and subsequently update the index file for future uses)\n",
    "sp = dc.spool(data_path).sort(\"time\").update()\n",
    "\n",
    "# print the contents of first 5 patches\n",
    "content_df = sp.get_contents()\n",
    "content_df.head()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get some metadata and define a sub spool (if needed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get sampling rate, channel spacing, and gauge length from the first patch\n",
    "patch_0 = sp[0]\n",
    "gauge_length = patch_0.attrs['gauge_length']\n",
    "print(\"Gauge length = \", gauge_length)\n",
    "channel_spacing = patch_0.attrs['d_distance']\n",
    "print(\"Channel spacing = \", channel_spacing)\n",
    "sampling_interval = patch_0.attrs['d_time']\n",
    "print(\"Sampling interval = \", sampling_interval)\n",
    "sampling_rate = 1/(sampling_interval / np.timedelta64(1, 's'))\n",
    "print(\"Sampling rate = \", sampling_rate)\n",
    "\n",
    "# select sub spool\n",
    "t_1 = '2023-03-22 03:00:00'\n",
    "t_2 = '2023-03-22 07:00:00'\n",
    "ch_start = 400\n",
    "ch_end = 1400\n",
    "d_1 = patch_0.coords['distance'][ch_start] # in meter\n",
    "d_2 = patch_0.coords['distance'][ch_end] # in meter\n",
    "sub_sp = sp.select(distance=(d_1, d_2), time=(t_1, t_2))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Apply the rolling mean function and save the results in [DASDAE](https://dascore.org/api/dascore/io/dasdae/core/DASDAEV1.html) format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the target sampling interval (cutoff_freq = 1/(2*d_t)) \n",
    "d_t = 10.0\n",
    "\n",
    "# determine window size in sec.\n",
    "window = d_t*s\n",
    "\n",
    "# determine step size in sec.\n",
    "step = d_t*s\n",
    "\n",
    "# define the scale to apply to the raw data\n",
    "scale_iDAS = float((116*sampling_rate/gauge_length)/1e9)\n",
    "\n",
    "# apply the rolling function on each patch in a for loop\n",
    "for i, patch in enumerate (sub_sp):\n",
    "    print (\"working on patch \", i)\n",
    "    # apply rolling mean function\n",
    "    rolling_mean_patch = patch.rolling(time=window, step=step, engine=\"numpy\").mean()\n",
    "    # scale data\n",
    "    new_scaled_patch = rolling_mean_patch.new(data=rolling_mean_patch.data*scale_iDAS) \n",
    "\n",
    "    # save the result to output folder\n",
    "    filename = _get_filename(new_scaled_patch.attrs['time_min'],\n",
    "                new_scaled_patch.attrs['time_max'])\n",
    "    filename = output_figure_folder + '/' + filename \n",
    "    new_scaled_patch.io.write(filename, \"dasdae\") \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the spool results and merge it\n",
    "rolling_spool = dc.spool(output_data_folder).chunk(time=None)\n",
    "# get the patch out of spool\n",
    "rolling_merged_patch = rolling_spool[0]\n",
    "# get the data out of patch\n",
    "rolling_merged_patch_data = rolling_merged_patch.data\n",
    "\n",
    "# define time axis\n",
    "n_samples = rolling_merged_patch_data.shape[0]\n",
    "num_sec = int(n_samples*d_t)\n",
    "time = np.linspace(0, num_sec, n_samples, dtype=np.float64, endpoint=False)\n",
    "# drop time associated with nan values\n",
    "time[np.isnan(rolling_merged_patch_data[:, 0])] = np.nan\n",
    "time_no_nans = time[~np.isnan(time)]\n",
    "\n",
    "# drop nan values of rolling mean results (same behavior as Pandas)\n",
    "rolling_merged_patch_no_nans = rolling_merged_patch.dropna(\"time\")\n",
    "\n",
    "# make sure time and rolling mean results have the same shape\n",
    "assert time_no_nans.shape[0] == rolling_merged_patch_no_nans.data.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,8))\n",
    "\n",
    "# define the channel of interest\n",
    "channel = 1330\n",
    "ch_inx = channel - ch_start\n",
    "\n",
    "plt.plot(time_no_nans, rolling_merged_patch_no_nans.data[:, ch_inx], label='channel: ' + str(channel))\n",
    "\n",
    "plt.legend(loc='best')\n",
    "plt.ylabel('Strain rate (1/sec)')\n",
    "plt.xlabel('Time (sec) \\n (4 hours, starting from 2023/03/22 03:00:00 UTC)')\n",
    "# plt.ylim(-3e-13, 3e-13)\n",
    "plt.title('Rolling mean results')\n",
    "plt.grid(True)\n",
    "\n",
    "file_name_lowfreq = '/rolling_mean_' + str(int(d_t*2)) + 'sec_sample_interval_channel' + str(channel) + '.jpeg'\n",
    "plt.savefig(output_figure_folder + file_name_lowfreq, dpi=600, format='jpeg')\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dascore",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}