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

Reworking Array Analysis for ObsPy #1002

Open
wants to merge 127 commits into
base: master
from

Conversation

@krischer
Copy link
Member

commented Mar 27, 2015

This is work in progress and will likely take a couple of months.

It attempts to unify and simplify the array processing in ObsPy and also adds a couple more methods.

The idea is to base it on an array class, e.g.

from obspy.signal.array_analysis import SeismicArray

arr = SeismicArray("test")
arr.add_inventory(...)

arr.beamforming(stream, ...)
@mbyt

This comment has been minimized.

Copy link
Member

commented Mar 27, 2015

Tons of changes here. Can we maintain such a huge functionality in the future?

@@ -272,6 +272,9 @@ def linspace2(val1, val2, N):
scatter = bmap.scatter(x, y, marker=marker, s=size, c=color,
zorder=10, cmap=colormap)

# Hack!!

This comment has been minimized.

Copy link
@QuLogic

QuLogic Mar 27, 2015

Member

aka, why aren't we using cartopy?

This comment has been minimized.

Copy link
@megies

megies May 22, 2015

Member

How can this work anyway, I don't see that ax.basemap is available normally, i.e. Axes are not aware of the Basemap instance, only the other way around. And I don't see where it gets attached manually.. or am I missing something..?

This comment has been minimized.

Copy link
@QuLogic

QuLogic May 23, 2015

Member

Uh, it's set right here on the line below...

This comment has been minimized.

Copy link
@megies

megies May 24, 2015

Member

euhm. indeed. nevermind..

from matplotlib import cm
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from obspy.core import UTCDateTime

This comment has been minimized.

Copy link
@QuLogic

QuLogic Mar 27, 2015

Member

Sorted incorrectly.

from scipy.integrate import cumtrapz

from obspy.core import Stream
from obspy.signal.headers import clibsignal
from obspy.signal.invsim import cosTaper
from obspy.signal.util import nextpow2, utlGeoKm
from obspy.core.util import AttribDict
import warnings

This comment has been minimized.

Copy link
@QuLogic

QuLogic Mar 27, 2015

Member

Also in the wrong place.

@QuLogic QuLogic added the .signal label Mar 27, 2015
"""
distances = []
geo = self.geometry
for station, coordinates in list(geo.items()):

This comment has been minimized.

Copy link
@QuLogic

QuLogic Apr 17, 2015

Member

It's not necessary to add list when you're just going to iterate through the items.

@Rothenhouser

This comment has been minimized.

Copy link
Contributor

commented Jul 8, 2015

So as far as I am concerned, this code is pretty usable by now and could probably be merged.
Check out the new version of the beamforming tutorial to get an idea of the changes.

@@ -9,7 +9,7 @@
:toctree: autogen
:nosignatures:

~array_analysis.array_processing
~array_analysis._covariance_array_processing

This comment has been minimized.

Copy link
@Rothenhouser

Rothenhouser Jul 8, 2015

Contributor

This probably doesn't need to be here any more, but I don't really know how the docs are built.

@QuLogic

This comment has been minimized.

Copy link
Member

commented Jul 8, 2015

This PR is non-mergeable due to conflicts. You will have to either rebase or attempt a merge. Then we'll see whether the tests pass, as they haven't been running for a long time.

@Rothenhouser

This comment has been minimized.

Copy link
Contributor

commented Jul 8, 2015

I meant: Have a look, leave a comment and then we can start polishing/merging. The basic functionality, insofar as I have been using it, is done but it definitely requires more testing and tests.


# Station names must match those set in the trace headers.
stn0 = Station("BW01", Latitude(48.108589), Longitude(11.582967),
elevation=.45,

This comment has been minimized.

Copy link
@QuLogic

QuLogic Jul 8, 2015

Member

Elevation is in metres. Is it really only 45 cm?

@@ -103,7 +103,7 @@ def test_array_dilation(self):
ts2[t, stat] = array_coords[stat, 1] * dilation[t]
ts3[t, stat] = array_coords[stat, 2] * dilation[t]

out = array_rotation_strain(subarray, ts1, ts2, ts3, Vp, Vs,
out = SeismicArray.array_rotation_strain(subarray, ts1, ts2, ts3, Vp, Vs,
array_coords, sigmau)

This comment has been minimized.

Copy link
@QuLogic

QuLogic Jul 8, 2015

Member

Needs re-alignment.

@@ -147,7 +147,7 @@ def test_array_horizontal_shear(self):
ts1[t, stat] = array_coords[stat, 1] * shear_strainh[t]
ts2[t, stat] = array_coords[stat, 0] * shear_strainh[t]

out = array_rotation_strain(subarray, ts1, ts2, ts3, Vp, Vs,
out = SeismicArray.array_rotation_strain(subarray, ts1, ts2, ts3, Vp, Vs,
array_coords, sigmau)

This comment has been minimized.

Copy link
@QuLogic

QuLogic Jul 8, 2015

Member

Needs re-alignment.

@@ -0,0 +1,79 @@
from unittest import TestCase

This comment has been minimized.

Copy link
@QuLogic

QuLogic Jul 8, 2015

Member

Missing __future__ imports.

from obspy.signal.array_analysis import SeismicArray, correct_with_3dplane
from obspy.core.inventory.station import Station
from obspy.core.inventory.network import Network
import numpy as np

This comment has been minimized.

Copy link
@QuLogic

QuLogic Jul 8, 2015

Member

NumPy should be imported before ObsPy.

self.testarray = SeismicArray()
self.testarray.add_inventory(self.simpleinv)

self.geo_exp = {'5pt.bl': {'absolute_height_in_km': 0.0,

This comment has been minimized.

Copy link
@QuLogic

QuLogic Jul 8, 2015

Member

This is the same data as self.simpleinv, right? I think it should be generated from it instead of duplicated here.

def test__get_geometry(self):
geo = self.testarray.geometry
self.assertEqual(geo, self.geo_exp)
# test for both inventories (or fake inventories) with and w/o channels

This comment has been minimized.

Copy link
@QuLogic

QuLogic Jul 8, 2015

Member

Is something supposed to be here following this comment?

@@ -4,6 +4,12 @@
unicode_literals)
from future.builtins import * # NOQA

from obspy import Trace, Stream, UTCDateTime

This comment has been minimized.

Copy link
@QuLogic

QuLogic Jul 8, 2015

Member

These imports almost all exist a little further down (and should be there, not here.)

@QuLogic

This comment has been minimized.

Copy link
Member

commented Jul 8, 2015

Unfortunately, the changes to obspy/signal/array_analysis.py are so large that GitHub will not show the diff. That will make it difficult to make comments on it, especially since that's the meat of the changes. I don't suppose there's any way to split these changes into smaller PRs?

@krischer krischer added this to the 0.11.0 milestone Aug 3, 2015
@QuLogic

This comment has been minimized.

Copy link
Member

commented Aug 10, 2015

@krischer So now that this is on 0.11.0, is it ready? It still needs some merge conflict fixing...

@krischer krischer modified the milestones: Future release, 0.11.0 Aug 12, 2015
@krischer

This comment has been minimized.

Copy link
Member Author

commented Aug 12, 2015

Yea this is really not ready yet...I moved it to a future milestone.

@Rothenhouser

This comment has been minimized.

Copy link
Contributor

commented Aug 26, 2015

Summary of changes so far compared to array analysis currently in obspy:

  • Use of a SeismicArray class to which an Inventory object, containing location information etc for all stations in the array, is attached. The array geometry is calculated from this inventory and (internally) passed to the array analysis methods (methods of the SeismicArray object).
  • The SeismicArray class has all the currently implemented methods for array processing, plus some new ones. Methods are available for Vespagrams, slowness whitened power, phase weighted stack, delay and sum, f-k analysis, three-component beamforming and array-derived rotation. The code of the methods present previously, e.g. the f-k analysis, is almost unchanged but it is wrapped in a 'front-end method' for the user which, most importantly, calculates the array geometry and passes it to the original method as an argument.
  • Using the methods requires data streams, which must of course be taken from stations/channels represented in the attached inventory. They are matched using the SEED ID in the trace headers. The location information from the inventory can then automatically be attached to the correct trace if the analysis method requires it. A method is provided to reduce the inventory to only contain stations or channels from which data were recorded.
  • The results from the beamforming are no longer just simple numpy arrays; these are instead returned as attributes of a BeamformerResult object which also contains some metadata such as what method created the results object, for what times the analysis was performed, and some of the processing parameters that were chosen. The Results object also contains some plotting methods to make outputs from different methods comparable. They can also be added (a bit rudimentary so far) to allow processing of e.g. one day at a time, writing to file, then combining for long-term plots later.
  • To get a quick overview of how the workflow has changes, have a look at the reworked tutorial for f-k analysis in the pull request. The same results are achieved more easily.

Caveats/Issues:

  • I used, and focussed on, the three-component beamforming and the f-k analysis and am fairly confident they are correct (this was my Master's thesis). The other methods should be correct but haven't been used or changed much by me.
  • Unit tests to prove the above are still missing.
  • Only one tutorial so far.
  • No citations for the methods used.
@megies

This comment has been minimized.

Copy link
Member

commented Apr 13, 2016

@Rothenhouser, note that we build docs for all PRs automatically now (see commit check with link next to Travis, Appveyor etc.), so it's easy to check if new pages render OK. (+DOCS)

See e.g. http://docs.obspy.org/pull-requests/1002/packages/autogen/obspy.signal.array_analysis.BeamformerResult.html, which has some broken links to functions.

@krischer krischer force-pushed the krischer:arraytools branch from 03bc689 to f916715 Jul 26, 2016
@krischer

This comment has been minimized.

Copy link
Member Author

commented Jul 26, 2016

Rebased and force pushed. This was nasty and I hope everything worked as expected - I still have a backup of the other branch in case I screwed up.

+DOCS

krischer and others added 19 commits Jan 13, 2017
vespagram routine
@ThomasLecocq

This comment has been minimized.

Copy link
Contributor

commented Jan 13, 2017

someone is on fire 😄 🔥 !!

@megies

This comment has been minimized.

Copy link
Member

commented Jan 16, 2017

Just for future reference, here's how to get the old exercise notebook (data is here) running (with a separate conda environment):

$ conda create -n mess2014_tmp python=2.7
$ source activate mess2014
$ conda install -c obspy obspy  # to install all dependencies with conda..
$ conda remove obspy
$ pip install https://github.com/jwassermann/obspy/archive/arraytools.zip
$ pip install https://github.com/obspy/mess2014/archive/master.zip
@megies

This comment has been minimized.

Copy link
Member

commented Jan 16, 2018

Hooray.. exactly one year since the last post in here.. 🎉 and we're sitting together discussing it again.. 🙈

@megies

This comment has been minimized.

Copy link
Member

commented Jan 17, 2018

@jwassermann @krischer I've put together a size-reduced test data set that is roughly 600k with one event for yellowknife array, both real and synthetic data, downsampled but still good, I think. This should be fine for both testing and tutorial purposes.

Of course, size could still be reduced further for testing, but we have much larger test files in some older test cases, so I think it's OK.

gist with jupyter notebook to create the data set from the MESS2014 data:
https://gist.github.com/58593bbde0e5773ae8426e5117dcf09f
(at the very bottom is a comparison of the raw + processed/downsampled + synthetic data for the first phase)

@krischer

This comment has been minimized.

Copy link
Member Author

commented Jan 17, 2018

Yea that is definitely small enough - just one note of caution: The synthetics from syngine will have no elevation difference as all receivers are always placed at the surface of the sphere.

@megies

This comment has been minimized.

Copy link
Member

commented Jan 17, 2018

The synthetics from syngine will have no elevation difference as all receivers are always placed at the surface of the sphere.

Yeah. That's why I chose Yellowknife. For subroutines that use 3D elevation correction we'll probably some different testing (maybe some simple made up synthetic test).

I thought it's maybe a good idea to throw in the axisem synthetics as they will likely give very crisp results for the array analysis results, dunno.

@megies

This comment has been minimized.

Copy link
Member

commented Aug 24, 2018

@jwassermann has worked on reintroducing routines that got scrapped during the refactoring.. see krischer#3

@trichter trichter referenced this pull request Aug 22, 2019
9 of 9 tasks complete
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
9 participants
You can’t perform that action at this time.