From 77228de54c10ea288850605fd72841d220bb7aee Mon Sep 17 00:00:00 2001 From: thurinj Date: Fri, 1 Mar 2024 14:04:05 +1100 Subject: [PATCH 01/31] Added generic pygmt plotting backend for attributes Implemented a full "generic" pygmt spider plot, with dynamic topography, coastlines, moment tensor and header option. The backend comes with a utility class where most individual functions are stored. --- mtuq/graphics/__init__.py | 2 +- mtuq/graphics/attrs.py | 314 ++++++++++++++++++++++++++++++++++++-- 2 files changed, 302 insertions(+), 14 deletions(-) diff --git a/mtuq/graphics/__init__.py b/mtuq/graphics/__init__.py index badf8d0b8..35ae6ee0a 100644 --- a/mtuq/graphics/__init__.py +++ b/mtuq/graphics/__init__.py @@ -5,7 +5,7 @@ from mtuq.graphics.attrs import\ plot_time_shifts, plot_amplitude_ratios, plot_log_amplitude_ratios,\ - _plot_attrs + _plot_attrs, plot_cross_corr, _pygmt_backend from mtuq.graphics.beachball import\ plot_beachball, plot_polarities diff --git a/mtuq/graphics/attrs.py b/mtuq/graphics/attrs.py index 45a0d9c86..30878b2fc 100644 --- a/mtuq/graphics/attrs.py +++ b/mtuq/graphics/attrs.py @@ -7,6 +7,8 @@ from os.path import join from mtuq.util import defaults, warn +from mtuq.graphics._pygmt import exists_pygmt +from mtuq.event import MomentTensor def plot_time_shifts(dirname, attrs, stations, origin, key='total_shift', @@ -22,15 +24,15 @@ def plot_time_shifts(dirname, attrs, stations, origin, key='total_shift', MTUQ distinguishes between the following different types of time shifts - + - `static_shift` is an initial user-supplied time shift applied during - data processing + data processing - `time_shift` is a subsequent cross-correlation time shift applied - during misfit evaluation + during misfit evaluation - `total_shift` is the total correction, or in other words the sum of - static and cross-correlation time shifts + static and cross-correlation time shifts .. rubric :: Required input arguments @@ -87,7 +89,7 @@ def plot_amplitude_ratios(dirname, attrs, stations, origin, **kwargs): """ defaults(kwargs, { - 'colormap': 'Reds', + 'colormap': 'inferno', 'label': '$A_{obs}/A_{syn}$', 'zero_centered': False, }) @@ -127,8 +129,8 @@ def plot_log_amplitude_ratios(dirname, attrs, stations, origin, **kwargs): def _plot_attrs(dirname, stations, origin, attrs, key, - components=['Z', 'R', 'T'], format='png', backend=None, - **kwargs): + components=['Z', 'R', 'T'], format='png', backend=None, + **kwargs): """ Reads the attribute given by `key` from the `attrs` data structure, and plots how this attribute varies @@ -162,27 +164,28 @@ def _plot_attrs(dirname, stations, origin, attrs, key, if backend is None: backend = _default_backend + elif backend == _pygmt_backend and not exists_pygmt(): + warn('PyGMT backend requested but PyGMT not found'); backend = _default_backend if not callable(backend): raise TypeError - os.makedirs(dirname, exist_ok=True) for component in components: values = [] - station_list = [] + active_stations_list = [] for _i, station in enumerate(stations): if component not in attrs[_i]: continue values += [attrs[_i][component][key]] - station_list += [stations[_i]] + active_stations_list += [stations[_i]] if len(values) > 0: filename = join(dirname, component+'.'+format) - backend(filename, values, station_list, origin, **kwargs) + backend(filename, values, active_stations_list, origin, stations_list = stations, **kwargs) # @@ -191,7 +194,7 @@ def _plot_attrs(dirname, stations, origin, attrs, key, def _default_backend(filename, values, stations, origin, colormap='coolwarm', zero_centered=True, colorbar=True, - label='', width=5., height=5.): + label='', width=5., height=5., **kwargs): """ Default backend for all other `mtuq.graphics.attrs` functions @@ -230,7 +233,7 @@ def _default_backend(filename, values, stations, origin, else: min_val = np.min(values) max_val = np.max(values) - + # plot stations im = pyplot.scatter( [station.longitude for station in stations], @@ -284,3 +287,288 @@ def _default_backend(filename, values, stations, origin, pyplot.close() +def _pygmt_backend(filename, values, active_stations, origin, + colormap='polar', zero_centered=True, display_topo=True, + label='', width=5, moment_tensor=None, process=None, + stations_list=None, station_labels=True, min_val=None, max_val=None, **kwargs): + """ + PyGMT backend for plotting station attributes with hillshading using the + Miller Cylindrical projection, with an azimuth of 0/90 and a normalization + of t1 for the hillshade intensity. + """ + import pygmt + + if not stations_list: + stations_list = active_stations + print('Complete station list not passed to pygmt plotting backend \nWill plot only active stations') + # Collection of longitudes and latitudes from all available stations + longitudes = [s.longitude for s in stations_list + [origin]] + latitudes = [s.latitude for s in stations_list + [origin]] + + # Calculate the region to display with a buffer around the stations + region, lat_buffer = PyGMTUtilities.calculate_plotting_region(stations_list, origin, buffer_percentage=0.1) + + # Setting up the figure + fig = pygmt.Figure() + + # Dynamically determine the grid resolution for topography based on the range of longitudes and latitudes + # (etopo topography file will be downloaded if not found) + resolution = PyGMTUtilities.get_resolution(max(longitudes) - min(longitudes), max(latitudes) - min(latitudes)) + grid = pygmt.datasets.load_earth_relief(region=region, resolution=resolution) + + # Define a grayscale colormap for topography + pygmt.makecpt(cmap='gray', series=[-7000, 7000]) + + # Calculate the gradient (hillshade) grid with azimuth 0/300 and normalization t1 + # + shade = pygmt.grdgradient(grid=grid, azimuth="0/300", normalize="t1") + # Plot the hillshade grid as an image + if display_topo: + fig.grdimage(grid=grid, shading=shade, projection=f'J{width}i', frame='a', cmap='gray', no_clip=True) + + # Overlay coastlines + PyGMTUtilities.draw_coastlines(fig) + + # Configure the colormap for station values + colormap, cmap_reverse_flag = PyGMTUtilities.configure_colormap(colormap) + if zero_centered: + pygmt.makecpt(cmap=colormap, series=[-np.max(np.abs(values))*1.01, np.max(np.abs(values))*1.01], reverse=cmap_reverse_flag) + elif min_val is not None and max_val is not None: + pygmt.makecpt(cmap=colormap, series=[min_val, max_val], continuous=True, reverse=cmap_reverse_flag) + else: + pygmt.makecpt(cmap=colormap, series=[np.min(values), np.max(values)], continuous=True, reverse=cmap_reverse_flag) + + + # Plotting lines from origin to stations + for station in stations_list: + if station in active_stations: + # Plot line for active station as colored line + value = values[active_stations.index(station)] if station in active_stations else 0 + fig.plot( + x=[origin.longitude, station.longitude], + y=[origin.latitude, station.latitude], + cmap=True, + zvalue=value, + pen="thick,+z,-" + ) + + # Plotting stations as triangles + fig.plot( + x=[station.longitude for station in active_stations], + y=[station.latitude for station in active_stations], + style='i0.8c', # Triangle + color=values, + cmap=True, + pen="0.5p,black" + ) + + # Plotting non-active stations as hollow triangles + non_active_stations = [station for station in stations_list if station not in active_stations] + if len(non_active_stations) > 0: + fig.plot( + x=[station.longitude for station in non_active_stations], + y=[station.latitude for station in non_active_stations], + style='i0.8c', # Triangle + color=None, # Hollow (white) triangle + pen="0.5p,black" # Outline color + ) + fig.plot( + x=[station.longitude for station in non_active_stations], + y=[station.latitude for station in non_active_stations], + style='i0.6c', # Triangle + color=None, # Hollow (white) triangle + pen="0.5p,white" # Outline color + ) + + # Plotting the origin as a star + fig.plot( + x=[origin.longitude], + y=[origin.latitude], + style='a0.6c', # Star, size 0.5 cm + color='yellow', + pen="0.5p,black" + ) + + if moment_tensor is not None: + # Normalize the moment tensor components to the desired exponent + + if type(moment_tensor) is MomentTensor: + moment_tensor = moment_tensor.as_vector() + + moment_tensor = np.array(moment_tensor)/np.linalg.norm(moment_tensor) + + moment_tensor_spec = { + 'mrr': moment_tensor[0], + 'mtt': moment_tensor[1], + 'mff': moment_tensor[2], + 'mrt': moment_tensor[3], + 'mrf': moment_tensor[4], + 'mtf': moment_tensor[5], + 'exponent': 21 # Merely for size control, as the MT is normalized prior to plotting + } + + # Plot the moment tensor as a beachball + fig.meca( + spec=moment_tensor_spec, + scale="1c", # Sets a fixed size for the beachball plot + longitude=origin.longitude, + latitude=origin.latitude, + depth=10, # Depth is required, even if not used, set to a small number + convention="mt", # Use GMT's mt convention + compressionfill="red", + extensionfill="white", + pen="black" + ) + + if station_labels is True: + # Plotting station labels + for station in stations_list: + fig.text( + x=station.longitude, + y=station.latitude, + text=station.station, + font="5p,Helvetica-Bold,black", + justify="LM", + offset="-0.45c/0.125c", + fill='white' + ) + + fig.colorbar(frame=f'+l"{PyGMTUtilities.prepare_latex_annotations(label)}"', position="JMR+o1.5c/0c+w7c/0.5c") + + fig.basemap(region=region, projection=f'J{width}i', frame=True) + + # Now starts the header text above the plot -- It is not a title and can be modified. + # Add an integer increment to the text_line_val bellow to add a new line above. + text_line_val = 1 + header_lines = PyGMTUtilities.get_header(label, origin, filename, process) + + for header_line in header_lines: + fig.text(x=-148, y=(max(latitudes) + lat_buffer)+(text_line_val)*0.25, text=header_line, font="14p,Helvetica-Bold,black", justify="MC", no_clip=True) + text_line_val += 1 + + # Saving the figure + fig.savefig(filename, crop=True, dpi=300) + +class PyGMTUtilities: + @staticmethod + def calculate_plotting_region(stations, origin, buffer_percentage=0.1): + longitudes = [station.longitude for station in stations] + [origin.longitude] + latitudes = [station.latitude for station in stations] + [origin.latitude] + + lon_buffer = (max(longitudes) - min(longitudes)) * buffer_percentage + lat_buffer = (max(latitudes) - min(latitudes)) * buffer_percentage + + region = [min(longitudes) - lon_buffer, max(longitudes) + lon_buffer, + min(latitudes) - lat_buffer, max(latitudes) + lat_buffer] + return region, lat_buffer + + + @staticmethod + def get_resolution(lon_range, lat_range): + """ + Determines the resolution based on the given longitude and latitude ranges. + + Args: + lon_range (float): The range of longitudes. + lat_range (float): The range of latitudes. + + Returns: + str: pygmt etopo grid resolution based on the given ranges. + """ + + if lon_range > 10 or lat_range > 10: + return '01m' + elif lon_range > 5 or lat_range > 5: + return '15s' + elif lon_range > 2 or lat_range > 2: + return '03s' + elif lon_range > 1 or lat_range > 1: + return '01s' + else: + return '05m' + + @staticmethod + def configure_colormap(colormap): + """ + Configures the colormap based on the given input - as conventions for matplotlib and pygmt can differ + + Args: + colormap (str): The name of the colormap. + + Returns: + tuple: A tuple containing the modified colormap name and a flag indicating + whether the colormap should be reversed. + """ + cmap_reverse_flag = True if colormap.endswith('_r') else False + colormap = colormap[:-2] if cmap_reverse_flag else colormap + return colormap, cmap_reverse_flag + + @staticmethod + def prepare_latex_annotations(label): + """ + Prepares LaTeX annotations for plotting. Uses HTML for compatibility with PyGMT/GMT. + + Args: + label (str): The LaTeX label to be prepared. + + Returns: + str: The prepared label. + + """ + if label.startswith('$') and label.endswith('$'): + # Convert LaTeX to HTML for compatibility with PyGMT/GMT + return f"{label[1:-1]}" + else: + return label + + @staticmethod + def get_header(label, origin, filename, process = None): + """ + Generates a header for a plot based on the provided parameters. + + Args: + label (str): The label for the plot. Defined in default kwargs. + origin (Origin): mtuq.event.Origin object. + filename (str): The filename of the plot. Defined by default the high-level function. Used to retrieve the component. + process (Process, optional): mtuq.process_data.ProcessData object for appropriate dataset. + + Returns: + list: A list containing two lines of the header. + """ + if process is not None: + # get type of waves used for the window + window_type = process.window_type + if window_type == 'surface_wave': + window_type = 'Surface wave' + elif window_type == 'body_wave': + window_type = 'Body wave' + + component = filename.split('/')[-1].split('.')[0] + origin_time = str(origin.time)[0:19] + origin_depth = origin.depth_in_m/1000 + + label = PyGMTUtilities.prepare_latex_annotations(label) + + # if window_type exists, define Rayleigh or Love wave + if process is not None: + if window_type == 'Surface wave' and component == 'Z' or window_type == 'Surface wave' and component == 'R': + # First line of the header defined as: label - Rayleigh wave (component) + header_line_1 = f"{label} - Rayleigh wave ({component})" + elif window_type == 'Surface wave' and component == 'T': + # First line of the header defined as: label - Love wave (component) + header_line_1 = f"{label} - Love wave ({component})" + elif window_type == 'Body wave': + # First line of the header defined as: label - (component) + header_line_1 = f"{label} - Body wave ({component})" + else: + # First line of the header defined as: label - (component) + header_line_1 = f"{label} - ({component})" + + header_line_2 = f"Event Time: {origin_time} UTC, Depth: {origin_depth:.1f} km" + + return [header_line_1, header_line_2] + + @staticmethod + def draw_coastlines(fig, area_thresh=100, water_color='paleturquoise', water_transparency=55): + fig.coast(shorelines=True, area_thresh=area_thresh) + fig.coast(shorelines=False, water=water_color, transparency=water_transparency, area_thresh=area_thresh) \ No newline at end of file From 2a413de3ac3f7954f7ee06cd266b4921e04aa577 Mon Sep 17 00:00:00 2001 From: thurinj Date: Fri, 1 Mar 2024 14:04:47 +1100 Subject: [PATCH 02/31] Added normalize-cross correlation attribute plotting Added a new plotting function for normalized cross-correlation, working with the default and pygmt backend. --- mtuq/graphics/attrs.py | 55 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/mtuq/graphics/attrs.py b/mtuq/graphics/attrs.py index 30878b2fc..179984a01 100644 --- a/mtuq/graphics/attrs.py +++ b/mtuq/graphics/attrs.py @@ -62,6 +62,61 @@ def plot_time_shifts(dirname, attrs, stations, origin, key='total_shift', _plot_attrs(dirname, stations, origin, attrs, key, **kwargs) +def plot_cross_corr(dirname, attrs, stations, origin, key='normalized_cc_max', + **kwargs): + + """ Plots how time shifts vary by location and component + + By default, total time shifts are plotted. To plot just static or + cross-correlation time shifts, use ``key='static_shift'`` or + ``key='time_shift'``, respectively + + .. note :: + + MTUQ distinguishes between the following different types of + time shifts + + - `static_shift` is an initial user-supplied time shift applied during + data processing + + - `time_shift` is a subsequent cross-correlation time shift applied + during misfit evaluation + + - `total_shift` is the total correction, or in other words the sum of + static and cross-correlation time shifts + + + .. rubric :: Required input arguments + + ``dirname`` (`str`): + Directory in which figures will be written + + ``attrs`` (`list` of `AttribDict`): + List returned by misfit function's `collect_attributes` method + + ``stations`` (`list` of `mtuq.Station` objects): + Used to plot station locations + + ``origin`` (`mtuq.Origin` object): + Used to plot origin location + + + .. rubric :: Optional input arguments + + For optional argument descriptions, + `see here `_ + + """ + defaults(kwargs, { + 'label': 'Maximum normalized CC', + 'zero_centered': False, + 'colormap': 'inferno', + 'min_val': 0.0, + 'max_val': 1.0, + }) + + _plot_attrs(dirname, stations, origin, attrs, key, **kwargs) + def plot_amplitude_ratios(dirname, attrs, stations, origin, **kwargs): """ Plots how Aobs/Asyn varies by location and component From 0b68624d1ed699acf83f2a46ed7aeea4699d2436 Mon Sep 17 00:00:00 2001 From: thurinj Date: Wed, 6 Mar 2024 09:59:58 +1100 Subject: [PATCH 03/31] Correction to automatic text placement Fixed text placement. Will now be automatically placed based of figure dimensions / latitude extent, and longitude. --- mtuq/graphics/attrs.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/mtuq/graphics/attrs.py b/mtuq/graphics/attrs.py index 179984a01..ada98d691 100644 --- a/mtuq/graphics/attrs.py +++ b/mtuq/graphics/attrs.py @@ -497,8 +497,14 @@ def _pygmt_backend(filename, values, active_stations, origin, text_line_val = 1 header_lines = PyGMTUtilities.get_header(label, origin, filename, process) + # Add the header text to the plot + # Text spacing is based on longitude range and latitude buffer size. + lon_mean = np.max(longitudes) - (np.max(longitudes) - np.min(longitudes)) / 2 + text_spacing = lat_buffer / 1.5 for header_line in header_lines: - fig.text(x=-148, y=(max(latitudes) + lat_buffer)+(text_line_val)*0.25, text=header_line, font="14p,Helvetica-Bold,black", justify="MC", no_clip=True) + fig.text(x=lon_mean, + y=max(latitudes) + lat_buffer + text_line_val*text_spacing, + text=header_line, font="14p,Helvetica-Bold,black", justify="MC", no_clip=True) text_line_val += 1 # Saving the figure From 0025507161fb694efa08d962d48845d178756fa4 Mon Sep 17 00:00:00 2001 From: thurinj Date: Wed, 6 Mar 2024 10:23:24 +1100 Subject: [PATCH 04/31] Correct docstrings. Updated cross_correlation docstring. --- mtuq/graphics/attrs.py | 25 +++++++++---------------- 1 file changed, 9 insertions(+), 16 deletions(-) diff --git a/mtuq/graphics/attrs.py b/mtuq/graphics/attrs.py index ada98d691..668b066a4 100644 --- a/mtuq/graphics/attrs.py +++ b/mtuq/graphics/attrs.py @@ -65,26 +65,19 @@ def plot_time_shifts(dirname, attrs, stations, origin, key='total_shift', def plot_cross_corr(dirname, attrs, stations, origin, key='normalized_cc_max', **kwargs): - """ Plots how time shifts vary by location and component + """ Plots how cross-correlation values vary by location and component - By default, total time shifts are plotted. To plot just static or - cross-correlation time shifts, use ``key='static_shift'`` or - ``key='time_shift'``, respectively + By default, maximum normalized cross-correlation values are plotted. To plot just + maximum cross-correlation values, use ``key='cc_max'`` .. note :: - MTUQ distinguishes between the following different types of - time shifts - - - `static_shift` is an initial user-supplied time shift applied during - data processing - - - `time_shift` is a subsequent cross-correlation time shift applied - during misfit evaluation - - - `total_shift` is the total correction, or in other words the sum of - static and cross-correlation time shifts - + MTUQ distinguishes between the following different types of + cross-correlation values + + - `cc_max` is the maximum cross-correlation value + + - `normalized_cc_max` is the maximum cross-correlation value normalized between 0 and 1 .. rubric :: Required input arguments From ae1dedcb5ba920147273dace3c9cf0d1f95d340c Mon Sep 17 00:00:00 2001 From: thurinj Date: Thu, 7 Mar 2024 11:26:59 +1100 Subject: [PATCH 05/31] Tweaked pygmt plot appearance Introduced histogram gaussian normalization for background topo map to avoid too-harsh contrasts. Change moment tensor coloring. --- mtuq/graphics/attrs.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/mtuq/graphics/attrs.py b/mtuq/graphics/attrs.py index 668b066a4..1412c0cfc 100644 --- a/mtuq/graphics/attrs.py +++ b/mtuq/graphics/attrs.py @@ -364,15 +364,17 @@ def _pygmt_backend(filename, values, active_stations, origin, resolution = PyGMTUtilities.get_resolution(max(longitudes) - min(longitudes), max(latitudes) - min(latitudes)) grid = pygmt.datasets.load_earth_relief(region=region, resolution=resolution) - # Define a grayscale colormap for topography - pygmt.makecpt(cmap='gray', series=[-7000, 7000]) - # Calculate the gradient (hillshade) grid with azimuth 0/300 and normalization t1 # shade = pygmt.grdgradient(grid=grid, azimuth="0/300", normalize="t1") + + # Define a grayscale colormap for topography + normal = pygmt.grdhisteq.equalize_grid(grid=grid, gaussian=True) + gray_cmap = pygmt.makecpt(cmap='gray', series=[np.min(normal.values), np.max((normal.values))]) + # Plot the hillshade grid as an image if display_topo: - fig.grdimage(grid=grid, shading=shade, projection=f'J{width}i', frame='a', cmap='gray', no_clip=True) + fig.grdimage(grid=normal, shading=shade, projection=f'J{width}i', frame='a', cmap=gray_cmap, no_clip=True) # Overlay coastlines PyGMTUtilities.draw_coastlines(fig) @@ -463,9 +465,9 @@ def _pygmt_backend(filename, values, active_stations, origin, latitude=origin.latitude, depth=10, # Depth is required, even if not used, set to a small number convention="mt", # Use GMT's mt convention - compressionfill="red", + compressionfill="gray15", extensionfill="white", - pen="black" + pen="0.5p,black" ) if station_labels is True: @@ -592,7 +594,7 @@ def get_header(label, origin, filename, process = None): if process is not None: # get type of waves used for the window window_type = process.window_type - if window_type == 'surface_wave': + if window_type == 'surface_wave' or window_type == 'group_velocity': window_type = 'Surface wave' elif window_type == 'body_wave': window_type = 'Body wave' From f2a54658512b8de00cce2663638e172f489f4444 Mon Sep 17 00:00:00 2001 From: thurinj Date: Thu, 7 Mar 2024 17:23:20 +1100 Subject: [PATCH 06/31] Populated new functions docstrings Added complete docstring for the pygmt functions. --- mtuq/graphics/attrs.py | 273 +++++++++++++++++++++++++++++++++++------ 1 file changed, 234 insertions(+), 39 deletions(-) diff --git a/mtuq/graphics/attrs.py b/mtuq/graphics/attrs.py index 1412c0cfc..7ed810040 100644 --- a/mtuq/graphics/attrs.py +++ b/mtuq/graphics/attrs.py @@ -22,17 +22,14 @@ def plot_time_shifts(dirname, attrs, stations, origin, key='total_shift', .. note :: - MTUQ distinguishes between the following different types of - time shifts - - - `static_shift` is an initial user-supplied time shift applied during - data processing + MTUQ distinguishes between the following different types of time shifts: + + - `static_shift` is an initial user-supplied time shift applied during data processing - - `time_shift` is a subsequent cross-correlation time shift applied - during misfit evaluation + - `time_shift` is a subsequent cross-correlation time shift applied during misfit evaluation - - `total_shift` is the total correction, or in other words the sum of - static and cross-correlation time shifts + - `total_shift` is the total correction, or in other words the sum of static and cross-correlation time shifts + .. rubric :: Required input arguments @@ -207,7 +204,18 @@ def _plot_attrs(dirname, stations, origin, attrs, key, for details. Otherwise, defaults to a generic matplotlib `backend `_. + + .. rubric :: Standard pygmt backend + + The standard `pygmt backend `_ is used to + plot station attributes over a hillshaded map. Default calls to front_end functions + can be supplemented with optional keyword arguments to customize the appearance of the plot. + .. code :: + + from mtuq.graphics.attrs import _pygmt_backend + plot_time_shifts('./SW/tshift', attributes_sw, stations, origin, + moment_tensor=best_mt, process=process_sw, backend=_pygmt_backend) """ if backend is None: @@ -244,7 +252,8 @@ def _default_backend(filename, values, stations, origin, colormap='coolwarm', zero_centered=True, colorbar=True, label='', width=5., height=5., **kwargs): - """ Default backend for all other `mtuq.graphics.attrs` functions + """ + Default backend for all frontend `mtuq.graphics.attrs` functions The frontend functions perform only data manipulation. All graphics library calls occur in the backend @@ -252,7 +261,7 @@ def _default_backend(filename, values, stations, origin, By isolating the graphics function calls in this way, users can completely interchange graphics libraries (matplotlib, GMT, PyGMT, and so on) - .. rubric:: Keyword arguments + .. rubric :: Keyword arguments ``colormap`` (`str`): Matplotlib color palette @@ -266,7 +275,6 @@ def _default_backend(filename, values, stations, origin, ``label`` (`str`): Optional colorbar label - """ fig = pyplot.figure(figsize=(width, height)) @@ -340,9 +348,79 @@ def _pygmt_backend(filename, values, active_stations, origin, label='', width=5, moment_tensor=None, process=None, stations_list=None, station_labels=True, min_val=None, max_val=None, **kwargs): """ - PyGMT backend for plotting station attributes with hillshading using the - Miller Cylindrical projection, with an azimuth of 0/90 and a normalization - of t1 for the hillshade intensity. + PyGMT backend for plotting station attributes over a hillshaded map. + + .. note :: + + This function requires the PyGMT library to be installed. + If called while pygmt is not installed, the default matplotlib backend will be used. + + The function accepts a number of optional keyword arguments to customize + the appearance of the plot. If passed to another backend, these arguments + will be ignored. + + .. rubric :: Required input arguments + + ``filename`` (`str`): + The name of the file to which the figure will be saved. This should be passed by the frontend function. + Expected filenames are in the format `component.png`, where `component` is the component of the data being plotted. + This might be revised in the future to allow for more flexible naming conventions. + + ``values`` (`list` of `float`): + List of values to be plotted for each station. The length of the list should match the number of stations. + + ``active_stations`` (`list` of `mtuq.Station` objects): + List of stations to be plotted, with a value entry for each station in the `values` list. + + ``origin`` (`mtuq.Origin` object): + Origin object used to plot the origin location. + + .. rubric :: Keyword arguments + + ``colormap`` (`str`): + GMT colormap name - see `GMT documentation ` + + ``zero_centered`` (`bool`): + Whether or not the colormap is centered on zero + + ``display_topo`` (`bool`): + Whether or not to display topography in the background -- will download the ETOPO1 topography file if not found + + ``label`` (`str`): + Optional colorbar text label -- supports LaTeX expressions + + ``width`` (`float`): + Width of the figure in inches -- default is 5 inches + + ``moment_tensor`` (`mtuq.event.MomentTensor` or `np.ndarray`): + Moment tensor to plot as a beachball -- will plot a star at the origin if not provided + Can be either a `mtuq.event.MomentTensor` object or a 1D numpy array with the six independent components + [m_rr, m_tt, m_pp, m_rt, m_rp, m_tp] + + ``process`` (`mtuq.ProcessData`): + ProcessData object used to determine the type of waves used for the window + + ``stations_list`` (`list` of `mtuq.Station` objects): + List of all stations available before data processing -- will plot only active stations if not provided + + ``station_labels`` (`bool`): + Whether or not to plot station labels -- default is True + + ``min_val`` (`float`): + Minimum value for the colorbar -- will be determined automatically if not provided + + ``max_val`` (`float`): + Maximum value for the colorbar -- will be determined automatically if not provided + + .. rubric :: Backend specific utility class + + The PyGMTUtilities class is a utility class designed to enhance and simplify the usage of PyGMT for plotting. + It includes methods for calculating plotting regions with buffers, configuring colormaps, preparing LaTeX + annotations for PyGMT, and generating standardized headers for plots. Documentation for the PyGMTUtilities + class can be found in the `PyGMTUtilities `_ module. + + + """ import pygmt @@ -506,8 +584,54 @@ def _pygmt_backend(filename, values, active_stations, origin, fig.savefig(filename, crop=True, dpi=300) class PyGMTUtilities: + """ + Utility class for PyGMT plotting backend. + + This class offers a set of static methods designed for enhancing and simplifying the usage of PyGMT + for plotting by handling plotting regions, color maps, LaTeX annotations, and plot headers. + + .. note :: + The class is designed to be used without instantiation due to its static methods. This approach + helps in organizing code related to the PyGMT plotting backend and avoids confusion with other plotting backends. + + Methods include calculating plotting regions with buffers, configuring colormaps, preparing LaTeX + annotations for PyGMT, and generating standardized headers for plots. + + Examples and more detailed method descriptions can be found in the documentation of each method. + """ + @staticmethod def calculate_plotting_region(stations, origin, buffer_percentage=0.1): + """ + Calculates the region for plotting, including a buffer area around specified stations and origin. + + .. rubric :: Parameters + + ``stations`` (`list` of `mtuq.Station` objects): + The stations to be included in the plot. + + ``origin`` (`mtuq.Origin` object): + The origin object is used to calculate the region for the plot in case the origin is outside the range of the stations. + + ``buffer_percentage`` (`float`, optional): + The percentage of the total longitude and latitude range to be added as a buffer around the specified region. + Defaults to 0.1 (10%). + + .. rubric :: Returns + + ``region`` (`list` of `float`), ``lat_buffer`` (`float`): + A tuple containing the calculated region as a list `[west, east, south, north]` and the latitude buffer value. + The latitude buffer is returned to later be used for adjusting text spacing in the plot header. + + .. rubric :: Example + + >>> region, lat_buffer = PyGMTUtilities.calculate_plotting_region(stations, origin) + >>> print(region) + [149.55, 151.45, -35.1, -32.9] + >>> print(lat_buffer) + 0.22 + """ + longitudes = [station.longitude for station in stations] + [origin.longitude] latitudes = [station.latitude for station in stations] + [origin.latitude] @@ -522,14 +646,35 @@ def calculate_plotting_region(stations, origin, buffer_percentage=0.1): @staticmethod def get_resolution(lon_range, lat_range): """ - Determines the resolution based on the given longitude and latitude ranges. + Determines the appropriate PyGMT etopo grid resolution based on longitude and latitude ranges. + + .. rubric :: Parameters + + ``lon_range`` (`float`): + The longitudinal range of the area of interest. + + ``lat_range`` (`float`): + The latitudinal range of the area of interest. + + .. rubric :: Returns + + ``resolution`` (`str`): + The resolution string for PyGMT, e.g., '01m', '15s', ..., based on the size of the specified area. + + .. note :: + The resolution is determined based on predefined thresholds for the ranges, aiming to balance + detail and performance for different scales of geographic areas - Args: - lon_range (float): The range of longitudes. - lat_range (float): The range of latitudes. + - If lon_range > 10 or lat_range > 10, the resolution is '01m'. + + - If lon_range > 5 or lat_range > 5, the resolution is '15s'. + + - If lon_range > 2 or lat_range > 2, the resolution is '03s'. + + - If lon_range > 1 or lat_range > 1, the resolution is '01s'. + + Otherwise, the resolution is '05m'. - Returns: - str: pygmt etopo grid resolution based on the given ranges. """ if lon_range > 10 or lat_range > 10: @@ -546,14 +691,31 @@ def get_resolution(lon_range, lat_range): @staticmethod def configure_colormap(colormap): """ - Configures the colormap based on the given input - as conventions for matplotlib and pygmt can differ + Adjusts the given colormap name for compatibility with PyGMT and matplotlib conventions. + + .. rubric :: Parameters + + ``colormap`` (`str`): + The name of the colormap to be used. If the colormap name ends with '_r', the colormap is + reversed, and the '_r' suffix is removed. + + .. rubric :: Returns + + ``colormap`` (`str`), ``cmap_reverse_flag`` (`bool`): + A tuple containing the adjusted colormap name and a boolean indicating whether the colormap should + be reversed. + + .. note :: - Args: - colormap (str): The name of the colormap. + The method accept only colormaps that are available in PyGMT. For a list of available colormaps, - Returns: - tuple: A tuple containing the modified colormap name and a flag indicating - whether the colormap should be reversed. + .. rubric :: Example + + >>> colormap, reverse = PyGMTUtilities.configure_colormap('viridis_r') + >>> print(colormap) + viridis + >>> print(reverse) + True """ cmap_reverse_flag = True if colormap.endswith('_r') else False colormap = colormap[:-2] if cmap_reverse_flag else colormap @@ -562,13 +724,18 @@ def configure_colormap(colormap): @staticmethod def prepare_latex_annotations(label): """ - Prepares LaTeX annotations for plotting. Uses HTML for compatibility with PyGMT/GMT. + Prepares LaTeX annotations for plotting. Uses HTML tags instead + of $•$ for compatibility with PyGMT/GMT. + + .. rubric :: Parameters + + ``label`` (`str`): + The LaTeX label to be prepared. - Args: - label (str): The LaTeX label to be prepared. + .. rubric :: Returns - Returns: - str: The prepared label. + ``str``: + The prepared label. """ if label.startswith('$') and label.endswith('$'): @@ -582,14 +749,24 @@ def get_header(label, origin, filename, process = None): """ Generates a header for a plot based on the provided parameters. - Args: - label (str): The label for the plot. Defined in default kwargs. - origin (Origin): mtuq.event.Origin object. - filename (str): The filename of the plot. Defined by default the high-level function. Used to retrieve the component. - process (Process, optional): mtuq.process_data.ProcessData object for appropriate dataset. + .. rubric :: Parameters + + ``label`` (`str`): + The label for the plot. Usually defined in the frontend function. + + ``origin`` (mtuq.Origin): + mtuq.event.Origin object, used to retrieve the event time and depth. - Returns: - list: A list containing two lines of the header. + ``filename`` (str): + The filename of the plot. Defined by default the high-level function. Used to retrieve the component. + + ``process`` (Process, optional): + mtuq.process_data.ProcessData object for appropriate dataset. + + .. rubric :: Returns + + ``list``: + A list containing two lines of the header. [Label - (component)], [Event Time: (time) UTC, Depth: (depth) km] """ if process is not None: # get type of waves used for the window @@ -626,5 +803,23 @@ def get_header(label, origin, filename, process = None): @staticmethod def draw_coastlines(fig, area_thresh=100, water_color='paleturquoise', water_transparency=55): + """ + Draws coastlines and fills water areas with a transparent blue shade. + + .. rubric :: Parameters + + ``fig`` (pygmt.Figure): + The PyGMT figure object to which the coastlines and water areas will be added. + + ``area_thresh`` (`int`, optional): + The minimum area of land to be displayed. Defaults to 100. + + ``water_color`` (`str`, optional): + The color of the water areas. Defaults to 'paleturquoise'. + + ``water_transparency`` (`int`, optional): + The transparency of the water areas. Defaults to 55. + + """ fig.coast(shorelines=True, area_thresh=area_thresh) fig.coast(shorelines=False, water=water_color, transparency=water_transparency, area_thresh=area_thresh) \ No newline at end of file From 2e31ccd9380d12dee8a9af51744a7a230ea4676c Mon Sep 17 00:00:00 2001 From: thurinj Date: Thu, 7 Mar 2024 17:24:10 +1100 Subject: [PATCH 07/31] Improved math expression extraction Changed the function so that it can handle all combination of text + latex math expressions. --- mtuq/graphics/attrs.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/mtuq/graphics/attrs.py b/mtuq/graphics/attrs.py index 7ed810040..3d649cd0d 100644 --- a/mtuq/graphics/attrs.py +++ b/mtuq/graphics/attrs.py @@ -738,12 +738,18 @@ def prepare_latex_annotations(label): The prepared label. """ - if label.startswith('$') and label.endswith('$'): - # Convert LaTeX to HTML for compatibility with PyGMT/GMT - return f"{label[1:-1]}" - else: + + parts = label.split('$') + if len(parts) == 1: # No '$' found return label - + new_text = '' + for i, part in enumerate(parts): + if i % 2 == 0: + new_text += part + else: + new_text += f"{part}" + return new_text + @staticmethod def get_header(label, origin, filename, process = None): """ From d7d1db936b86d1cf6281c4fcbd2987b18f9a56f3 Mon Sep 17 00:00:00 2001 From: thurinj Date: Thu, 7 Mar 2024 17:25:29 +1100 Subject: [PATCH 08/31] Updated docs sources Added a few key functions for the attributes pygmt plotting backend. Compilation was tested with the building script and html renders as expected. --- docs/library/autogen.rst | 3 +++ docs/library/index.rst | 3 ++- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/docs/library/autogen.rst b/docs/library/autogen.rst index 426fa47c0..bc738ef10 100644 --- a/docs/library/autogen.rst +++ b/docs/library/autogen.rst @@ -44,8 +44,11 @@ autogen mtuq.graphics.plot_time_shifts mtuq.graphics.plot_amplitude_ratios mtuq.graphics.plot_log_amplitude_ratios + mtuq.graphics.plot_cross_corr mtuq.graphics._plot_attrs mtuq.graphics.attrs._default_backend + mtuq.graphics.attrs._pygmt_backend + mtuq.graphics.attrs.PyGMTUtilities mtuq.grid.DeviatoricGridRandom mtuq.grid.DeviatoricGridSemiregular mtuq.grid.DoubleCoupleGridRandom diff --git a/docs/library/index.rst b/docs/library/index.rst index 1f7c7022e..a796d05a2 100644 --- a/docs/library/index.rst +++ b/docs/library/index.rst @@ -93,12 +93,13 @@ Depth and hypocenter visualization ============================================================================================================ ============================================================================================================ -Time shift and amplitude ratio visualization +Station attributes visualization -------------------------------------------- ============================================================================================================ ============================================================================================================ `mtuq.graphics.plot_time_shifts `_ Plots time shifts by location and component `mtuq.graphics.plot_amplitude_ratios `_ Plots amplitude ratios by location and component +`mtuq.graphics.plot_cross_corr `_ Plots normalized cross-correlation by location and component ============================================================================================================ ============================================================================================================ From e58d5ee024e090f1bf13643f4eddbee1316e7e62 Mon Sep 17 00:00:00 2001 From: Amanda McPherson Date: Tue, 12 Mar 2024 14:06:18 -0800 Subject: [PATCH 09/31] Permutation and naming fix --- mtuq/greens_tensor/SPECFEM3D.py | 18 +++++++++--------- mtuq/io/clients/SPECFEM3D_SAC.py | 6 +++--- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/mtuq/greens_tensor/SPECFEM3D.py b/mtuq/greens_tensor/SPECFEM3D.py index 27be1d81d..acc474c37 100644 --- a/mtuq/greens_tensor/SPECFEM3D.py +++ b/mtuq/greens_tensor/SPECFEM3D.py @@ -85,18 +85,18 @@ def _precompute_force(self): for _i, component in enumerate(self.components): if component=='Z': - array[_i, _j+0, :] = self.select(channel="Z.Fe")[0].data - array[_i, _j+1, :] = self.select(channel="Z.Fn")[0].data - array[_i, _j+2, :] = self.select(channel="Z.Fz")[0].data + array[_i, _j+0, :] = self.select(channel="Z.Fz")[0].data + array[_i, _j+1, :] = self.select(channel="Z.Fs")[0].data + array[_i, _j+2, :] = self.select(channel="Z.Fe")[0].data elif component=='R': - array[_i, _j+0, :] = self.select(channel="R.Fe")[0].data - array[_i, _j+1, :] = self.select(channel="R.Fn")[0].data - array[_i, _j+2, :] = self.select(channel="R.Fz")[0].data + array[_i, _j+0, :] = self.select(channel="R.Fz")[0].data + array[_i, _j+1, :] = self.select(channel="R.Fs")[0].data + array[_i, _j+2, :] = self.select(channel="R.Fe")[0].data elif component=='T': - array[_i, _j+0, :] = self.select(channel="T.Fe")[0].data - array[_i, _j+1, :] = self.select(channel="T.Fn")[0].data - array[_i, _j+2, :] = self.select(channel="T.Fz")[0].data + array[_i, _j+0, :] = self.select(channel="T.Fz")[0].data + array[_i, _j+1, :] = self.select(channel="T.Fs")[0].data + array[_i, _j+2, :] = self.select(channel="T.Fe")[0].data diff --git a/mtuq/io/clients/SPECFEM3D_SAC.py b/mtuq/io/clients/SPECFEM3D_SAC.py index 27cb83292..919529e48 100644 --- a/mtuq/io/clients/SPECFEM3D_SAC.py +++ b/mtuq/io/clients/SPECFEM3D_SAC.py @@ -35,9 +35,9 @@ 'R.Fe', 'T.Fe', 'Z.Fe', - 'R.Fn', - 'T.Fn', - 'Z.Fn', + 'R.Fs', + 'T.Fs', + 'Z.Fs', 'R.Fz', 'T.Fz', 'Z.Fz', From e51f65b985f84c7ad3734cfd582a8dce95d6e7c7 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Tue, 6 Feb 2024 00:01:11 -0700 Subject: [PATCH 10/31] Updated URLs --- docs/install/issues.rst | 2 +- tests/check_urls.bash | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/install/issues.rst b/docs/install/issues.rst index 14fc4664e..8ee8368a3 100644 --- a/docs/install/issues.rst +++ b/docs/install/issues.rst @@ -37,7 +37,7 @@ Older versions of the conda package manager can be very slow. For a significant conda update -n base conda -For reference, the largest potential speed up comes from the new `mamba `_ dependency solver, which was `adopted `_ in the 23.10 release. +For reference, the largest potential speed up comes from the new `mamba `_ dependency solver, which was `adopted `_ in the 23.10 release. diff --git a/tests/check_urls.bash b/tests/check_urls.bash index 9053377f3..9bef6b62d 100755 --- a/tests/check_urls.bash +++ b/tests/check_urls.bash @@ -20,6 +20,8 @@ URLS="\ https://docs.xarray.dev/en/stable/generated/xarray.DataArray.html\ https://github.com/uafgeotools/mtuq/blob/master/docs/user_guide/05/code/gallery_mt.py https://github.com/uafgeotools/mtuq/blob/master/docs/user_guide/05/code/gallery_force.py + https://conda.org/blog/2023-11-06-conda-23-10-0-release\ + https://www.anaconda.com/blog/a-faster-conda-for-a-growing-community\ " From ef33c26737cb5943d89126e7b628e157993e2fdb Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Tue, 6 Feb 2024 13:19:22 -0700 Subject: [PATCH 11/31] New Greens function docs pages --- docs/user_guide/03.rst | 39 ++++++++++++++++------------ docs/user_guide/03/receiver_side.rst | 15 +++++++++++ docs/user_guide/03/source_side.rst | 15 +++++++++++ 3 files changed, 53 insertions(+), 16 deletions(-) create mode 100644 docs/user_guide/03/receiver_side.rst create mode 100644 docs/user_guide/03/source_side.rst diff --git a/docs/user_guide/03.rst b/docs/user_guide/03.rst index d66655088..4918a826b 100644 --- a/docs/user_guide/03.rst +++ b/docs/user_guide/03.rst @@ -11,6 +11,24 @@ Role of Green's functions By combining Green's functions, it is possible to obtain synthetics from any moment tensor or force source. Generating synthetics in this way is useful because it provides cost savings compared with on-the-fly wavefield simulations. Synthetics can then be compared with data to determine a best-fitting source. +Green's function conventions +---------------------------- + +MTUQ consistently uses an `up-south-east` `convention `_ for internally storing Green's functions (as well as moment tensors and forces). + +A major goal of MTUQ is to avoid exposing users to unnecessary complexity. For Green's functions, MTUQ tries to accomplish this by understanding external conventions and converting to a common internal format. + +Under the hood, MTUQ deals with a variety of externel conventions related to + +- the symmetry of the medium (for example, 1D media require fewer independent Green's functions than 3D media) + +- the choice of local Cartesian basis convention (for example, some authors employ `up-south-east`, others `north-east-down`; see `ObsPy documentation `_ for more information) + +- the type of medium under consideration (for example, acoustic media require fewer independent Green's functions than elastic media) + +A better sense of what's involved can be obtained by browsing the `source code `_ and references therein. + + `GreensTensor` objects ---------------------- @@ -81,7 +99,7 @@ Computing Green's functions from 3D Earth models MTUQ currently supports 3D Green's functions from SPECFEM3D/3D_GLOBE. -To generate a complete Green's function database for a given hypocenter and depth, six SPECFEM3D wavefield simulations are required. Output must be saved as/converted to SAC files at individual stations using the following filename convention, which comes from `GRD_CMT3D `_. Place all SAC files corresponding to a single hypocenter and depth in the same directory as follows: +To generate a complete Green's function database for a given hypocenter and depth, six SPECFEM3D/3D_GLOBE wavefield simulations are required. Output must be saved as/converted to SAC files at individual stations using the following filename convention, which comes from `GRD_CMT3D `_. Place all SAC files corresponding to a single hypocenter and depth in the same directory as follows: .. code :: @@ -114,27 +132,16 @@ To open an SPECFEM3D database client in MTUQ: db = open_db(path_SPECFEM3D_output_directory, format="SPECFEM3D") -Once opened, a SPECFEM3D database client can be used to generate `GreensTensor `_ objects as follows: +Once opened, a SPECFEM3D/3D_GLOBE database client can be used to generate `GreensTensor `_ objects as follows: .. code:: greens_tensors = db.get_greens_tensors(stations, origin) +For more information, please see: -Green's function conventions ----------------------------- - -A variety of Green's function conventions exist. Figuring out which are used in a particular application can be challenging because it depends on - -- the type of medium under consideration (for example, acoustic media require fewer independent Green's functions than elastic media) - -- the symmetry of the medium (for example, 1D media require fewer independent Green's functions than 3D media) - -- the choice of local Cartesian basis conventions (for example, some authors employ `up-south-east`, others `north-east-down`; see `ObsPy documentation `_ for more information) - -A major goal is to avoid exposing users to unnecessary basis complexity. MTUQ accomplishes this by understanding external formats and converting to a common internal format that works for both 1D and 3D media. - -For internally storing moment tensors, forces, and Green's functions, MTUQ consistently uses an `up-south-east` Cartesian convention. +`Detailed guide to source-side 3D Green's functions (under construction) `_ +`Detailed guide to receiver-side 3D Green's functions (under construction) `_ diff --git a/docs/user_guide/03/receiver_side.rst b/docs/user_guide/03/receiver_side.rst new file mode 100644 index 000000000..f1a4be8a7 --- /dev/null +++ b/docs/user_guide/03/receiver_side.rst @@ -0,0 +1,15 @@ + +.. warning:: + + This page is currently just as stub. + + To fill in missing documentation, feel free to submit a pull request. + + +Working example +=============== + +A working of example using receiver-side 3D Green's functions from SPECFEM3D in a moment tensor inversion: + +`test_greens_SPECFEM3D_SGT.py `_ + diff --git a/docs/user_guide/03/source_side.rst b/docs/user_guide/03/source_side.rst new file mode 100644 index 000000000..56fdfc872 --- /dev/null +++ b/docs/user_guide/03/source_side.rst @@ -0,0 +1,15 @@ + +.. warning:: + + This page is currently just as stub. + + To fill in missing documentation, feel free to submit a pull request. + + +Working example +=============== + +A working of example using source-side 3D Green's functions from SPECFEM3D_GLOBE in a moment tensor inversion: + +`test_greens_SPECFEM3D_SAC.py `_ + From dc5116b9c92f36303a2f7890bce93825c95540f2 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Tue, 6 Feb 2024 17:33:58 -0700 Subject: [PATCH 12/31] Filled in source-side Greens function notes --- docs/user_guide/03.rst | 33 ++++------- docs/user_guide/03/receiver_side.rst | 4 +- docs/user_guide/03/source_side.rst | 85 +++++++++++++++++++++++++++- 3 files changed, 98 insertions(+), 24 deletions(-) diff --git a/docs/user_guide/03.rst b/docs/user_guide/03.rst index 4918a826b..ce045e12f 100644 --- a/docs/user_guide/03.rst +++ b/docs/user_guide/03.rst @@ -18,15 +18,7 @@ MTUQ consistently uses an `up-south-east` `convention `_ for more information) - -- the type of medium under consideration (for example, acoustic media require fewer independent Green's functions than elastic media) - -A better sense of what's involved can be obtained by browsing the `source code `_ and references therein. +Under the hood, MTUQ deals with a variety of external conventions related to (1) the symmetry of the medium (for example, 1D media require fewer independent Green's functions than 3D media); (2) the choice of local Cartesian basis convention (for example, some authors employ `up-south-east`, others `north-east-down`). A sense of what's involved can be obtained by browsing the `source code `_ and references therein. @@ -35,8 +27,7 @@ A better sense of what's involved can be obtained by browsing the `source code < `GreensTensor `_ objects provide access to a set of Green's functions describing the medium response between a given station and origin. -Methods built into the `GreensTensor` class allow data processing and synthetics generation. In particular, the `get_synthetics `_ method accepts a `MomentTensor` or `Force` and returns corresponding synthetics for the given station. - +Methods built into the `GreensTensor` class allow data processing and synthetics generation. In particular, the `get_synthetics `_ method accepts a `MomentTensor` or `Force` and returns corresponding synthetics. @@ -52,7 +43,7 @@ To download AK135 Green's functions and generate MTUQ `GreensTensor` objects: from mtuq import download_greens_functions greens_tensors = download_greens_functions(stations, origins, model="ak135f_2s") -A limitation of syngine is that Green's functions can be downloaded only on a station-by-station basis, not over an entire subset of Earth. An alternative that avoids this limitation is to compute your own Green's functions. +A limitation of syngine is that Green's functions can be downloaded only on a station-by-station basis, not over an entire area or volume. An alternative that avoids this limitation is to compute your own Green's functions. @@ -62,7 +53,7 @@ Computing Green's functions from 1D Earth models MTUQ supports the following 1D Green's functions formats: AxiSEM (preferred) and FK. -`AxiSEM `_ performs spectral element wave simulations in radially-symmetric Earth models. To generate synthetics in a format MTUQ natively supports, follow the instructions for in the `AxiSEM user manual `_ under "Output wavefields in NetCDF formated needed by instaseis." AxiSEM NetCDF files can be used to retrieve vertical, radial, and transverse displacement in units of m*(N-m)^-1. +`AxiSEM `_ performs spectral element wave simulations in radially-symmetric Earth models. To generate synthetics in a format MTUQ natively supports, follow the instructions for in the `AxiSEM user manual `_ under "Output wavefields in NetCDF format needed by instaseis." AxiSEM NetCDF files can be used to retrieve vertical, radial, and transverse displacement in units of m*(N-m)^-1. To open an AxiSEM database client in MTUQ: @@ -99,11 +90,11 @@ Computing Green's functions from 3D Earth models MTUQ currently supports 3D Green's functions from SPECFEM3D/3D_GLOBE. -To generate a complete Green's function database for a given hypocenter and depth, six SPECFEM3D/3D_GLOBE wavefield simulations are required. Output must be saved as/converted to SAC files at individual stations using the following filename convention, which comes from `GRD_CMT3D `_. Place all SAC files corresponding to a single hypocenter and depth in the same directory as follows: +To generate a full set of Green's functions for a given hypocenter and depth, six SPECFEM3D/3D_GLOBE wavefield simulations are required. Output must be saved as/converted to SAC files at individual stations using the following filename convention, which comes from `GRD_CMT3D `_. Place all SAC files corresponding to a single hypocenter and depth in the same directory as follows: .. code :: - {basedir}/{event_id}/ + {event_id}/ {net}.{sta}.{loc}.Z.Mrr.sac {net}.{sta}.{loc}.Z.Mtt.sac {net}.{sta}.{loc}.Z.Mpp.sac @@ -124,24 +115,24 @@ To generate a complete Green's function database for a given hypocenter and dept {net}.{sta}.{loc}.T.Mtp.sac -To open an SPECFEM3D database client in MTUQ: +To open a SPECFEM3D/3D_GLOBE database client: -.. code :: +.. code:: from mtuq import open_db db = open_db(path_SPECFEM3D_output_directory, format="SPECFEM3D") -Once opened, a SPECFEM3D/3D_GLOBE database client can be used to generate `GreensTensor `_ objects as follows: +Once opened, the database client can be used to generate `GreensTensor `_ objects as follows: .. code:: greens_tensors = db.get_greens_tensors(stations, origin) -For more information, please see: +For more information, also see: -`Detailed guide to source-side 3D Green's functions (under construction) `_ +`Source-side SPECFEM3D/3D_GLOBE Green's functions `_ -`Detailed guide to receiver-side 3D Green's functions (under construction) `_ +`Receiver-side SPECFEM3D/3D_GLOBE Green's functions (under construction) `_ diff --git a/docs/user_guide/03/receiver_side.rst b/docs/user_guide/03/receiver_side.rst index f1a4be8a7..deaa470d3 100644 --- a/docs/user_guide/03/receiver_side.rst +++ b/docs/user_guide/03/receiver_side.rst @@ -6,8 +6,8 @@ To fill in missing documentation, feel free to submit a pull request. -Working example -=============== +Notes on receiver-side 3D Green's functions +=========================================== A working of example using receiver-side 3D Green's functions from SPECFEM3D in a moment tensor inversion: diff --git a/docs/user_guide/03/source_side.rst b/docs/user_guide/03/source_side.rst index 56fdfc872..605db9db1 100644 --- a/docs/user_guide/03/source_side.rst +++ b/docs/user_guide/03/source_side.rst @@ -6,8 +6,91 @@ To fill in missing documentation, feel free to submit a pull request. +Notes on source-side 3D Green's functions +========================================= + +**File format** + +- Currently, MTUQ reads source-side 3D Green's functions from SAC files + +- SAC output format is natively supported by SPECFEM3D_GLOBE, but not SPECFEM3D Cartesian (output from the latter can be manually converted) + + +**Units convention** + +- MTUQ treats each individual SAC file as the response in meters to a unit (1 Newton-meter) force couple + +- Users must ensure that source-side Green's functions are normalized according to the above units + + +**Basis convention** + +- For moment tensors, MTUQ uses an `up-south-east` convention, identical to one the used by SPECFEM3D_GLOBE + + +**Origin time convention** + +- For origin time, MTUQ uses a centroid convention (`more details `_), so that `t=0` in the `GreensTensor` time discretization corresponds to mean source excitation time + +- MTUQ uses the `REF_TIME` header from the SPECFEM3D_GLOBE SAC output files, which gives the peak excitation of the source relative to the simulation start time + +- MTUQ ignores the origin time given in the CMTSOLUTION file and `ORIGIN_TIME` header + + +**Depth searches (experimental)** + +- Only depth searches are possible with source-side 3D Green's functions (no other hypocenter parameters) + +- The current `depth search `_ implementation is somewhat crude and experimental; consider local modifcations to suit your needs + +- To allow depth searches, create subdirectories for each centroid depth, as below + +.. code :: + + {event_id}/ + {depth_in_km}/ + {net}.{sta}.{loc}.Z.Mrr.sac + {net}.{sta}.{loc}.Z.Mtt.sac + {net}.{sta}.{loc}.Z.Mpp.sac + {net}.{sta}.{loc}.Z.Mrt.sac + {net}.{sta}.{loc}.Z.Mrp.sac + {net}.{sta}.{loc}.Z.Mtp.sac + {net}.{sta}.{loc}.R.Mrr.sac + {net}.{sta}.{loc}.R.Mtt.sac + {net}.{sta}.{loc}.R.Mpp.sac + {net}.{sta}.{loc}.R.Mrt.sac + {net}.{sta}.{loc}.R.Mrp.sac + {net}.{sta}.{loc}.R.Mtp.sac + {net}.{sta}.{loc}.T.Mrr.sac + {net}.{sta}.{loc}.T.Mtt.sac + {net}.{sta}.{loc}.T.Mpp.sac + {net}.{sta}.{loc}.T.Mrt.sac + {net}.{sta}.{loc}.T.Mrp.sac + {net}.{sta}.{loc}.T.Mtp.sac + {depth_in_km}/ + {net}.{sta}.{loc}.Z.Mrr.sac + {net}.{sta}.{loc}.Z.Mtt.sac + {net}.{sta}.{loc}.Z.Mpp.sac + {net}.{sta}.{loc}.Z.Mrt.sac + {net}.{sta}.{loc}.Z.Mrp.sac + {net}.{sta}.{loc}.Z.Mtp.sac + {net}.{sta}.{loc}.R.Mrr.sac + {net}.{sta}.{loc}.R.Mtt.sac + {net}.{sta}.{loc}.R.Mpp.sac + {net}.{sta}.{loc}.R.Mrt.sac + {net}.{sta}.{loc}.R.Mrp.sac + {net}.{sta}.{loc}.R.Mtp.sac + {net}.{sta}.{loc}.T.Mrr.sac + {net}.{sta}.{loc}.T.Mtt.sac + {net}.{sta}.{loc}.T.Mpp.sac + {net}.{sta}.{loc}.T.Mrt.sac + {net}.{sta}.{loc}.T.Mrp.sac + {net}.{sta}.{loc}.T.Mtp.sac + ... + + Working example -=============== +--------------- A working of example using source-side 3D Green's functions from SPECFEM3D_GLOBE in a moment tensor inversion: From bcebc6b5a8fbb1bd3de77c952467fa942b14e17d Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Fri, 9 Feb 2024 20:24:01 -0700 Subject: [PATCH 13/31] Revised comments/docstrings --- docs/install/issues.rst | 2 +- docs/user_guide/03.rst | 6 ++++-- docs/user_guide/03/receiver_side.rst | 2 +- docs/user_guide/03/source_side.rst | 2 +- mtuq/greens_tensor/base.py | 2 +- tests/test_time_shifts.py | 9 ++++----- 6 files changed, 12 insertions(+), 11 deletions(-) diff --git a/docs/install/issues.rst b/docs/install/issues.rst index 8ee8368a3..5ac380216 100644 --- a/docs/install/issues.rst +++ b/docs/install/issues.rst @@ -31,7 +31,7 @@ We note that some versions of GMT and ObsPy do not plot `full moment tensors `_ performs spectral element wave simulations in radially-symmetric Earth models. To generate synthetics in a format MTUQ natively supports, follow the instructions for in the `AxiSEM user manual `_ under "Output wavefields in NetCDF format needed by instaseis." AxiSEM NetCDF files can be used to retrieve vertical, radial, and transverse displacement in units of m*(N-m)^-1. +`AxiSEM `_ performs spectral element wave simulations in radially-symmetric Earth models. AxiSEM NetCDF files can be used to retrieve vertical, radial, and transverse displacement in units of m*(N-m)^-1. + +To generate AxiSEM synthetics in a format MTUQ supports, follow the instructions in the `AxiSEM user manual `_ under "Output wavefields in NetCDF format needed by instaseis." To open an AxiSEM database client in MTUQ: diff --git a/docs/user_guide/03/receiver_side.rst b/docs/user_guide/03/receiver_side.rst index deaa470d3..10994583b 100644 --- a/docs/user_guide/03/receiver_side.rst +++ b/docs/user_guide/03/receiver_side.rst @@ -1,7 +1,7 @@ .. warning:: - This page is currently just as stub. + This page is currently just a stub. To fill in missing documentation, feel free to submit a pull request. diff --git a/docs/user_guide/03/source_side.rst b/docs/user_guide/03/source_side.rst index 605db9db1..df27bdaf9 100644 --- a/docs/user_guide/03/source_side.rst +++ b/docs/user_guide/03/source_side.rst @@ -1,7 +1,7 @@ .. warning:: - This page is currently just as stub. + This page is currently just a stub. To fill in missing documentation, feel free to submit a pull request. diff --git a/mtuq/greens_tensor/base.py b/mtuq/greens_tensor/base.py index 34c5cfea1..3ec868d4a 100644 --- a/mtuq/greens_tensor/base.py +++ b/mtuq/greens_tensor/base.py @@ -471,7 +471,7 @@ def sort_by_azimuth(self, reverse=False): def sort_by_function(self, function, reverse=False): - """ Sorts in-place using the python built-in `sort` + """ Sorts in-place by user-supplied function """ self.sort(key=function, reverse=reverse) diff --git a/tests/test_time_shifts.py b/tests/test_time_shifts.py index 0a6b31374..2c5db6a25 100644 --- a/tests/test_time_shifts.py +++ b/tests/test_time_shifts.py @@ -161,13 +161,16 @@ def __call__(self, traces): return super(ProcessData, self).__call__( traces, station=station, origin=origin) + + # + # plotting functions + # def _plot_dat(axis, t, data, attrs, pathspec='-k'): stream = data[0] trace = data[0][0] stats = trace.stats axis.plot(t, trace.data, pathspec) - def _plot_syn(axis, t, data, attrs, pathspec='-r'): stream = data[0] trace = data[0][0] @@ -176,7 +179,6 @@ def _plot_syn(axis, t, data, attrs, pathspec='-r'): idx2 = attrs.idx_stop axis.plot(t, trace.data[idx1:idx2], pathspec) - def _annotate(axis, attrs): text = 'static_shift: %.1f' % attrs.static_shift pyplot.text(0.02, 0.85, text, transform=axis.transAxes) @@ -187,15 +189,12 @@ def _annotate(axis, attrs): text = 'total_shift: %.1f' % attrs.total_shift pyplot.text(0.02, 0.55, text, transform=axis.transAxes) - def _get_synthetics(greens, mt): return greens.get_synthetics(mt, components=['Z']) - def _get_attrs(data, greens, misfit): return misfit.collect_attributes(data, greens, mt)[0]['Z'] - def _get_time_sampling(data): stream = data[0] t1 = float(stream[0].stats.starttime) From a138cf7c2db3613af51de59501460181b22bd945 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Mon, 11 Mar 2024 16:26:40 -0600 Subject: [PATCH 14/31] Added CPS GreensTensor subclass --- mtuq/greens_tensor/CPS.py | 109 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 109 insertions(+) create mode 100644 mtuq/greens_tensor/CPS.py diff --git a/mtuq/greens_tensor/CPS.py b/mtuq/greens_tensor/CPS.py new file mode 100644 index 000000000..ba229c6e0 --- /dev/null +++ b/mtuq/greens_tensor/CPS.py @@ -0,0 +1,109 @@ + +import obspy +import numpy as np + +from mtuq.greens_tensor.base import GreensTensor as GreensTensorBase + + + +print('WARNING: CPS Greens functions are not fully tested yet') + + + +class GreensTensor(GreensTensorBase): + """ + FK Green's tensor object + + Overloads base class with machinery for working with CPS-style + Green's functions + + """ + def __init__(self, *args, **kwargs): + super(GreensTensor, self).__init__(*args, **kwargs) + + if 'type:greens' not in self.tags: + self.tags += ['type:greens'] + + if 'type:velocity' not in self.tags: + self.tags += ['type:velocity'] + + if 'units:m' not in self.tags: + self.tags += ['units:m'] + + def _precompute(self): + """ Computes NumPy arrays used by get_synthetics + """ + if self.include_mt: + self._precompute_mt() + + if self.include_force: + self._precompute_force() + + + def _precompute_mt(self): + """ Recombines CPS time series so they can be used in straightforward + linear combination with Mrr,Mtt,Mpp,Mrt,Mrp,Mtp + """ + + array = self._array + phi = np.deg2rad(self.azimuth) + _j = 0 + + # The formulas below were obtained by reverse engineering FK + + for _i, component in enumerate(self.components): + if component=='Z': + ZSS = self.select(channel="ZSS")[0].data + ZDS = self.select(channel="ZDS")[0].data + ZDD = self.select(channel="ZDD")[0].data + ZEP = self.select(channel="ZEP")[0].data + array[_i, _j+0, :] = ZSS/2. * np.cos(2*phi) - ZDD/6. + ZEP/3. + array[_i, _j+1, :] = -ZSS/2. * np.cos(2*phi) - ZDD/6. + ZEP/3. + array[_i, _j+2, :] = ZDD/3. + ZEP/3. + array[_i, _j+3, :] = ZSS * np.sin(2*phi) + array[_i, _j+4, :] = ZDS * np.cos(phi) + array[_i, _j+5, :] = ZDS * np.sin(phi) + + elif component=='R': + RSS = self.select(channel="RSS")[0].data + RDS = self.select(channel="RDS")[0].data + RDD = self.select(channel="RDD")[0].data + REP = self.select(channel="REP")[0].data + array[_i, _j+0, :] = RSS/2. * np.cos(2*phi) - RDD/6. + REP/3. + array[_i, _j+1, :] = -RSS/2. * np.cos(2*phi) - RDD/6. + REP/3. + array[_i, _j+2, :] = RDD/3. + REP/3. + array[_i, _j+3, :] = RSS * np.sin(2*phi) + array[_i, _j+4, :] = RDS * np.cos(phi) + array[_i, _j+5, :] = RDS * np.sin(phi) + + elif component=='T': + TSS = self.select(channel="TSS")[0].data + TDS = self.select(channel="TDS")[0].data + array[_i, _j+0, :] = TSS/2. * np.sin(2*phi) + array[_i, _j+1, :] = -TSS/2. * np.sin(2*phi) + array[_i, _j+2, :] = 0. + array[_i, _j+3, :] = -TSS * np.cos(2*phi) + array[_i, _j+4, :] = TDS * np.sin(phi) + array[_i, _j+5, :] = -TDS * np.cos(phi) + + else: + raise ValueError + + # + # CPS uses a north-east-down basis convention, while mtuq uses an + # up-south-east basis convention, so a permutation is necessary + # + array_copy = array.copy() + array[:, 0, :] = array_copy[:, 2, :] + array[:, 1, :] = array_copy[:, 0, :] + array[:, 2, :] = array_copy[:, 1, :] + array[:, 3, :] = array_copy[:, 4, :] + array[:, 4, :] = -array_copy[:, 5, :] + array[:, 5, :] = -array_copy[:, 3, :] + + + def _precompute_force(self): + raise NotImplementedError() + + + From 9f162f3dba3d3b8c9a52e214d9e8ff496a7af7e1 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Wed, 13 Mar 2024 14:35:21 -0600 Subject: [PATCH 15/31] Typo fixes --- docs/README | 2 +- mtuq/misfit/polarity.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/README b/docs/README index 732d3375c..5e1130515 100644 --- a/docs/README +++ b/docs/README @@ -14,7 +14,7 @@ HOW TO UPDATE DOCS >> ./build.bash -3. If build fails with "Encountered unknown tag 'do', then try the following: +3. If build fails with "Encountered unknown tag 'do'", try the following: Modify the sphinx autosummary source code by adding 'jinja2.ext.do' to the jinja Environment() extensions list diff --git a/mtuq/misfit/polarity.py b/mtuq/misfit/polarity.py index 401284de3..08021f4ed 100644 --- a/mtuq/misfit/polarity.py +++ b/mtuq/misfit/polarity.py @@ -58,9 +58,9 @@ class PolarityMisfit(object): .. rubric:: Other input arguments that may be required, depending on the above ``taup_model`` (`str`): Name of built-in ObsPy TauP model or path to custom - ObsPy TauP model, required for `type=taup` + ObsPy TauP model, required for `method=taup` - ``FK_database`` (`str`): Path to FK database, required for for `type=FK_metadata`. + ``FK_database`` (`str`): Path to FK database, required for for `method=FK_metadata`. .. note:: From 0af5838416c2794c79af3ac75e89d1ca0cfb38fa Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Fri, 15 Mar 2024 15:49:49 -0600 Subject: [PATCH 16/31] Revised documentation --- docs/user_guide/03.rst | 62 ++++++------ docs/user_guide/03/receiver_side.rst | 9 +- docs/user_guide/03/source_side.rst | 120 ++++++++++++++++++++---- docs/user_guide/05/gallery_force.rst | 4 + docs/user_guide/05/gallery_mt.rst | 6 ++ docs/user_guide/06/trace_attributes.rst | 2 +- 6 files changed, 149 insertions(+), 54 deletions(-) diff --git a/docs/user_guide/03.rst b/docs/user_guide/03.rst index 405acc3d1..90b119556 100644 --- a/docs/user_guide/03.rst +++ b/docs/user_guide/03.rst @@ -10,24 +10,21 @@ Role of Green's functions By combining Green's functions, it is possible to obtain synthetics from any moment tensor or force source. Generating synthetics in this way is useful because it provides cost savings compared with on-the-fly wavefield simulations. Synthetics can then be compared with data to determine a best-fitting source. +`GreensTensor` objects +---------------------- -Green's function conventions ----------------------------- - -MTUQ consistently uses an `up-south-east` `convention `_ for internally storing Green's functions (as well as moment tensors and forces). - -A major goal of MTUQ is to avoid exposing users to unnecessary complexity. For Green's functions, MTUQ tries to accomplish this by understanding external conventions and converting to a common internal format. +`GreensTensor `_ objects provide access to a set of Green's functions describing the medium response between a given hypocenter and station. -Under the hood, MTUQ deals with a variety of external conventions related to (1) the symmetry of the medium (for example, 1D media require fewer independent Green's functions than 3D media); (2) the choice of local Cartesian basis convention (for example, some authors employ `up-south-east`, others `north-east-down`). A sense of what's involved can be obtained by browsing the `source code `_ and references therein. +Methods built into the `GreensTensor` class allow data processing and synthetics generation. In particular, the `get_synthetics `_ method accepts a `MomentTensor` or `Force` and returns corresponding synthetics. -`GreensTensor` objects ----------------------- +Green's function conventions +---------------------------- -`GreensTensor `_ objects provide access to a set of Green's functions describing the medium response between a given station and origin. +A major goal of MTUQ is to avoid exposing users to unnecessary complexity. For Green's functions, MTUQ tries to accomplish this by understanding external conventions and converting to a common internal format. Specifically, MTUQ uses an `up-south-east` `convention `_ for internally storing impulse responses corresponding to individual force couples. (Moment tensors and forces are internally represented using the same `up-south-east` basis.) -Methods built into the `GreensTensor` class allow data processing and synthetics generation. In particular, the `get_synthetics `_ method accepts a `MomentTensor` or `Force` and returns corresponding synthetics. +Under the hood, MTUQ deals with a variety of external conventions related to (1) the symmetry of the medium (for example, 1D media require fewer independent Green's functions than 3D media); (2) the choice of local Cartesian basis convention (for example, some authors employ `up-south-east`, others `north-east-down`). A sense of what's involved can be obtained by browsing the `source code `_ and references therein. @@ -97,24 +94,25 @@ To generate a full set of Green's functions for a given hypocenter and depth, si .. code :: {event_id}/ - {net}.{sta}.{loc}.Z.Mrr.sac - {net}.{sta}.{loc}.Z.Mtt.sac - {net}.{sta}.{loc}.Z.Mpp.sac - {net}.{sta}.{loc}.Z.Mrt.sac - {net}.{sta}.{loc}.Z.Mrp.sac - {net}.{sta}.{loc}.Z.Mtp.sac - {net}.{sta}.{loc}.R.Mrr.sac - {net}.{sta}.{loc}.R.Mtt.sac - {net}.{sta}.{loc}.R.Mpp.sac - {net}.{sta}.{loc}.R.Mrt.sac - {net}.{sta}.{loc}.R.Mrp.sac - {net}.{sta}.{loc}.R.Mtp.sac - {net}.{sta}.{loc}.T.Mrr.sac - {net}.{sta}.{loc}.T.Mtt.sac - {net}.{sta}.{loc}.T.Mpp.sac - {net}.{sta}.{loc}.T.Mrt.sac - {net}.{sta}.{loc}.T.Mrp.sac - {net}.{sta}.{loc}.T.Mtp.sac + {depth_in_m}/ + {net}.{sta}.{loc}.Z.Mrr.sac + {net}.{sta}.{loc}.Z.Mtt.sac + {net}.{sta}.{loc}.Z.Mpp.sac + {net}.{sta}.{loc}.Z.Mrt.sac + {net}.{sta}.{loc}.Z.Mrp.sac + {net}.{sta}.{loc}.Z.Mtp.sac + {net}.{sta}.{loc}.R.Mrr.sac + {net}.{sta}.{loc}.R.Mtt.sac + {net}.{sta}.{loc}.R.Mpp.sac + {net}.{sta}.{loc}.R.Mrt.sac + {net}.{sta}.{loc}.R.Mrp.sac + {net}.{sta}.{loc}.R.Mtp.sac + {net}.{sta}.{loc}.T.Mrr.sac + {net}.{sta}.{loc}.T.Mtt.sac + {net}.{sta}.{loc}.T.Mpp.sac + {net}.{sta}.{loc}.T.Mrt.sac + {net}.{sta}.{loc}.T.Mrp.sac + {net}.{sta}.{loc}.T.Mtp.sac To open a SPECFEM3D/3D_GLOBE database client: @@ -132,9 +130,9 @@ Once opened, the database client can be used to generate `GreensTensor `_ +`Source-side Green's function details (under construction) `_ -`Receiver-side SPECFEM3D/3D_GLOBE Green's functions (under construction) `_ +`Receiver-side Green's function details (under construction) `_ diff --git a/docs/user_guide/03/receiver_side.rst b/docs/user_guide/03/receiver_side.rst index 10994583b..c13fdc2c6 100644 --- a/docs/user_guide/03/receiver_side.rst +++ b/docs/user_guide/03/receiver_side.rst @@ -1,13 +1,12 @@ .. warning:: - This page is currently just a stub. + This page is still under construction. To help improve the + documentation, feel free to submit a pull request. - To fill in missing documentation, feel free to submit a pull request. - -Notes on receiver-side 3D Green's functions -=========================================== +Receiver-side 3D Green's functions +================================== A working of example using receiver-side 3D Green's functions from SPECFEM3D in a moment tensor inversion: diff --git a/docs/user_guide/03/source_side.rst b/docs/user_guide/03/source_side.rst index df27bdaf9..b43c17f12 100644 --- a/docs/user_guide/03/source_side.rst +++ b/docs/user_guide/03/source_side.rst @@ -1,49 +1,137 @@ .. warning:: - This page is currently just a stub. + This page is still under construction. To help improve the + documentation, feel free to submit a pull request. - To fill in missing documentation, feel free to submit a pull request. +Source-side 3D Green's functions +================================ +Generating source-side 3D Green's functions using SPECFEM3D/3D_GLOBE +-------------------------------------------------------------------- -Notes on source-side 3D Green's functions -========================================= +In principle, any 3D solver can be used to generate Green's functions, as long as the `requirements `_ below are satisfied. + +So far, however, the source-side Green's function machinery has been tested using only SPECFEM3D/3D_GLOBE. To convert SPECFEM3D/3D_GLOBE output to MTUQ-compliant Green's functions, the following steps are necessary. + + +**Generate SAC binary files** + +SAC binary output format is natively supported by SPECFEM3D_GLOBE (in the parameter file, set `OUTPUT_SEISMOS_SAC_BINARY = .true.`). + +Unfortunately, SAC binary output format is not natively supported by SPECFEM3D, so it becomes necessary to manually convert SPECFEM3D output to SAC binary format. + + +**Convert to SI units** + +MTUQ uses the fully SI convention described below. In contrast, SPECFEM3D/3D_GLOBE use a mixed SI and CGS convention, in which moment tensor elements are input in terms of dynes and centimeters, and seismograms are output in meters. As a result, it is necessary to scale SPECFEM3D/3D_GLOBE seismograms by 10^7 prior to using them as MTUQ Green's functions. + + +**Additional amplitude scaling** + +In addition to converting to SI units, it is also necessary to account for any scaling factors in the SPECFEM3D/3D_GLOBE input files. Such scaling factors can enter, for example, through the `M_rr,M_tt,M_pp,M_rt,M_rp,M_tp` values in the moment tensor input file or through the `scaling factor `_ in the force input file. + + +**Rotate to vertical, radial, transverse components** + +Conveniently, SPECFEM3D_GLOBE can be made to automatically rotate output seismograms into vertical, radial and transverse components (set `ROTATE_SEISMOGRAMS_RT = .true.` in the parameter file). + +No modifications are necessary on account of moment tensor basis convention, since MTUQ's `up-south-east` convention matches SPECFEM3D/3D_GLOBE's. + + + + +Requirements for MTUQ source-side 3D Green's functions +------------------------------------------------------ **File format** -- Currently, MTUQ reads source-side 3D Green's functions from SAC files +Individual Green's functions must be written to SAC binary files. + +A total 18 SAC binary files are required to represent the response between a given hypocenter and station (corresponding to 6 moment tensor elements times 3 directions of motion). -- SAC output format is natively supported by SPECFEM3D_GLOBE, but not SPECFEM3D Cartesian (output from the latter can be manually converted) **Units convention** -- MTUQ treats each individual SAC file as the response in meters to a unit (1 Newton-meter) force couple +For a moment tensor inversion, each SAC binary file must give the response in meters to a 1 Newton-meter force couple. + +For a force inversion, each SAC binary file must give the response in meters to a 1 Newton force. + +In both cases, MTUQ uses a fully SI units convention (compare with SPECFEM3D/3D_GLOBE notes, above). -- Users must ensure that source-side Green's functions are normalized according to the above units **Basis convention** -- For moment tensors, MTUQ uses an `up-south-east` convention, identical to one the used by SPECFEM3D_GLOBE +MTUQ uses an `up-south-east` basis convention in which `r` denotes up, `t` denotes south, and `p` denotes east. + +Green's functions must be rotated into into vertical `Z`, radial `R` and transverse `T` components relative to the source-receiver backazimuth. + +Place all seismograms for the same hypocenter in a single directory as follows: + +.. code :: + + {event_id}/ + {depth_in_km}/ + {net}.{sta}.{loc}.Z.Mrr.sac + {net}.{sta}.{loc}.Z.Mtt.sac + {net}.{sta}.{loc}.Z.Mpp.sac + {net}.{sta}.{loc}.Z.Mrt.sac + {net}.{sta}.{loc}.Z.Mrp.sac + {net}.{sta}.{loc}.Z.Mtp.sac + {net}.{sta}.{loc}.R.Mrr.sac + {net}.{sta}.{loc}.R.Mtt.sac + {net}.{sta}.{loc}.R.Mpp.sac + {net}.{sta}.{loc}.R.Mrt.sac + {net}.{sta}.{loc}.R.Mrp.sac + {net}.{sta}.{loc}.R.Mtp.sac + {net}.{sta}.{loc}.T.Mrr.sac + {net}.{sta}.{loc}.T.Mtt.sac + {net}.{sta}.{loc}.T.Mpp.sac + {net}.{sta}.{loc}.T.Mrt.sac + {net}.{sta}.{loc}.T.Mrp.sac + {net}.{sta}.{loc}.T.Mtp.sac + ... + +The corresponding convention for force responses is: + +.. code :: + + {event_id}/ + {depth_in_km}/ + {net}.{sta}.{loc}.Z.Fr.sac + {net}.{sta}.{loc}.Z.Ft.sac + {net}.{sta}.{loc}.Z.Fp.sac + {net}.{sta}.{loc}.R.Fr.sac + {net}.{sta}.{loc}.R.Ft.sac + {net}.{sta}.{loc}.R.Fp.sac + {net}.{sta}.{loc}.T.Fr.sac + {net}.{sta}.{loc}.T.Fp.sac + {net}.{sta}.{loc}.T.Fp.sac + ... **Origin time convention** -- For origin time, MTUQ uses a centroid convention (`more details `_), so that `t=0` in the `GreensTensor` time discretization corresponds to mean source excitation time +For origin time, MTUQ uses a centroid convention (`more details `_), so that `t=0` in the `GreensTensor` time discretization corresponds to mean source excitation time. + +MTUQ uses the begin time (`B`) and end time (`E`) headers from the SAC binary files to align the Green's functions relative to centroid origin time. + +Currently, these are the only SAC headers used in reading Green's functions. -- MTUQ uses the `REF_TIME` header from the SPECFEM3D_GLOBE SAC output files, which gives the peak excitation of the source relative to the simulation start time +(Note that different `SAC headers `_ are required in reading `observed data `_.) -- MTUQ ignores the origin time given in the CMTSOLUTION file and `ORIGIN_TIME` header -**Depth searches (experimental)** +Hypocenter searches (experimental) +---------------------------------- -- Only depth searches are possible with source-side 3D Green's functions (no other hypocenter parameters) +Currently, only searches over source depth are possible with source-side 3D Green's functions (no other hypocenter parameters). -- The current `depth search `_ implementation is somewhat crude and experimental; consider local modifcations to suit your needs +The current `depth search `_ implementation is especially crude and experimental (consider local modifications to suit your needs). -- To allow depth searches, create subdirectories for each centroid depth, as below +To allow depth searches, create subdirectories for each centroid depth as follows: .. code :: diff --git a/docs/user_guide/05/gallery_force.rst b/docs/user_guide/05/gallery_force.rst index deaa33576..1383e0757 100644 --- a/docs/user_guide/05/gallery_force.rst +++ b/docs/user_guide/05/gallery_force.rst @@ -33,6 +33,10 @@ For a grid of regulary-spaced forces, `ds` may look something like: * origin_idx (origin_idx) int64 0 +.. note:: + + Force grids are implemented using parameters `F0, phi, h`, which are related to `r, phi, theta` spherical coordinates (`physics convention `_) by `F0 = r`, `phi = phi`, `h = cos(theta)`. In addition, `F0, phi, h` are related to geographic directions by these `formulas `_. + Misfit values """"""""""""" diff --git a/docs/user_guide/05/gallery_mt.rst b/docs/user_guide/05/gallery_mt.rst index 8b57c8d84..ff9880e58 100644 --- a/docs/user_guide/05/gallery_mt.rst +++ b/docs/user_guide/05/gallery_mt.rst @@ -37,6 +37,12 @@ For a grid of regulary-spaced moment tensors, `ds` may look something like: +.. note:: + + Moment tensor grids are implemented using the `rho, v, w, kappa, sigma, h` parameterization of `Tape2012 `_ and `Tape2015 `_, from which `formulas `_ converting to parameterizations can be derived. + + + Misfit values """"""""""""" diff --git a/docs/user_guide/06/trace_attributes.rst b/docs/user_guide/06/trace_attributes.rst index 41014330e..5e267d2f9 100644 --- a/docs/user_guide/06/trace_attributes.rst +++ b/docs/user_guide/06/trace_attributes.rst @@ -2,7 +2,7 @@ Trace attributes ================ -At various points during an inversion, waveform differences, phase shifts, and other values are calculated from observed and synthetic seismic traces. Such `trace attribute` quantities provide important information about how data misfit varies by geographic location and seismic component. +Waveform differences, time-shift corrections, and other values are calculated during an inversion on a trace-by-trace basis. Such `trace attributes` provide important information about how data misfit varies by geographic location and seismic component. Collecting trace attributes From 3d33d54be4732a4b731767719187dd38756a4741 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Fri, 15 Mar 2024 15:53:48 -0600 Subject: [PATCH 17/31] Updated check_urls.bash --- tests/check_urls.bash | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/check_urls.bash b/tests/check_urls.bash index 9bef6b62d..f712d96b5 100755 --- a/tests/check_urls.bash +++ b/tests/check_urls.bash @@ -10,8 +10,9 @@ URLS="\ https://www.eas.slu.edu/People/LZhu/home.html\ https://github.com/geodynamics/axisem\ https://github.com/Liang-Ding/seisgen\ - http://ds.iris.edu/ds/products/syngine\ - http://ds.iris.edu/ds/products/syngine/#models\ + https://ds.iris.edu/ds/products/syngine\ + https://ds.iris.edu/ds/products/syngine/#models\ + https://ds.iris.edu/files/sac-manual/manual/file_format.html\ https://instaseis.net\ https://docs.obspy.org/tutorial/index.html\ https://docs.obspy.org/packages/autogen/obspy.core.stream.Stream.html\ From 55ea83a0903d34bdc97f9a8c8e9511d63cf7e7a2 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Thu, 21 Mar 2024 16:37:56 -0600 Subject: [PATCH 18/31] Tweaked get_synthetics() method --- mtuq/greens_tensor/base.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/mtuq/greens_tensor/base.py b/mtuq/greens_tensor/base.py index 3ec868d4a..baace937d 100644 --- a/mtuq/greens_tensor/base.py +++ b/mtuq/greens_tensor/base.py @@ -6,6 +6,7 @@ from mtuq.event import Origin from mtuq.station import Station from mtuq.dataset import Dataset +from mtuq.util import Null from mtuq.util.signal import check_time_sampling from obspy.core import Stream, Trace from obspy.geodetics import gps2dist_azimuth @@ -176,6 +177,9 @@ def get_synthetics(self, source, components=None, stats=None, inplace=False): List containing zero or more of the following components: ``Z``, ``R``, ``T``. (Defaults to ``['Z', 'R', 'T']``.) + ``stats`` (`obspy.Trace.Stats` object): + ObsPy Stats object that will be attached to the synthetics + """ if components is None: @@ -300,9 +304,17 @@ def get_synthetics(self, source, components=None, stats=None, mode='apply', **kw ``components`` (`list`): List containing zero or more of the following components: ``Z``, ``R``, ``T``. (Defaults to ``['Z', 'R', 'T']``.) + + ``stats`` (`obspy.Trace.Stats` object): + ObsPy Stats object that will be attached to the synthetics """ if mode=='map': + if components is None: + components = [None for _ in range(len(self))] + if stats is None: + stats = [None for _ in range(len(self))] + synthetics = Dataset() for _i, tensor in enumerate(self): synthetics.append(tensor.get_synthetics( From fb8ca068cfd2594caa6cf6a1e08fe0ee47c5f4a3 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Thu, 21 Mar 2024 17:13:57 -0600 Subject: [PATCH 19/31] Improved misfit.waveform._stats --- mtuq/misfit/waveform/_stats.py | 61 ++++++++++++++++++++++++++++------ 1 file changed, 50 insertions(+), 11 deletions(-) diff --git a/mtuq/misfit/waveform/_stats.py b/mtuq/misfit/waveform/_stats.py index 067f503f9..b31db475d 100644 --- a/mtuq/misfit/waveform/_stats.py +++ b/mtuq/misfit/waveform/_stats.py @@ -7,36 +7,62 @@ from mtuq.util.signal import get_components -def calculate_norm_data(data, norm, components): - # error checking - assert norm in ('L1', 'L2') +def calculate_norm_data(data, norm, components, apply_weights=True): + """ Calculates the norm of the given waveform data + """ + + # Here we try to mimic the implementation in + # mtuq/misfit/waveform/level0.py + + # Mimicking the misfit function implementation helps ensure consistency + # between the misfit values (i.e. residual norms) and data norms + + # The idea here is to use the regular misfit function machinery, but set + # synthetics to zero, so that instead of the residual norm || d - s ||, + # we get the data norm || d - 0 || norm_data = 0. - for _j, d in enumerate(data): + + for _j, stream in enumerate(data): _components, indices = list_intersect_with_indices( - get_components(d), components) + get_components(stream), _flatten(components)) if not indices: continue # time sampling scheme - npts = d[0].data.size - dt = d[0].stats.delta + npts = stream[0].data.size + dt = stream[0].stats.delta for _k in indices: - r = d[_k].data + d = stream[_k].data if norm=='L1': - norm_data += np.sum(np.abs(r))*dt + value = np.sum(np.abs(d))*dt elif norm=='L2': - norm_data += np.sum(r**2)*dt + value = np.sum(d**2)*dt + + elif norm=='hybrid': + value = np.sqrt(np.sum(r**2))*dt + + # optionally, applies user-supplied weights attached during + # process_data() + if apply_weights: + try: + value *= d[_k].weight + except: + pass + + norm_data += value return norm_data def estimate_sigma(data, greens, best_source, norm, components, time_shift_min, time_shift_max): + """ Makes an a posteriori estimate of the data error standard deviation + """ # error checking assert norm in ('L1', 'L2') @@ -81,7 +107,6 @@ def estimate_sigma(data, greens, best_source, norm, components, stop = start + npts for _k in indices: - # subtract data from shifted synthetics r = s[_k].data[start:stop] - d[_k].data @@ -95,3 +120,17 @@ def estimate_sigma(data, greens, best_source, norm, components, return np.mean(residuals)**0.5 + + +def _flatten(lists): + # TODO - better implementation? + + # example: + # converts ['ZR','T'] to ['Z','R','T'] + + flattened = [] + for list in lists: + for element in list: + flattened.append(element) + return flattened + From ffe5058c8ca5ac77bd6a65f60a7155bb30c7db3b Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Thu, 21 Mar 2024 18:16:30 -0600 Subject: [PATCH 20/31] Normalization option for waveform misfit functions --- mtuq/misfit/waveform/__init__.py | 14 ++++++++------ mtuq/misfit/waveform/level0.py | 12 +++++++++++- mtuq/misfit/waveform/level1.py | 10 +++++++++- mtuq/misfit/waveform/level2.py | 12 +++++++++++- 4 files changed, 39 insertions(+), 9 deletions(-) diff --git a/mtuq/misfit/waveform/__init__.py b/mtuq/misfit/waveform/__init__.py index cfb2de468..5b55640b8 100644 --- a/mtuq/misfit/waveform/__init__.py +++ b/mtuq/misfit/waveform/__init__.py @@ -162,7 +162,7 @@ def __init__(self, def __call__(self, data, greens, sources, progress_handle=Null(), - set_attributes=False, optimization_level=None): + normalize=False, set_attributes=False, optimization_level=None): """ Evaluates misfit on given data """ if optimization_level is None: @@ -195,17 +195,19 @@ def __call__(self, data, greens, sources, progress_handle=Null(), return level0.misfit( data, greens, sources, self.norm, self.time_shift_groups, self.time_shift_min, self.time_shift_max, progress_handle, - set_attributes) + normalize=normalize, set_attributes=set_attributes) if optimization_level==1: return level1.misfit( data, greens, sources, self.norm, self.time_shift_groups, - self.time_shift_min, self.time_shift_max, progress_handle) + self.time_shift_min, self.time_shift_max, progress_handle, + normalize=normalize) if optimization_level==2: return level2.misfit( data, greens, sources, self.norm, self.time_shift_groups, - self.time_shift_min, self.time_shift_max, progress_handle) + self.time_shift_min, self.time_shift_max, progress_handle, + normalize=normalize) def collect_attributes(self, data, greens, source): @@ -229,7 +231,7 @@ def collect_attributes(self, data, greens, source): _ = level0.misfit( data, greens, iterable(source), self.norm, self.time_shift_groups, self.time_shift_min, self.time_shift_max, msg_handle=Null(), - set_attributes=True) + normalize=normalize, set_attributes=True) # collects attributes attrs = [] @@ -265,7 +267,7 @@ def collect_synthetics(self, data, greens, source): _ = level0.misfit( data, greens, iterable(source), self.norm, self.time_shift_groups, self.time_shift_min, self.time_shift_max, msg_handle=Null(), - set_attributes=True) + normalize=normalize, set_attributes=True) return deepcopy(synthetics) diff --git a/mtuq/misfit/waveform/level0.py b/mtuq/misfit/waveform/level0.py index e22123c82..2469fa2e7 100644 --- a/mtuq/misfit/waveform/level0.py +++ b/mtuq/misfit/waveform/level0.py @@ -6,13 +6,15 @@ import numpy as np +from mtuq.misfit.waveform._stats import _flatten, calculate_norm_data from mtuq.util import AttribDict from mtuq.util.math import isclose, list_intersect_with_indices from mtuq.util.signal import get_components def misfit(data, greens, sources, norm, time_shift_groups, - time_shift_min, time_shift_max, msg_handle, set_attributes=False): + time_shift_min, time_shift_max, msg_handle, + normalize=False, set_attributes=False): """ Waveform misfit function (non-optimized pure Python version) @@ -20,6 +22,10 @@ def misfit(data, greens, sources, norm, time_shift_groups, """ values = np.zeros((len(sources), 1)) + if normalize: + components = _flatten(time_shift_groups) + norm_data = calculate_norm_data(data, norm, components) + # # initialize Green's function machinery # @@ -164,6 +170,10 @@ def misfit(data, greens, sources, norm, time_shift_groups, s[_k].attrs.amplitude_ratio = d_max/s_max s[_k].attrs.log_amplitude_ratio = np.log(d_max/s_max) + + if normalize: + values /= norm_data + return values diff --git a/mtuq/misfit/waveform/level1.py b/mtuq/misfit/waveform/level1.py index 4a0f70813..fdab3ddf7 100644 --- a/mtuq/misfit/waveform/level1.py +++ b/mtuq/misfit/waveform/level1.py @@ -5,13 +5,14 @@ """ import numpy as np +from mtuq.misfit.waveform._stats import _flatten, calculate_norm_data from mtuq.util.math import correlate, isclose, list_intersect_with_indices from mtuq.util.signal import get_components, get_time_sampling def misfit(data, greens, sources, norm, time_shift_groups, - time_shift_min, time_shift_max, msg_handle): + time_shift_min, time_shift_max, msg_handle, normalize=False): """ Waveform misfit function (fast pure Python version) @@ -20,6 +21,10 @@ def misfit(data, greens, sources, norm, time_shift_groups, helpers = [] values = np.zeros((len(sources), 1)) + if normalize: + components = _flatten(time_shift_groups) + norm_data = calculate_norm_data(data, norm, components) + # # initialize Green's function machinery # @@ -84,6 +89,9 @@ def misfit(data, greens, sources, norm, time_shift_groups, except: values[_i] += value + if normalize: + values /= norm_data + return values diff --git a/mtuq/misfit/waveform/level2.py b/mtuq/misfit/waveform/level2.py index 229120fac..c99839229 100644 --- a/mtuq/misfit/waveform/level2.py +++ b/mtuq/misfit/waveform/level2.py @@ -7,6 +7,7 @@ import numpy as np import time from copy import deepcopy +from mtuq.misfit.waveform._stats import _flatten, calculate_norm_data from mtuq.misfit.waveform.level1 import correlate from mtuq.util.math import to_mij, to_rtp from mtuq.util.signal import get_components, get_time_sampling @@ -14,12 +15,18 @@ def misfit(data, greens, sources, norm, time_shift_groups, - time_shift_min, time_shift_max, msg_handle, debug_level=0): + time_shift_min, time_shift_max, msg_handle, debug_level=0, + normalize=False): """ Data misfit function (fast Python/C version) See ``mtuq/misfit/waveform/__init__.py`` for more information """ + + if normalize: + components = _flatten(time_shift_groups) + norm_data = calculate_norm_data(data, norm, components) + # # collect metadata # @@ -86,6 +93,9 @@ def misfit(data, greens, sources, norm, time_shift_groups, print(' Elapsed time (C extension) (s): %f' % \ (time.time() - start_time)) + if normalize: + results /= norm_data + return results From 37fa4d34f0b21190d422566707cfa0643b58d401 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Thu, 21 Mar 2024 18:20:20 -0600 Subject: [PATCH 21/31] Tweaked utilities --- mtuq/util/__init__.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/mtuq/util/__init__.py b/mtuq/util/__init__.py index b9595f89a..9f66cf6bd 100644 --- a/mtuq/util/__init__.py +++ b/mtuq/util/__init__.py @@ -287,6 +287,9 @@ def __init__(self, *args, **kwargs): def __call__(self, *args, **kwargs): return self + def __getitem__(self, *args, **kwargs): + return None + def __nonzero__(self): return False From 378113fc2672692baa78a387e9717e15c5a59f91 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Thu, 21 Mar 2024 18:22:20 -0600 Subject: [PATCH 22/31] Docs revisions --- docs/user_guide/03.rst | 2 +- docs/user_guide/03/source_side.rst | 61 ++++++++++++++++++------------ 2 files changed, 38 insertions(+), 25 deletions(-) diff --git a/docs/user_guide/03.rst b/docs/user_guide/03.rst index 90b119556..0a73c94bd 100644 --- a/docs/user_guide/03.rst +++ b/docs/user_guide/03.rst @@ -94,7 +94,7 @@ To generate a full set of Green's functions for a given hypocenter and depth, si .. code :: {event_id}/ - {depth_in_m}/ + {depth_in_km}/ {net}.{sta}.{loc}.Z.Mrr.sac {net}.{sta}.{loc}.Z.Mtt.sac {net}.{sta}.{loc}.Z.Mpp.sac diff --git a/docs/user_guide/03/source_side.rst b/docs/user_guide/03/source_side.rst index b43c17f12..bef15939b 100644 --- a/docs/user_guide/03/source_side.rst +++ b/docs/user_guide/03/source_side.rst @@ -4,35 +4,40 @@ This page is still under construction. To help improve the documentation, feel free to submit a pull request. + Source-side 3D Green's functions ================================ Generating source-side 3D Green's functions using SPECFEM3D/3D_GLOBE -------------------------------------------------------------------- -In principle, any 3D solver can be used to generate Green's functions, as long as the `requirements `_ below are satisfied. +In principle, any 3D solver can be used to generate source-side Green's functions, as long as the `requirements `_ below are satisfied. -So far, however, the source-side Green's function machinery has been tested using only SPECFEM3D/3D_GLOBE. To convert SPECFEM3D/3D_GLOBE output to MTUQ-compliant Green's functions, the following steps are necessary. +So far, however, the machinery has only been tested using SPECFEM3D/3D_GLOBE. To convert SPECFEM3D/3D_GLOBE output to MTUQ-compliant Green's functions, the following steps are necessary. -**Generate SAC binary files** +Generate SAC binary files +......................... SAC binary output format is natively supported by SPECFEM3D_GLOBE (in the parameter file, set `OUTPUT_SEISMOS_SAC_BINARY = .true.`). -Unfortunately, SAC binary output format is not natively supported by SPECFEM3D, so it becomes necessary to manually convert SPECFEM3D output to SAC binary format. +Unfortunately, SAC binary output format is not natively supported by SPECFEM3D, so it is necessary to manually convert SPECFEM3D output. -**Convert to SI units** +Convert to SI units +................... -MTUQ uses the fully SI convention described below. In contrast, SPECFEM3D/3D_GLOBE use a mixed SI and CGS convention, in which moment tensor elements are input in terms of dynes and centimeters, and seismograms are output in meters. As a result, it is necessary to scale SPECFEM3D/3D_GLOBE seismograms by 10^7 prior to using them as MTUQ Green's functions. +SPECFEM3D/3D_GLOBE use a mixed SI and CGS convention, in which moment tensor elements are input in terms of dynes and centimeters, and seismograms are output in meters. In contrast, MTUQ uses a fully SI :ref:`convention`. As a result, it is necessary to scale SPECFEM3D/3D_GLOBE seismograms by 10^7 prior to using them as MTUQ Green's functions. -**Additional amplitude scaling** +Other amplitude scaling +....................... In addition to converting to SI units, it is also necessary to account for any scaling factors in the SPECFEM3D/3D_GLOBE input files. Such scaling factors can enter, for example, through the `M_rr,M_tt,M_pp,M_rt,M_rp,M_tp` values in the moment tensor input file or through the `scaling factor `_ in the force input file. -**Rotate to vertical, radial, transverse components** +Rotate traces +............. Conveniently, SPECFEM3D_GLOBE can be made to automatically rotate output seismograms into vertical, radial and transverse components (set `ROTATE_SEISMOGRAMS_RT = .true.` in the parameter file). @@ -44,15 +49,16 @@ No modifications are necessary on account of moment tensor basis convention, sin Requirements for MTUQ source-side 3D Green's functions ------------------------------------------------------ -**File format** +File format +........... Individual Green's functions must be written to SAC binary files. -A total 18 SAC binary files are required to represent the response between a given hypocenter and station (corresponding to 6 moment tensor elements times 3 directions of motion). - +A total of 18 SAC binary files are required to represent the response between a given hypocenter and station (corresponding to 6 moment tensor elements times 3 directions of motion). -**Units convention** +Units convention +................ For a moment tensor inversion, each SAC binary file must give the response in meters to a 1 Newton-meter force couple. @@ -62,9 +68,10 @@ In both cases, MTUQ uses a fully SI units convention (compare with SPECFEM3D/3D_ -**Basis convention** +Basis convention +................ -MTUQ uses an `up-south-east` basis convention in which `r` denotes up, `t` denotes south, and `p` denotes east. +MTUQ uses an `up-south-east` basis convention in which `r` denotes vertically upward, `t` denotes south, and `p` denotes east. Green's functions must be rotated into into vertical `Z`, radial `R` and transverse `T` components relative to the source-receiver backazimuth. @@ -96,23 +103,29 @@ Place all seismograms for the same hypocenter in a single directory as follows: The corresponding convention for force responses is: +.. warning:: + + Not yet tested; subject to change without notice. + + .. code :: {event_id}/ {depth_in_km}/ - {net}.{sta}.{loc}.Z.Fr.sac - {net}.{sta}.{loc}.Z.Ft.sac - {net}.{sta}.{loc}.Z.Fp.sac - {net}.{sta}.{loc}.R.Fr.sac - {net}.{sta}.{loc}.R.Ft.sac - {net}.{sta}.{loc}.R.Fp.sac - {net}.{sta}.{loc}.T.Fr.sac - {net}.{sta}.{loc}.T.Fp.sac - {net}.{sta}.{loc}.T.Fp.sac + {net}.{sta}.{loc}.Z.Fz.sac + {net}.{sta}.{loc}.Z.Fs.sac + {net}.{sta}.{loc}.Z.Fe.sac + {net}.{sta}.{loc}.R.Fz.sac + {net}.{sta}.{loc}.R.Fs.sac + {net}.{sta}.{loc}.R.Fe.sac + {net}.{sta}.{loc}.T.Fz.sac + {net}.{sta}.{loc}.T.Fs.sac + {net}.{sta}.{loc}.T.Fe.sac ... -**Origin time convention** +Origin time convention +...................... For origin time, MTUQ uses a centroid convention (`more details `_), so that `t=0` in the `GreensTensor` time discretization corresponds to mean source excitation time. From 5693463d047b2cfaef899833c12eb8c3795422a5 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Thu, 21 Mar 2024 18:25:55 -0600 Subject: [PATCH 23/31] Tweaked examples --- setup/code_generator.py | 1 + tests/test_greens_SPECFEM3D_SAC.py | 1 + tests/test_greens_SPECFEM3D_SGT.py | 1 + tests/test_grid_search_mt.py | 1 + 4 files changed, 4 insertions(+) diff --git a/setup/code_generator.py b/setup/code_generator.py index 1a9e76009..1198ce31b 100755 --- a/setup/code_generator.py +++ b/setup/code_generator.py @@ -1588,6 +1588,7 @@ def isclose(a, b, atol=1.e6, rtol=1.e-6): print('%s: %.e <= %.1e + %.1e * %.1e' %\\ ('passed' if _bool else 'failed', abs(_a-_b), atol, rtol, abs(_b))) + print('') return np.all( diff --git a/tests/test_greens_SPECFEM3D_SAC.py b/tests/test_greens_SPECFEM3D_SAC.py index 14aa3c091..7b2734270 100644 --- a/tests/test_greens_SPECFEM3D_SAC.py +++ b/tests/test_greens_SPECFEM3D_SAC.py @@ -169,6 +169,7 @@ def isclose(a, b, atol=1.e6, rtol=1.e-6): print('%s: %.e <= %.1e + %.1e * %.1e' %\ ('passed' if _bool else 'failed', abs(_a-_b), atol, rtol, abs(_b))) + print('') return np.all( diff --git a/tests/test_greens_SPECFEM3D_SGT.py b/tests/test_greens_SPECFEM3D_SGT.py index f7584642f..2bfd7c963 100644 --- a/tests/test_greens_SPECFEM3D_SGT.py +++ b/tests/test_greens_SPECFEM3D_SGT.py @@ -170,6 +170,7 @@ def isclose(a, b, atol=1.e6, rtol=1.e-6): print('%s: %.e <= %.1e + %.1e * %.1e' %\ ('passed' if _bool else 'failed', abs(_a-_b), atol, rtol, abs(_b))) + print('') return np.all( diff --git a/tests/test_grid_search_mt.py b/tests/test_grid_search_mt.py index c2c4fda1b..74ba6951f 100644 --- a/tests/test_grid_search_mt.py +++ b/tests/test_grid_search_mt.py @@ -179,6 +179,7 @@ def isclose(a, b, atol=1.e6, rtol=1.e-6): print('%s: %.e <= %.1e + %.1e * %.1e' %\ ('passed' if _bool else 'failed', abs(_a-_b), atol, rtol, abs(_b))) + print('') return np.all( From fe49ce6445f12c4dc9cc49a0e390787398b82da9 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Thu, 21 Mar 2024 23:41:03 -0600 Subject: [PATCH 24/31] Fixed GridSearch.Force.py output filenames --- examples/GridSearch.Force.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/GridSearch.Force.py b/examples/GridSearch.Force.py index 15cc1f649..43489c2ab 100755 --- a/examples/GridSearch.Force.py +++ b/examples/GridSearch.Force.py @@ -189,17 +189,17 @@ print('Generating figures...\n') - plot_data_greens1(event_id+'FMT_waveforms1.png', + plot_data_greens1(event_id+'_force_waveforms.png', data_sw, greens_sw, process_sw, misfit_sw, stations, origin, best_force, direction_dict) - plot_misfit_force(event_id+'DC_misfit_force.png', results) + plot_misfit_force(event_id+'_force_misfit.png', results) print('Saving results...\n') # save best-fitting source - save_json(event_id+'Force_solution.json', merged_dict) + save_json(event_id+'_force_solution.json', merged_dict) # save misfit surface - results.save(event_id+'DC_misfit.nc') + results.save(event_id+'_force_misfit.nc') - print('\nFinished\n') \ No newline at end of file + print('\nFinished\n') From 9c53730e65ef83f4deb18a1a55efd642749a98ee Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Thu, 21 Mar 2024 23:49:14 -0600 Subject: [PATCH 25/31] Fixed misfit function normalize argument --- mtuq/misfit/waveform/__init__.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/mtuq/misfit/waveform/__init__.py b/mtuq/misfit/waveform/__init__.py index 5b55640b8..4070ad68e 100644 --- a/mtuq/misfit/waveform/__init__.py +++ b/mtuq/misfit/waveform/__init__.py @@ -162,7 +162,7 @@ def __init__(self, def __call__(self, data, greens, sources, progress_handle=Null(), - normalize=False, set_attributes=False, optimization_level=None): + normalize=None, set_attributes=False, optimization_level=None): """ Evaluates misfit on given data """ if optimization_level is None: @@ -210,7 +210,7 @@ def __call__(self, data, greens, sources, progress_handle=Null(), normalize=normalize) - def collect_attributes(self, data, greens, source): + def collect_attributes(self, data, greens, source, normalize=False): """ Collects misfit, time shifts and other attributes corresponding to each trace """ @@ -247,7 +247,7 @@ def collect_attributes(self, data, greens, source): return deepcopy(attrs) - def collect_synthetics(self, data, greens, source): + def collect_synthetics(self, data, greens, source, normalize=False): """ Collects synthetics with misfit, time shifts and other attributes attached """ # checks that dataset is nonempty @@ -267,7 +267,7 @@ def collect_synthetics(self, data, greens, source): _ = level0.misfit( data, greens, iterable(source), self.norm, self.time_shift_groups, self.time_shift_min, self.time_shift_max, msg_handle=Null(), - normalize=normalize, set_attributes=True) + normalize=False, set_attributes=True) return deepcopy(synthetics) From 30496a16b6533574ad1758be91a90f39c9d9d6b2 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Fri, 22 Mar 2024 18:19:34 -0600 Subject: [PATCH 26/31] Graphics updates --- mtuq/graphics/header.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/mtuq/graphics/header.py b/mtuq/graphics/header.py index 14c305223..c56231e3a 100644 --- a/mtuq/graphics/header.py +++ b/mtuq/graphics/header.py @@ -15,7 +15,7 @@ class Base(object): - """ Base class for storing and writing text to a matplotlib figure + """ Base class for writing headers to matplotlib figures """ def __init__(self): raise NotImplementedError("Must be implemented by subclass") @@ -24,6 +24,9 @@ def __init__(self): def _get_axis(self, height, fig=None): """ Returns matplotlib axes of given height along top of figure """ + if hasattr(self, '_axis'): + return self._axis + if fig is None: fig = pyplot.gcf() width, figure_height = fig.get_size_inches() @@ -35,7 +38,9 @@ def _get_axis(self, height, fig=None): x0 = 0. y0 = 1.-height/figure_height - ax = fig.add_axes([x0, y0, 1., height/figure_height]) + self._axis = fig.add_axes([x0, y0, 1., height/figure_height]) + ax = self._axis + ax.set_xlim([0., width]) ax.set_ylim([0., height]) @@ -54,9 +59,10 @@ def write(self, *args, **kwargs): raise NotImplementedError("Must be implemented by subclass") - class TextHeader(Base): - """ Prints header text from a list ((xp, yp, text), ...) + """ Generic text header + + Prints header text from a list ((xp, yp, text), ...) """ def __init__(self, items): # validates items From eb7e5adaf835e2cc2ce183a7a904be471a7e35c5 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Fri, 22 Mar 2024 18:20:13 -0600 Subject: [PATCH 27/31] New saturation option (turned off by default) --- mtuq/graphics/uq/_gmt.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/mtuq/graphics/uq/_gmt.py b/mtuq/graphics/uq/_gmt.py index d07cba7ee..01042709d 100644 --- a/mtuq/graphics/uq/_gmt.py +++ b/mtuq/graphics/uq/_gmt.py @@ -69,7 +69,8 @@ def _plot_latlon_gmt(filename, lon, lat, values, best_latlon=None, lune_array=No def _call(shell_script, filename, data, supplemental_data=None, title='', colormap='viridis', flip_cpt=False, colorbar_type=1, - colorbar_label='', colorbar_limits=None, marker_coords=None, marker_type=0): + colorbar_label='', colorbar_limits=None, marker_coords=None, marker_type=0, + saturation=0.): # # Common wrapper for all GMT plotting functions involving 2D surfaces @@ -95,6 +96,9 @@ def _call(shell_script, filename, data, supplemental_data=None, except: minval, maxval, exp = _parse_limits(data[:,-1]) + if saturation > 0.: + maxval = maxval - saturation*(maxval-minval) + data[:,-1] /= 10.**exp cpt_step=(maxval-minval)/20. From c6da78b4663a7438a264d9f8b91ff60cf8b5ccc7 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Fri, 22 Mar 2024 18:37:13 -0600 Subject: [PATCH 28/31] Backward compatible keyword argument changes --- mtuq/graphics/header.py | 180 ++++++++++++++++++++-------------------- 1 file changed, 88 insertions(+), 92 deletions(-) diff --git a/mtuq/graphics/header.py b/mtuq/graphics/header.py index c56231e3a..44a0cf239 100644 --- a/mtuq/graphics/header.py +++ b/mtuq/graphics/header.py @@ -87,28 +87,15 @@ def write(self, height, width, margin_left, margin_top): _write_text(text, xp, yp, ax, **kwargs) +class SourceHeader(Base): + """ Base class for moment tensor and force headers -class MomentTensorHeader(Base): - """ Writes UAF-style moment tensor inversion summary to the top of a - matplotlib figure + (Added to reduce duplication, still somewhat of an afterthought) """ - def __init__(self, process_bw, process_sw, misfit_bw, misfit_sw, - best_misfit_bw, best_misfit_sw, model, solver, mt, lune_dict, origin, - event_name=None): - - if not event_name: - # YYYY-MM-DDTHH:MM:SS.??????Z - event_name = '%s' % origin.time - - # trim of fraction of second - event_name = event_name[:-8] - - self.event_name = event_name - self.origin = origin - self.magnitude = mt.magnitude() - depth_in_m = origin.depth_in_m - depth_in_km = origin.depth_in_m/1000. + def parse_origin(self): + depth_in_m = self.origin.depth_in_m + depth_in_km = self.origin.depth_in_m/1000. if depth_in_m < 1000.: self.depth_str = '%.0f m' % depth_in_m elif depth_in_km <= 100.: @@ -117,46 +104,84 @@ def __init__(self, process_bw, process_sw, misfit_bw, misfit_sw, self.depth_str = '%.0f km' % depth_in_km - self.model = model - self.solver = solver - - self.mt = mt - self.lune_dict = lune_dict - - self.process_bw = process_bw - self.process_sw = process_sw - - self.norm = misfit_sw.norm - - self.best_misfit_bw = best_misfit_bw - self.best_misfit_sw = best_misfit_sw + def parse_misfit(self): + # TODO - keep track of body and surface wave norms + self.norm = self.misfit_sw.norm self.best_misfit = self.best_misfit_bw + self.best_misfit_sw - if not process_bw: + + def parse_data_processing(self): + if not self.process_bw: pass - if not process_sw: + if not self.process_sw: raise Excpetion() - if process_sw.freq_max > 1.: + if self.process_sw.freq_max > 1.: units = 'Hz' else: units = 's' - if process_bw and units=='Hz': + if self.process_bw and units=='Hz': self.passband_bw = '%.1f - %.1f Hz' %\ - (process_bw.freq_min, process_bw.freq_max) + (self.process_bw.freq_min, self.process_bw.freq_max) - elif process_bw and units=='s': + elif self.process_bw and units=='s': self.passband_bw = '%.1f - %.1f s' %\ - (process_bw.freq_max**-1, process_bw.freq_min**-1) + (self.process_bw.freq_max**-1, self.process_bw.freq_min**-1) - if process_sw and units=='Hz': + if self.process_sw and units=='Hz': self.passband_sw = '%.1f - %.1f Hz' %\ - (process_sw.freq_min, process_sw.freq_max) + (self.process_sw.freq_min, self.process_sw.freq_max) - elif process_sw and units=='s': + elif self.process_sw and units=='s': self.passband_sw = '%.1f - %.1f s' %\ - (process_sw.freq_max**-1, process_sw.freq_min**-1) + (self.process_sw.freq_max**-1, self.process_sw.freq_min**-1) + + + +class MomentTensorHeader(SourceHeader): + """ Writes moment tensor inversion summary to the top of a + matplotlib figure + """ + def __init__(self, process_bw, process_sw, misfit_bw, misfit_sw, + best_misfit_bw, best_misfit_sw, model, solver, mt, lune_dict, origin, + data=None, synthetics=None, mt_grid=None, event_name=None): + + if not event_name: + # YYYY-MM-DDTHH:MM:SS.??????Z + event_name = '%s' % origin.time + + # trim fraction of second + event_name = event_name[:-8] + + self.event_name = event_name + + # required arguments + self.process_bw = process_bw + self.process_sw = process_sw + self.misfit_bw = misfit_bw + self.misfit_sw = misfit_sw + self.best_misfit_bw = best_misfit_bw + self.best_misfit_sw = best_misfit_sw + self.model = model + self.solver = solver + self.mt = mt + self.lune_dict = lune_dict + self.origin = origin + + # optional arguments, reserved for possible future use + # (or for use by subclasses) + self.data = data + self.synthetics = synthetics + self.mt_grid = mt_grid + + # moment tensor-derived attributes + self.magnitude = mt.magnitude() + + + self.parse_origin() + self.parse_misfit() + self.parse_data_processing() def display_source(self, ax, height, width, offset): @@ -174,7 +199,6 @@ def display_source(self, ax, height, width, offset): # beach(self.mt, xy=(xp, yp), width=diameter, # linewidth=0.5, facecolor='gray')) - # # Instead, we must use this workaround # @@ -252,73 +276,46 @@ def write(self, height, width, margin_left, margin_top): _write_text(line, px, py, ax, fontsize=14) -class ForceHeader(Base): + +class ForceHeader(SourceHeader): """ Writes force inversion summary to the top of a matplotlib figure """ def __init__(self, process_bw, process_sw, misfit_bw, misfit_sw, best_misfit_bw, best_misfit_sw, model, solver, force, force_dict, origin, - event_name=None): + data=None, synthetics=None, force_grid=None, event_name=None): if not event_name: # YYYY-MM-DDTHH:MM:SS.??????Z event_name = '%s' % origin.time - # trim of fraction of second + # trim fraction of second event_name = event_name[:-8] self.event_name = event_name - self.origin = origin - - depth_in_m = origin.depth_in_m - depth_in_km = origin.depth_in_m/1000. - if depth_in_m < 1000.: - self.depth_str = '%.0f m' % depth_in_m - elif depth_in_km <= 100.: - self.depth_str = '%.1f km' % depth_in_km - else: - self.depth_str = '%.0f km' % depth_in_km - - self.model = model - self.solver = solver - - self.force = force - self.force_dict = force_dict + # required arguments self.process_bw = process_bw self.process_sw = process_sw - - self.norm = misfit_sw.norm - + self.misfit_bw = misfit_bw + self.misfit_sw = misfit_sw self.best_misfit_bw = best_misfit_bw self.best_misfit_sw = best_misfit_sw - self.best_misfit = self.best_misfit_bw + self.best_misfit_sw - - if not process_bw: - pass - if not process_sw: - raise Excpetion() - - if process_sw.freq_max > 1.: - units = 'Hz' - else: - units = 's' - - if process_bw and units=='Hz': - self.passband_bw = '%.1f - %.1f Hz' %\ - (process_bw.freq_min, process_bw.freq_max) - - elif process_bw and units=='s': - self.passband_bw = '%.1f - %.1f s' %\ - (process_bw.freq_max**-1, process_bw.freq_min**-1) + self.model = model + self.solver = solver + self.force = force + self.force_dict = force_dict + self.origin = origin - if process_sw and units=='Hz': - self.passband_sw = '%.1f - %.1f Hz' %\ - (process_sw.freq_min, process_sw.freq_max) + # optional arguments, reserved for possible future use + # (or for use by subclasses) + self.data = data + self.synthetics = synthetics + self.force_grid = force_grid - elif process_sw and units=='s': - self.passband_sw = '%.1f - %.1f s' %\ - (process_sw.freq_max**-1, process_sw.freq_min**-1) + self.parse_origin() + self.parse_misfit() + self.parse_data_processing() def write(self, height, width, margin_left, margin_top): @@ -400,7 +397,6 @@ def display_source(self, ax, height, width, offset): - def _lat_lon(origin): if origin.latitude >= 0: latlon = '%.2f%s%s' % (+origin.latitude, u'\N{DEGREE SIGN}', 'N') From 7df2f47b24d8fa49b20ef96d9ee81cc3b63b03ae Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Tue, 26 Mar 2024 23:40:03 -0600 Subject: [PATCH 29/31] Optional GMT colormap clipping --- mtuq/graphics/uq/_gmt.py | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/mtuq/graphics/uq/_gmt.py b/mtuq/graphics/uq/_gmt.py index 01042709d..73b250c17 100644 --- a/mtuq/graphics/uq/_gmt.py +++ b/mtuq/graphics/uq/_gmt.py @@ -69,8 +69,8 @@ def _plot_latlon_gmt(filename, lon, lat, values, best_latlon=None, lune_array=No def _call(shell_script, filename, data, supplemental_data=None, title='', colormap='viridis', flip_cpt=False, colorbar_type=1, - colorbar_label='', colorbar_limits=None, marker_coords=None, marker_type=0, - saturation=0.): + colorbar_label='', cpt_limits=None, cpt_steps=20, + cpt_clip=[0., 1.], marker_coords=None, marker_type=0): # # Common wrapper for all GMT plotting functions involving 2D surfaces @@ -89,18 +89,28 @@ def _call(shell_script, filename, data, supplemental_data=None, # parse colorbar limits try: - minval, maxval = colorbar_limits + minval, maxval = cpt_limits exp = _exponent((minval, maxval)) minval /= 10.**exp maxval /= 10.**exp except: minval, maxval, exp = _parse_limits(data[:,-1]) - if saturation > 0.: - maxval = maxval - saturation*(maxval-minval) + # optionally, clip colorbar limits to bring out more detail + try: + assert 0. <= cpt_clip[0] < cpt_clip[1] + minval = minval + cpt_clip[0]*(maxval - minval) + except: + pass + try: + assert cpt_clip[0] < cpt_clip[1] <= 1. + maxval = minval + cpt_clip[1]*(maxval - minval) + except: + pass + data[:,-1] /= 10.**exp - cpt_step=(maxval-minval)/20. + cpt_step=(maxval-minval)/cpt_steps # write values to be plotted as ASCII table ascii_file_1 = _safename('tmp_'+filename+'_ascii1.txt') From 33e8c94b807bde9d91a108fdf96cd47edc5f1c48 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Tue, 26 Mar 2024 23:40:46 -0600 Subject: [PATCH 30/31] Updated force header --- mtuq/graphics/header.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/mtuq/graphics/header.py b/mtuq/graphics/header.py index 44a0cf239..2f3398fa8 100644 --- a/mtuq/graphics/header.py +++ b/mtuq/graphics/header.py @@ -145,7 +145,7 @@ class MomentTensorHeader(SourceHeader): """ def __init__(self, process_bw, process_sw, misfit_bw, misfit_sw, best_misfit_bw, best_misfit_sw, model, solver, mt, lune_dict, origin, - data=None, synthetics=None, mt_grid=None, event_name=None): + data_bw=None, data_sw=None, mt_grid=None, event_name=None): if not event_name: # YYYY-MM-DDTHH:MM:SS.??????Z @@ -171,8 +171,8 @@ def __init__(self, process_bw, process_sw, misfit_bw, misfit_sw, # optional arguments, reserved for possible future use # (or for use by subclasses) - self.data = data - self.synthetics = synthetics + self.data_bw = data_bw + self.data_sw = data_sw self.mt_grid = mt_grid # moment tensor-derived attributes @@ -283,7 +283,7 @@ class ForceHeader(SourceHeader): def __init__(self, process_bw, process_sw, misfit_bw, misfit_sw, best_misfit_bw, best_misfit_sw, model, solver, force, force_dict, origin, - data=None, synthetics=None, force_grid=None, event_name=None): + data_bw=None, data_sw=None, force_grid=None, event_name=None): if not event_name: # YYYY-MM-DDTHH:MM:SS.??????Z @@ -309,8 +309,8 @@ def __init__(self, process_bw, process_sw, misfit_bw, misfit_sw, # optional arguments, reserved for possible future use # (or for use by subclasses) - self.data = data - self.synthetics = synthetics + self.data_bw = data_bw + self.data_sw = data_sw self.force_grid = force_grid self.parse_origin() From dbb2fc15b4aa04907582d08f22aebe8f7805c637 Mon Sep 17 00:00:00 2001 From: Ryan Modrak Date: Tue, 26 Mar 2024 23:44:54 -0600 Subject: [PATCH 31/31] Fixed GridSearch.Force.py comments --- examples/GridSearch.Force.py | 18 +++--------------- 1 file changed, 3 insertions(+), 15 deletions(-) diff --git a/examples/GridSearch.Force.py b/examples/GridSearch.Force.py index 43489c2ab..8decd7906 100755 --- a/examples/GridSearch.Force.py +++ b/examples/GridSearch.Force.py @@ -18,19 +18,7 @@ if __name__=='__main__': # - # Carries out grid search over 64,000 double couple moment tensors - # - # USAGE - # mpirun -n python GridSearch.DoubleCouple.py - # - # For a simpler example, see SerialGridSearch.DoubleCouple.py, - # which runs the same inversion in serial - # - - - # - # We will investigate the source process of an Mw~4 earthquake using data - # from a regional seismic array + # Carries out grid search over force orientation and magnitude # path_data= fullpath('data/examples/20210809074550/*[ZRT].sac') @@ -40,7 +28,7 @@ # - # We are only using surface waves in this example. Check out the DC or FMT examples for multi-mode inversions. + # We are only using surface waves in this example. Check out the DC or FMT examples for multi-mode inversions # process_sw = ProcessData( @@ -77,7 +65,7 @@ # - # Next, we specify the moment tensor grid and source-time function + # Next, we specify the force grid and source-time function # grid = ForceGridRegular(magnitudes_in_N=10**np.arange(9,12.1,0.1), npts_per_axis=90)