examples GridSearch.DoubleCouple.py - beachball.png station azimuths are rotated

[mikehagerty]

I'm trying to follow along with last year's virtual workshop on youtube.
On day 1 at [1:13:52] Felix shows the beachball.png produced by the GridSearch.DoubleCouple example.

My results - waveform fits, misfit contours, solution s,d,r - are identical to what he shows, except
the station triangles/names plotted on the beachball look rotated by 180 deg.
e.g., the beachball is the same but the station azimuths are totally different.
For instance, you have TRF plotting at 353 deg (towards north) and mine is rotated to the south.

I thought maybe this could be an upper/lower hemisphere focal mechanism thing, since the upper hemisphere plot is rotated 180deg from the lower hemisphere plot, but as I say, only the station locations (azimuths) are rotated, not the fault planes, so that doesn't make sense.

After some debugging ... I think you might have an error in your example.
For instance, here is the tmp station file that polar1 plots:

cat tmp.20090407201255351DC_beachball.png.sta
BIGB 345.527769 152.458534
PMR 64.454566 128.320995
KASH 338.666196 118.526754
TUPA 157.288772 105.902894
LSUM 170.675982 103.422575
MPEN 206.794148 103.171812
DEVL 175.365655 101.057717
BLAK 136.074440 100.572742
LSKI 200.106721 99.145813
NSKI 223.868000 98.472389
AVAL 169.692917 98.364355
PERI 129.875238 97.862826
SOLD 213.907196 97.457392
SWD 173.872713 96.070030
HEAD 173.413687 95.478760
DIV 97.918007 93.427289
TRF 353.014642 93.206074
EYAK 113.279988 92.876460
PAX 50.896740 92.049359
BMR 98.847989 92.010388

TRF has an azimuth of 353 but a takeoff angle of 93 which is > 90. so in fact,
in the lower hemisphere plot, it is an upgoing ray and should plot 180 deg.
away from the azimuth(or towards the south) as it does in my plot (see attached).

I don't think I have a way to know the takeoff angles you used, so maybe in a different
velocity model it was < 90. for TRF ?

Thanks!
-Mike

Versions:
python: 3.11.0
mtuq: 0.2.0
pygmt: 0.8.0

[rmodrak]

Mike,

Thanks for your careful attention. Your report matches my recollection of an issue from 2021/early 2022, when we began switching to PyGMT and encountered many bugs and regressions.

A couple things that seem relevant:

• We still have separate GMT and PyGMT implementations of plot_beachball, which makes the code fragile
• Until PyGMT becomes more mature, however, it may be difficult overcome this fragility
• A particular stumbling block, it seems, is that PyGMT does not yet implement pspolar-- the function actually responsible for plotting station locations relative to beachball first motions
• For the moment, we have only a crude workaround that involves a Python wrapper over the GMT pspolar executable

From experimenting just now with different MTUQ versions:

• I can confirm-- station locations changed by 180 degrees in mid 2022
• Since mid 2022, station locations are being plotted correctly, I believe, and have undergone no further changes
• However, I have not used polarity measurements much myself, so would like to ask other developers--

Could others please comment on/confirm the above?

thanks,
Ryan

[rmodrak]

Now zeroing in a bit more, I believe the 180 degree rotation occurred starting with the following commit:
2c66213

@thurinj-- I think you may have a better grasp. Could you please double check, as well as comment about @mikehagerty 's overall issue, which seems very thoughtful and plausible? I just noticed that the MTUQ version Mike mentioned (0.2.0) is the up to date version, which I misunderstood in the last post.

[mikehagerty]

Hi.
Thanks for your replies.
In the past I've tried to use pyGMT in my own applications but I've run into bugs and abandoned it.
However, I did debug the beachball plots a bit -- I made it use _plot_beachball_pygmt vs _plot_beachball_gmt --
and the results were identical, so no issue there.