pygmt.Figure.meca offset not working

[rdcorona]

I am using pygmt to graph focal mechanisms with magnitude scale.
However, when I try to give it an offset so that they don't pile up the function doesn't seem to work.

I'm using :

spec = dict (strike=strikes[0], dip=dips[0], rake=rakes[0], magnitude=magnitude [0])

fig.meca (spec, scale=f ".{0.002 * (2 ** magnitude[0])}c", longitude = lons [0], latitude = lats [0], depth = 25.0, convention = "aki", plot_longitude = lons [0], plot_latitude = lats[0] -10, offset=True) 

where all the objects are lists with real (float) values

System information

# Name                    Version                   Build  Channel
pygmt                     0.3.1              pyhd8ed1ab_0    conda-forge

Do I make something wrong?

[welcome]

👋 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.

[core-man]

Thanks, @rdcorona. I think you find a bug in the meca method. I can re-produce this potential bug using Focal mechanisms gallery example. The earthquake is at longitude=-124.3 and latitude=48.1 and the beachball should be plotted at plot_longitude=-124.3 and plot_latitude=47.5. However, the beachball is stilled plotted at the epicenter.

import pygmt
fig = pygmt.Figure()
fig.coast(region=[-125, -122, 47, 49], projection="M6c", land="grey", water="lightblue", shorelines=True, frame="a")
focal_mechanism = dict(strike=330, dip=30, rake=90, magnitude=3)
fig.meca(focal_mechanism, scale="1c", longitude=-124.3, latitude=48.1, depth=12.0, plot_longitude=-124.3, plot_latitude=47.5, offset=True)
fig.show()

[seisman]

I can reproduce the issue, and now know why it fails. Unfortunately, there is no workaround now.

[rdcorona]

Hello seisman.
So, we'll have to wait for the 0.4.0 release.🤷🏻‍♂️

Thanks for the reply

[seisman]

Currently, offsetting beachballs only works for files.

I'll see what I can do to fix it.

[rdcorona]

Ok, meanwhile I'll try it with a dataFrame and an ASCII file

Thanks again

[weiji14]

I can reproduce the issue, and now know why it fails. Unfortunately, there is no workaround now.

Still a problem on PyGMT v0.5.0 unfortunately. Could you elaborate more on why it fails? I tried debugging this while handling #1611, and couldn't see any obvious issue on the Python side, e.g. the plot_longitude and plot_latitude variables seem to be appended to spec ok here:

pygmt/pygmt/src/meca.py

Lines 330 to 337 in 9e8a7bd

	# Add in plotting options, if given, otherwise add 0s
	for arg in plot_longitude, plot_latitude:
	if arg is None:
	spec.append(0)
	else:
	if "A" not in kwargs:
	kwargs["A"] = True
	spec.append(arg)

Only thing I can think of is that this is a problem with the spec array to GMT virtualfile convention at

pygmt/pygmt/src/meca.py

Line 454 in 9e8a7bd

file_context = lib.virtualfile_from_matrix(np.atleast_2d(spec))

so should this be marked as an upstream GMT API bug then?

[seisman]

Ah I see, maybe we could use lib.virtualfile_from_data as with #949?

Yes, ideally we should use lib.virtualfile_from_data.

Just to be clear, do you mean that plot_longitude, plot_latitude and event_name is to be passed in as 3 columns or 1 column using GMT_Put_Strings?

We can only call GMT_Put_Strings once, so we have to combine the three parameters into one string vector and pass it to GMT_Put_Strings. Actually, we already did it in the virtualfile_from_vector function. The virtualfile_from_vector function doesn't work for meca, because plot_longitude and plot_latitude are numeric vectors but must be passed as strings.

pygmt/pygmt/clib/session.py

Lines 1159 to 1172 in 4b80f3c

	# Use put_strings for last column(s) with string type data
	# Have to use modifier "GMT_IS_DUPLICATE" to duplicate the strings
	string_arrays = arrays[columns:]
	if string_arrays:
	if len(string_arrays) == 1:
	strings = string_arrays[0]
	elif len(string_arrays) > 1:
	strings = np.apply_along_axis(
	func1d=" ".join, axis=0, arr=string_arrays
	)
	strings = np.asanyarray(a=strings, dtype=str)
	self.put_strings(
	dataset, family="GMT_IS_VECTOR|GMT_IS_DUPLICATE", strings=strings
	)

.

[seisman]

I don't know what happened to GitHub.

My reply to your comment #1129 (comment) is displayed before your comment at here: #1129 (comment)

[weiji14]

I don't know what happened to GitHub.

My reply to your comment #1129 (comment) is displayed before your comment at here: #1129 (comment)

Yeah I see the same thing, must be time going backwards 😆

Ah I see, maybe we could use lib.virtualfile_from_data as with #949?

Yes, ideally we should use lib.virtualfile_from_data.

Just to be clear, do you mean that plot_longitude, plot_latitude and event_name is to be passed in as 3 columns or 1 column using GMT_Put_Strings?

We can only call GMT_Put_Strings once, so we have to combine the three parameters into one string vector and pass it to GMT_Put_Strings. Actually, we already did it in the virtualfile_from_vector function. The virtualfile_from_vector function doesn't work for meca, because plot_longitude and plot_latitude are numeric vectors but must be passed as strings.

Ok, I think I got a clear idea on what to do now. We'll need to convert the plot_longitude and plot_latitude vectors from numerical arrays to string arrays, and then lib.virtualfile_from_vectors (via lib.virtualfile_from_data) should do the concatenation magic (joining event_name if it is present) and call GMT_Put_Strings only once on the concatenated string column.

Probably do 2 separate PRs:

1. Refactor from using lib.virtualfile_from_matrix to lib.virtualfile_from_data (which will enable other input data structures like xarray.Dataset arrays) and possibly allow for text labels to work on non-file inputs (fixing Wrap meca #516 (comment)).
2. Do the actual bugfix so that offset works.

[seisman]

Ok, I think I got a clear idea on what to do now. We'll need to convert the plot_longitude and plot_latitude vectors from numerical arrays to string arrays, and then lib.virtualfile_from_vectors (via lib.virtualfile_from_data) should do the concatenation magic (joining event_name if it is present) and call GMT_Put_Strings only once on the concatenated string column.

Probably do 2 separate PRs. 1) Refactor from using lib.virtualfile_from_matrix to lib.virtualfile_from_data (which will enable other input data structures like xarray.Dataset arrays) and 2) Do the actual bugfix so that offset works.

Yes!

[seisman]

The input parameters of the meca module are a little complicated. I will take the Aki & Richard convention (-Sa) as an example. The input parameters are:

longitude latitude depth strike dip rake magnitude [plot_longitude plot_latitude] [event_name]

As you can see, plot_longitude & plot_latitude and event_name are optional. Thus users can give input parameters in many different ways:

1. don't give plot longitude&latitude and event_name
2. give plot longitude&latitude but not event_name
3. give event_name but not plot longitude&latitude

Due to these complications, the first few columns (7 columns for the -Sa case) are read as numeric values (S->data in the GMT source code), thus data passed as vectors (virtualfile_from_vector) or a matrix (virtualfile_from_matrix) can be read correctly. However, because GMT can't know exactly the number of trailing columns and the event_name parameter can contain whitespace, all the trailing columns are read as trailing text (S->text in the GMT source code).

See https://github.com/GenericMappingTools/gmt/blob/3e68b8de02774f2eda0068ea356586d61770f070/src/seis/psmeca.c#L775 and https://github.com/GenericMappingTools/gmt/blob/3e68b8de02774f2eda0068ea356586d61770f070/src/seis/psmeca.c#L824 if you want to read the GMT source code.

Due to the above reasons, we have to pass the trailing columns (plot longitude&latitude and event_name) using GMT_Put_Strings, not GMT_Put_Matrix or GMT_Put_Vectors. Actually, we once had similar issues with the text function (related to #483, #520).

Thus, what we need to do is:

• Replace lib.virtualfile_from_matrix with lib.virtualfile_from_vector
• Build a string list/vector for plot_longitude plot_latitude event_name and pass it to lib.virtual_from_vector.

[weiji14]

Thus, what we need to do is:

• Replace lib.virtualfile_from_matrix with lib.virtualfile_from_vector

Ah I see, maybe we could use lib.virtualfile_from_data as with #949? Since it basically passes things to lib.virtualfile_from_vectors anyways:

pygmt/pygmt/clib/session.py

Lines 1427 to 1431 in 4b80f3c

	# Note: virtualfile_from_matrix is not used because a matrix can be
	# converted to vectors instead, and using vectors allows for better
	# handling of string type inputs (e.g. for datetime data types)
	"matrix": self.virtualfile_from_vectors,
	"vectors": self.virtualfile_from_vectors,

though it sounds like we'll need to be careful and ensure that it doesn't get into this case:

pygmt/pygmt/clib/session.py

Lines 1454 to 1458 in 4b80f3c

	# Just use virtualfile_from_matrix for 2D numpy.ndarray
	# which are signed integer (i), unsigned integer (u) or
	# floating point (f) types
	assert data.ndim == 2 and data.dtype.kind in "iuf"
	_virtualfile_from = self.virtualfile_from_matrix

so some very careful refactoring and unit test writing to do!

• Build a string list/vector for plot_longitude plot_latitude event_name and pass it to lib.virtual_from_vector.

Just to be clear, do you mean that plot_longitude, plot_latitude and event_name is to be passed in as 3 columns or 1 column using GMT_Put_Strings?

[weiji14]

Ok, I think I got a clear idea on what to do now. We'll need to convert the plot_longitude and plot_latitude vectors from numerical arrays to string arrays, and then lib.virtualfile_from_vectors (via lib.virtualfile_from_data) should do the concatenation magic (joining event_name if it is present) and call GMT_Put_Strings only once on the concatenated string column.
Probably do 2 separate PRs. 1) Refactor from using lib.virtualfile_from_matrix to lib.virtualfile_from_data (which will enable other input data structures like xarray.Dataset arrays) and 2) Do the actual bugfix so that offset works.

Yes!

On third thought, maybe 3 separate PRs (because pygmt.meca is a beast to refactor...):

1. Refactor meca to use virtualfile_from_data (doing at Refactor meca to use virtualfile_from_data #1613)
2. Handle text labels and angle/font/justify (c.f. Figure.meca: Allow specifying the label angle, font, justification, and offset #1466)
3. Fix this offset parameter bug (pygmt.Figure.meca offset not working #1129).

[weiji14]

Apparently the beachball offset bug wasn't properly fixed in #1784 for PyGMT v0.7.0, but there's a workaround at #2016 (comment) (use str types for plot_longitude and plot_latitude inputs), and we're looking at fixing things properly in #2202 for PyGMT v0.8.0.