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.