Wrap meca

[tjnewton]

Description of proposed changes
Add meca functionality to plot focal mechanisms. Started during SciPy2020 sprint.

Adds meca function to base_plotting.py to add focal mechanism plotting to pygmt.

# MWE:
import pygmt

fig = pygmt.Figure()
fig.meca(focal_mechanism=[239.384, 34.556, 12., 180, 18, -88, 0, 72, -90, 5.5, 0, 0, 0],
         region=[239, 240, 34, 35.2], convention='c2c', projection='m4c')
fig.show(method='external')

Live documentation preview at https://pygmt-git-fork-tjnewton-add-meca.gmt.now.sh/api/generated/pygmt.Figure.meca.html

Reminders

• Run make format and make check to make sure the code follows the style guide.
• Add tests for new features or tests that would have caught the bug that you're fixing.
• Add new public functions/methods/classes to doc/api/index.rst.
• Write detailed docstrings for all functions/methods.
• If adding new functionality, add an example to docstrings or tutorials.

[welcome]

💖 Thanks for opening this pull request! 💖

Please make sure you read our contributing guidelines and abide by our code of conduct.

A few things to keep in mind:

• If you need help writing tests, take a look at the existing ones for inspiration. If you don't know where to start, let us know and we'll walk you through it.
• All new features should be documented. It helps to write the docstrings for your functions/classes before writing the code. This will help you think about your code design and results in better code.
• No matter what, we are really grateful that you put in the effort to do this! 🎉

[liamtoney]

@tjnewton and I will work on this more tomorrow, adding tests and whatnot.

[rmodrak]

My question for you guys, as end-users, is: In what formats would you want to provide MT/FMT info? As a filename? List of params? Dictionary? Something with Pandas?

Hi Liam, It sounds like you are looking to replace the classic psmeca command syntax, in which values are read from an ascii table, with something more Pythonic?   Then perhaps one (or both) of the following  conventions from obspy:

• a dictionary with keys 'm_rr', 'm_tt', 'm_pp', 'm_rt', 'm_rp', 'm_tp'
• a list of moment tensor elements in the order [m_rr, m_tt, m_pp, m_rt, m_rp, m_tp]

[liamtoney]

Hi Ryan, thanks for reminding us that we should try to align with existing libraries! Tyler and I were discussing two options — the original file input as well as an array input with the same info for non-file usage — but I think your suggestion is much more Pythonic. So we can add a third option for a keyword (i.e., dict) based option.

---

@tjnewton I propose the following sketched-out method definition, similar to what we had but with a third option that involves providing a dictionary for spec (same as focal_mechanism above, shortening here for demo only) along with lon, lat, and depth required. This is more verbose source code but I think will make for a nicer user experience, such as:

fig.meca(
    lon=239.384,
    lat=34.556,
    depth=12.0,
    spec=dict(
        strike1=180,
        dip1=18,
        rake1=-88,
        strike2=0,
        dip2=72,
        rake2=-90,
        mantissa=5.5,
        exponent=0,
    ),
    region=[239, 240, 34, 35.2],
    **other_kwargs
)

All names/details are drafts of course.

def meca(
    self,
    lon=None,
    lat=None,
    depth=None,
    spec=None,
    convention=None,
    plot_lon=None,
    plot_lat=None,
    text=None,
    **kwargs
):

    if type(spec) is dict:  # The user opted for the nice keyword-based approach
        if (
            lat is None or lon is None or depth is None
        ):  # They need to give the location!
            raise Error('Location not fully specified.')

        # These are constants related to GMT and could be defined elsewhere (top of file?)
        AKI_PARAMS = [
            'strike',
            'dip',
            'rake',
            'magnitude'
        ]
        GCMT_PARAMS = [
            'strike1',
            'dip1',
            'rake1',
            'strike2',
            'dip2',
            'rake2',
            'mantissa',
            'exponent',
        ]

        # Aki and Richards
        if set(spec.keys()) == set(AKI_PARAMS):
            convention = 'a'
            foc_params = AKI_PARAMS

        # GCMT
        elif set(spec.keys()) == set(GCMT_PARAMS):
            convention = 'c'
            foc_params = GCMT_PARAMS

        # Other cases....
        elif ...:
            pass  # Same strategy as above

        else:
            raise Error('Focal mech params not understood')

        # Construct the vector (note that order matters here, hence the list comprehension!)
        spec = [lon, lat, depth] + [spec[key] for key in foc_params]

        # Add in plotting options, if given, otherwise add 0s as required by GMT
        for arg in plot_lon, plot_lat, text:
            if arg is None:
                spec.append(0)
            else:
                spec.append(arg)

    elif convention is None:
        raise Error('We need a convention to know how to interpret the input!')

    kind = data_kind(spec)
    if kind == 'file':
        pass  # Use file as-is (dummy file context)...
    elif kind == 'matrix':
        pass  # Convert to virtual file and then use as-is...
    else:
        raise Error('File type not understood')

    # Build arg str and call meca...

[liamtoney]

I also wonder if we could replace the conventions that require mantissa and/or exponent with just the values themselves. It seems odd to specify them separately, perhaps a consideration for CLI GMT only that we needn't worry about in Python. I.e. for the GCMT convention we could replace those two with just moment to make things nicer:

moment = 5000

mantissa = 10  # Hard-coded
exponent = np.log10(moment)

mantissa ** exponent  # == 4999.999999999999

[liamtoney]

@tjnewton thanks for cleaning up my bad test. 😉 Just one small tweak with the file test — we can use GMTTempFile for this

from pygmt.helpers import GMTTempFile

Thanks!

[liamtoney]

I think this is a really good application for the GMTTempFile, docstring w/ usage example below:

pygmt/pygmt/helpers/tempfile.py

Line 27 in 2884fd0

class GMTTempFile:

[liamtoney]

I think the usage example would be:

from pygmt.helpers import GMTTempFile
with GMTTempFile() as temp:
    with open(temp.name, mode="w") as temp_file:
        temp_file.write(...)
    some_gmt_module(temp.name, ...)

Then you don't need to use the os.remove() since we have the context manager. LMK if that has issues!

[tjnewton]

Oooo thanks Liam, that's clever! Testing this test is a little slow on my end because make test on my machine freezes on ../pygmt/tests/test_grdview.py::test_grdview_grid_dataarray every time starting today -_-

[tjnewton]

Finally passed after a restart but the grdview tests are very slow now. Looks like test_meca is passing now.

[weiji14]

Sorry about the slow test (we'll need to streamline this somehow at some point). A workaround is to only run the meca tests locally. Use pytest -v --mpl pygmt/tests/test_meca.py.

[seisman]

Including the as temp that Liam mentioned in the second line will throw a pylint warning for an unused variable.

After using with GMTTempFile() as temp:, you should also remove the double quotes around "temp.name".

[tjnewton]

I see, so the temp.file object needed to be passed to meca rather than a temporary file name string. Should be fixed now.

[seisman]

Yes, I think this PR is good to merge. Thanks for your contribution!

[liamtoney]

temp.name is a string with the system filename.

[tjnewton]

thanks @liamtoney I apparently need some rest

[seisman]

@tjnewton Again, thanks for your contributions to PyGMT. Feel free to open another PR to add yourself to the AUTHORS list after this PR is merged.

[tjnewton]

@tjnewton Again, thanks for your contributions to PyGMT. Feel free to open another PR to add yourself to the AUTHORS list after this PR is merged.

Will do, thanks! I plan to continue contributing to the project.

[weiji14]

Thanks @tjnewton, really appreciate the work you've done here! 😄

@liamtoney, could you squash and merge in this PR since you've worked on it too? Will need to shorten the commit message to capture the gist of what was achieved in this PR.

[welcome]

🎉🎉🎉 Congrats on merging your first pull request and welcome to the team! 🎉🎉🎉

Please open a new pull request to add yourself to the AUTHORS.md file. We hope that this was a good experience for you. Let us know if there is any way that the contributing process could be improved.

[liamtoney]

Congrats @tjnewton for your first contribution to PyGMT! Definitely let us know how we could improve this process and make it more accessible. I think a lot of seismologists will be excited about this implementation!