With the newly added chinook/examples scripts, we are able to meaningfully compare the computational cost of MTUQ and legacy CAPUAF inversions.

Currently, a grid search over 50,000 moment tensors (chinook/examples/CapStyleGridSearch.DoubleCouple.py) takes about 30 seconds on a single chinook node. This excludes the initial I/O walltime, which varies depending on cluster usage. A single chinook node has 3 Intel Xeon E5-2620 CPUs, each with 8 cores.

Although the current 30 second walltime represents an acceptable performance level, significant speedup may still be possible. In particular, accelerating the get_synthetics and misfit functions seems promising.

If not carried out thoughtfully, though, acceleration may have negative side effects. So for now we will proceed outside of the main development branch.

Hi all,

I am doing a synthetic test in MTUQ (use synthetics as observed data) and I see that MTUQ uses the ceil method for rounding numbers . In my synthetic test, all the stations are at a distance of 201 km and the MTUQ calculation shows similar results. I added a print command below the line I am highlighting in the link and I got these results in meters for the receiver-station distance:

200999.75896759419 200999.7623664013 201000.0049258777 201000.00492587863 201000.00740668134 ...

Which rounded using the ceil method gives a distance of 201 or 202. My FK GFs are calculated by a distance of 201 km which is the actual distance. Therefore, when MTUQ tries to find the GFs sometimes looks for 201 or 202 and fails with the latest one because there are no GFs at a distance of 202.

May I ask why not use np.round instead of np.ceil where you are forcing to the closest largest integer?

Thanks

Here's what I was aiming for:

(F0, theta, h)
Uniform parameterization for forces, where F0 is magnitude of force in Newtons, theta is azimuth in degrees counterclockwise from east, and h is cos(colatitude)

function to_xyz
converts from (F0,theta,h) to east-north-up

function to_rtp
converts from (F0,theta,h) to up-south-east

After running MTUQ only with surface waves, I got this error:

mtuq/mtuq/graphics/waveforms.py", line 103, inplot_waveforms1   normalize, trace_labels, max_amplitude, total_misfit)NameError: name 'trace_labels' is not defined

I guess this happens because when plotting the waveforms, there is a subroutine expecting a parameter named trace_labels. It seems that the parameter in this latest version of MTUQ should be named trace_label_writer

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

We could do this while still keeping backward compatibility

Hi @rmodrak

With @ammcpherson noticed that the example test_greens_SPECFEM3D_SAC.py is not working after the last merge. I noticed the same discrepancy after running my scripts using the latest version of MTUQ. @ammcpherson noticed that this concerns the default units expected in the SPECFEM3D-derived GFs. After the last merge, the Specfem3D-derived GFs are expected in velocity rather than displacement. The easiest solution is to change these lines of code back to displacement, but it is important to verify that this solution is compatible with the new HMC feature where SPECFEM3D_SGT-derived GFs are read.

Thanks!

I tried to install mtuq but it fails because it cannot build instaseis, because gfortran is missing. Should I install gfortran globablly so any conda env can access? What about instaseis? Should this be install in the mtuq conda env or by itself?

I'm not able to install gfortran. I'm running on a relatively new Mac OSX 13.5.2 (Ventura) with Apple M2 Max chip. Is there a way to bypass instaseis and use other Greens functions, say fk?

MTUQ has four options for handling the arrival times (pick_type):

- ``'taup'``
  calculates P, S arrival times from Tau-P model (uses `obspy.taup.TauPyModel.get_arrival_times`)

- ``'FK_metadata'``
  reads P, S arrival times from FK metadata

- ``'SAC_metadata'``
  reads P, S arrival times from SAC metadata fields `t5`, `t6`

- ``'user_supplied'``
  reads P, S arrival times from columns 8, 10 of `capuaf_file`

At present MTUQ cannot accommodate the case when we want to provide (user_supplied) arrival time information at only selected stations. If we update the arrival time in weight files for only a few stations, leaving others unchanged, it will take the arrival times for other stations as 0. It will be good to have a fallback_option (say 'taup'). In such a case, one can use user_supplied arrival times when available in weight file, and falling back to taup if not supplied.

Note that MTUQ installation is already automatically tested. Every time there is a pull request, GitHub’s continuous integration runs something like git clone https://github.com/uafgeotools/mtuq.git && cd mtuq && pip install -e ..

As a result, MTUQ and all dependencies specified in setup.py get installed somewhere on GitHub server. If the installation procedure fails, a prominent red mark is displayed on mtuq’s GitHub page and an email notification may be sent out.

Currently, mtuq's GitHub CI tests all dependencies except

• GMT (not a PyPI package)
• mpi4py (not included in setup.py, because GitHub CI does not provide a parallel environment)
• instaseis (not included in setup.py, because GitHub CI does not correctly test Fortran extension models)

To remove these exclusions, we could explore conda based CI strategies (e.g. https://autobencoder.com/2020-08-24-conda-actions/).

Note that instaseis itself uses conda directives in its Travis CI tests. Also note the very involved Fortran compilation options .

Station location codes can be important, especially for data sets using GSN stations, which have multiple sensors and location codes. It would be good to ensure that location codes are tracked and also displayed on the waveform fits plot and wherever else the station label is used. Many location codes will be blank; if they're not blank, they're often 00 or 10. Not urgent.

pygmt.exceptions.GMTCLibError: Module 'meca' failed with status code 72:
meca [ERROR]: Option -c is not a recognized common option

Since waveform modeling for small events, particularly in volcanic regions, can be difficult would it be possible to introduce amplitude ratios (P/SV and P/SH) for 1D models into the estimation in addition to just the polarities?

I wanted to push the recent changes through quickly to fix the conda installation procedure. But in the process, it appears the SCOPED badge was broken? Perhaps we could discuss next month, when everyone is not so overcommitted.

Reqest ability to modify GFs output path. Writing files within the mtuq install folder becomes a problem on a multi-user system where not everyone has write permissions to that path.

To contribute changes to the upstream repository:

• if you haven't yet, create your own fork of the upstream (i.e. uafgeotools) repository
• submit changes via a pull request from your fork
• after submitting a pull request, double check that the continuous integration tests pass by looking to see if a green checkmark appears next to the pull request in the GitHub web interface.  If a red X appears instead, make any necessary fixes
• Carl, Julien, or I will review the pull request

To obtain the latest changes from the upstream repository:

• first, make sure there are no uncommitted changes in your local fork (the git stash command can be useful for this purpose)
• then, within your local fork

>>> git remote set-url upstream [email protected]:uafgeotools/mtuq.git 
>>> git fetch upstream master
>>> git rebase upstream/master

I've been getting errors when running the MTUQ container on TACC's frontera through apptainer. The errors have been indeterminant, however have always happened after the third "about 75 percent finished" message. See below for the std out, but I have also gotten error like malloc(): invalid size (unsorted) , double free or corruption (out) , corrupted size vs. prev_size in fastbins

The sif image was freshly pulled and it is the newest version.

c202-001[clx](423)$ APPTAINERENV_SYNGINE_CACHE=syngine_output ibrun apptainer run mtuq_ubuntu20.04.sif python3 /home/scoped/mtuq/examples/DetailedAnalysis.py
TACC:  Starting up job 5796560 
TACC:  Starting parallel tasks... 
  about 0 percent finished
  about 25 percent finished
  about 50 percent finished
  about 75 percent finished
  about 0 percent finished
  about 25 percent finished
  about 50 percent finished
  about 75 percent finished
  about 0 percent finished
  about 25 percent finished
  about 50 percent finished
  about 75 percent finished
free(): invalid pointer

@qnishida pointed me to the following Instaseis issue #82 which shows a polarity issue on the Radial component synthetics for a point force source. The issue was also raised in Instaseis issue #77, and never got addressed.

I made a quick test, based on MTUQ objects, to replicate @liamtoney results. By modifying one of MTUQ example scripts so that data import and green's function download is the following:

data = read(path_data, format='sac', 
            event_id=event_id,
            station_id_list=station_id_list,
            tags=['units:cm', 'type:velocity']) 
data.sort_by_distance()

#Here, we create a source is located westward of the receiver, so that the direction of propagation is East
origin = Origin({
'time': '2009-04-07T20:12:55.000000Z',
'latitude': data[0].station.latitude,
'longitude':data[0].station.longitude-0.01, 
'depth_in_m': 0,
})

greens = download_greens_tensors(stations, origin, model, include_mt=False, include_force=True)

The force's response for a station, eastward of the source can be computed with the following:

greens[0].get_synthetics(Force([0,0,1]), components=['Z','R','T'])

for which the synthetics are the following:

which effectively reproduces the results shown in Instaseis' issue #77. Provided we expect the force to be positive away from the source, we should have a positive motion in both Z and R.

MTUQ can be patched for this by modifying the sign of the green's function in _precompute_force of line 121 to 123 to:

R0 = self.select(channel="R0")[0].data * -1
R1 = self.select(channel="R1")[0].data * -1
R2 = self.select(channel="R2")[0].data * -1

Provided we implement this fix, I also have a functional test ready, that will download synthetic seismograms from Syngine and check that E = -R. If the test fails, it means that Instaseis has been patched and that we can revert this fix.

The waveform fits plots for moment tensors has a nice beachball at the upper left. When the user is trying a force source, we need to display a "force sphere" showing a symbol (dot/circle/symbol) for the best-fitting force direction. W and E meridians need to be indicated somehow. Thank you!

There seems to be a point in MTUQ where, while the data is being read in (for example), the event depth is read in from the SAC headers. Later, when the data processing class is mapped onto the data, a np.ceil operation is applied to this event depth.

This causes problems when using local FK databases if said databases contain a boundary at that depth. For example, if I have an event at depth 8850 meters, and MTUQ pulls that depth from the SAC header and rounds up to 9000 meters, then MTUQ will fail if I am using a FK database with a boundary at 9000 meters depth. The error it throws is:

 File "run_mtuq.py", line 181, in <module>
    data_bw = data.map(process_bw)
  File "/home/ammcpherson/REPOSITORIES/mtuq/mtuq/dataset.py", line 143, in map
    processed += [function(stream, *args)]
 File "/home/ammcpherson/REPOSITORIES/mtuq/mtuq/process_data.py", line 497, in __call__
    sac_headers = obspy.read('%s/%s_%s/%s.grn.0' %
  File "/home/ammcpherson/.conda/envs/mtuq/lib/python3.8/site-packages/decorator.py", line 232, in fun
    return caller(func, *(extras + args), **kw)
  File "/home/ammcpherson/.conda/envs/mtuq/lib/python3.8/site-packages/obspy/core/util/decorator.py", line 300, in _map_example_filename
    return func(*args, **kwargs)
  File "/home/ammcpherson/.conda/envs/mtuq/lib/python3.8/site-packages/obspy/core/stream.py", line 212, in read
    st = _generic_reader(pathname_or_url, _read, **kwargs)
  File "/home/ammcpherson/.conda/envs/mtuq/lib/python3.8/site-packages/obspy/core/util/base.py", line 700, in _generic_reader
    raise IOError(2, "No such file or directory", pathname)
FileNotFoundError: [Errno 2] No such file or directory: '/store/wf/FK_synthetics/scak/scak_9/69.grn.0'

I'm not sure how to address this, but currently it means there are events that are simply unusable in MTUQ using local FK databases. I have run into at least two of them so far. Perhaps if the function instead took the event depth from the Origin object we set up earlier in the code, so that the user can fudge the depth as needed to account for these boundaries?

Thanks.

Hi,
I am trying to apply statics in MTUQ by modifying 'apply_statics=True' in mtuq/process_data.py as well as the weight file. The code is running well when I apply 'download_greens_tensors' (downloads Green's tensors from syngine). But the problem lies with 'get_greens_tensors' (extracts green's tensors from FK database). I am attaching the snapshot of the error with this message. Hope you will resolve the issue.

Thanks
Rinku

Working our way down the list of priorities, faster handling of sources is now a goal.

Some things we might want to try:

• C extension for to_mij?

• as_array method for Grid, UnstructuredGrid?

• finish implementing MomentTensor.change_convention

• use MomentTensor.as_array in GreensTensor.get_synthetics

• new CompositeSource class

Looking ahead to hypocenter grid searches or random walks, some notes on instaseis performance:

• major instaseis bottleneck (gets called 10 times per station)

• the cost of numerically differentiating displacement to get strain doesn't seem to matter

due to changes in conda itself, these installation instructions no longer work:
https://uafgeotools.github.io/mtuq/install/env_conda.html

Here are details about how conda installation fails:
https://github.com/rmodrak/mtuq/actions/runs/3247505451/jobs/5327561521

it seems we encounter this only because we load NumPy libraries in Python/C extensions

I think a relatively quick fix should be possible...

To allow support for other file formats beyond SAC, please try to adopt the following

• the initial mtuq.read() call handles all file format dependent metadata parsing
• after that, metadata are stored in a common internal structure that works the same regardless of the original file format
• in particular, this means not using observed data SAC headers at any point after the initial mtuq.read() call
• the whole motivation is that data processing and inversion functions should work the same, regardless of the original file format

It seems that MTUQ is not implemented for single-station plotting, although the grid-search seems to work well.

I have this error and no files are generated (NetCDF, JSON, ps, and png):

Generating figures...

Traceback (most recent call last):
File "/home/jovyan/scoped/mtuq_workshop_2022/GridSearch.DoubleCouple_options.py", line 235, in
plot_data_greens2(event_id+'DC_waveforms.png',
File "/home/jovyan/scoped/mtuq/mtuq/graphics/waveforms.py", line 314, in plot_data_greens2
plot_waveforms2(filename,
File "/home/jovyan/scoped/mtuq/mtuq/graphics/waveforms.py", line 134, in plot_waveforms2
fig, axes = _initialize(
File "/home/jovyan/scoped/mtuq/mtuq/graphics/waveforms.py", line 357, in _initialize
_hide_axes(axes)
File "/home/jovyan/scoped/mtuq/mtuq/graphics/waveforms.py", line 643, in _hide_axes
for col in row:
TypeError: 'AxesSubplot' object is not iterable

Hello everyone,
I want to plot each value of misfit against omega (angular distance between two MTs) for both regular as well as random points. Is there any code available in MTUQ for this? I have not found such function inside MTUQ. Hope you will look into the matter.

Thank you.
Rinku

MTUQ includes an option to read a CAPUAF-style text file with user-supplied data weights and phase picks.

Jonas K (@kintner92) discovered that, while fully implemented in the pure Python misfit functions, data weights are not fully implemented for the Python/C misfit function. More precisely, user-supplied weights get cast to 0s or 1s, so that rather than a data weighting matrix per se, only a binary mask is implemented.

A bugfix is required so that any nonnegative scalar can be used as a weight (not just 0s and 1s).

It looks like mtuq continues to crashes when the event lies at the model interface.
To avoid this, we generally change the inversion depth by +-1 km. In original version of CAP this was handled directly during inversion, however, in mtuq this requires changing the SAC headers.
I believe this capability should be available in mtuq.

A bug in newest version (1.16.0) of NumPy causes failed ObsPy imports. A workaround is to specify numpy==1.15.4 in MTUQ's setup.py file.

There are currently three implementations of the misfit function:

• simple
• fast1
• fast2

The simple version is already very clean--this easy-to-read version exists exists to understand what the code is doing and to check the correctness of the other two versions.

A test comparing the three misfit functions is a high-level priority. This test should be automatically generated from code_generator.py and include the following cases

• single and multiple stations
• with and without time shifts
• with and without independently varying time shifts
• uniform and nonuniform weighting

Without sacrificing speed, the two fast implementations could probably be cleaned up a lot, but this is only a mid-level priority.

The current depth plot implementation is lacking the option to display force amplitude (as inputs are all directed toward MT grid format).

Which of these makes the most sense?

1. pymc
   Use pymc in combination with a pure Python or C extension misfit function. I started looking into this-- my sense is it would be hard to get good enough performance. NumPy and pymc functions are simple to work with but very complex under the hood. Understanding these implementation details, which seem important for good performance, may be just about as time consuming as creating one's own C/C++ extension modules.

2. C++
   Improve the existing extension module by adding a C++ class similar to misfit.level1.Helper. MPI parallelization over sources is not applicable to MCMC. OPENMP paralleziation over stations might be possible. Use numerical libraries (e.g. https://github.com/kthohr/mcmc) to implement Metropolis-Hastings MCMC in lune or Mij parameters. Hamiltonian MCMC in Mij parameters should also be possible, since get_synthetics is linear w.r.t. these parameters.

3. Bayesian Earthquake Analysis Tool
   find some way to interface with BEAT, which seems to have very good MCMC impelemenations. Unfortunately, BEAT lies outside the ObsPy ecosystem.

travis-ci.org has been discontinued, so github.com/uafgeotools/mtuq is now without automatic testing of new commits.

It appears GitHub itself is now offering such automatic testing, without reliance on third parties. The way you set up this GitHub testing may be different, however, so it would be good to try it out in a fork or dummy repo first.

This is the first time I noticed these changes. If anyone else has been following more closely, feel free to comment.

Multiple lines with the 'tmp.'+filename+... in the function def _plot_beachball_gmt(...) in the script beachball.py may raise errors when passing the absolute filename.

GMT/PyGMT provide graphics capabilities missing from matplotlib and other native Python packages. Interfacing with these external graphics libraries, however, causes all sorts of problems (non-backward compatible changes between GMT versions; missing or immature PyGMT implementations; dependency conflicts etc.).

In response, I’m working on isolating the graphics library calls as much as possible within backends. Ultimately, we can use these backends to

• adopt PyGMT without breaking the existing GMT plotting functionality.
• fall back gracefully from graphics library errors

The idea is that, if a problem arises, we only need to adjust the backend, not the entire function stack. For an example, see _plot_lune, which contains the the _plot_lune_gmt backend. When the time comes, we can add _plot_lune_pygmt without having to reimplement the whole function stack.

Such changes will pay off long term through more improved installation and user experience.

Time discretization

• MTUQ uses a centroid origin time convention, in which t=0 corresponds to the mean source excitation (usually equivalent to the peak source excitation)

• CAPUAF uses a rupture origin time convention, in which t=0 corresponds to the initial source excitation

• MTUQ still works even if the time discretization varies from station to station, or if data and Green's functions have different time discretization.  Certain code optimizations may be lost, however, making run times much longer.  Also such cases are not yet well tested.

Time shifts

• Systematic difference between CAP and MTUQ time shifts arise from the way CAP implements source-time function convolution. CAP's "conv" function results in systematic magnitude-dependent shifts between centroid origin times and arrival times.

• The amount of the time shift is half the sum of the earthquake's rupture time and rise time, as given by relations in the CAP Perl wrapper, t_offset = (cap_rupture_time(Mw) + cap_rise_time(Mw))/2

• Shifting a processed MTUQ trace relative to a processed CAP trace does not completely remove the effect of the idiosyncratic time shift;  the two traces will still be tapered differently at the end points.

Distances

• MTUQ uses obspy.geodetics to compute station-event distances, with results that differ slightly from CAPUAF. 

• Distances shown in CAPUAF figure labels are sometimes off by an integer compared with the actual values CAPUAF uses in computing synthetics.

• CAPUAF rounds down when working with FK Green's functions, while MTUQ rounds to the nearest integer.

• MTUQ also supports AxiSEM and SPECFEM3D Green's functions, in which case GLL interpolation is used rather than rounding.

Cross correlations

• MTUQ provides the option to pad synthetics prior to cross correlating them with data.  CAPUAF does not, so cross correlation values become less accurate at large lag times.

It would be nice to have some indication in the waveform fits plot of any constraint that was used. Consider the case where any parameter (Mw, strike, dip, rake, gamma, delta) is fixed. Then a label could be added such as
delta 0*
or
delta (0)
If that were the only constrained parameter, then it would be a deviatoric inversion. Or
gamma 0* delta 0*
or
gamma (0) delta (0)
would indicate a double couple inversion.
Or Mw (6.21) or Mw 6.21* would indicated a fixed-magnitude constraint. And so on. This would allow us to distinguish between cases where parameters were constrained and where they were unconstrained. At present, the label
gamma 0 delta 0
could either be a DC-constrained inversion or a FMT inversion that happened to find a double couple as the best-fitting moment tensor.

It looks like the current version of MTUQ does not have the capability of plotting depth against misfit function for random grid search. Although I am able to achieve the MT solution using DoubleCoupleGridRandom funtion in GridSearch.DoubleCouple+Magnitude+Depth.py .
I believe there should be an option in the function plot_misfit_depth for random distribution also.

Thank you.

Julien pointed out that when many hypocenters or depths are considered, Green's function data processing becomes a bottleneck.

One idea would be to start by implementing this method:

I recently ran into an issue where, whilst attempting to use a Green's Function database generated in SPECFEM3D_GLOBE, the following error was thrown:

Reading Greens functions...

/home/ammcpherson/REPOSITORIES/mtuq/mtuq/greens_tensor/SPECFEM3D.py:26: UserWarning: SPECFEM3D modules have not yet been tested.
  warnings.warn("SPECFEM3D modules have not yet been tested.")
Processing Greens functions...

Traceback (most recent call last):
  File "3D_run_mtuq.py", line 237, in <module>
    greens_bw = greens.map(process_bw)
  File "/home/ammcpherson/REPOSITORIES/mtuq/mtuq/greens_tensor/base.py", line 363, in map
    [function(tensor, *args)]
  File "/home/ammcpherson/REPOSITORIES/mtuq/mtuq/process_data.py", line 599, in __call__
    cut(trace, starttime, endtime)
  File "/home/ammcpherson/REPOSITORIES/mtuq/mtuq/util/signal.py", line 22, in cut
    raise Exception('The chosen window ends after the trace.  Consider '
Exception: The chosen window ends after the trace.  Consider using an earlier window, or to automatically pad the end of the trace with zeros, use mtuq.util.signal.resample instead

I would receive this error regardless of if I was using taup or FK metadata to calculate P and S wave arrivals. After doing a little digging, I found that the windows for the data and the windows for the Green's Functions had extremely different values for the same station:

As it turns out, I had mistyped the time in my Origin object as '2022-11-30' instead of '2020-11-30'. However, this does raise the question of why does MTUQ use the origin time in SAC headers for the data to create a window for a trace, but uses the origin time in the Origin object to create a window for the synthetics/Green's Functions for SPECFEM3D databases?

The mtuq.event.Force() is currently lacking a method that outputs the force vector as a dictionary.

It would be convenient to have a similar method that in [mtuq.event.Force()] (

so that it is easier to build an output dictionary to save it as a JSON file (as in DetailedAnalysis.py)

Please consider some additional annotations for the MTUQ waveform fits plots. This excerpt from Figure S2 of Alvizuri and Tape (2018), generated with the older code, includes some items that we found useful.

It shows the length of the P window and surface wave windows (1s -- not used here; 500 s), the number of stations used (N = 34), the number of P window measurements (Np = 0), and number of surface wave measurements (Ns = 102). Ideally it would list the number of first-motion polarities used in the header, but this information is listed as the 4th text row on the left. So "1 (0.74)" represents the measured polarity (1 = up) and the theoretical amplitude factor for this particular mechanisms (0.74). A theoretical value near 0 would represent a near-nodal station. Thanks.

On the waveform fits plots, the station distance is listed in km. When stations are >3000 km or so, it makes sense to list the station distance as an arc-distance in degrees.

Hello,

We would like the capability to read in SPECFEM3D-generated, source-side force Green's Functions. Currently, MTUQ is able to read in SPECFEM3D source-side MT Green's Functions and compute synthetics from them. MTUQ also has the capability to use AxiSEM force Green's Functions to compute synthetics.

The end goal would be something similar to what @SeismoFelix has established here, in which we would like to compute source-side GFs in a chosen velocity model, and use them in an inversion with MTUQ.

Thank you,
Amanda

I am facing problem with attaching first motion polarity in MTUQ. Can you please provide a solution.

@SeismoFelix raised an interesting point concerning the current implementation of the double-couple misfit map.
It currently relies on the built-in min() method in xarray.DataArray format, in which the results are stored.
What this does is that it aggregates the minimum misfit values along all dimensions of the misfit space.

It would be nice to add an option that first isolates the misfit cube (in strike, dip, and slip dimensions) containing the minimum misfit, and then return the three slices that cross the minimum misfit value.

This is the function that needs to be updated:

potentially by adding an alternative function to

After the installation, I tried to run an example but there where an error when importing obspy

 python SerialGridSearch.DoubleCouple.py 

Traceback (most recent call last):
  File "SerialGridSearch.DoubleCouple.py", line 6, in <module>
    from mtuq import read, open_db, download_greens_tensors
  File "/Users/felix/Documents/INVESTIGACION/software/mtuq/mtuq/__init__.py", line 18, in <module>
    from mtuq.dataset import Dataset
  File "/Users/felix/Documents/INVESTIGACION/software/mtuq/mtuq/dataset.py", line 2, in <module>
    import obspy
  File "/Users/felix/opt/anaconda3/envs/mtuq/lib/python3.8/site-packages/obspy/__init__.py", line 39, in <module>
    from obspy.core.utcdatetime import UTCDateTime  # NOQA
  File "/Users/felix/opt/anaconda3/envs/mtuq/lib/python3.8/site-packages/obspy/core/__init__.py", line 124, in <module>
    from obspy.core.utcdatetime import UTCDateTime  # NOQA
  File "/Users/felix/opt/anaconda3/envs/mtuq/lib/python3.8/site-packages/obspy/core/utcdatetime.py", line 27, in <module>
    from obspy.core.util.deprecation_helpers import ObsPyDeprecationWarning
  File "/Users/felix/opt/anaconda3/envs/mtuq/lib/python3.8/site-packages/obspy/core/util/__init__.py", line 27, in <module>
    from obspy.core.util.base import (ALL_MODULES, DEFAULT_MODULES,
  File "/Users/felix/opt/anaconda3/envs/mtuq/lib/python3.8/site-packages/obspy/core/util/base.py", line 36, in <module>
    from obspy.core.util.misc import to_int_or_zero, buffered_load_entry_point
  File "/Users/felix/opt/anaconda3/envs/mtuq/lib/python3.8/site-packages/obspy/core/util/misc.py", line 214, in <module>
    loadtxt(np.array([0]), ndmin=1)
  File "/Users/felix/opt/anaconda3/envs/mtuq/lib/python3.8/site-packages/numpy/lib/npyio.py", line 1086, in loadtxt
    ncols = len(usecols or split_line(first_line))
  File "/Users/felix/opt/anaconda3/envs/mtuq/lib/python3.8/site-packages/numpy/lib/npyio.py", line 977, in split_line
    line = line.split(comment, 1)[0]
AttributeError: 'numpy.int64' object has no attribute 'split'

When I try to use obspy in python 3.7 I have no problems. This only happens with mtuq environment (obspy in python 3.8).

I've just finished following along with the 2022 MTUQ Virtual Workshop.
I was able to reproduce all of the plots shown.

One thing I noticed is that, although my waveform fits and mechanisms were very similar, often
my L2 misfit (as reported on the waveform plot) was ~1.4 - 2 times larger than
what was shown in the workshop video.

I'm just curious - did the misfit calculation change (e.g., by factor of sqrt(2)) since the
workshop ?

I know if I had enough patience I could dig through the commit history of the
misfit module and try to figure it out, but I thought I'd just ask instead.

Thanks!
-Mike

TODO list
(For new TODO items, update this comment rather than add a new comment.)

• finish implementing graphics.waveforms.plot_time_shifts

• option to display best-fitting magnitude on force plots

• display force orientation in upper left corner of waveform comparison plot by modifying mtuq.graphics.header.ForceHeader

• why are isotropic beachballs from gmt_beachball slightly off center?

• can we safely remove gmt.conf and gmt.history files when plotting scripts finish, e.g. what if multiple instances of the script are running?

• mtuq.graphics.uq.dc: finish implementing double couple plotting functions

• additional colorbar placement options (e.g., bottom center, upper right)

• for lune, vw, and force plots, write likelihood colorbar values in terms of PI

• support all GMT file formats in gmt_beachball

• plot_beachball: option to display station locations on beachball

• add backend keyword argument to plotting utilities (accepts either matplotlib or GMT)

• increase white space around vw figures to avoid marker getting clipped at edges

COMPLETED
(move items from TODOLIST to here as they are finished)

• optional colorbar labels
• implemented plot_misfit_depth
• option to display best-fitting beachballs on lune plots
• export panoply, hot color maps from GMT to matplotlib
• more flexible keyword arguments, e.g. add_colorbar -> colorbar_type
• use intelligent color limits for misfit and likelihood colorbars