Update monitoring post-process scrips and tutorials (#308)

* Add a new ConfigParameter_monitoring class to the parameters for dv/v and attenuation * update attenuation notebook - adding timerange to stackstore * update test scripts * update setuptools version * remove setuptool deps Signed-off-by: Yiyu Ni <[email protected]> * replace kernelspec Signed-off-by: Yiyu Ni <[email protected]> * replace kernelspec Signed-off-by: Yiyu Ni <[email protected]> --------- Signed-off-by: Yiyu Ni <[email protected]> Co-authored-by: Yiyu Ni <[email protected]>
noisepy · Apr 13, 2024 · 9a8684c · 9a8684c
1 parent 23e0dbf
commit 9a8684c
Show file tree

Hide file tree

Showing 8 changed files with 1,706 additions and 1,461 deletions.
diff --git a/src/noisepy/monitoring/esyn_utils.py → src/noisepy/monitoring/attenuation_utils.py b/src/noisepy/monitoring/esyn_utils.py → src/noisepy/monitoring/attenuation_utils.py
diff --git a/src/noisepy/monitoring/monitoring_methods.py b/src/noisepy/monitoring/monitoring_methods.py
diff --git a/src/noisepy/monitoring/monitoring_utils.py b/src/noisepy/monitoring/monitoring_utils.py
diff --git a/src/noisepy/monitoring/esyn_plotting.py → ...oisepy/monitoring/plotting_attenuation.py b/src/noisepy/monitoring/esyn_plotting.py → ...oisepy/monitoring/plotting_attenuation.py
diff --git a/tests/test_monitoring_utils.py b/tests/test_monitoring_utils.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 
-from noisepy.monitoring.monitoring_utils import (
+from noisepy.monitoring.monitoring_methods import (
     dtw_dvv,
     mwcs_dvv,
     stretching,

diff --git a/tests/test_stretching.py b/tests/test_stretching.py
@@ -4,7 +4,7 @@
 import pytest
 from obspy.signal.invsim import cosine_taper
 
-from noisepy.monitoring.monitoring_utils import (
+from noisepy.monitoring.monitoring_methods import (
     mwcs_dvv,
     stretching,
     stretching_vect,

diff --git a/tutorials/monitoring/attenuation_singlestation.ipynb b/tutorials/monitoring/attenuation_singlestation.ipynb
@@ -35,10 +35,14 @@
     "import matplotlib.pyplot as plt\n",
     "\n",
     "from obspy.signal.filter import bandpass\n",
-    "from noisepy.monitoring.esyn_utils import *\n",
-    "#from noisepy.monitoring.esyn_plotting import plot_waveforms\n",
+    "from noisepy.monitoring.attenuation_utils import *\n",
+    "from noisepy.monitoring.monitoring_utils import *\n",
+    "#from noisepy.monitoring.plotting_attenuation import plot_waveforms\n",
     "from noisepy.seis.noise_module import mad\n",
-    "from noisepy.seis.io.asdfstore import ASDFStackStore"
+    "from noisepy.seis.io.asdfstore import ASDFStackStore\n",
+    "\n",
+    "from datetimerange import DateTimeRange\n",
+    "from datetime import datetime, timezone"
    ]
   },
   {
@@ -64,6 +68,18 @@
     "wave_store = ASDFStackStore(data_path)\n"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c38c5726",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "start = datetime(2019, 2, 1, tzinfo=timezone.utc)\n",
+    "end = datetime(2019, 2, 2, tzinfo=timezone.utc)\n",
+    "timerange = DateTimeRange(start, end)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -78,7 +94,7 @@
     "stations = set(pair[0] for pair in pairs)\n",
     "print('Stations: ', stations)\n",
     "\n",
-    "sta_stacks = wave_store.read_bulk(None, pairs) # no timestamp used in ASDFStackStore\n",
+    "sta_stacks = wave_store.read_bulk(timerange, pairs) # no timestamp used in ASDFStackStore\n",
     "stacks = sta_stacks[0][1]\n",
     "print(\"Target pair: \",sta_stacks[0])\n",
     "target_pair=sta_stacks[0][0][0].network+\".\"+sta_stacks[0][0][0].name"
@@ -179,6 +195,27 @@
     "--> normalized MS envelope is referred to as the observed energy density Eobs "
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "45721a3d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "config_monito = ConfigParameters_monitoring() # default config parameters which can be customized\n",
+    "\n",
+    "# --- parameters for measuring attenuation ---\n",
+    "config_monito.smooth_winlen = 5.0  # smoothing window length\n",
+    "config_monito.cvel = 2.6  # Rayleigh wave velocities over the freqency bands\n",
+    "config_monito.atten_tbeg = 2.0\n",
+    "config_monito.atten_tend = 10.0 # or it will be determined by a ratio of MAD value\n",
+    "config_monito.ratio = 5.0 # ratio for determining noise level by Mean absolute deviation (MAD)\n",
+    "\n",
+    "# basic parameters\n",
+    "config_monito.freq = [0.5, 1.0, 2.0, 4.0]  # targeted frequency band for waveform monitoring\n",
+    "nfreq = len(config_monito.freq) - 1\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -210,8 +247,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "freq = [0.6, 0.8, 2,4]  # targeted frequency band for waveform monitoring\n",
-    "nfreq = len(freq) - 1\n",
     "\n",
     "MSE=np.ndarray((fnum,num_cmp,nfreq+1,npts)) # filtered two-side averaged stack CF\n",
     "\n",
@@ -223,14 +258,14 @@
     "        print(fname[aa],ccomp)\n",
     "        \n",
     "        for fb in range(nfreq):\n",
-    "            fmin=freq[fb]\n",
-    "            fmax=freq[fb+1]\n",
+    "            fmin=config_monito.freq[fb]\n",
+    "            fmax=config_monito.freq[fb+1]\n",
     "            tt = np.arange(0, npts) * dt\n",
     "            data = stackf[aa][ncmp][1]\n",
     "            dafbp[fb] = bandpass(data, fmin, fmax, int(1 / dt), corners=4, zerophase=True)\n",
     "        \n",
     "        MSE[aa][ncmp]=[stackf[aa][ncmp][0],dafbp[0],dafbp[1],dafbp[2]] \n",
-    "        plot_filtered_waveforms(freq, stackf[aa][ncmp][0], dafbp, fname[aa], ccomp)\n",
+    "        plot_filtered_waveforms(config_monito.freq, stackf[aa][ncmp][0], dafbp, fname[aa], ccomp)\n",
     "\n"
    ]
   },
@@ -282,16 +317,17 @@
     "msv[:][:][:][:]=0.\n",
     "msv_mean[:][:][:]=0.\n",
     "\n",
-    "winlen=[8,4,2]   # smoothing window lengths corresponding to the frequency bands\n",
+    "# smoothing window lengths corresponding to the frequency bands\n",
+    "winlen=np.full(nfreq, config_monito.smooth_winlen)\n",
     "for aa in range(fnum):\n",
     "\n",
     "    for ncmp in  range(len(comp_arr)):\n",
     "        ccomp=comp_arr[ncmp]\n",
     "        msv[aa][ncmp][0]=MSE[aa][ncmp][0][:]\n",
     "        for fb in range(nfreq):\n",
     "            data=MSE[aa][ncmp][fb+1][:]\n",
-    "            fmin=freq[fb]\n",
-    "            fmax=freq[fb+1]\n",
+    "            fmin=config_monito.freq[fb]\n",
+    "            fmax=config_monito.freq[fb+1]\n",
     "            \n",
     "            para = { 'winlen':winlen[fb], 'dt':dt , 'npts': len(data)}\n",
     "            msv[aa][ncmp][fb+1]=get_smooth(data, para)\n",
@@ -301,13 +337,13 @@
     "    # get average waveforms from components\n",
     "    msv_mean[aa][0]=msv[aa][0][0][:]\n",
     "    for fb in range(nfreq):\n",
-    "        fmin=freq[fb]\n",
-    "        fmax=freq[fb+1]\n",
+    "        fmin=config_monito.freq[fb]\n",
+    "        fmax=config_monito.freq[fb+1]\n",
     "        for ncmp in range(len(comp_arr)):\n",
     "            msv_mean[aa][fb+1]+=msv[aa][ncmp][fb+1][:]\n",
     "        msv_mean[aa][fb+1]=msv_mean[aa][fb+1]/len(comp_arr)\n",
     "        \n",
-    "    plot_envelope(comp_arr,freq,msv[aa],msv_mean[aa],fname[aa],vdist[aa])\n"
+    "    plot_envelope(comp_arr, config_monito.freq, msv[aa], msv_mean[aa], fname[aa], vdist[aa])\n"
    ]
   },
   {
@@ -351,14 +387,14 @@
     "fmsv_mean=np.ndarray((fnum,nfreq+1,indx+1))\n",
     "\n",
     "# noise level setting\n",
-    "ratio=3\n",
+    "ratio=config_monito.ratio\n",
     "level=np.ndarray((fnum,nfreq,1))\n",
     "twinbe=np.ndarray((fnum,nfreq,2))\n",
     "\n",
     "for aa in range (fnum):\n",
     "    for fb in range(nfreq):\n",
-    "        fmin=freq[fb]\n",
-    "        fmax=freq[fb+1]            # stack positive and negative lags  \n",
+    "        fmin=config_monito.freq[fb]\n",
+    "        fmax=config_monito.freq[fb+1]            # stack positive and negative lags  \n",
     "        sym=get_symmetric(msv_mean[aa][fb+1],indx)\n",
     "        data_sym[fb]=sym\n",
     "        Val_mad=mad(sym)\n",
@@ -371,7 +407,7 @@
     "                print(aa,fb,pt,sym[pt],level[aa][fb],twinbe[aa][fb])\n",
     "                break\n",
     "    fmsv_mean[aa]=[msv[aa][0][0][indx:],data_sym[0],data_sym[1],data_sym[2]]\n",
-    "    plot_fmsv_waveforms(freq,fmsv_mean[aa],fname[aa],level[aa],twinbe[aa])\n",
+    "    plot_fmsv_waveforms(config_monito.freq,fmsv_mean[aa],fname[aa],level[aa],twinbe[aa])\n",
     "\n"
    ]
   },
@@ -439,7 +475,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cvel=[2.6, 2.0, 1.8]    # Rayleigh wave velocities over the freqency bands\n",
+    "cvel=np.full(nfreq, config_monito.cvel)    # Rayleigh wave velocities over the freqency bands\n",
     "mfpx=np.zeros(1)        # mean_free_path search array\n",
     "intby=np.zeros(30)      # intrinsic_b search array\n"
    ]
@@ -459,8 +495,8 @@
     "SSR[:][:][:]=0.\n",
     "\n",
     "for fb in range(nfreq):\n",
-    "    fmin=freq[fb]\n",
-    "    fmax=freq[fb+1]\n",
+    "    fmin=config_monito.freq[fb]\n",
+    "    fmax=config_monito.freq[fb+1]\n",
     "    c=cvel[fb]\n",
     "    SSR_final[:][:]=0.\n",
     "    vdist[:]=0.000001  # To avoid zero value at denominator\n",
@@ -527,17 +563,17 @@
     "aa=0\n",
     "r=np.take(vdist[aa],0)+0.000001\n",
     "for fb in range(nfreq): \n",
-    "    fmin=freq[fb]\n",
-    "    fmax=freq[fb+1]\n",
+    "    fmin=config_monito.freq[fb]\n",
+    "    fmax=config_monito.freq[fb+1]\n",
     "\n",
     "    # parameters for getting optimal value from the sum of squared residuals (SSR) between Eobs and Esyn \n",
     "    para={ 'fb':fb, 'fmin':fmin, 'fmax':fmax, 'vdist':vdist, 'npts':npts, 'dt':dt, 'cvel':c, 'filenum':aa, \\\n",
     "        'mfp':mfpx, 'intb':intby,'twin':twinbe, 'fmsv':fmsv_mean, 'SSR':SSR , 'sta':target_pair}\n",
     "    # call function get_optimal\n",
     "    result_intb[fb], result_mfp[fb], Eobs[fb], Esyn[fb] = get_optimal_Esyn(fnum,para)\n",
     "    \n",
-    "    plot_fitting_result(result_mfp[fb],result_intb[fb],fmsv_mean[aa][0][:],\n",
-    "                        Eobs[fb],Esyn[fb],target_pair,vdist[aa],twinbe[aa][fb],fmin,fmax)\n"
+    "    plot_fitting_result(result_mfp[fb], result_intb[fb], fmsv_mean[aa][0][:],\n",
+    "                        Eobs[fb], Esyn[fb], target_pair, vdist[aa], twinbe[aa][fb], fmin, fmax)\n"
    ]
   },
   {
@@ -551,7 +587,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "noisepy_test",
+   "display_name": ".env",
    "language": "python",
    "name": "python3"
   },
@@ -565,7 +601,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.13"
+   "version": "3.10.14"
   }
  },
  "nbformat": 4,

diff --git a/tutorials/monitoring/monitoring_methods_tutorial.ipynb b/tutorials/monitoring/monitoring_methods_tutorial.ipynb
@@ -143,7 +143,7 @@
     "import os, glob\n",
     "import numpy as np\n",
     "import h5py\n",
-    "from noisepy.monitoring.monitoring_utils import (\n",
+    "from noisepy.monitoring.monitoring_methods import (\n",
     "    wcc_dvv, \n",
     "    stretching, \n",
     "    dtw_dvv, \n",
@@ -1383,7 +1383,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": ".env",
    "language": "python",
    "name": "python3"
   },