Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update monitoring post-process scrips and tutorials #308

Merged
merged 8 commits into from
Apr 13, 2024
Merged

Update monitoring post-process scrips and tutorials #308

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
42f5030
Add a new ConfigParameter_monitoring class to the parameters for dv/v…
kuanfufeng Mar 18, 2024
c752b28
update attenuation notebook - adding timerange to stackstore
kuanfufeng Apr 12, 2024
035fc97
update test scripts
kuanfufeng Apr 12, 2024
3328f09
Merge branch 'main' into kuanfu_monito2
kuanfufeng Apr 12, 2024
8c29936
update setuptools version
kuanfufeng Apr 12, 2024
4d511be
remove setuptool deps
niyiyu Apr 12, 2024
a766ea6
replace kernelspec
niyiyu Apr 13, 2024
8c37278
replace kernelspec
niyiyu Apr 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
0 src/noisepy/monitoring/esyn_utils.py → src/noisepy/monitoring/attenuation_utils.py
View file
Open in desktop
File renamed without changes.
1,462 changes: 1,462 additions & 0 deletions src/noisepy/monitoring/monitoring_methods.py
View file
Open in desktop

Large diffs are not rendered by default.

1,605 changes: 176 additions & 1,429 deletions src/noisepy/monitoring/monitoring_utils.py
View file
Open in desktop

Large diffs are not rendered by default.

0 src/noisepy/monitoring/esyn_plotting.py → ...oisepy/monitoring/plotting_attenuation.py
View file
Open in desktop
File renamed without changes.
2 changes: 1 addition & 1 deletion tests/test_monitoring_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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,
Expand Down
2 changes: 1 addition & 1 deletion tests/test_stretching.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
94 changes: 65 additions & 29 deletions tutorials/monitoring/attenuation_singlestation.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
]
},
{
Expand All @@ -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,
Expand All @@ -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"
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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",
Expand All @@ -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"
]
},
Expand Down Expand Up @@ -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",
Expand All @@ -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"
]
},
{
Expand Down Expand Up @@ -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",
Expand All @@ -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"
]
},
Expand Down Expand Up @@ -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"
]
Expand All @@ -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",
Expand Down Expand Up @@ -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"
]
},
{
Expand All @@ -551,9 +587,9 @@
],
"metadata": {
"kernelspec": {
"display_name": "noisepy_test",
"display_name": "noisepy_dev",
"language": "python",
"name": "python3"
"name": "noisepy_dev"
},
"language_info": {
"codemirror_mode": {
Expand All @@ -565,7 +601,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.13"
"version": "3.10.14"
}
},
"nbformat": 4,
Expand Down
2 changes: 1 addition & 1 deletion tutorials/monitoring/monitoring_methods_tutorial.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down
Loading