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.

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

PPSD: Add earthquake models from Clinton & Heaton (2002) #2455

Open
wants to merge 5 commits into
base: master
from

Conversation

@FMassin
Copy link

commented Sep 13, 2019

This edit adds earthquake models within PPSD plot as in Clinton & Heaton (2002). This is very tentative and probably not in a correct form. I am expecting to receive (lot of) suggestions on how to make this correct...

Clinton, J. F., & Heaton, T. H. (2002). Potential Advantages of a Strong-motion Velocity Meter over a Strong-motion Accelerometer. Seismological Research Letters, 73(3), 332–342. http://doi.org/10.1785/gssrl.73.3.332

PR Checklist

  • Correct base branch selected? master for new features, maintenance_... for bug fixes
  • This PR is not directly related to an existing issue (which has no PR yet).
  • If the PR is making changes to documentation, docs pages can be built automatically.
    Just remove the space in the following string after the + sign: "+ DOCS"
  • If any network modules should be tested for the PR, add them as a comma separated list
    (e.g. clients.fdsn,clients.arclink) after the colon in the following magic string: "+TESTS:"
    (you can also add "ALL" to just simply run all tests across all modules)
  • All tests still pass.
  • Any new features or fixed regressions are be covered via new tests.
  • Any new or changed features have are fully documented.
  • Significant changes have been added to CHANGELOG.txt .
  • First time contributors have added your name to CONTRIBUTORS.txt .
@megies

This comment has been minimized.

Copy link
Member

commented Sep 13, 2019

Can you post an example picture what it looks like here in comments?

@megies megies added the .signal label Sep 13, 2019
@FMassin

This comment has been minimized.

Copy link
Author

commented Sep 13, 2019

It seems like my coding style is complete garbage, sorry... Will fix it on Monday.

@aringler-usgs

This comment has been minimized.

Copy link
Member

commented Sep 13, 2019

Your plot has units of (m/s^2))^2/Hz. However, I think the earthquake spectra don't have density units. So you would need to integrate the PSDs using a band-width. In the Clinton and Heaton paper, you can see an example of the integration in Figure 5.

@FMassin FMassin closed this Sep 13, 2019
@FMassin FMassin reopened this Sep 13, 2019
@FMassin

This comment has been minimized.

Copy link
Author

commented Sep 13, 2019

you would need to integrate the PSDs using a band-width.

Absolutely correct, now with Eq1 from Cauzzi & Clinton (2013):
image

@megies

This comment has been minimized.

Copy link
Member

commented Sep 16, 2019

+DOCS to trigger automated docs build

@megies megies changed the title Update spectral_estimation.py with earthquake models from Clinton & Heaton (2002) PPSD: Add earthquake models from Clinton & Heaton (2002) Sep 16, 2019
Copy link
Author

left a comment

+DOCS to trigger automated docs build

I am not sure what that means.

@megies

This comment has been minimized.

Copy link
Member

commented Sep 23, 2019

I am not sure what that means.

@FMassin It's a magic command in the comment, which will trigger an automated docs build, see below in the PR status, where you can follow this link to check how your changes look in the documentation: http://docs.obspy.org/pull-requests/2455

Screenshot from 2019-09-23 11-10-58

http://docs.obspy.org/pull-requests/2455/packages/autogen/obspy.signal.spectral_estimation.PPSD.plot.html

  • filename extension .bib missing, so our build scripts did not parse your bib file
  • second citation misses bib file
@megies megies added this to the 1.2.0 milestone Sep 23, 2019
@megies megies added this to In Progress in Release 1.2.0 Sep 23, 2019
@FMassin

This comment has been minimized.

Copy link
Author

commented Sep 24, 2019

Sorry for all these mistakes.

  • filename extension .bib missing, so our build scripts did not parse your bib file
  • second citation misses bib file
@megies

This comment has been minimized.

Copy link
Member

commented Sep 24, 2019

Sorry for all these mistakes.

No problem at all.

OK but I regret to alter the original data set.

Oh, OK, I thought it was kind of a glitch since this was the only dataset that had a name as a third key. If you want to keep that one, maybe a work around could be using empty strings on the other events. It's just a bit awkward if the shape of the keys changes across the dictionary.

earthquakes = {
    (1.5, 10, ""): [[7.0700000e-01, 1.4140000e+00, 2.8280000e+00, 5.6600000e+00, 
                     1.1300000e+01, 2.2600000e+01],
                    [1.1940722e-06, 5.9652811e-06, 4.2022150e-05, 2.1857010e-04, 
                     6.4011669e-04, 1.0630361e-03]],
    ...
    (.., .., "Bolivia"): ...,
    ...}
for (magnitude, distance, name), (frequencies, power) in earthquakes.items():
    ...
    if name:
        ...
@megies

This comment has been minimized.

Copy link
Member

commented Sep 24, 2019

@FMassin some more comments for improvements but this is close to being ready. I know we're kind of nit picky sometimes but it will help to maintain the code in the future.

@FMassin

This comment has been minimized.

Copy link
Author

commented Sep 25, 2019

You guys are providing a tool than spares tons of efforts to a lot of people for free. You are not picky by asking people to some respect some simple rules. It would be pure craziness not to.

(.., .., "Bolivia"): ...,

Really, it's OK, removing this one is better than adding a blank key element I think.

text = ax.tests[index]
assert all(line.get_ydata() == power)
assert all(line.get_xdata() == periods)
assert text.get_text() == native_str('$^{M3.5}_{100km}$')

This comment has been minimized.

Copy link
@megies

megies Sep 27, 2019

Member

Sorry if my example snippet was misleading there, please use the TestCases methods, like..

self.assertEqual(line.get_ydata(), power)
self.assertEqual(line.get_xdata(), periods)
self.assertEqual(text.get_text(), native_str('$^{M3.5}_{100km}$'))

The string at the end has to be modified as well i think '$^{M%.1f}_{10km}$' % magnitude

This comment has been minimized.

Copy link
@FMassin

FMassin Sep 27, 2019

Author

Done.

@megies

This comment has been minimized.

Copy link
Member

commented Sep 27, 2019

@FMassin test looks good for the most part :-)
There's a typo in it, in your above test run, your newly added test wasn't executed actually ("collected 14 items / 14 deselected" -- no test casess actually run), no idea why. You can see an overview of test report from the CI runners here (in addition to goin to the CI pages and scrolling down in the output to see the stdout of obspy-runtests): http://tests.obspy.org/?pr=2455 and then e.g. http://tests.obspy.org/105541/#1

@megies

This comment has been minimized.

Copy link
Member

commented Sep 27, 2019

@aringler-usgs this is mostly ready from the coding point of view, you seem to be familiar with the papers used here, would you be willing to give this a last check from the science point of view? Thanks!

@megies megies requested a review from aringler-usgs Sep 27, 2019
@@ -550,16 +550,17 @@ def test_earthquake_models(self):
fig = ppsd.plot(show_earthquakes=(1, 4, 5, 15), show=False)
ax = fig.axes[0]
for index, magnitude in enumerate([1.5, 2.5, 3.5]):
frequencies, power = earthquake_models[(magnitude, 10)]
frequencies, power = earthquake_models[(magnitude, 10))

This comment has been minimized.

Copy link
@megies

megies Sep 27, 2019

Member

needs reverting

This comment has been minimized.

Copy link
@FMassin

FMassin Sep 27, 2019

Author

Actually I managed to get pytest to work and there are more issues as there are other lines ... I am fixing it.

This comment has been minimized.

Copy link
@FMassin

FMassin Sep 27, 2019

Author

Changed approach but done? See below

@FMassin

This comment has been minimized.

Copy link
Author

commented Sep 27, 2019

I had to add a label, scan all lines in plot and to simplify array comparisons... now works on my laptop:

 ❯ ~/anaconda3/bin/pytest /Users/fmassin/anaconda3/lib/python3.5/site-packages/obspy/signal/tests/test_spectral_estimation.py -k test_earthquake_models
================================================================================ test session starts ================================================================================
platform darwin -- Python 3.5.5, pytest-3.8.1, py-1.8.0, pluggy-0.11.0
rootdir: /Users/fmassin, inifile:
plugins: remotedata-0.3.1, openfiles-0.4.0, doctestplus-0.3.0, arraydiff-0.3
collected 15 items / 14 deselected

../../../../anaconda3/lib/python3.5/site-packages/obspy/signal/tests/test_spectral_estimation.py .                                                                            [100%]

====================================================================== 1 passed, 14 deselected in 9.19 seconds ======================================================================
@FMassin

This comment has been minimized.

Copy link
Author

commented Sep 27, 2019

@aringler-usgs this is mostly ready from the coding point of view, you seem to be familiar with the papers used here, would you be willing to give this a last check from the science point of view? Thanks!

I will also have the actually authors checking plots and data included in code.

@@ -1958,7 +1958,7 @@ def plot(self, filename=None, show_coverage=True, show_histogram=True,
xdata = periods
if xaxis_frequency:
xdata = frequencies
ax.plot(xdata, ydata, '0.4', linewidth=2)
ax.plot(xdata, ydata, '0.4', linewidth=2, label=key)

This comment has been minimized.

Copy link
@megies

megies Sep 27, 2019

Member

I don't really like this. You're adding something real to the plot that just stays hidden because we (at the moment) never call .legend() on the plot, and that's only because it's used during testing. I don't really understand why you need this, but there has to be another way to do it..

This comment has been minimized.

Copy link
@FMassin

FMassin Sep 27, 2019

Author

OK I also admit, it is kind of cheating...

Copy link
Member

left a comment

This is a nice and useful contribution. There is only one issue I see. The units need to be dealt with.

Clinton's earthquake spectra are in units of (m/s/s)^2. A psd is in units of (m/s/s)^2/Hz. So you need to do an integration to remove the density. The second figure is from John Clinton's AGU poster (2010). Sorry, I can't seem to find a copy of Clinton and Cauzzi (2013).

Two ideas (neither of which are great):

Add a method for doing integrations. This would be very useful (I find myself often needing such a thing).

When including the earthquake spectra you could do the integration.

The integration should look like the following:

# Convert from dB back to ground acceleration
power = 10**(power/10.)
# Now in (m/s/s)^2/Hz
power = power*freq*(2**(1/2)-2**(-1/2))
# Now in one-octave integrated bands

If you wanted to do half-octave you would use 2**(1/4)-2**(1/4)

Let me know if you need more clarification.

Figure_1

Screen Shot 2019-09-27 at 9 33 48 AM

@FMassin

This comment has been minimized.

Copy link
Author

commented Sep 27, 2019

@aringler-usgs: I though the purpose of eq1 in Cauzzi & Clinton 2013 was to address this issue, that is why I followed them and used:

ydata = acceleration / (periods ** (-.5))
ydata = 20 * np.log10(ydata / 2)

I will revise with them on Monday.

@aringler-usgs

This comment has been minimized.

Copy link
Member

commented Sep 27, 2019

@FMassin is there any chance you could send me a copy of the paper? I can't seem to find my copy of it

I will double-check your comparisons.

@FMassin

This comment has been minimized.

Copy link
Author

commented Sep 27, 2019

@aringler-usgs : done (on your usgs email).

@aringler-usgs

This comment has been minimized.

Copy link
Member

commented Sep 27, 2019

@FMassin Thanks for the paper. Your earthquake spectra are in units of m/s/s. However, the PSDs for BW KW1 (the four yellow curves) are in units of (m/s/s)^2/Hz. So you want to convert the BW KW1 curves to be in m/s/s. Which could be done using equation (1) in Clinton and Cauzzi (2013).

You would also want the y-axis label to have the /Hz removed.

I apologize if I am not explaining this clearly. Do let me know if you disagree with what I am saying.

periods = 1.0 / frequencies
# Eq.1 Cauzzi & Clinton (2013)
ydata = accelerations / (periods ** (-.5))
ydata = 20 * np.log10(ydata / 2)

This comment has been minimized.

Copy link
@FMassin

FMassin Oct 4, 2019

Author

This conversion from m/s/s to dB has been checked by @ccauzzi and @jclinton.

This comment has been minimized.

Copy link
@aringler-usgs

aringler-usgs Oct 4, 2019

Member

Okay, yes. You are converting spectra to densities. Usually, folks go the other way, but this works. I think that was my confusion.

Would you be willing to put a comment or a note? I apologize for making this such an issue, but it has been done incorrectly in various papers, which has produced confusion.

These curves will be very useful!

This comment has been minimized.

Copy link
@FMassin

FMassin Oct 6, 2019

Author

@aringler-usgs : thanks a lot, but could you suggest the kind of comment that should be added? And where exactly?

This comment has been minimized.

Copy link
@aringler-usgs

aringler-usgs Oct 7, 2019

Member

Maybe just change your current comment to the following?

# Eq.1 from Clinton and Cauzzi (2013) converts power to density

That would be sufficient for me. This helps someone know there is a unit change here.

This comment has been minimized.

Copy link
@FMassin

FMassin Oct 7, 2019

Author

OK thanks a lot.

@FMassin

This comment has been minimized.

Copy link
Author

commented Oct 4, 2019

I had this reviewed by @jclinton and @ccauzzi, and they confirm that the plot looks correct, and that the usage of equation (1) in Clinton and Cauzzi (2013) is actually correctly implemented (https://github.com/obspy/obspy/pull/2455/files#r331446758).

I would like to improve the doc and the distance selection and then, if no further corrections are required, I would be happy to finalise this.

More images:
image

@FMassin FMassin force-pushed the FMassin:patch-1 branch 2 times, most recently from cec7bf8 to 047a635 Oct 4, 2019
FMassin and others added 3 commits Sep 13, 2019
This edit adds earthquake spectral acceleration models from Clinton & Heaton (2002) converted into decibels within PPSD plot following Cauzzi & Clinton (2013) equation 1.

Clinton, J. F., & Heaton, T. H. (2002). Potential Advantages of a Strong-motion Velocity Meter over a Strong-motion Accelerometer. Seismological Research Letters, 73(3), 332–342. http://doi.org/10.1785/gssrl.73.3.332
Cauzzi, C., & Clinton, J. (2013). A High- and Low-Noise Model for High-Quality Strong-Motion Accelerometer Stations. dx.doi.org (Vol. 29, pp. 85–102).  Earthquake Engineering Research Institute. http://doi.org/10.1193/1.4000108
@megies megies force-pushed the FMassin:patch-1 branch from c185823 to a87d72a Oct 9, 2019
@megies
megies approved these changes Oct 9, 2019
Copy link
Member

left a comment

Thanks for the contribution @FMassin and thanks for the review and comment @aringler-usgs. Will merge this as soon as CI looks good.

@megies megies moved this from In Progress to Waiting on CI in Release 1.2.0 Oct 9, 2019
@FMassin

This comment has been minimized.

Copy link
Author

commented Oct 9, 2019

Thanks you @megies and @aringler-usgs !

@megies

This comment has been minimized.

Copy link
Member

commented Oct 9, 2019

Looks like there can be weird issues with the labeling @FMassin , see http://tests.obspy.org/105891/#1

I really only just noticed this strange super-/subscript combo. Any reason to not just make it 'M%.1f\n%dkm' instead of '$^{M%.1f}_{%dkm}$'?

megies and others added 2 commits Oct 9, 2019
@FMassin

This comment has been minimized.

Copy link
Author

commented Oct 11, 2019

@megies : fine, with smaller font (text spills over the next line otherwise). Fixed.

@FMassin

This comment has been minimized.

Copy link
Author

commented Oct 17, 2019

@megies is there anything else that I should do?

@megies

This comment has been minimized.

Copy link
Member

commented Oct 18, 2019

@FMassin think you need to adjust tests, see http://tests.obspy.org/?pr=2455 e.g. http://tests.obspy.org/105972/#1

Those numerous Darwin fails are probably related to that Mac/Anaconda issue mentioned on gitter, no clue, just ignore them

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
Release 1.2.0
Waiting on CI
5 participants
You can’t perform that action at this time.