eqcorrscan / EQcorrscan

Earthquake detection and analysis in Python.
https://eqcorrscan.readthedocs.io/en/latest/
Other
167 stars 86 forks source link

Problem with low sample rates #515

Closed vSahakian closed 2 years ago

vSahakian commented 2 years ago

I'm using eqcorrscan for a slightly obscure application, on 1 minute sampled turbidity data on an oceanographic cabled array.

When generating templates of turbidity current events, I notice that the shortproc and processing functions are forcing a resample on line 570 of pre_processing.py (https://github.com/eqcorrscan/EQcorrscan/blob/72e769dc5aed2ac1a4e4bc0dd96fe7df5de1e88e/eqcorrscan/utils/pre_processing.py#L570), which is undesirable for my purposes (I'd rather it not resample).

It appears to be because the sampling rate of the data ("tr.stats.sampling_rate") is very low (1/60 Hz), the precision of this float is such that tr.stats.sampling_rate turns out to be 0.01666666753590107, while it should be 0.016666666666666666 (which is what I'm specifying 1/60 as the sampling rate "samp_rate" in the template_gen function).

Is it possible to modify this to handle low sample rates?

I'm working on Python 3.9.12.

calum-chamberlain commented 2 years ago

Thanks for raising this! This is an unusual use case, but hopefully the code can be useful for you.

Do all your traces have the same (imprecise) sampling-rate? If not then you should make sure they are all set to the same value if you expect them to be. Then I would suggest that you set samp_rate to that value. I don't think that the precision of samp_rate is changed anywhere (I hope it isn't), so that should compare as equal. Perhaps you have tried that though and this issue still happens?

It would be good to expose an option to not resample (e.g. by setting samp_rate=None to disable resampling), but in the matched-filter steps the sample-rate of the continuous data are compared to the sample rate of the template to ensure that they are the same and if they are not it will resample. If this is a floating point accuracy bug that gives different results for different epochs then that will be an issue. Keen to know what you have found in your experience and happy to try and help.

vSahakian commented 2 years ago

Thanks for such a fast response!

I dug in a bit more and realized that the problem is in my conversion to a SAC file, so I can use the SAC method in template_gen. I convert csv files to miniseed, then convert and write to SAC (with a pick time), and subsequently template_gen reads these back in. Somewhere in either the writing or reading of the SAC file, the sampling rate gets converted from delta=60.0 to 0.01666666753590107 (instead of 1/60 or 0.01666666666666666).

I agree that forcing it to resample would not be ideal. What are your thoughts on either using np.isclose on line 570 of pre_processing.py (instead of a logical testing equality - the default could be atol=0 which is the same as identical values), or otherwise adding an option to read in mseed files to template_gen? If easier, I could make suggested code changes.

Thank you!

calum-chamberlain commented 2 years ago

You can pass any Stream to most of the template creation functions/methods, but we don't support passing filenames - you have to read the file using obspy first, so you can pass miniseed to those (so long as you read the file first). You would use the from_meta_file method and pass a stream. In general I would avoid the SAC methods!

Also, generally it is better to use the object oriented API (Template and Tribe objects). The functions in template_gen are only there for legacy purposes. Similar to template_gen you should use the from_meta_file method.

Using np.isclose is a nice idea, we could do that, but I'm loath to introduce something that could be used to hide other issues (e.g. file conversion issues).

vSahakian commented 2 years ago

Perfect, the from_meta_file method works! Thank you so much!

calum-chamberlain commented 2 years ago

I really need to get back to finishing the local file tutorial in based on #468! If there is anything else you would like to see in a tutorial then let me know. I'm going to close this, but feel free to reply.