problems with cross-correlation

[Evalb44]

Hello everyone
First of all, thank you very much for this complete and understandable code, and above all for making it accessible to a whole community for the progress of science.

I am currently using this code only for the part concerning the correlation and estimation of seismic velocity changes from ambient noise. Data recovery and H/V processing to determine the specific frequency range on which to work is carried out beforehand.

I process data from three seismometers (RaspberryShake RS3D). The autocorrelations work very well for all three. However, when it comes to cross-correlations between the different stations and also between different components, that poses a problem.
Running the code for the cross-correlation works well, in fact the plots on the graph are generated and the H5 files in the corr folders are also present. The problem comes in the execution of the monitor code:

(
try:
ipynb_path
except NameError: ipynb_path = os.getcwd()

os.chdir(ipynb_path)

ncpus = os.cpu_count()//6

ip = get_ipython()
P:\PRODUCTION\MESURES_CONTINUES\Scripts_python\SEISMIC_survey\TEST.py")
ip.run_cell(f"!mpiexec -n {ncpus} python -m mpi4py P:\PRODUCTION\MESURES_CONTINUES\Scripts_python\SEISMIC_survey\monitor1.py")
)

the only files created are solely from autocorrelations (for example, there are only files such as DV-AM-AM.R343A-R343A.EHE-EHE.npz and not those from cross-correlation).
The problem would therefore only come from this dV part, but I'm having trouble pinpointing the specific problem. Especially as the functions in the Monitor class are not designed to be modified from the outset.
The params.yaml file also seems correct.

Here it is :

Project wide parameters

proj_dir : 'P:\PRODUCTION\MESURES_CONTINUES\Scripts_python\SEISMIC_survey\codes_evanna\data'
save_comps_separately: True
log_subdir : 'log'

log_level: 'WARNING'
fig_subdir : 'figures'

parameters that are network specific

net:
# list of stations used in the project
# type: list of strings or string, wildcards allowed
network : 'AM'
station : ['R343A','RA97C','RA63A']
component : 'Z'

parameters for correlation (emperical Green's function creation)

co:
# subdirectory of 'proj_dir' to store correlation
# type: string
subdir : 'corr'
# times sequences to read for cliping or muting on stream basis
# These should be long enough for the reference (e.g. the standard
# deviation) to be rather independent of the parts to remove
# type: string
read_start : '2024-04-01 00:00:01.0'
read_end : '2024-04-18 00:00:00.0'
# type: float [seconds]
# The longer the faster, but higher RAM usage.
# Note that this is also the length of the correlation batches
# that will be written (i.e., length that will be
# kept in memory before writing to disk)
read_len : 86398
read_inc : 86400

# New sampling rate in Hz. Note that it will try to decimate
# if possible (i.e., there is an integer factor from the
# native sampling_rate)
sampling_rate : 25
# Remove the instrument response, will take substantially more time
remove_response : False

# Method to combine different traces
combination_method : 'betweenStations'

# If you want only specific combinations to be computed enter them here
# In the form [Net0-Net0.Stat0-Stat1]
# This option will only be consider if combination_method == 'betweenStations'
# Comment or set == None if not in use
xcombinations : None

# preprocessing of the original length time series
# these function work on an obspy.Stream object given as first argument
# and return an obspy.Stream object.
preProcessing : [
                {'function':'seismic.correlate.preprocessing_stream.stream_filter',
                  'args':{'ftype':'bandpass',
                          'filter_option':{'freqmin':0.01, #0.01
                                           'freqmax':12.49}}}
                 #{'function':'seismic.correlate.preprocessing_stream.stream_mute',
                 # 'args':{'taper_len':100,
                 #         'mute_method':'std_factor',
                 #         'mute_value':3}}
                ]
# subdivision of the read sequences for correlation
# if this is set the stream processing will happen on the hourly subdivision. This has the
# advantage that data that already exists will not need to be preprocessed again
# On the other hand, computing a whole new database might be slower
# Recommended to be set to True if:
# a) You update your database and a lot of the data is already available (up to a magnitude faster)
# b) corr_len is close to read_len
# Is automatically set to False if you are computing a completely new db
preprocess_subdiv: True
# type: presence of this key
subdivision:
    # type: float [seconds]
    corr_inc : 3600
    corr_len : 3600
    # recombine these subdivisions
    # unused at the time
    # type: boolean
    recombine_subdivision : True
    # delete
    # type: booblean
    delete_subdivision : False

# parameters for correlation preprocessing
# Standard functions reside in seismic.correlate.preprocessing_td and preprocessing_fd, respectively
corr_args : {'TDpreProcessing':[
                              {'function':'seismic.correlate.preprocessing_td.taper',
                                'args':{'type':'cosine_taper', 'p': 0.02}},
                                {'function':'seismic.correlate.preprocessing_td.detrend',
                                'args':{'type':'linear'}},
                               {'function':'seismic.correlate.preprocessing_td.TDfilter',
                               'args':{'type':'bandpass','freqmin':1.2,'freqmax':1.5}},
                                #{'function':'seismic.correlate.preprocessing_td.mute',
                                # 'args':{'taper_len':100.,
                                       # 'threshold':1000, absolute threshold
                                #         'std_factor':3,
                                #         'filter':{'type':'bandpass','freqmin':1.2,'freqmax':1.5},
                                #         'extend_gaps':True}},
                              #  {'function':'seismic.correlate.preprocessing_td.clip',
                              #   'args':{'std_factor':2.5}},
                               {'function':'seismic.correlate.preprocessing_td.signBitNormalization',
                                'args': {}}
                               ],
              # Standard functions reside in seismic.correlate.preprocessing_fd
             'FDpreProcessing':[
                                {'function':'seismic.correlate.preprocessing_fd.spectralWhitening',
                                 'args':{'joint_norm':False}},
                                {'function':'seismic.correlate.preprocessing_fd.FDfilter',
                                 'args':{'flimit':[1.33, 2, 4, 6]}}
                                #  {'function':seismic.correlate.preprocessing_fd.FDsignBitNormalization,
                                # 'args':{}}
                                ],
             'lengthToSave':100,
             'center_correlation':True,      # make sure zero correlation time is in the center
             'normalize_correlation':True,
             'combinations':[]
            }

parameters for the estimation of time differences

dv:
# subfolder for storage of time difference results
subdir : 'vel_change'

# Plotting
plot_vel_change : True

### Definition of calender time windows for the time difference measurements
start_date : '2024-04-01 00:00:00.0'   # %Y-%m-%dT%H:%M:%S.%fZ'
end_date : '2024-04-18 00:00:00.0'
win_len : 3600                         # length of window in which EGFs are stacked
date_inc : 3600                        # increment of measurements

### Frequencies
freq_min : 1.2
freq_max : 1.5

### Definition of lapse time window
tw_start : 1.5     # lapse time of first sample [s]
tw_len : 5.0      # length of window [s] Can be None if the whole (rest) of the coda should be used
sides : 'both'   # options are left (for acausal), right (causal), both, or single (for active source experiments where the first sample is the trigger time)
compute_tt : True  # Computes the travel time and adds it to tw_start (tt is 0 if not xstations). If true a rayleigh wave velocity has to be provided
rayleigh_wave_velocity : 1  # rayleigh wave velocity in km/s, will be ignored if compute_tt=False


### Range to try stretching
stretch_range : 0.03
stretch_steps : 1001
# Return the similarity matrix for each dv/v estimate
# required for post processing. If False, saves some RAM and disk space
return_sim_mat : True

#### Reference trace extraction
#  win_inc : Length in days of each reference time window to be used for trace extraction
# If == 0, only one trace will be extracted
# If > 0, mutliple reference traces will be used for the dv calculation
# Can also be a list if the length of the windows should vary
# See seismic.correlate.stream.CorrBulk.extract_multi_trace for details on arguments
dt_ref : {'win_inc' : 0, 'method': 'mean', 'percentile': 50}

# preprocessing on the correlation bulk or corr_stream before stretch estimation
preprocessing: [
                #{'function': 'pop_at_utcs', 'args': {'utcs': np.array([UTCDateTime()])},
                #{'function': 'select_time', 'args': {'start': (0, 0, 0), 'end': (23, 0, 0), exclude=False}, # include only certain time of day (h, minute, second)
                {'function': 'smooth', 'args': {'wsize': 4, 'wtype': 'hanning', 'axis': 1}}
]

# postprocessing of dv objects before saving and plotting
postprocessing: [
                  # {'function': 'smooth_sim_mat', 'args': {'win_len': 7, exclude_corr_below: 0}}
]

parameters to compute the waveform coherence

wfc:

subfolder for storage of time difference results

subdir : 'wfc'

### Definition of calender time windows for the time difference measurements
start_date : '2024-04-01 00:00:00.0'   # %Y-%m-%dT%H:%M:%S.%fZ'
end_date : '2024-04-18 00:00:00.0'
win_len : 3600                         # length of window in which EGFs are stacked
date_inc : 3600                        # increment of measurements

### Frequencies
# can be lists of same length
freq_min : [0.0625, 0.09, 0.125, 0.1875, 0.25, 0.375, 0.5, 0.75, 1, 1.5, 2, 3, 4, 5, 6, 8]
freq_max : [0.125, 0.18, 0.25, 0.375, 0.5, 0.75, 1, 1.5, 2, 3, 4, 6, 8, 10, 12, 12.49]

### Definition of lapse time window
# can be lists of same length or tw_start: List and tw_len: single value (will be applied to all)
tw_start : [0, 1.25, 2.5, 3.75, 5, 6.25, 7.5, 8.75, 10, 11.25, 12.5, 13.75, 15, 16.25, 17.5, 17.75, 20]     # lapse time of first sample [s]
tw_len : 5       # length of window [s]

#### Reference trace extraction
#  win_inc : Length in days of each reference time window to be used for trace extraction
# If == 0, only one trace will be extracted
# If > 0, mutliple reference traces will be used for the dv calculation
# Can also be a list if the length of the windows should vary
# See seismic.correlate.stream.CorrBulk.extract_multi_trace for details on arguments
dt_ref : {'win_inc' : 0, 'method': 'mean', 'percentile': 50}

# preprocessing on the correlation bulk before stretch estimation
preprocessing: [
                {'function': 'smooth', 'args': {'wsize': 4, 'wtype': 'hanning', 'axis': 1}}
]

### SAVING
# save components separately or only their average?
save_comps: False

Could the problem be more in the data itself?
Thanks in advance for any help.
Best regards

[PeterMakus]

Hi @Evalb44
Thank you for getting in touch and for the flowers :)
To me it seems that the problem is in your parameter file. Specifically here:

parameters that are network specific
net:
# list of stations used in the project
# type: list of strings or string, wildcards allowed
network : 'AM'
station : ['R343A','RA97C','RA63A']
component : 'Z'

With this configuration, it will only use data from the vertical/Z channel and, thus, no cross-components can be computed

If you change component to allow all (i.e. to, '*' or '?'), it should work.

I will close the discussion, if this doesn't solve the issue, feel free to comment again and it will be reopened.

Cheers,
Peter

[Evalb44]

Hello,
thank you very much for your reply.
However, it still doesn't work if you change the Z to '*' or '?'....
As I said before, it still works fine for the first part of the generation of the cross correlation files and graphs but when I run the monitor1.py script as shared in the first post it displays :

try:
ipynb_path
except NameError: ipynb_path = os.getcwd()

os.chdir(ipynb_path)

ncpus = os.cpu_count()//6

ip = get_ipython()

ip.run_cell(f"!mpiexec -n {ncpus} python -m mpi4py P:\PRODUCTION\MESURES_CONTINUES\Scripts_python\SEISMIC_survey\monitor1.py")

0it [00:00, ?it/s]
0it [00:00, ?it/s]

0it [00:00, ?it/s]
0it [00:00, ?it/s]

0it [00:00, ?it/s]
0it [00:00, ?it/s]

0it [00:00, ?it/s]
0it [00:00, ?it/s]

0it [00:00, ?it/s]
0it [00:00, ?it/s]
Out[14]: <ExecutionResult object at 1182a7fda90, execution_count=None error_before_exec=None error_in_exec=None info=<ExecutionInfo object at 11780184690, raw_cell="!mpiexec -n 5 python -m mpi4py P:\PRODUCTION\MESUR.." store_history=False silent=False shell_futures=True cell_id=None> result=None>

However, 'mpiexec' is correctly installed on the system, and I've already run MPI commands to generate the correlations, for example. The path is the right one, it generates the files for the autocorrelations correctly. I don't really understand this error...
have you ever seen this error?
it's as if it didn't take into account or couldn't find the files from the correlate script.

[PeterMakus]

Hi again,
This indicates that SeisMIC could not find any correlation files (*.h5 files) that fit the search criteria. MPI seems to work fine.

Could you check the content of the following folder:
P:\PRODUCTION\MESURES_CONTINUES\Scripts_python\SEISMIC_survey\codes_evanna\data\corr
and send me the list of files that are in the folder?

[Evalb44]

here are the files found in the :

[PeterMakus]

Hmm ok. This seems fine actually.

After going through the source code, I found something that could cause a bug.

Could you follow these three steps:

1. Do not change anything in your parameter file but the logging level, which should be changed to "DEBUG". Then, post the output here.

2. Could you make the following change to your parameters before executing the monitoring (as a workaround):

parameters that are network specific
net:
# list of stations used in the project
# type: list of strings or string, wildcards allowed
network : 'AM'
station : '*'
component : '*'

3. If this does not work, could you try the following: Instead of using the compute_velocity_change_bulk()``-method, could you use the compute_velocity_change()``` method of the Monitor class on one of your files and see if this works.

I am sorry for the inconvenience, I am trying to isolate the source of this.

[Evalb44]

Hello,
So I tested the steps:
Step 1: It doesn't find any combinations just by running the code with the 'DEBUG' logging level. I have this:

2024-06-18 09:27:02,848 - INFO - Found correlation data for the following station and network combinations []
2024-06-18 09:27:02,848 - INFO - Found correlation data for the following station and network combinations []
2024-06-18 09:27:02,848 - INFO - Correlation data has been found for the following station and network combinations []
2024-06-18 09:27:02,848 - INFO - Correlation data has been found for the following station and network combinations []
2024-06-18 09:27:02,848 - INFO - Correlation data has been found for the following station and network combinations []

0it [00:00, ?it/s]
0it [00:00, ?it/s]

0it [00:00, ?it/s]
0it [00:00, ?it/s]

0it [00:00, ?it/s]
0it [00:00, ?it/s]

0it [00:00, ?it/s]
0it [00:00, ?it/s]

0it [00:00, ?it/s]
0it [00:00, ?it/s]
Out[3] : <ExecutionResult object at 23116082b50, execution_count=None error_before_exec=None error_in_exec=None info=<ExecutionInfo object at 230b27bc1d0, raw_cell="! mpiexec -n 5 python -m mpi4py P:\PRODUCTION\MESUR..." store_history=False silent=False shell_futures=True cell_id=None> result=None>

Step 2: By modifying the parameters as you indicated (station: '' and component: ''), and running the monitor code, it works ! :

runfile('P:/PRODUCTION/MESURES_CONTINUES/Scripts_python/SEISMIC_survey/monitor1.py', wdir='P:/PRODUCTION/MESURES_CONTINUES/Scripts_python/SEISMIC_survey')
2024-06-18 10:03:00,540 - INFO - Found correlation data for the following station and network combinations ['AM-AM.R343A-RA63A', 'AM-AM.R343A-RA63A', 'AM-AM.R343A-RA63A', 'AM-AM.R343A-RA63A', 'AM-AM.R343A-RA63A', 'AM-AM.R343A-RA63A', 'AM-AM.R343A-RA63A', 'AM-AM.R343A-RA63A', 'AM-AM.R343A-RA63A', 'AM-AM.R343A-RA97C', 'AM-AM.R343A-RA97C', 'AM-AM.R343A-RA97C', 'AM-AM.R343A-RA97C', 'AM-AM.R343A-RA97C', 'AM-AM.R343A-RA97C', 'AM-AM.R343A-RA97C', 'AM-AM.R343A-RA97C', 'AM-AM.R343A-RA97C', 'AM-AM.RA63A-RA97C', 'AM-AM.RA63A-RA97C', 'AM-AM.RA63A-RA97C', 'AM-AM.RA63A-RA97C', 'AM-AM.RA63A-RA97C', 'AM-AM.RA63A-RA97C', 'AM-AM.RA63A-RA97C', 'AM-AM.RA63A-RA97C', 'AM-AM.RA63A-RA97C']
0%| | 0/27 [00:00<?, ?it/s]2024-06-18 10:03:00,894 - INFO - Computing velocity change for file: P:\PRODUCTION\MESURES_CONTINUES\Scripts_python\SEISMIC_survey\codes_evanna\data\corr\AM-AM.R343A-RA63A.EHE-EHE.h5 and channel:EHE-EHE
saving to AM-AM_R343A-RA63A_EHE-EHE

Windows fatal exception: code 0xc0000374

Thread 0x00005058 (most recent call first):
File "D:\Projets\Py_Envs\py_seisnoise\Lib\threading.py", line 324 in wait
File "D:\Projets\Py_Envs\py_seisnoise\Lib\threading.py", line 622 in wait
File "D:\Projets\Py_Envs\py_seisnoise\Lib\site-packages\tqdm_monitor.py", line 60 in run
File "D:\Projets\Py_Envs\py_seisnoise\Lib\threading.py", line 1038 in _bootstrap_inner
File "D:\Projets\Py_Envs\py_seisnoise\Lib\threading.py", line 995 in _bootstrap

Main thread:
Current thread 0x000058b8 (most recent call first):
File "D:\Projets\Py_Envs\py_seisnoise\Lib\logging_init_.py", line 624 in formatTime
File "D:\Projets\Py_Envs\py_seisnoise\Lib\logging_init_.py", line 689 in format
File "D:\Projets\Py_Envs\py_seisnoise\Lib\logging_init_.py", line 953 in format
File "D:\Projets\Py_Envs\py_seisnoise\Lib\logging_init_.py", line 1110 in emit
File "D:\Projets\Py_Envs\py_seisnoise\Lib\logging_init_.py", line 1230 in emit
File "D:\Projets\Py_Envs\py_seisnoise\Lib\logging_init_.py", line 978 in handle
File "D:\Projets\Py_Envs\py_seisnoise\Lib\logging_init_.py", line 1706 in callHandlers
File "D:\Projets\Py_Envs\py_seisnoise\Lib\logging_init_.py", line 1644 in handle
File "D:\Projets\Py_Envs\py_seisnoise\Lib\logging_init_.py", line 1634 in log
File "D:\Projets\Py_Envs\py_seisnoise\Lib\logging_init.py", line 1489 in info
File "D:\Projets\Py_Envs\py_seisnoise\Lib\site-packages\seismic\monitor\monitor.py", line 190 in compute_velocity_change
File "D:\Projets\Py_Envs\py_seisnoise\Lib\site-packages\seismic\monitor\monitor.py", line 356 in compute_velocity_change_bulk
File "p:\production\mesures_continues\scripts_python\seismic_survey\monitor1.py", line 19 in
File "D:\Projets\Py_Envs\py_seisnoise\Lib\site-packages\spyder_kernels\py3compat.py", line 356 in compat_exec
File "D:\Projets\Py_Envs\py_seisnoise\Lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 473 in exec_code
File "D:\Projets\Py_Envs\py_seisnoise\Lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 615 in _exec_file
File "D:\Projets\Py_Envs\py_seisnoise\Lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 528 in runfile
File "C:\Users\blancharde\AppData\Local\Temp\2\ipykernel_12104\1848473341.py", line 1 in

As you can see, it finds the correlation data this time for each pair of stations. It calculates the change in seismic velocity for the different files. However, the program does not complete the traceback. If you only run the monitor script, the kernel inevitably reboots after only one or two files have been generated. Perhaps this is due to a memory problem? But the files are not that large.

Cheers,
Evanna

[Evalb44]

for example, when I run all the code, I only have these files :

[PeterMakus]

Hey Evanna
Ok,
So great to hear it works when modifying the parameters, so there was indeed a bug in SeisMIC, which I will fix in a future version that I am just working on. Thank you for finding it :)

Concerning your second issue, this seems to be related to some multicore threading that is internally invoked by one of the dependencies, so not something that happens internally in SeisMIC's source code. We have to do some trial and error here, unfortunately.

Try the following steps:

1. Set the logging level back to WARNING and try again
2. Delete your old correlation files (h5 files) and recompute
3. reinstall hdf5 pip install --upgrade --force-reinstall h5py
4. reinstall mpi and mpi4py (as above)
5. try reinstalling the whole SeisMIC environment.

If none of these helps, it probably points towards an issue in one of your system files (e.g., your system-wide MPI installation).

I would also recommend googling the error messages you posted here. As I said, this is not an error that is raised by SeisMIC itself, so I can only be of limited help here.

[Evalb44]

Hello,
I'd like to apologize for not replying sooner, I was actually trying to fix this problem. Indeed, I think it is a problem in the cores, at least the number of cores that is used in the process and as I work on a shared server, this is surely the problem, according to the tests conducted.
In any case, I'd like to thank you very much for your help from the start. I'm very happy if it has helped you to understand some of the things that went right or wrong in SeisMIC.
Good continuation and maybe I'll see you soon!

Best regards,
Evanna