thomaslecocq / seismorms Goto Github PK

I assume that for the clock plots, the 'after lockdown' subplot also takes into account data after there is a potential lifting of the lockdown situation. In these cases, such as in China, where the noise level has risen back to pre-lockdown levels, I am just wondering if the comparison between before and after can be a little bit biased by a long end-of-lockdown period.

Maybe we a separate variable, like the bans dict to set end of lockdown dates?

Please say if this isn't clear!

In the underlying plotting code, seismosocialdistancing.py, there is now a matplotlib warning raised when plotting the clock plots. Matplotlib version is "matplotlib-base 3.3.3", and Pandas is 1.1.4, both from conda-forge. This isn't a problem for the code itself, it's the imported packages, and has been raised as an issue there for a few months. However, it would prove un-nerving to users unfamiliar to these warnings, and can be alleviated.

The warning is "UserWarning: FixedFormatter should only be used together with FixedLocator".
THe error can be alleviated by a dummy call - I added the following after line 605:
ax.set_yticks(ax.get_yticks())

@claudiodsf @FMassin @megies

should I include this ipynb that shows the correctness of the computations: ?? I mean we know it's Parsaval's identity, but still:

http://msnoise.org/tmp/Amplitude%20&%20Spectrum.html

And shows advantages:

1. computing the PSD doesn't require pre_filtering the whole trace
2. PSD doens't need to remove_response (ok this could be done by just tr.data / total_sensitivity, right?)
3. RMS derived from PSD (calculated once) can be extracted for any frequency band, even wide/narrow

PS: in the example above, you'll see a ~33% noise reduction when it snows in Brussels (2 to 5cm max) and all the trafic is stopped.. Kinda matches what we see today :-)

Pandas version should be > 0.23.4 to support clock plots.

Somehow it might be an issue when using (outdated?) anaconda version... Anyhow, in case it might help somebody else, I had the following error for execution of clock plots while using pandas-0.23.4:

AttributeError: 'DataFrame' object has no attribute 'droplevel'

and I manage to fix it with:

~/anaconda3/bin/pip  install --upgrade  pandas

Which gave me pandas-0.25.3 and works like a charm.

For disk space saving, could we have an option to only download the mseed if the corresponding .npz file isn't already existing? As far as I can see the .npz file doesn't actually depend on user-defined parameters so I don't think these ever really will change (unless the waveform data changes).
That way, the user could delete the .mseed files at the end of every run. This might be quite useful if regularly updating a station or running for different time periods.

As discussed in the network channel - some stations only plot the most recent period of continuous data after a data gap, despite older data being available (green bars show as available and .npz files exist). This does not seem to be consistent between stations (e.g. station RAEBE or R2E51 do not plot all data, but station R5982 plots fine with a big data gap).

The clock plots are amazing and my expectations may be wrong, but:

1. I would expect the three graphs to all have the same y-axis scale, but the maximum amplitude in the top graph is 0.226nm and in the bottom graphs, it is nearly 40nm.
2. The scale in the polar graphs goes down to -10nm, where I would expect the centre of the graph to be zero.
3. The colours are different between the upper graph and the lower two graphs, making it hard to compare, so Wedensday is red in the top graph and green in the lower two.

I was having an issue with Line 10 in Step 5, which was resulting in this (shortened here) error:

ValueError: Loading PPSD results saved with ObsPy versions < 1.2 requires setting the allow_pickle parameter of PPSD.load_npz to True (needs numpy>=1.10).

So I changed Line 10 in Step 5 to:
ppsds[mseedid].add_npz(fn, allow_pickle=True)

@ThomasLecocq, nice work!!, however if we had data with gaps due to the transmission the script will stop at each gap, can we full with zeros if we find those gaps.

Hi,

in SeismoSocialDistancing.ipynb

days = [ 'Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday']

should be defined in stack_wday_time(df)

I am trying to run the code SeismoRMS for Power spectral density. It is successfully running in your given portal. But while running on my computer, it is giving some error. The difference is your data and resp is taken online, and I am giving data from the folder in my computer. Could you help me to solve this?

Error is coming while plotting. Without this, I am not able to proceed.

Exception Traceback (most recent call last)
Input In [13], in <cell line: 1>()
----> 1 [ppsd.plot(max_percentage=10) for mseedid, ppsd in ppsds.items()]
2 [ppsd.plot_temporal(0.10) for mseedid, ppsd in ppsds.items()]
3 [ppsd.plot_spectrogram(clim=(-160,-100)) for mseedid, ppsd in ppsds.items()]

Input In [13], in (.0)
----> 1 [ppsd.plot(max_percentage=10) for mseedid, ppsd in ppsds.items()]
2 [ppsd.plot_temporal(0.10) for mseedid, ppsd in ppsds.items()]
3 [ppsd.plot_spectrogram(clim=(-160,-100)) for mseedid, ppsd in ppsds.items()]

File ~\anaconda3\lib\site-packages\obspy\signal\spectral_estimation.py:1894, in PPSD.plot(self, filename, show_coverage, show_histogram, show_percentiles, percentiles, show_noise_models, grid, show, max_percentage, period_lim, show_mode, show_mean, cmap, cumulative, cumulative_number_of_colors, xaxis_frequency, show_earthquakes)
1821 """
1822 Plot the 2D histogram of the current PPSD.
1823 If a filename is specified the plot is saved to this file, otherwise
(...)
1891 in Hertz as opposed to the default of period in seconds.
1892 """
1893 import matplotlib.pyplot as plt
-> 1894 self.__check_histogram()
1895 fig = plt.figure()
1896 fig.ppsd = AttribDict()

File ~\anaconda3\lib\site-packages\obspy\signal\spectral_estimation.py:820, in PPSD.__check_histogram(self)
818 else:
819 msg = 'No data accumulated'
--> 820 raise Exception(msg)

Exception: No data accumulated

Nice work Thomas (and others) ...

Here is something I noticed about how you cache the downloaded files.

In STEP3 of SeismoSD.ipynb You currently don't use cached data from "today" but I think this runs the risk that tomorrow you would consider this file to be just fine even if it is the incomplete data from running the notebook in the middle of the day. One simple check might be to check the creation date of the cached file though that is not really a check that the file is not corrupted.

The same argument applies to the downloaded data from step 3 and presumably the npz files in Step 4 too. There is no check to see if the npz file is out of date compared to the mseed file which would help, I think.

Hi, great that you share this! When I try to run the example code I get the following error:

c = Client("ODC") FDSNException: No FDSN services could be discovered at 'http://www.orfeus-eu.org'. This could be due to a temporary service outage or an invalid FDSN service address.

What would be alternatives here?

Save processing times by splitting up the daylist into nproc chunks.
I have a recipe for this for running in a standalone python script, but haven't tested this using a Jupyter Notebook.

I know @FMassin that this is going to hurt your ego... 😄