Plot Event far-field radiation pattern / refactor core/event.py into folder structure

[MMesch]

This addition to obspy allows to compute and plot P/S-wave farfield radiation patterns for (hopefully) arbitray moment tensors. You can run the test or try:

from obspy.imaging.radpattern import plot_3drpattern
mt = [0.91, -0.89, -0.02, 1.78, -1.55, 0.47]
plot_3drpattern(mt,kind='p_sphere')

it is neither well tested nor optimized. Some improvements could be:

• add colour and other things to quiver plots
• making sure that axes are well aligned
• plotting beachball projections on the sides of the matplotlib 3d plot
• compute radiation pattern for different a set of taup phases and stations
• add nearfield radiation
• add single force

[megies]

Looks nice!

compute radiation pattern for different a set of taup phases and stations

That could be very nice for many people, I think..

[claudiodsf]

That's very nice!
It only seems strange to me to have farfield_p and farfield_s functions in obspy.imaging.
Maybe it's time for an obspy.source package?

[MMesch]

Agree. Things that could/should be improved are then:

• better function names and organisation. Inside the imaging package? Inside a class?
• add colour and other things to quiver plots. This requires non-standard mlab3d. Adapting the vector length to the radiation pattern would be nice as well. So far it can only take a direction.
• making sure that axes are well aligned with other moment tensor conventions
• plotting beachball projections on the sides of the matplotlib 3d plot. Should be possible but bugs (see commented lines in the p_sphere section of the code). Keypoint is to use sth. like art3d.patch_collection_2d_to_3d(bball). Does anyone know how this works?
• compute radiation pattern for a set of taup phases and stations. I guess this could be combined with a method to plot all of it on a map.
• add nearfield radiation pattern. Possible but more complicated than farfield.
• add single force. Should be straight forward for the farfield.
• tetrahedral spherical grid to reduce plotting time. Does anyone know how to get this in mlab3d?

[krischer]

I'm playing around with the quiver plots - not perfect yet but its getting better :-)

[megies]

It only seems strange to me to have farfield_p and farfield_s functions in obspy.imaging.

I agree.

Maybe it's time for an obspy.source package?

But the place to go is already there: obspy/core/event.py

EDIT: There could even be methods attached to FocalMechanism and MomentTensor making use of these routines -- and to encourage more use of these classes.. (we should migrate a lot of beach ball code to use those objects anyway, compare e.g. #1038)

[krischer]

Hmm...I did not realize that matplotlibs 3D plotting lacks basic features like figuring out what should be drawn in the front if one plots two objects...I don't think its possible to make a good 3D quiver plot with matplotlib. Even the "nicer" version of me is almost impossible to interpret because the orthographic projection messes with the mind and is just very confusing.

@MMesch I would be fine to have your original mayavi plots in there as well. Just make sure to import mayavi within the plotting functions so it does not become a hard ObsPy dependency.

Other opinions on this?

[megies]

Just make sure to import mayavi within the plotting functions so it does not become a hard ObsPy dependency.

Other opinions on this?

Yeah, why not. If it's optional visualizations of more special stuff and done as soft dependency.. 👍

[MMesch]

OK. If someone finds out how to plot the beachball patch_collection in the 3D plot this would be great. I didn't get very far but it should be possible (see http://stackoverflow.com/questions/18228966/how-can-matplotlib-2d-patches-be-transformed-to-3d-with-arbitrary-normals )

[claudiodsf]

@megies
But obspy/core/event.py contains: Classes for handling event metadata

Here we have a possible new class of functions for seismic source representation, which seems something different to me.

[MMesch]

I managed to plot the beachball together with the surface plot but matplotlib is just not capable of really doing 3D. It can't tell which patch is above the other for the beachball (therefore facecolor='none') and it cannot say when the beachball is behind the surface. 3D plotting in matplotlib is too limited and I'm giving up on this...

[QuLogic]

It doesn't look like any released version of Mayavi2 supports Python 3 😦 .

FYI @carltape.

[carltape]

Cool stuff in this thread! I've got some cool stuff in the link below, but most of it is mathematica (which very few seismologists use). I agree, it would be great to have some open-source alternatives.
http://www.giseis.alaska.edu/input/carl/research/beachball_gallery.html

[celsoa]

Here's a nice spherical harmonics example with Mayavi, this may be what @QuLogic was referring to
http://docs.enthought.com/mayavi/mayavi/auto/example_spherical_harmonics.html#example-spherical-harmonics

[megies]

@megies
But obspy/core/event.py contains: Classes for handling event metadata

Here we have a possible new class of functions for seismic source representation, which seems something different to me.

It's just that we're having a lot of pain currently (and I mean really a lot), cleaning up and sanitizing the module structure (which has exploded over the last years), esp. to clean up the top levels of our submodule hierarchy (see #999, #1039, ...).
So we should be very careful not to open Pandora's box again. Any proposals for new top-level submodules immediately ring three or four alarm bells in my head. In my opinion, something concerned with seismic sources belongs somewhere close to our event related stuff. The Classes for handling event metadata in event.py is just some info stub somebody put there when we started that file (with QuakeML support in mind) and can be adapted. I agree that event.py is a bit stuffed, but we can easily make it into a subfolder and separate things in there logically, that wouldn't break user code.

[MMesch]

Here are plots with mayavi that I did before. Another alternative would be to just output a vtk file if there is no mayavi/vtk and plot otherwise.

http://pythology.blogspot.fr/2015/11/earthquake-moment-tensor-radiation.html

[MMesch]

[megies]

Looks great! 👍

[ThomasLecocq]

👍 amazing look in 3D vtk (mayavi or paraview, whatever) !!

[MMesch]

@megies @claudiodsf I open an event folder now in obspy.core that combines the event.py file and radpattern.py . Is this what you were thinking? I think the radiation pattern functions should somehow be grouped. I am not sure whether a submodule or a class is best. Any thoughts?

[megies]

@MMesch, yes, that's what I think makes most sense, to group these event/origin/source related functions close to where all the event stuff is defined. Since up to now it was a plain py file, making it a folder and having everything in there feels natural to me.
But I'm open for any reasoning that suggests something else, of course.

What do you think about renaming the old event.py which is now event/event.py to event/base.py?

I think the radiation pattern functions should somehow be grouped. I am not sure whether a submodule or a class is best.

I think having them grouped in that submodule event/radpattern.py is fine. The pure plotting functions should in principle be put back into imaging..? Also, I'd propose adding methods to FocalMechanism/MomentTensor and Event to plot/calculate the radiation pattern -- as a high-level, preferred usage option for users. Most users should not need to lookup those functions and use those methods instead. I think a map plot would also be nice, but that probably would not work well for very short distances with taup?

[MMesch]

There is direct vtk file output now. Renamed event.py -> base.py . Still need to move the plotting routines, and add a taup plotting routine. Normalization of the vector field needs to be checked as well.

[MMesch]

vtk output for paraview is implemented. I am really having trouble installing mayavi with anaconda. It wants to downgrade numpy which leads to some errors in matplotlib afterwards. Did anyone succesfully install the newest numpy matplotlib and mayavi with anaconda?

[MMesch]

This is the new plot. I hope the orientation is correct. It seems to be OK because as far as I know a beachball is the projection of the lower hemisphere on a plane seen from above (which inverts left and right if seen from below in the 3d view). Please tell me if something looks weird :)

[MMesch]

I'm done for the moment. We could add some kind of map in this plot if latlon keyword is provided. Connecting it to the event object is then very quick and simple.

I thought about the maps: it seems that a function to plot 2d greatcircle paths between stations and one or several sources is still missing in obspy. I think we should make it a separate branch/issue to add a function that gives 2d greatcircle paths and 3d paths using taup. Once we have such a function, it is simple to combine it with the event radiation pattern and also plotting routines in 2D and 3D.

I therefore suggest to finish/test/include the radiation pattern functions as they are right now and to look after the taup map in the next step.

Matthias

[MMesch]

Mesh orientation switches now automatically between null and other axes

[MMesch]

And yet more options. Pressure stretches and tension squeezes the sphere. I personally find this more intuitive.

[MMesch]

I have pulled the Event and Catalog classes out of base.py and attached a plot function to the Event class. Not sure about the best way to better organize this huge class hierarchy @megies @krischer but base.py is really too big and the class hierarchy too deep to be easily accessible in my opinion. If it has to be as deep as it is now (I guess adapted to quakeml) then all classes that are unlikely to be changed should probably be hidden in one file and the others separate from them in different files.

You plot the event data now with:

from obspy import read_events
cat = read_events('.../obspy/obspy/io/quakeml/tests/data/quakeml_1.2_focalmechanism.xml')
ev = cat[0]
ev.plot()

[megies]

OK, green light by Travis, finally..

Since nobody objected to the core/event/ structure in here I'm gonna merge now.

[krischer]

🎉

[claudiodsf]

👍

[MMesch]

yey!

[QuLogic]

@MMesch I think one comment regarding the usage of RTP vs. USE was not addressed; can you please follow up with a correction in another PR?

[QuLogic]

I'm re-adding this comment so that @MMesch can address it:

Isn't RTP the same as USE? Why is this name used here, but not anywhere else (e.g., in MoPaD)?

[MMesch]

yes good point. Technically they are the same, but I think I added USE for regional models (cartesian definition), just to be complete. If you prefer we can remove it!

[QuLogic]

It'd be preferable to use consistent naming between this and MoPaD, which doesn't use RTP at all. I'm fine with adding aliases; I can never remember which one is which anyway.

The above statement is incorrect though, since USE is implemented as it's the same as RTP?

[MMesch]

How can I open a good pull request for this that includes the stuff that @megies merged? I can go through the doc once more then...

[QuLogic]

This branch is already merged to master; create a new branch starting from the latest master. There's no need to use the existing radpattern branch when creating a new PR.

[MMesch]

OK!

[krischer]

I just discovered this here which also has some very fancy plots and is based on this very PR:

https://github.com/FMassin/SeismicSource/blob/master/README.ipynb

Any chance for also getting those into main ObsPy?

@FMassin, @MMesch

[MMesch]

looks great! I think adapting the obspy plot to look similar would be very good. Essentially this means different colors and a smaller surface representation of the radiation pattern? Then all of them should be put into one plot. Also we can benchmark the radiation pattern values with this independent code!

[krischer]

Also maybe center the colormap around 0?