low_pass_dascore.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Low-frequency Processing of DAS Data\n",
    "\n",
    "This Jupyter Notebook is created for low-frequency processing of 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 and it is inspired by Dr. Ge Jin's [low frequency processing](https://github.com/DASDAE/DASLowFreqProcessing). \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``` 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.quarto.pub/dascore/).\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",
    "import scipy\n",
    "\n",
    "from time import time\n",
    "from lf_das import LFProc, get_edge_effect_time, get_patch_time, waterfall_plot\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/a'\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()"
   ]
  },
  {
   "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['distance_step']\n",
    "print(\"Channel spacing = \", channel_spacing)\n",
    "sampling_interval = patch_0.attrs['time_step']\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 a 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",
    "# or:\n",
    "# d_1 = -115 # in meter\n",
    "# d_2 = 2000 # in meter\n",
    "sub_sp = sp.select(distance=(d_1, d_2), time=(t_1, t_2)) \n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get low-pass filter parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define memory size that you'd like to dedicate to low-frequency processing \n",
    "memory_size = 10000 # in MB\n",
    "patch_length = get_patch_time(memory_size=memory_size, sampling_rate=sampling_rate, num_ch=ch_end-ch_start)\n",
    "print('patch_length = ', patch_length, str(' sec.'))\n",
    "\n",
    "# define the target sampling interval in seconds\n",
    "d_t = 10.0 # so, cutoff_freq = Nyq_new = 1/(2*d_t) = 0.05 hz\n",
    "\n",
    "# define the desired tolerance for getting the edge time (smaller tolerance results a longer eliminated edges in each patch and higher accuracy. 1e-3 is recommended.)\n",
    "tolerance = 1e-3\n",
    "edge_buffer = get_edge_effect_time(sampling_interval=1/sampling_rate, total_T=patch_length, tol=tolerance, freq=1/d_t)\n",
    "print('edge_buffer = ', edge_buffer, str(' sec.'))\n",
    "\n",
    "# pass the spool to the LFProc class\n",
    "lfp = LFProc(sub_sp)\n",
    "lfp.update_processing_parameter(output_sample_interval=d_t, process_patch_size=int(patch_length/d_t), edge_buff_size=int(np.ceil(edge_buffer/d_t)))\n",
    "\n",
    "# set the output folder - Caution: If you set delete_existing=True, you will remove all contents in the output_data_folder directory. If you set delete_existing=False, you need to have a empty output_data_folder directory to proceed.\n",
    "lfp.set_output_folder(output_data_folder, delete_existing=False)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Do low-frequency processing a portion of the spool"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tic = time()\n",
    "# do lowpass processing on (t_1,t_2) time range\n",
    "t_1 = np.datetime64('2023-03-22T03:00:00')\n",
    "t_2 = np.datetime64('2023-03-22T07:00:00')\n",
    "lfp.process_time_range(t_1,t_2)\n",
    "toc = time()\n",
    "print(f'processing time (sec): {toc-tic}')\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Merge the results before visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp_result = dc.spool(output_data_folder)\n",
    "sp_result = sp_result.chunk(time=None) \n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize the low-pass filtered results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### a) Seismogram using Matplotlib "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the channel of interest\n",
    "channel = 1330\n",
    "ch_inx = channel - ch_start\n",
    "filtered_data = sp_result[0].data[:, ch_inx]\n",
    "n_samples = filtered_data.shape[0]\n",
    "num_sec = int(n_samples*d_t)\n",
    "t_filt = np.linspace(0, num_sec, n_samples, endpoint=False)\n",
    "\n",
    "# define the scale to apply to the raw data\n",
    "scale_iDAS = float((116*sampling_rate/gauge_length)/1e9) # to strain rate\n",
    "\n",
    "# define the channel range whose mean value you want to subtract from the DAS trace\n",
    "ch_start_demean = 400\n",
    "ch_end_demean = 600\n",
    "demeaned_scaled = (filtered_data - np.mean(sp_result[0].data[:,ch_start_demean-ch_start:ch_end_demean-ch_end], axis=1)) * scale_iDAS\n",
    "\n",
    "# plot the result\n",
    "plt.figure(figsize=(12,8))\n",
    "plt.plot(t_filt, demeaned_scaled, label='Low-freq. Silixa - channel: ' + str(channel))\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.title('Filtered signal using DASDAE - Stage 8')\n",
    "plt.grid(True)\n",
    "\n",
    "# save the figure into the output_figure_folder directory\n",
    "file_name_lowfreq = '/dascore_lowpass_' + str(int(d_t*2)) + 'sec_filter_channel' + str(channel) + '.jpeg'\n",
    "plt.savefig(output_figure_folder + file_name_lowfreq, dpi=600, format='jpeg')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Apply a median filter to remove noise and microseismic events\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define window length (number of samples) for the median filter\n",
    "window_size = 9\n",
    "median_filtered_signal = scipy.ndimage.median_filter(demeaned_scaled, size=window_size)       \n",
    "\n",
    "plt.figure(figsize=(12,8))\n",
    "plt.plot(t_filt, median_filtered_signal, label='Low-freq. Silixa - channel: ' + str(channel))\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(-6e-6, 6e-6)\n",
    "plt.title('Filtered signal using DASDAE - Stage 8')\n",
    "plt.grid(True)\n",
    "\n",
    "# save the figure into the output_figure_folder directory\n",
    "file_name_lowfreq = '/dascore_lowpass_' + str(int(d_t*2)) + 'sec_filtered_median_filtered_channel' + str(channel) + '.jpeg'\n",
    "plt.savefig(output_figure_folder + file_name_lowfreq, dpi=600, format='jpeg')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### b) Waterfall plot using Matplotlib "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the scale to apply to the raw data\n",
    "scale_iDAS = float((116*sampling_rate/gauge_length)/1e9) # to strain rate\n",
    "\n",
    "# get the filtered results from the result spool\n",
    "filtered_data = sp_result[0].data\n",
    "\n",
    "# define the channel range whose mean you want to subtract from the DAS trace\n",
    "ch_start_demean = 400\n",
    "ch_end_demean = 600\n",
    "mean_array = np.mean(filtered_data[:,ch_start_demean-ch_start:ch_end_demean-ch_start], axis=1)  # Calculate the mean along axis 1\n",
    "mean_array = mean_array.reshape(-1, 1)  # Reshape mean_array to match the shape of full_array\n",
    "\n",
    "# demean and scale the results\n",
    "demeaned_scaled_data = (filtered_data - mean_array) * scale_iDAS\n",
    "\n",
    "# define the channles and seconds you want to plot\n",
    "min_sec = 0\n",
    "max_sec = 13990\n",
    "min_ch = 400\n",
    "max_ch = 1355\n",
    "fig_title = \"Silixa's low-freq. DAS - Stage 8\"\n",
    "waterfall_plot(demeaned_scaled_data.T, min_sec, max_sec, min_ch-ch_start, max_ch-ch_start, 1185, 1/d_t, fig_title, output_figure_folder, \"low_freq_raster_plot\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Apply a median filter to remove noise and microseismic events\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define window length (number of samples) for the median filter\n",
    "window_size = 5 # means 5*d_t seconds\n",
    "median_filtered_signal = scipy.ndimage.median_filter(demeaned_scaled_data, size=window_size)  \n",
    "\n",
    "# define the channles and seconds you want to plot\n",
    "min_sec = 0\n",
    "max_sec = 13990\n",
    "min_ch = 400\n",
    "max_ch = 1355\n",
    "fig_title = \"Silixa's low-freq. DAS - Median filtered - Stage 8\"\n",
    "waterfall_plot(median_filtered_signal.T, min_sec, max_sec, min_ch-ch_start, max_ch-ch_start,1185, 1/d_t, fig_title, output_figure_folder, \"low_freq_raster_plot_median_filtered\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### c) Waterfall plot using DASCore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the scale to apply to the raw data\n",
    "scale_iDAS = float((116*sampling_rate/gauge_length)/1e9) # to strain rate\n",
    "\n",
    "# select a sub spool for visualization\n",
    "min_ch = 400\n",
    "max_ch = 1355\n",
    "d_1 = sp_result[0].coords['distance'][min_ch-ch_start] # in meter\n",
    "d_2 = sp_result[0].coords['distance'][max_ch-ch_start] # in meter\n",
    "sub_sp_result = sp_result.select(distance=(d_1, d_2))\n",
    "\n",
    "# get the filtered data from the saved result\n",
    "sub_sp_result_merged = sub_sp_result.chunk(time=None)\n",
    "filtered_data = sub_sp_result_merged[0].data\n",
    "\n",
    "# define the channel range whose mean you want to subtract from the DAS trace\n",
    "ch_start_demean = 400\n",
    "ch_end_demean = 600\n",
    "mean_array = np.mean(filtered_data[:,ch_start_demean-ch_start:ch_end_demean-ch_start], axis=1)  # Calculate the mean along axis 1\n",
    "mean_array = mean_array.reshape(-1, 1)  # Reshape mean_array to match the shape of full_array\n",
    "\n",
    "# demean and scale the results\n",
    "demeaned_scaled_data = (filtered_data - mean_array) * scale_iDAS\n",
    "\n",
    "# make a new patch containing the demeaned_scaled_data data\n",
    "filtered_demeaned_scaled = sub_sp_result_merged[0].new(data=demeaned_scaled_data)\n",
    "\n",
    "# plot\n",
    "filtered_demeaned_scaled.viz.waterfall(scale=0.01) # scale acts as a clipper for data\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
}