Strange bug with psmeca and Oblique Mercator projection

Q&A
1

Hi, I ran into a strange issue while plotting focal mechanisms. Some of my data are moment tensors generated from strike/dip/rake data (pure double couple). I am using -Sm to plot the moment tensors instead of -Sa because I am combining different datasets, some of which have non-double-couple components. When I plot these in using -JOa and psmeca -Sm, I see that a subset of events will have the nodal planes plotted using the -T0 option rotated around a vertical axis as compared to the focal mechanism (which always seems to plot correctly). The behavior is sensitive to small variations in the moment tensor. Notably, I haven’t run into this effect when plotting moment tensors from GCMT, which almost always have a significant non-double-couple component.

Here’s a test script that replicates the behavior using moment tensors generated from slightly different strikes, and an image showing how the rotation amount changes with obliquity of the projection. The labels are the value of strike that went into generating each moment tensor.

As always, thanks for any help.

gmt --version = 6.4.0_a2eafda_2022.02.21

#!/usr/bin/env bash

cat <<-EOF > testfile.dat
90 20 40 -0.566585 -1.09383 1.66042 0.929833 0.243698 0.473602 21 s330.1
90 21 40 -0.566585 -1.0923 1.65888 0.930343 0.242165 0.47871 21 s330.2
90 22 40 -0.566585 -1.09077 1.65684 0.930854 0.240632 0.483309 21 s330.3
90 23 40 -0.566585 -1.08872 1.65531 0.931365 0.2391 0.488418 21 s330.4
91 20 40 -0.566585 -1.08719 1.65377 0.931876 0.237567 0.493016 21 s330.5
91 21 40 -0.566585 -1.08566 1.65173 0.931876 0.235524 0.497614 21 s330.6
91 22 40 -0.566585 -1.08361 1.6502 0.932387 0.233991 0.502723 21 s330.7
91 23 40 -0.566585 -1.08208 1.64815 0.932898 0.232458 0.507321 21 s330.8
92 20 40 -0.566585 -1.08004 1.64662 0.933409 0.230925 0.511919 21 s330.9
92 21 40 -0.566585 -1.0785 1.64458 0.93392 0.229393 0.517028 21 s331.
92 22 40 -0.566585 -1.07646 1.64304 0.93392 0.227349 0.521626 21 s331.1
92 23 40 -0.566585 -1.07493 1.641 0.934431 0.225816 0.526224 21 s331.2
93 20 40 -0.566585 -1.07288 1.63947 0.934942 0.224284 0.531333 21 s331.3
93 21 40 -0.566585 -1.07084 1.63742 0.935452 0.222751 0.535931 21 s331.4
93 22 40 -0.566585 -1.06931 1.63538 0.935452 0.221218 0.540529 21 s331.5
93 23 40 -0.566585 -1.06726 1.63385 0.935963 0.219686 0.545127 21 s331.6
EOF

gmt psmeca testfile.dat -Ggray -T0/1p,black -Sm0.25i -L0.25p,black -R-300/300/-300/300+uk -JOa91.5/21.5/90/2.5i -Bbltr+t"-JOa91.5/21.5/90/2.5i" --FONT_TITLE=10p -Bg2 -K > cmt_test.ps
gmt psmeca testfile.dat -Ggray -T0/1p,black -Sm0.25i -L0.25p,black -R-300/300/-300/300+uk -JOa91.5/21.5/110/2.5i -Bbltr+t"-JOa91.5/21.5/110/2.5i" --FONT_TITLE=10p -Bg2 -X3i -K -O >> cmt_test.ps
gmt psmeca testfile.dat -Ggray -T0/1p,black -Sm0.25i -L0.25p,black -R-300/300/-300/300+uk -JOa91.5/21.5/130/2.5i -Bbltr+t"-JOa91.5/21.5/130/2.5i" --FONT_TITLE=10p -Bg2 -X3i -O >> cmt_test.ps
gmt psconvert cmt_test.ps -Tf -A+m0.5i

image
image2174×810 106 KB

2

Hi, I have figured out the bug in psmeca.c

For some focal mechanisms, meca.NP1.str and meca.NP2.str have delaz added twice. This happens when meca.NP1 and meca.NP2 are recalculated from the principal axes, which have already had delaz added, and then delaz is added again a bit later.

I have fixed this by adding a check for whether we have already added delaz. I’m not very good with git so I’m unsure how to propose changes to the code directly. Here’s a diff on my fixed psmeca.c and the old version - I’m happy to propose the fix in a better more git-friendly way with a bit of guidance.

diff …/src/seis/psmeca.c …/…/gmt-old/src/seis/psmeca.c

673c673
< 	bool transparence_old = false, not_defined = false, has_text, added_delaz = false;
---
> 	bool transparence_old = false, not_defined = false, has_text;
984d983
< 				gmt_geo_to_xy (GMT, in[GMT_X], in[GMT_Y], &plot_x, &plot_y);
986,987c985
< 				/* Keep track of whether we have corrected for delaz */
< 				added_delaz=false;
---
> 				gmt_geo_to_xy (GMT, in[GMT_X], in[GMT_Y], &plot_x, &plot_y);
1051d1048
< 						added_delaz = true;
1069,1072c1066,1067
< 					if (! added_delaz) {
< 						meca.NP1.str = meca_zero_360(meca.NP1.str + delaz);
< 						meca.NP2.str = meca_zero_360(meca.NP2.str + delaz);
< 					}
---
> 					meca.NP1.str = meca_zero_360(meca.NP1.str + delaz);
> 					meca.NP2.str = meca_zero_360(meca.NP2.str + delaz);
1087,1090c1082,1083
< 					if (! added_delaz) {
< 						meca.NP1.str = meca_zero_360(meca.NP1.str + delaz);
< 						meca.NP2.str = meca_zero_360(meca.NP2.str + delaz);
< 					}
---
> 					meca.NP1.str = meca_zero_360(meca.NP1.str + delaz);
> 					meca.NP2.str = meca_zero_360(meca.NP2.str + delaz);

And the working test case showing the fix:

image
image2280×856 111 KB

3

Ping @moderators

4

I’ve asked @maxrjones to work with Kyle on this.

5

Hi @kyleedwardbradley, thanks for the detailed bug report and proposed fix!

The most helpful path forward would be to propose the fix through a “pull request” on GitHub, which is basically a request to incorporate your changes into the GMT code base (https://github.com/GenericMappingTools/gmt). I can provide some instructions for making a pull request. To start with, you would need an account on GitHub (https://github.com/). Do you already have an account or would you be willing to create one to propose this bug fix?

6

Hi Meghan, I poked around and found the Contributors Guide and followed the instructions to create a fork and make a pull request. I also added an issue and linked it to the PR. Hopefully I did things correctly!

Kyle

7

Great, thanks for opening the pull request! We’ll give it a review soon!