ray path plotting routine #1236

MMesch · 2016-02-05T23:57:42Z

This is probably another longer duration pull request. It adds functions for greatcircle path plotting in 2D and 3D to obspy, based on geographiclib and taup. It is already working but needs to be structured nicely. Here is an example, plotted again with mayavi (green = P, blue = Pdiff, red = PKP), but plotting could be done with matplotlib and vtk file output again. Right now it works for one event and and inventory, but I think now it is easy to expand it to a catalog and an inventory together.

+DOCS

QuLogic · 2016-02-06T03:32:51Z

Perhaps you have missed it (likely because of our slow release rate), but there is already a plotting routine in TauPyModel (though not in 3D) as well as a function to get ray paths via geographic coordinates.

Expanding the existing stuff to 3D and to handle an entire Catalog/Inventory would certainly be interesting.

MMesch · 2016-02-06T11:09:23Z

ok I use the taup function now and put the subroutines directly into the inventory class. I'll add some non-mayavi plotting options later...

krischer · 2016-02-06T11:44:10Z

Looks great :-)

In terms of code structure I think it would be much better to but all plotting code into obspy.imaging and just add thin wrappers to the inventory class. If we start to add fancy plotting to all core objects it will really bloat the core module. The same might be good for the moment tensor plotting things in the other PR.

class Inventory(...):
    ...

    def plot_rays(self, *args, **kwargs):
        from obspy.imaging.inventory import plot_rays
        plot_rays(inventory=self, *args, **kwargs)

The plot_rays() function should probably also take event objects. The optional event coordinates are just there for convenience.

plot_rays(inventory, event=None, event_latitude=None,
          event_longitude=None, event_depth_in_m=None):
    ...

MMesch · 2016-02-06T12:06:34Z

OK agree. I just don't see that these things are really 'imaging' in the sense of tomography. They are more visualization routines. That's why it seems strange to put them in this module.

krischer · 2016-02-06T12:11:27Z

Haha...i did not realize that confusion. obspy.imaging actually is the visualization module of ObsPy ;-) I'm so used to the name I never saw that the name is ambiguous in seismology...

MMesch · 2016-02-06T12:15:35Z

ok I'll try to see it as visualization then :)

MMesch · 2016-02-06T17:59:18Z

looks pretty cool now if I may say so :)

P

PKiKP

PPP

Pdiff

PcP

PKiKP bright colorscheme

megies · 2016-02-07T12:50:50Z

Bright color scheme is much better to see for me, at least in the examples above..

QuLogic · 2016-02-08T00:24:58Z

It looks like you're creating your own colour scheme by hand. That should not be necessary. Since we have a hard dependency on matplotlib, you can just pull the colour map from there. Even better would be to use the existing colour set.

Also, it seems you have not configured your git client with your email address, is there a reason you have not done so? The first commit here is attributed to <matthias@matthias-Latitude-E6520.(none)> and the rest to <MMesch@users.noreply.github.com>. It is okay if you do not wish to expose your email address, but it would be best if a consistent address were used.

MMesch · 2016-02-08T16:52:18Z

ok. I use the taup colors now but I have to change the lightness and saturation because it is essential to see something in the 3d plots.

sorry about my email address, I'll try not to change it anymore!

krischer · 2016-02-08T23:20:41Z

I pushed some small things to make it a bit easier to use. Looks really cool :-)

A couple notes:

Diffracted phases need to be interpolated as they only store two points along the whole diffracted part of the phase which results in a straight line in xyz which is wrong. This probably has to go into the geographic ray calculation routine.
The colors don't seem to work for me - I'll look into that.
Mayavi as a whole seems very unstable - it keeps crashing randomly for me.

krischer · 2016-02-09T10:06:38Z

Considering that Mayavi is very unstable for me I played around a bit and maybe we can do this completely in WebGL. The nice thing is that the only dependency is a web browser which everyone has and it would even work in the notebooks. Also the render quality seems to be higher (I cannot get mayavi to render with antialiasing) and it is probably a bit easier to add interactivity.

Super unpolished toy example (12 MB of data...so you might have to wait a second): http://www.geophysik.uni-muenchen.de/~krischer/raypaths/

claudiodsf · 2016-02-09T10:13:31Z

@krischer Wow! That is impressive! (except for the earthquake at north pole :). Which tool do you use to generate WebGL?

krischer · 2016-02-09T10:16:08Z

Yea...the coordinate system for the rays (or the textures) has to be rotated so its not fully correct as of now ;-)

This uses three.js: http://threejs.org/

MMesch · 2016-02-09T10:42:29Z

https://youtu.be/_4YiHAO70KA

MMesch · 2016-02-09T10:44:03Z

@krischer great thing!!! This will certainly attract some people to the lmu server :)

MMesch · 2016-02-09T10:44:15Z

also hopefully to obspy

MMesch · 2016-02-09T10:45:46Z

my mayavi is absolutely stable btw. Maybe it is an issue of your graphic cards driver?

ThomasLecocq · 2016-02-09T10:46:40Z

@MMesch WOOOOW amazing!

MMesch · 2016-02-09T10:47:10Z

Another though (sorry about all the mini posts):

The Pdiff issue should be solved in the taup module directly in my opinion

krischer · 2016-02-09T10:55:10Z

Great movie :-)

@krischer great thing!!! This will certainly attract some people to the lmu server :)

I'll give this a real shot then. I'll try to mimik your plot and the website will just be output of inv.plot_rays(). A really cool side effect is that it works on Python 3, has zero dependencies, and could also be embedded in the notebooks.

my mayavi is absolutely stable btw. Maybe it is an issue of your graphic cards driver?

Hmm...who knows :-( I have a macbook with the same graphic card as ten million other people and it usually works fine enough.

The Pdiff issue should be solved in the taup module directly in my opinion

I'm not sure - it is not wrong in the spherical coordinate systems. The plotting just forces a cartesian interpolation which is wrong. But it has to happen before the geographical routines take over, otherwise we have to deal with circular longitudes and wrap around effects...

MMesch · 2016-02-09T12:33:05Z

The right thing would be to do a spherical interpolation in taup_geo before it is passed to geographiclib then. EDIT: should we open an issue for this?

Also the complete mayavi scene can be composed as pure vtk files. It is then independent from any instabilities and can be displayed e.g. with paraview.

I'm not sure what limitations you might encounter with pure OpenGL. I guess since we use very basic geometries only (lines, spheres, cones), it shouldn't be to difficult to rewrite it. It would probably be more stable but we wouldn't have lights or the whole graphical user interface of mayavi.

webGL on the other hand really adds more fuctionality, so I think it is worth it!

MMesch · 2016-02-09T12:36:06Z

One more thing:

the taup_geo routine seems to be somewhat inaccurate for PKIKP phases. E.g. in the video that I have posted, rays that arrive at stations that are close to the antipode do not exactly converge at the station longitude latitude.

I am not sure why this is the case. Either a taup problem, or maybe the ellipticity is not treated correctly ?

megies · 2016-02-09T12:57:41Z

Hmm.. a general comment regarding the ray path plots:

To me those plots seem incredibly overcrowded, I wouldn't know where to start trying to interpret those.. I am not a tomo-guy but I guess this is mostly about knowing where and which rays hit some zone of interest (e.g. CMB). Wouldn't it make more sense to do e.g. a 2D (or 3D) ortho map-plot and plot into it the pierce points (e.g. pierce points at CMB but plotting coastlines over it), labelled by e.g. "{event_id}:{station_code}", or do a heatmap/hexbin of pierce points (to visualize sensitivity of the data set regarding some region and layer of interest)?

MMesch · 2016-02-09T13:41:24Z

Don't worry, these plots are using the maximum of phases/stations/events
possible to put into a single plot. Less paths are way better to see and
look really great (I'm preparing a video).
I think an ortho map would be really nice to have because it's easier to
work with it. It is less eye candy/fancy though. The plot right now is more
for educational purposes and visual impression. There is no reason why we
shouldn't have both.

megies · 2016-02-09T14:06:48Z

There is no reason why we shouldn't have both.

Yeah, that's absolutely right of course.

QuLogic · 2016-02-09T16:30:10Z

I'm not saying it'll be any easier, but maybe you should try VisPy too? It supports many backends including WebGL and the notebook without need for any JS (hopefully). Unfortunately, it's not quite as mature as matplotlib, so you definitely need to use a git clone.

megies · 2017-10-11T15:34:07Z

I've added a small test for the path resampling, rebased and force pushed.

I remove the vtk output, @MMesch if you offer to add a test for it I would consider keeping it.

Otherwise I'll improve the test "test_compute_ray_paths" maybe in the evening because right now it's hardly testing anything, just asserting that a list with one item is returned.

megies · 2017-10-11T16:06:44Z

@MMesch next time please take some time to implement a real test. This is really not cutting it:
https://github.com/obspy/obspy/pull/1236/files#diff-dd572eba2d259313f2468ec5ae6d42d9R56

MMesch · 2017-10-11T17:32:00Z

Haha yes. The real test used to be an image comparison. But I agree that it is better not to have any mayavi/vtk stuff on travis... Am 11.10.2017 18:06 schrieb "Tobias Megies" <notifications@github.com>:

…

@MMesch <https://github.com/mmesch> next time please take some time to implement a real test. This is really not cutting it: https://github.com/obspy/obspy/pull/1236/files#diff- dd572eba2d259313f2468ec5ae6d42d9R56 [image: screenshot from 2017-10-11 18-06-14] <https://user-images.githubusercontent.com/1842780/31452334-eba8a958-aeae-11e7-85b4-28bf77f8ac99.png> — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1236 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ACejq_wqCpTSOnf1Uhpo6P0QFwirL3OAks5srOebgaJpZM4HUuJw> .

megies · 2017-10-12T12:33:11Z

obspy/taup/ray_paths.py

+                                        radii * np.sin(thetas) * np.sin(phis),
+                                        radii * np.cos(thetas)])
+
+                value = 0.


@MMesch can you elaborate on the significance of always adding 0.0 to the output tuple?

Also it seems like bad design to include spaces for padding in the textual headers of the paths.. this should be done in the plotting part.. (coordinates in data array for path omitted below)

[(array([[...]]), u'P', 0.0, u' NET.STA', u' 2017-10-12 | M7.0')]```

I'll also include the event id to the output tuple

in the long run this routine could also use the information in the Event.arrivals to figure out what rays really exist in the data set.. but I guess in many cases such detailed info is not included in the Event data..

MMesch · 2017-10-12T14:21:25Z

Also it seems like bad design to include spaces for padding in the textual headers of the paths.. this should be done in the plotting part..

agree that's better

can you elaborate on the significance of always adding 0.0 to the output tuple

originally, I wanted to leave the routine flexible to add a value such as the amplitude from the radiation pattern to the ray path. This can be removed or organized in a better way if you have a better idea!

add event and origin id to output, remove leading spaces from labels of rays (station code and basic event info string)

megies · 2017-10-12T14:58:50Z

I hope this should do it now.. added some proper tests and sanitized output of get_ray_paths some

outdated

megies · 2017-10-12T16:01:38Z

Looks like CI is gonna give the nod.. @krischer @MMesch I think this is good to go when Travis completes.

MMesch · 2017-10-12T20:47:26Z

Great, thanks for your help @megies. The question is whether we put the mayavi/vtk moment tensor also in a separate repository - to some kind of vtk based visualisation library based on obspy.. Am 12.10.2017 18:01 schrieb "Tobias Megies" <notifications@github.com>:

…

Looks like CI is gonna give the nod.. @krischer <https://github.com/krischer> @MMesch <https://github.com/mmesch> I think this is good to go when Travis completes. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1236 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ACejq1bHVXJE8uJTzNqd9q2yAgJ-f3veks5srjfogaJpZM4HUuJw> .

megies · 2017-10-12T21:32:55Z

The question is whether we put the mayavi/vtk moment tensor also in a separate repository

Sure, go ahead.

string containing both

megies · 2017-10-13T07:10:50Z

Gonna push one additional minor change and will merge this once CI goes green then.

krischer

Minor comments. I can also do it if you are currently not working on it.

krischer · 2017-10-13T08:14:03Z

obspy/core/util/testing.py

@@ -203,10 +203,21 @@ def compare_images(expected, actual, tol):
            "shape %s." % (str(actual_image.shape),
                           str(expected_image.shape)))

+    if actual_image.shape[2] == 4:
+        has_alpha_channel = True


Do we need to keep this alpha channel handling? If it is not going to be used it might become stale.

Yeah, guess you're right, I'll remove it.

krischer · 2017-10-13T08:14:28Z

obspy/imaging/tests/data/IU_stations.txt

@@ -0,0 +1,82 @@
+#Network | Station | Latitude | Longitude | Elevation | SiteName | StartTime | EndTime 


This file can also be removed I think

currently not used in any test

Done

megies · 2017-10-13T08:48:11Z

Should be ready now.. gonna merge when CI goes green.

megies · 2017-10-13T08:57:12Z

I'll merge this locally and squash everything into one commit, to avoid polluting our master history with all the binary stuff. The history and full diff across individual commits can then still be looked up here and at https://github.com/mmesch/obspy/tree/gc_paths

(locally squashed into a single commit)

megies · 2017-10-13T09:32:35Z

This was merged locally after squashing history, see 6b43db0.

QuLogic added enhancement feature request .taup .imaging Issues with our plotting functionalities labels Feb 6, 2016

megies added 6 commits October 11, 2017 17:30

remove images again related to 3d ray path mayavi plotting

b5b7434

move ray path helpers from imaging to taup

c678d16

remove obsolete file

9259a0d

refactoring

687e215

taup: add unit test for ray path resampling

72113d3

taup: remove vtk output of paths for now

bc72ed1

megies force-pushed the gc_paths branch from 061fd80 to bc72ed1 Compare October 11, 2017 15:31

megies reviewed Oct 12, 2017

View reviewed changes

megies added 2 commits October 12, 2017 16:56

adjust get_ray_paths routine

1c209b1

add event and origin id to output, remove leading spaces from labels of rays (station code and basic event info string)

add proper test for get_ray_paths routine

b221b04

taup ray path helper: store magnitude and origin time instead of a label

220817e

string containing both

krischer previously requested changes Oct 13, 2017

View reviewed changes

megies added 2 commits October 13, 2017 10:19

remove obsolete test data file, not used anymore

9fcc3b5

image testing: remove handling of alpha channels again

9589f18

currently not used in any test

megies added a commit that referenced this pull request Oct 13, 2017

Merge pull request #1236 from MMesch/gc_paths

6b43db0

(locally squashed into a single commit)

megies closed this Oct 13, 2017

		@@ -0,0 +1,82 @@
		#Network \| Station \| Latitude \| Longitude \| Elevation \| SiteName \| StartTime \| EndTime

ray path plotting routine #1236

ray path plotting routine #1236

Conversation

MMesch commented Feb 5, 2016 • edited by megies

QuLogic commented Feb 6, 2016

MMesch commented Feb 6, 2016

krischer commented Feb 6, 2016

MMesch commented Feb 6, 2016

krischer commented Feb 6, 2016

MMesch commented Feb 6, 2016

MMesch commented Feb 6, 2016

megies commented Feb 7, 2016

QuLogic commented Feb 8, 2016

MMesch commented Feb 8, 2016

krischer commented Feb 8, 2016

krischer commented Feb 9, 2016

claudiodsf commented Feb 9, 2016

krischer commented Feb 9, 2016

MMesch commented Feb 9, 2016

MMesch commented Feb 9, 2016

MMesch commented Feb 9, 2016

MMesch commented Feb 9, 2016

ThomasLecocq commented Feb 9, 2016

MMesch commented Feb 9, 2016

krischer commented Feb 9, 2016

MMesch commented Feb 9, 2016

MMesch commented Feb 9, 2016

megies commented Feb 9, 2016

MMesch commented Feb 9, 2016

megies commented Feb 9, 2016

QuLogic commented Feb 9, 2016

megies commented Oct 11, 2017

megies commented Oct 11, 2017

MMesch commented Oct 11, 2017 via email

megies Oct 12, 2017

Choose a reason for hiding this comment

megies Oct 12, 2017

Choose a reason for hiding this comment

megies Oct 12, 2017

Choose a reason for hiding this comment

megies Oct 12, 2017

Choose a reason for hiding this comment

MMesch commented Oct 12, 2017 • edited

megies commented Oct 12, 2017

megies commented Oct 12, 2017

MMesch commented Oct 12, 2017 via email

megies commented Oct 12, 2017

megies commented Oct 13, 2017

krischer left a comment

Choose a reason for hiding this comment

krischer Oct 13, 2017

Choose a reason for hiding this comment

megies Oct 13, 2017

Choose a reason for hiding this comment

krischer Oct 13, 2017

Choose a reason for hiding this comment

megies commented Oct 13, 2017

megies commented Oct 13, 2017

megies commented Oct 13, 2017

MMesch commented Feb 5, 2016 •

edited by megies

MMesch commented Oct 12, 2017 •

edited