diff --git a/doc/api/index.rst b/doc/api/index.rst index 53334b04384..c4f199b3dbc 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -35,6 +35,7 @@ Plotting data and laying out the map: Figure.image Figure.shift_origin Figure.text + Figure.meca Color palette table generation: diff --git a/examples/gallery/plot/meca.py b/examples/gallery/plot/meca.py new file mode 100644 index 00000000000..17631a4f51d --- /dev/null +++ b/examples/gallery/plot/meca.py @@ -0,0 +1,33 @@ +""" +Focal mechanisms +---------------- + +The :meth:`pygmt.Figure.meca` method can plot focal mechanisms, or beachballs. +We can specify the focal mechanism nodal planes or moment tensor components as +a dict using the ``spec`` argument (or they can be specified as a 1d or 2d array, +or within a specified file). The size of plotted beachballs can be specified +using the ``scale`` argument. +""" + +import pygmt + +fig = pygmt.Figure() + +# generate a basemap near Washington state showing coastlines, land, and water +fig.coast( + region=[-125, -122, 47, 49], + projection="M6c", + land="grey", + water="lightblue", + shorelines=True, + resolution="f", + frame="a", +) + +# store focal mechanisms parameters in a dict +focal_mechanism = dict(strike=330, dip=30, rake=90, magnitude=3) + +# pass the focal mechanism data to meca in addition to the scale and event location +fig.meca(focal_mechanism, scale="1c", longitude=-124.3, latitude=48.1, depth=12.0) + +fig.show() diff --git a/pygmt/base_plotting.py b/pygmt/base_plotting.py index a9ca7dd26b8..f978e9c9cf5 100644 --- a/pygmt/base_plotting.py +++ b/pygmt/base_plotting.py @@ -8,7 +8,7 @@ import pandas as pd from .clib import Session -from .exceptions import GMTInvalidInput +from .exceptions import GMTError, GMTInvalidInput from .helpers import ( build_arg_string, dummy_context, @@ -994,3 +994,384 @@ def text( arg_str = " ".join([fname, build_arg_string(kwargs)]) lib.call_module("text", arg_str) + + @fmt_docstring + @use_alias(R="region", J="projection", B="frame", C="offset") + @kwargs_to_strings(R="sequence",) + def meca( + self, + spec, + scale, + longitude=None, + latitude=None, + depth=None, + convention=None, + component="full", + plot_longitude=None, + plot_latitude=None, + **kwargs, + ): + """ + Plot focal mechanisms. + + Full option list at :gmt-docs:`supplements/seis/meca.html` + + Note + ---- + Currently, labeling of beachballs with text strings is only + supported via providing a file to `spec` as input. + + {aliases} + + Parameters + ---------- + spec: dict, 1D array, 2D array, pd.DataFrame, or str + Either a filename containing focal mechanism parameters as columns, + a 1- or 2-D array with the same, or a dictionary. If a filename or + array, `convention` is required so we know how to interpret the + columns/entries. If a dictionary, the following combinations of + keys are supported; these determine the convention. Dictionary + may contain values for a single focal mechanism or lists of + values for many focal mechanisms. A Pandas DataFrame may + optionally contain columns latitude, longitude, depth, + plot_longitude, + and/or plot_latitude instead of passing them to the meca method. + + - ``"aki"`` — *strike, dip, rake, magnitude* + - ``"gcmt"`` — *strike1, dip1, rake1, strike2, dip2, rake2, + mantissa, exponent* + - ``"mt"`` — *mrr, mtt, mff, mrt, mrf, mtf, exponent* + - ``"partial"`` — *strike1, dip1, strike2, fault_type, magnitude* + - ``"principal_axis"`` — *t_exponent, t_azimuth, t_plunge, + n_exponent, n_azimuth, n_plunge, p_exponent, p_azimuth, p_plunge, + exponent* + + scale: str + Adjusts the scaling of the radius of the beachball, which is + proportional to the magnitude. Scale defines the size for + magnitude = 5 (i.e. scalar seismic moment M0 = 4.0E23 dynes-cm) + longitude: int, float, list, or 1d numpy array + Longitude(s) of event location. Ignored if `spec` is not a + dictionary. List must be the length of the number of events. + Ignored if `spec` is a DataFrame and contains a 'longitude' column. + latitude: int, float, list, or 1d numpy array + Latitude(s) of event location. Ignored if `spec` is not a + dictionary. List must be the length of the number of events. + Ignored if `spec` is a DataFrame and contains a 'latitude' column. + depth: int, float, list, or 1d numpy array + Depth(s) of event location in kilometers. Ignored if `spec` is + not a dictionary. List must be the length of the number of events. + Ignored if `spec` is a DataFrame and contains a 'depth' column. + convention: str + ``"aki"`` (Aki & Richards), ``"gcmt"`` (global CMT), ``"mt"`` + (seismic moment tensor), ``"partial"`` (partial focal mechanism), + or ``"principal_axis"`` (principal axis). Ignored if `spec` is a + dictionary or dataframe. + component: str + The component of the seismic moment tensor to plot. ``"full"`` (the + full seismic moment tensor), ``"dc"`` (the closest double couple + with zero trace and zero determinant), ``"deviatoric"`` (zero + trace) + plot_longitude: int, float, list, or 1d numpy array + Longitude(s) at which to place beachball, only used if `spec` is a + dictionary. List must be the length of the number of events. + Ignored if `spec` is a DataFrame and contains a 'plot_longitude' + column. + plot_latitude: int, float, list, or 1d numpy array + Latitude(s) at which to place beachball, only used if `spec` is a + dictionary. List must be the length of the number of events. + Ignored if `spec` is a DataFrame and contains a 'plot_latitude' + column. + offset: bool or str + Offsets beachballs to the longitude, latitude specified in + the last two columns of the input file or array, + or by `plot_longitude` and `plot_latitude` if provided. A small + circle is plotted at the initial location and a line connects + the beachball to the circle. Specify pen and optionally append + ``+ssize`` to change the line style and/or size of the circle. + {J} + {R} + {B} + """ + + # pylint warnings that need to be fixed + # pylint: disable=too-many-locals + # pylint: disable=too-many-nested-blocks + # pylint: disable=too-many-branches + # pylint: disable=no-self-use + # pylint: disable=too-many-statements + + def set_pointer(data_pointers, spec): + """Set optional parameter pointers based on DataFrame or dict, if + those parameters are present in the DataFrame or dict.""" + for param in list(data_pointers.keys()): + if param in spec: + # set pointer based on param name + data_pointers[param] = spec[param] + + def update_pointers(data_pointers): + """Updates variables based on the location of data, as the + following data can be passed as parameters or it can be + contained in `spec`.""" + # update all pointers + longitude = data_pointers["longitude"] + latitude = data_pointers["latitude"] + depth = data_pointers["depth"] + plot_longitude = data_pointers["plot_longitude"] + plot_latitude = data_pointers["plot_latitude"] + return (longitude, latitude, depth, plot_longitude, plot_latitude) + + # Check the spec and parse the data according to the specified + # convention + if isinstance(spec, (dict, pd.DataFrame)): + # dicts and DataFrames are handed similarly but not identically + if ( + longitude is None or latitude is None or depth is None + ) and not isinstance(spec, (dict, pd.DataFrame)): + raise GMTError("Location not fully specified.") + + param_conventions = { + "AKI": ["strike", "dip", "rake", "magnitude"], + "GCMT": ["strike1", "dip1", "dip2", "rake2", "mantissa", "exponent"], + "MT": ["mrr", "mtt", "mff", "mrt", "mrf", "mtf", "exponent"], + "PARTIAL": ["strike1", "dip1", "strike2", "fault_type", "magnitude"], + "PRINCIPAL_AXIS": [ + "t_exponent", + "t_azimuth", + "t_plunge", + "n_exponent", + "n_azimuth", + "n_plunge", + "p_exponent", + "p_azimuth", + "p_plunge", + "exponent", + ], + } + + # to keep track of where optional parameters exist + data_pointers = { + "longitude": longitude, + "latitude": latitude, + "depth": depth, + "plot_longitude": plot_longitude, + "plot_latitude": plot_latitude, + } + + # make a DataFrame copy to check convention if it contains + # other parameters + if isinstance(spec, (dict, pd.DataFrame)): + # check if a copy is necessary + copy = False + drop_list = [] + for pointer in data_pointers: + if pointer in spec: + copy = True + drop_list.append(pointer) + if copy: + spec_conv = spec.copy() + # delete optional parameters from copy for convention check + for item in drop_list: + del spec_conv[item] + else: + spec_conv = spec + + # set convention and focal parameters based on spec convention + convention_assigned = False + for conv in param_conventions: + if set(spec_conv.keys()) == set(param_conventions[conv]): + convention = conv.lower() + foc_params = param_conventions[conv] + convention_assigned = True + break + if not convention_assigned: + raise GMTError( + "Parameters in spec dictionary do not match known " "conventions." + ) + + # create a dict type pointer for easier to read code + if isinstance(spec, dict): + dict_type_pointer = list(spec.values())[0] + elif isinstance(spec, pd.DataFrame): + # use df.values as pointer for DataFrame behavior + dict_type_pointer = spec.values + + # assemble the 1D array for the case of floats and ints as values + if isinstance(dict_type_pointer, (int, float)): + # update pointers + set_pointer(data_pointers, spec) + # look for optional parameters in the right place + ( + longitude, + latitude, + depth, + plot_longitude, + plot_latitude, + ) = update_pointers(data_pointers) + + # Construct the array (order matters) + spec = [longitude, latitude, depth] + [spec[key] for key in foc_params] + + # 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 "C" not in kwargs: + kwargs["C"] = True + spec.append(arg) + + # or assemble the 2D array for the case of lists as values + elif isinstance(dict_type_pointer, list): + # update pointers + set_pointer(data_pointers, spec) + # look for optional parameters in the right place + ( + longitude, + latitude, + depth, + plot_longitude, + plot_latitude, + ) = update_pointers(data_pointers) + + # before constructing the 2D array lets check that each key + # of the dict has the same quantity of values to avoid bugs + list_length = len(list(spec.values())[0]) + for value in list(spec.values()): + if len(value) != list_length: + raise GMTError( + "Unequal number of focal mechanism " + "parameters supplied in 'spec'." + ) + # lets also check the inputs for longitude, latitude, + # and depth if it is a list or array + if ( + isinstance(longitude, (list, np.ndarray)) + or isinstance(latitude, (list, np.ndarray)) + or isinstance(depth, (list, np.ndarray)) + ): + if (len(longitude) != len(latitude)) or ( + len(longitude) != len(depth) + ): + raise GMTError( + "Unequal number of focal mechanism " + "locations supplied." + ) + + # values are ok, so build the 2D array + spec_array = [] + for index in range(list_length): + # Construct the array one row at a time (note that order + # matters here, hence the list comprehension!) + row = [longitude[index], latitude[index], depth[index]] + [ + spec[key][index] for key in foc_params + ] + + # Add in plotting options, if given, otherwise add 0s as + # required by GMT + for arg in plot_longitude, plot_latitude: + if arg is None: + row.append(0) + else: + if "C" not in kwargs: + kwargs["C"] = True + row.append(arg[index]) + spec_array.append(row) + spec = spec_array + + # or assemble the array for the case of pd.DataFrames + elif isinstance(dict_type_pointer, np.ndarray): + # update pointers + set_pointer(data_pointers, spec) + # look for optional parameters in the right place + ( + longitude, + latitude, + depth, + plot_longitude, + plot_latitude, + ) = update_pointers(data_pointers) + + # lets also check the inputs for longitude, latitude, and depth + # just in case the user entered different length lists + if ( + isinstance(longitude, (list, np.ndarray)) + or isinstance(latitude, (list, np.ndarray)) + or isinstance(depth, (list, np.ndarray)) + ): + if (len(longitude) != len(latitude)) or ( + len(longitude) != len(depth) + ): + raise GMTError( + "Unequal number of focal mechanism locations supplied." + ) + + # values are ok, so build the 2D array in the correct order + spec_array = [] + for index in range(len(spec)): + # Construct the array one row at a time (note that order + # matters here, hence the list comprehension!) + row = [longitude[index], latitude[index], depth[index]] + [ + spec[key][index] for key in foc_params + ] + + # Add in plotting options, if given, otherwise add 0s as + # required by GMT + for arg in plot_longitude, plot_latitude: + if arg is None: + row.append(0) + else: + if "C" not in kwargs: + kwargs["C"] = True + row.append(arg[index]) + spec_array.append(row) + spec = spec_array + + else: + raise GMTError( + "Parameter 'spec' contains values of an unsupported type." + ) + + # Add condition and scale to kwargs + if convention == "aki": + data_format = "a" + elif convention == "gcmt": + data_format = "c" + elif convention == "mt": + # Check which component of mechanism the user wants plotted + if component == "deviatoric": + data_format = "z" + elif component == "dc": + data_format = "d" + else: # component == 'full' + data_format = "m" + elif convention == "partial": + data_format = "p" + elif convention == "principal_axis": + # Check which component of mechanism the user wants plotted + if component == "deviatoric": + data_format = "t" + elif component == "dc": + data_format = "y" + else: # component == 'full' + data_format = "x" + # Support old-school GMT format options + elif convention in ["a", "c", "m", "d", "z", "p", "x", "y", "t"]: + data_format = convention + else: + raise GMTError("Convention not recognized.") + + # Assemble -S flag + kwargs["S"] = data_format + scale + + kind = data_kind(spec) + with Session() as lib: + if kind == "matrix": + file_context = lib.virtualfile_from_matrix(np.atleast_2d(spec)) + elif kind == "file": + file_context = dummy_context(spec) + else: + raise GMTInvalidInput("Unrecognized data type: {}".format(type(spec))) + with file_context as fname: + arg_str = " ".join([fname, build_arg_string(kwargs)]) + lib.call_module("meca", arg_str) diff --git a/pygmt/tests/baseline/test_meca_loc_array.png b/pygmt/tests/baseline/test_meca_loc_array.png new file mode 100644 index 00000000000..73c9675c1d7 Binary files /dev/null and b/pygmt/tests/baseline/test_meca_loc_array.png differ diff --git a/pygmt/tests/baseline/test_meca_spec_1d_array.png b/pygmt/tests/baseline/test_meca_spec_1d_array.png new file mode 100644 index 00000000000..e6c1f0c25c2 Binary files /dev/null and b/pygmt/tests/baseline/test_meca_spec_1d_array.png differ diff --git a/pygmt/tests/baseline/test_meca_spec_2d_array.png b/pygmt/tests/baseline/test_meca_spec_2d_array.png new file mode 100644 index 00000000000..51c62146df3 Binary files /dev/null and b/pygmt/tests/baseline/test_meca_spec_2d_array.png differ diff --git a/pygmt/tests/baseline/test_meca_spec_dataframe.png b/pygmt/tests/baseline/test_meca_spec_dataframe.png new file mode 100644 index 00000000000..45a2efda1b4 Binary files /dev/null and b/pygmt/tests/baseline/test_meca_spec_dataframe.png differ diff --git a/pygmt/tests/baseline/test_meca_spec_dict_list.png b/pygmt/tests/baseline/test_meca_spec_dict_list.png new file mode 100644 index 00000000000..c4b54723c7b Binary files /dev/null and b/pygmt/tests/baseline/test_meca_spec_dict_list.png differ diff --git a/pygmt/tests/baseline/test_meca_spec_dictionary.png b/pygmt/tests/baseline/test_meca_spec_dictionary.png new file mode 100644 index 00000000000..9663b3a9c1d Binary files /dev/null and b/pygmt/tests/baseline/test_meca_spec_dictionary.png differ diff --git a/pygmt/tests/baseline/test_meca_spec_file.png b/pygmt/tests/baseline/test_meca_spec_file.png new file mode 100644 index 00000000000..d3a57997080 Binary files /dev/null and b/pygmt/tests/baseline/test_meca_spec_file.png differ diff --git a/pygmt/tests/test_meca.py b/pygmt/tests/test_meca.py new file mode 100644 index 00000000000..4776be1e5b1 --- /dev/null +++ b/pygmt/tests/test_meca.py @@ -0,0 +1,239 @@ +""" +Tests for meca +""" +import os +import pandas as pd +import numpy as np +import pytest +from pygmt.helpers import GMTTempFile + +from .. import Figure + + +TEST_DATA_DIR = os.path.join(os.path.dirname(__file__), "data") + + +@pytest.mark.mpl_image_compare +def test_meca_spec_dictionary(): + """ + Test supplying a dictionary containing a single focal mechanism to the + `spec` argument. + """ + + fig = Figure() + + # Right lateral strike slip focal mechanism + fig.meca( + dict(strike=0, dip=90, rake=0, magnitude=5), + longitude=0, + latitude=5, + depth=0, + scale="2.5c", + region=[-1, 1, 4, 6], + projection="M14c", + frame=2, + ) + + return fig + + +@pytest.mark.mpl_image_compare +def test_meca_spec_dict_list(): + """ + Test supplying a dictionary containing a list of focal mechanism to the + `spec` argument. + """ + + fig = Figure() + + # supply focal mechanisms as a dict of lists + focal_mechanisms = dict( + strike=[330, 350], dip=[30, 50], rake=[90, 90], magnitude=[3, 2] + ) + + fig.meca( + focal_mechanisms, + longitude=[-124.3, -124.4], + latitude=[48.1, 48.2], + depth=[12.0, 11.0], + region=[-125, -122, 47, 49], + scale="2c", + projection="M14c", + ) + + return fig + + +@pytest.mark.mpl_image_compare +def test_meca_spec_dataframe(): + """ + Test supplying a pandas DataFrame containing focal mechanisms and + locations to the `spec` argument. + """ + + fig = Figure() + + # supply focal mechanisms to meca as a dataframe + focal_mechanisms = dict( + strike=[324, 353], + dip=[20.6, 40], + rake=[83, 90], + magnitude=[3.4, 2.9], + longitude=[-124, -124.4], + latitude=[48.1, 48.2], + depth=[12, 11.0], + ) + spec_dataframe = pd.DataFrame(data=focal_mechanisms) + + fig.meca(spec_dataframe, region=[-125, -122, 47, 49], scale="2c", projection="M14c") + + return fig + + +@pytest.mark.mpl_image_compare +def test_meca_spec_1d_array(): + """ + Test supplying a 1D numpy array containing focal mechanisms and + locations to the `spec` argument. + """ + + fig = Figure() + + # supply focal mechanisms to meca as a 1D numpy array, here we are using + # the Harvard CMT zero trace convention but the focal mechanism + # parameters may be specified any of the available conventions. Since we + # are not using a dict or dataframe the convention and component should + # be specified. + focal_mechanism = [ + -127.40, # longitude + 40.87, # latitude + 12, # depth + -3.19, # mrr + 0.16, # mtt + 3.03, # mff + -1.02, # mrt + -3.93, # mrf + -0.02, # mtf + 23, # exponent + 0, # plot_lon, 0 to plot at event location + 0, # plot_lat, 0 to plot at event location + ] + focal_mech_array = np.asarray(focal_mechanism) + + fig.meca( + focal_mech_array, + convention="mt", + component="full", + region=[-128, -127, 40, 41], + scale="2c", + projection="M14c", + ) + + return fig + + +@pytest.mark.mpl_image_compare +def test_meca_spec_2d_array(): + """ + Test supplying a 2D numpy array containing focal mechanisms and + locations to the `spec` argument. + """ + + fig = Figure() + + # supply focal mechanisms to meca as a 2D numpy array, here we are using + # the GCMT convention but the focal mechanism parameters may be + # specified any of the available conventions. Since we are not using a + # dict or dataframe the convention and component should be specified. + focal_mechanisms = [ + [ + -127.40, # longitude + 40.87, # latitude + 12, # depth + 170, # strike1 + 20, # dip1 + -110, # rake1 + 11, # strike2 + 71, # dip2 + -83, # rake2 + 5.1, # mantissa + 23, # exponent + 0, # plot_lon, 0 means we want to plot at the event location + 0, # plot_lat + ], + [-127.50, 40.88, 12.0, 168, 40, -115, 20, 54, -70, 4.0, 23, 0, 0], + ] + focal_mechs_array = np.asarray(focal_mechanisms) + + fig.meca( + focal_mechs_array, + convention="gcmt", + region=[-128, -127, 40, 41], + scale="2c", + projection="M14c", + ) + + return fig + + +@pytest.mark.mpl_image_compare +def test_meca_spec_file(): + """ + Test supplying a file containing focal mechanisms and locations to the + `spec` argument. + """ + + fig = Figure() + + focal_mechanism = [-127.43, 40.81, 12, -3.19, 1.16, 3.93, -1.02, -3.93, -1.02, 23] + + # writes temp file to pass to gmt + with GMTTempFile() as temp: + with open(temp.name, mode="w") as temp_file: + temp_file.write(" ".join([str(x) for x in focal_mechanism])) + # supply focal mechanisms to meca as a file + fig.meca( + temp.name, + convention="mt", + component="full", + region=[-128, -127, 40, 41], + scale="2c", + projection="M14c", + ) + + return fig + + +@pytest.mark.mpl_image_compare +def test_meca_loc_array(): + """ + Test supplying lists and np.ndarrays as the event location (longitude, + latitude, and depth). + """ + + fig = Figure() + + # specify focal mechanisms + focal_mechanisms = dict( + strike=[327, 350], dip=[41, 50], rake=[68, 90], magnitude=[3, 2] + ) + + # longitude, latitude, and depth may be specified as an int, float, + # list, or 1d numpy array + longitude = np.array([-123.3, -124.4]) + latitude = np.array([48.4, 48.2]) + depth = [12.0, 11.0] # to test mixed data types as inputs + + scale = "2c" + + fig.meca( + focal_mechanisms, + scale, + longitude, + latitude, + depth, + region=[-125, -122, 47, 49], + projection="M14c", + ) + + return fig