psmeca plotting full moment tensors

[carltape]

Thank you for the sustained support of GMT. I have made some input files for testing plotting routines for full moment tensors. Examples, including figures and input files, can be found here:
https://github.com/carltape/compearth/blob/master/momenttensor/plot/gmt/beachball_plots/
This includes the test script lune_beachballs_min.csh

gmtset PAPER_MEDIA Custom_3.2ix7.7i PAGE_ORIENTATION portrait MEASURE_UNIT inch BASEMAP_TYPE plain TICK_LENGTH 0.0c FRAME_PEN 2p TICK_PEN 2p
psbasemap -JXi/6.72i -R-30/30/-90/90 -Ba10f5g10:" ":/a10f5g10:" ":wesn -G200 -K -V -X0.5 -Y0.5 > lune_ngeo_iplot2_lplot1_kplot1_min.ps
psmeca beachballs_ipts1_iref1_lune_psmeca -JXi/6.72i -R-30/30/-90/90 -Sm0.45/8p -L0.5p/0/0/0 -G255/0/0 -N -O -V >> lune_ngeo_iplot2_lplot1_kplot1_min.ps

which was run using GMT 4.5.3 and using a modified version of utilmeca.c. The output figure is included, and it is also attached here as
lune_ngeo_iplot2_lplot1_kplot1_min_psmeca_custom.pdf

The example, partially converted for GMT 5, is:

gmt gmtset PS_MEDIA Custom_3.2ix7.7i PS_PAGE_ORIENTATION portrait PROJ_LENGTH_UNIT inch MAP_FRAME_TYPE plain MAP_TICK_LENGTH 0.0c MAP_FRAME_PEN 2p MAP_TICK_PEN 2p
gmt psbasemap -JXi/6.72i -R-30/30/-90/90 -Ba10f5g10:" ":/a10f5g10:" ":wesn -G200 -K -V -X0.5 -Y0.5 > test.ps
gmt psmeca beachballs_ipts1_iref1_lune_psmeca -JXi/6.72i -R-30/30/-90/90 -Sm0.45/8p -L0.5p,0/0/0 -W0.5p,0/0/0 -G255/0/0 -N -O -V >> test.ps

This generates the file test.ps, which I have renamed and attached as
lune_ngeo_iplot2_lplot1_kplot1_min_psmeca_gmt5.4.5.pdf

Comparison between the two pdf files reveals two main differences.

1. 14 beachballs are plotted incorrectly in GMT 5.4.5:
   6 are all-white and should be all-red (top).
   6 are all-red and should be all-white (bottom).
   1 is all-white and should have some red (upper middle).
   1 has its colors flipped (lower middle).

For the all-red and all-white beachballs, this means that the eigenvalues of the moment tensors are either all positive (red) or all negative (white). This would seem to be a straightforward check. I'm not sure what fix is needed for the other two beachballs.

2. For the beachballs plotted in GMT 5.4.5, the nodal lines are not plotted with the pen specification. (Except for the one double couple, at center.) This issue is much less significant compared with the first issue. It's possible that I need to fix the syntax, but I did not see how. (Plus, I don't see why the nodal lines would appear for one beachball but not the others.)

The correct plot was generated using a modified version of utilmeca.c, which can be found at the link at the top. It seems like this would be a better starting point than the current version of utilmeca.c. I have not tried to compile GMT 5.4.5 with the modified utilmeca.c to see how it works. Finally, in the repository above I have provided five additional sets of beachballs that would be useful for more comprehensive testing.

System information

• Operating system: Mac OS 10.13.6
• Version of GMT: 5.4.5

[welcome]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. We appreciate that you took the time to contribute!

Please make sure you read our Contributing Guide and abide by our Code of Conduct.

[PaulWessel]

Thanks @carltape, we will have a look at all this as soon as we can. Perhaps @seisman could have a look?

[seisman]

The difference between the old utilmeca.c and our current one is quite large. I'll give a try.

[PaulWessel]

Thanks, hopefully not breaking existing test scripts.

[PaulWessel]

Hi @seisman, any luck testing Carl's revised code? Would be nice to get into GMT 6 beta. We are aiming for June 15.

[seisman]

I'll try it again this weekend.

[carltape]

@seisman I am trying to test the 3-line script (at the top of this issue) in GMT 6.0.0. I am getting some errors and was wondering if you could send me (or post here) your version of the test script, so that I can try. Thanks.

[seisman]

@carltape I modified your script to GMT5 syntax, and included it in GMT's test directory. The test script is available from https://github.com/GenericMappingTools/gmt/blob/6.0/test/seis/full_moment_tensor.sh.

[carltape]

@seisman Thank you! That was enough for me to update my own script here Huge thanks for tackling this issue -- everyone who plots full moment tensors will benefit. There are likely some lingering issues with psmeca, but this was the biggest one.

Rather than open a new issue, can I ask why psmeca requires entries in the 11 and 12 columns? Is that true for all input files for psmeca? I did not see this in the 6.0.0 documentation here. We're not using the -C option, so I don't see why those columns are needed. They were not needed in 5.4.4. (You can see in my script that I created an awk command, but it'd be a pain to do this for all existing psmeca input files.)

[seisman]

Yes, I agree those columns should be optional. I remember that I tried to fix the issue but realized that the current GMT I/O mechanism is too complicated to understand and I had no clue how to fix it.

I've opened a new issue in #2016. Feel free to comment there.

[carltape]

At the start of this issue, I listed input files for six moment tensors, along with a test script. The target, correct output figure is here. The figure produced by GMT 6.0.0 is here.

5/6 input files seem to plot correctly. The problem is associated with the lambda2=0 moment tensors, shown in the sixth set of moment tensors. These produce the errors like
Warning: bigisotest failed for record 9, please report the issue to us!

The input moment tensors are extremely simple, as they only have two non-zero entries.

I documented the testing in this README. Somehow the capabilities in the Dreger-modified utilmeca.c are not being reproduced in GMT 6.0.0. Please let me know if you have questions. Thanks for your help.

[seisman]

@carltape Thanks for your extensive tests.

I'll look into your test script. As I understand it, I should run lune_beachballs_gmt601.csh and get the figure lune_hammer_iplot2_lplot1_kplot1_gmt601.pdf, right? I can't lot the 5th set, because the file ../dfiles//beachballs_ipts1_iref5_lune_psmeca is not found.

[carltape]

@seisman Sorry about the missing input files! They were sitting in my local repo, not checked in. They should be here now.

[seisman]

Great! Now I can reproduce your figure. I'll try debugging this weekend.

[carltape]

@seisman Thanks for the quick update. I confirm here that it works for me: the lambda2=0 beachballs in the 6th set look good. I renamed some of my scripts and output files to better reflect the gmt versions I was testing. The "pr2018" script names will handle the variable number of columns for the psmeca input file (which is not possibly in 6.0.0) -- these scripts can be used with 6.0.0pr2067 as well -- for that version of GMT you get the correct output. I'll plan to update my own tests for a future official release of GMT. Thanks again.