Wrap coupe

[konalindsey]

Description of the desired feature

I am mapping cross sections of subduction zones, and need to plot focal mechanisms. I tried using Session.call_module, but ran into various errors... I probably am not using it correctly, but a real pyGMT function of pscoupe would be preferable.

Are you willing to help implement and maintain this feature? Yes/No

[welcome[bot]]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. You might also want to take a look at our contributing guidelines and code of conduct.

[yvonnefroehlich]

Thanks @konalindsey for trying out PyGMT :smiley: !

Maybe the method pygmt.Figure.meca is helpful for you?

There is also a gallery example for this method: https://www.pygmt.org/latest/gallery/seismology/meca.html#sphx-glr-gallery-seismology-meca-py

Please note: This method was fully rewritten after release v0.7.0 (see https://github.com/GenericMappingTools/pygmt/issues/2016#issuecomment-1189135191). Despite this, there are still some known issues and limitations with this method (e.g. https://github.com/GenericMappingTools/pygmt/issues/2016).

---

Edit Maybe I understood you / the question wrong, and you want to plot a cross section of a focal mechanisms, as the title of this issue refers to pscoupe module of GMT (https://docs.generic-mapping-tools.org/latest/supplements/seis/pscoupe.html)? Sorry for that!

[konalindsey]

Hi @yvonnefroehlich. Yes, meca works great! But I need to make cross sections of the focal mechanisms which would correspond to pscoupe in GMT.

[yvonnefroehlich]

I have no experience with pscoupe. So, probably the PyGMT maintainers can help you better and faster :wink:, regarding both wrapping pscoupe and using pscoupe in PyGMT via Session.call_module.

I tried using Session.call_module, but ran into various errors... I probably am not using it correctly,

I think it could be helpful if you post your code or a minimal example showing your issues when using pscoupe via Session.call_module.

[yvonnefroehlich]

Ping at @seisman, @weiji14, @maxrjones, @willschlitzer, and @michaelgrund for help on this, please.

[seisman]

I don't think I have time to wrap the coupe module recently, but here is a minimal script showing how to call a module that is not wrapped in PyGMT yet.

The example data is available at https://raw.githubusercontent.com/GenericMappingTools/gmt-for-geodesy/main/5_seismology/meca.dat

import pygmt

spec = "meca.dat"
fig = pygmt.Figure()
with pygmt.clib.Session() as lib:
    file_context = lib.virtualfile_from_data(check_kind="vector", data=spec)
    with file_context as fname:
        lib.call_module(module="coupe", args=f"{fname} -Sa1c -Aa110/33/120/33+r -JX15c/-5c -Baf")
fig.show()

[yvonnefroehlich]

Thanks @seisman for your answer and especially for the code example! 😃

[willschlitzer]

@GenericMappingTools/pygmt-maintainers Is there a remote file that would be applicable as a table for this? The example in the GMT docs does not point to an input file.

[maxrjones]

I'm not sure how easy it would be to create a simpler example from this, but @GCMT_1976-2017_meca.gmt contains focal mechanism data and is used in https://docs.generic-mapping-tools.org/latest/animations/anim14.html

[willschlitzer]

Thanks! I'm mostly looking for a remote file so I can write a test; I'll give that a go at some point!

[weiji14]

Another request from the forum at https://forum.generic-mapping-tools.org/t/how-to-plot-focal-mechanisme-cross-section-in-pygmt/3209. Looks like coupe is very in demand! Is anyone interested in implementing this in PyGMT?

We've got some draft instructions on how to wrap a new module at https://github.com/GenericMappingTools/pygmt/pull/1687, specifically at https://github.com/GenericMappingTools/pygmt/blob/wrap-module-instructions/doc/contributing.md#initial-feature-implementation. Let us know if anyone is interested in helping out with this and the team will be happy to provide guidance.

[willschlitzer]

I'm planning on trying to wrap it in September/October (my computer is still enroute from the U.K.); happy to help if anyone wants to give it a go!

[willschlitzer]

I was unsuccessful in getting coupe to work. I pushed my edits so far to the branch wrap/coupe if anyone wants to take a look at it/wrap it and submit a PR. I'm giving up on it for the time being.

[Josesx506]

fig = pygmt.Figure()
# Use the session helper to import the 'clip' module from GMT
with pygmt.helpers.GMTTempFile() as temp_file:
    with open(temp_file.name, 'w') as f:
        for i in range(len(mts_array)):
            f.write(f'{mts_array[i,0]} {mts_array[i,1]} {mts_array[i,2]} {mts_array[i,3]} {mts_array[i,4]} {mts_array[i,5]} {mts_array[i,6]}\n')
    with pygmt.clib.Session() as session:
        session.call_module('coupe', f'{temp_file.name} -R0/5/-0.4/3.2 -JX8.72c/-3c -Sa1.5 -Aa-90.82/30.174/-90.76/30.174+d90+w0.7/0/0 -Gred -Baf255')
fig.show()

The format of mts_array is longitude latitude depth strike dip rake Mw Each row is a new event. e.g for 20 events, the array shape is (20,7)

The 2 coordinates in the -Aa argument are the start and end long-lat locations of the profile, w is the width of profile in kilometers and d90 makes sure the moment tensors are viewed vertically from the side

[weiji14]

Found this blog post (in Japanese) from @yasuit21 who implemented some code to wrap coupe (see also code at https://gist.github.com/yasuit21/db1d676d784769b80365a3528d1ab547). Might be a good reference for someone who wants to bring this into PyGMT.

[DianKusumawati]

Found this blog post (in Japanese) from @yasuit21 who implemented some code to wrap coupe (see also code at https://gist.github.com/yasuit21/db1d676d784769b80365a3528d1ab547). Might be a good reference for someone who wants to bring this into PyGMT.

Thank you for the guidance. I have checked on the website and tried to implement it (the coupe.py). However, I encountered a problem related to:

decorators.py module_func.doc = docstring.format(**filler_text) KeyError: 'projection'

I tried to trace the script. As I checked on this website (https://github.com/GenericMappingTools/pygmt/blob/wrap-module-instructions/doc/contributing.md#add-missing-aliases), I tried to guess that the issue is possibly related to missing aliases in my case (?) (please correct me if I'm wrong). However, I am still puzzled about how to solve the problem.

May I ask for help regarding this problem?

Below is the screenshot of the error. The script that I use is: Crosssec_tes.py and coupe.py from https://gist.github.com/yasuit21/db1d676d784769b80365a3528d1ab547 meca.py from https://github.com/GenericMappingTools/pygmt/blob/v0.10.0/pygmt/src/meca.py

Your guidance will be very valuable to me. Thank you.

[seisman]

It seems the error comes from meca_edit.py. Please attach the script so that we can see what's inside.

[DianKusumawati]

It seems the error comes from meca_edit.py. Please attach the script so that we can see what's inside.

for the meca_edit.py script, I copied it entirely from the script in https://github.com/GenericMappingTools/pygmt/blob/v0.10.0/pygmt/src/meca.py

I renamed it to meca_edit.py because there is already a meca.py script in the Pygmt installation on my machine, and I noticed there is a slight difference with the meca script posted in Git Hub.

Below is the attachment of the script (meca_edit.py as well as with the other script that I use) in a zip file.

Thank you

script.zip

[seisman]

The script works well for me after some necessary minor changes, and I can't reproduce your issue.

[DianKusumawati]

The script works well for me after some necessary minor changes, and I can't reproduce your issue.

After I checked on the PyGMT version requirement (python=3.10.12 pygmt=0.10.0 numpy=1.25.2 pandas=2.0.3) on the website from @yasuit21, I was able to produce the figure in the example script. My previous PyGMT version was 0.5.0.

Kindly apologize for the missed version.

Thank you.