A comprehensive tutorial for plotting focal mechanism "beach-balls" using the PyGMT package for Python.
(Resulting map of this demo)
Background
To better understand the faulting motions that create earthquakes, it can be useful to graphically display focal mechanisms. A focal mechanism is the orientation and type of slip of a fault, and it is represented as a “beach-ball” symbol.
(Image is sourced from the United States Geological Survey earthquake hazards webpage on focal mechanisms)
Within a “beach-ball,” the dark quarters represent compressional forces and the white quarters represent tensional forces. A “beach-ball” has two possible orientations of fault movement, each along with one of the planes dividing the “beach-ball” into fourths, and determining the correct orientation requires additional information. When “beach-balls” present four equally sized quadrants in a “top-down” view like the one in the above diagram, strike-slip faulting is represented. Rotated “beach-balls” with a white quarter in the center represent normal faulting, while “beach-balls” with a dark quarter in the center represent thrust faulting.
Please keep in mind that with the exception of the folder creation script, all code blocks shown in this demo are intended to be run as part of a complete script, which is included at the bottom of the page.
Getting Started
Create Project Folders
To follow this demo as written, project folders need to be created on the Desktop. A diagram of the desired directory tree for Windows 10 is shown below:
(Image created with Lucidchart)
These folders can be created either by running the script provided below from anywhere after changing the main_dir file path to reflect your username or by manually creating a Focal_Mechanism_Demo folder and within it Data, Methods, and Results folders.
'''
Decription:
This script creates folders for the focal_mechanism_demo.py script if they don't already exist
'''
import os
# main directory path
main_dir = r'C:\Users\USER\Desktop\Focal_Mechanism_Demo'
# path to the directory holding the project data files
data_folder = os.path.join(main_dir, 'Data')
# path to the directory holding the project Python scripts
methods_folder = os.path.join(main_dir, 'Methods')
# path to the directory holding the map generated by the scripts
results_folder = os.path.join(main_dir, 'Results')
directories = [main_dir, data_folder, methods_folder, results_folder]
# Iterates through the list of directories and creates them if they don't already exist
for directory in directories:
os.makedirs(directory, exist_ok = True)
Data files used and created by the demo script will be located in the Data folder, and the focal mechanism map will save to the Results folder. The Methods folder should contain the demo scripts.
Retrieve Focal Mechanism Data
Focal mechanism data will be acquired by downloading it from the SCEDC focal mechanism catalog. This demo uses the following search parameters:
Date range: 2019/01/01 to 2020/01/01
Magnitude range: 5 to 10
Latitude range: 35.5 to 36.1
Longitude range: -118 to -117.15
Depth range: -5 to 100
(Example of entering the focal mechanism parameters into the SCEDC focal mechanism catalog search)
After clicking the search button, the results are presented as a text page. Select the data and column headers as shown below, and copy them into a text file. Then save the text file to the Data folder and name it focal_mechanism_data.
(Example of selecting and copying the focal mechanism data and column headers)
Creating The Project Files
Conditioning The Focal Mechanism Data
The focal mechanism data as retrieved from the SCEDC does not have a uniform delimiter separating the columns and has an empty line between the column header and the column values. In order to read the data into a data frame, it needs to be conditioned so that the aforementioned issues are resolved.
The following demo script function removes the empty line and standardizes the delimiter to be tabs, then saves the conditioned data as a CSV file:
# Conditions focal mechanism files downloaded from the Southern California Earthquake Data Center so that they can be used
def Condition_Focal_Mechanism_Data(self):
import os
import re
print('Conditioning focal mechanisms...')
input_file = os.path.join(self.main_dir, 'Data', 'focal_mechanism_data.txt')
# reads in the file by line rather than the default of by character
with open(input_file, 'r') as f:
lines = f.readlines()
# holds the conditioned lines, namely the removal of empty lines
conditioned_lines = []
# iterates over each line and adds them to a list as long as they are not empty
for line in lines:
if line == '\n':
continue
conditioned_lines.append(line)
# empty placeholder for the conditioned lines list to be rejoined into a string
conditioned_data = ''
# rejoins the conditioned lines into a string so that it can be further conditioned
conditioned_data = conditioned_data.join(conditioned_lines)
# regular expression pattern for finding 1 or more white spaces (key) and the replacment value for matches (value)
replacements = {' +':'\t'}
# subsitutes text using the patterns in "replacements"
for key, value in replacements.items():
data = re.sub(key, value, conditioned_data)
output_file = os.path.join(self.main_dir, 'Data', 'focal_mechanism_data.csv')
with open(output_file, 'w') as f:
f.write(data)
Filtering Focal Mechanism Data Into The “Aki” Format
After conditioning the focal mechanism data, it needs to be filtered so that the resulting file can be read by the PyGMT meca function, which plots focal mechanism “beach-balls.” As the data provided by the SCEDC is in terms of strike, drip, rake, and magnitude, the focal mechanism “beach-balls” will be determined by way of the Aki and Richards (1980) convention. In this demo, earthquakes greater than or equal to M7.0 will be offset to showcase how offsetting is accomplished. The meca function requires a file to be passed directly to the spec parameter when offsetting, so a CSV file with specific columns and column order pertaining to the convention needs to be created. Additionally, the column headers must be removed, so that header names don’t plot as text.
The following demo script function splits the focal mechanism data into two data frames at M7.0 earthquakes, then converts the magnitudes so that they are exponentially scaled; this is how the “beach-balls” will be sized when they are plotted. Next, it adds offset coordinates and text label columns to the data frame of M7.0 and greater earthquakes so they will be offset from their point of origin when plotted. Finally, the two data frames are saved as CSV files:
< magnitude_filter)] # focal mechanisms with earthquake magnitude greater than or equal to the magnitude filter df_offset_mag_filtered_focal_mecha =df_focal_mechanisms[(df_focal_mechanisms.MAG >= magnitude_filter)] # reads the unmodifed magnitude of the focal mechanisms that will be offset into a list so it can be later used to label them fm_label = [] for i in df_offset_mag_filtered_focal_mecha.MAG: fm_label.append(i) fm_dataframes = [df_mag_filtered_focal_mecha, df_offset_mag_filtered_focal_mecha] # iterates through each focal mechanism dataframe and scales them expoentially so that they will be plotted as exponentially larger sizes, and resets the index in place while not including the old index as a column which is done for tidyness for dataframe in fm_dataframes: dataframe.reset_index(inplace=True, drop=True) dataframe['MAG'] = 0.1 * (2 ** dataframe.MAG) # new dataframes to hold only the columns used for the "aki" focal mechanism format fm_aki_format = pd.DataFrame() fm_aki_format_offset = pd.DataFrame() aki_dataframes = [fm_aki_format, fm_aki_format_offset] # column names of the aki format, with keys being columns from the focal mechanism dataframes and values being the new column names to be used aki_format_attributes = {'LON':'longitude', 'LAT':'latitude', 'DEPTH':'depth', 'STRIKE':'strike', 'DIP':'dip', 'RAKE':'rake', 'MAG':'magnitude'} # add the desired columns for the 'aki' focal mechanism format to the new dataframes in the desired order for aki_data, fm_data in zip(aki_dataframes, fm_dataframes): for key, value in aki_format_attributes.items(): aki_data[value] = fm_data[key] # offset coordinates (lon:lat) used for offsetting the selected focal mechanisms offset_coordinates = {-117.70:35.78} fm_aki_format_offset['offset_longitude'] = offset_coordinates.keys() fm_aki_format_offset['offset_latitude'] = offset_coordinates.values() # adds label data for plotting above each focal mechanism beachball fm_aki_format_offset['label_mag'] = fm_label # path to the "small" EQ focal mechanisms focal_mechanisms_lesser_output = os.path.join(self.main_dir, 'Data', 'focal_mechanism_data_less_than_M{}.csv').format(magnitude_filter) # path to the "large" EQ focal mechanisms focal_mechanisms_greater_output = os.path.join(self.main_dir, 'Data', 'focal_mechanism_data_greater_than_or_equal_to_M{}.csv').format(magnitude_filter) focal_mechanism_outputs = [focal_mechanisms_lesser_output, focal_mechanisms_greater_output] # writes the aki format dataframes to csv files, leaving out the header. This is important because if the header is included, fig.meca will plot the header text for dataframe, path in zip(aki_dataframes, focal_mechanism_outputs): dataframe.to_csv(path, sep='\t', index=False, header=None) return fm_depths, fm_mags, magnitude_filter ">
# Creates a focal mechanism dataframe in the "aki" format from the focal mechanism data
def Filter_AKI_Format_Focal_Mechanism_Data(self):
import os
import pandas as pd
print('Writing focal mechanisms to CSV...')
focal_mechanisms = os.path.join(self.main_dir, 'Data', 'focal_mechanism_data.csv')
df_focal_mechanisms = pd.read_csv(focal_mechanisms , sep='\t')
# saves the unified depth column as a pandas series so that it can be recalled later for coloring the focal mechanisms by depth
fm_depths = df_focal_mechanisms['DEPTH']
# saves the unified depth magnitude column as a pandas series so that it can be recalled later for creating the legend
fm_mags = df_focal_mechanisms['MAG']
# magnitude of earthquakes used to filter out the focal mechanisms that will be offset
magnitude_filter = 7
# focal mechanisms with earthquake magnitude less than the magnitude filter
df_mag_filtered_focal_mecha = df_focal_mechanisms[(df_focal_mechanisms.MAG < magnitude_filter)]
# focal mechanisms with earthquake magnitude greater than or equal to the magnitude filter
df_offset_mag_filtered_focal_mecha =df_focal_mechanisms[(df_focal_mechanisms.MAG >= magnitude_filter)]
# reads the unmodifed magnitude of the focal mechanisms that will be offset into a list so it can be later used to label them
fm_label = []
for i in df_offset_mag_filtered_focal_mecha.MAG:
fm_label.append(i)
fm_dataframes = [df_mag_filtered_focal_mecha, df_offset_mag_filtered_focal_mecha]
# iterates through each focal mechanism dataframe and scales them expoentially so that they will be plotted as exponentially larger sizes, and resets the index in place while not including the old index as a column which is done for tidyness
for dataframe in fm_dataframes:
dataframe.reset_index(inplace=True, drop=True)
dataframe['MAG'] = 0.1 * (2 ** dataframe.MAG)
# new dataframes to hold only the columns used for the "aki" focal mechanism format
fm_aki_format = pd.DataFrame()
fm_aki_format_offset = pd.DataFrame()
aki_dataframes = [fm_aki_format, fm_aki_format_offset]
# column names of the aki format, with keys being columns from the focal mechanism dataframes and values being the new column names to be used
aki_format_attributes = {'LON':'longitude', 'LAT':'latitude', 'DEPTH':'depth', 'STRIKE':'strike', 'DIP':'dip', 'RAKE':'rake', 'MAG':'magnitude'}
# add the desired columns for the 'aki' focal mechanism format to the new dataframes in the desired order
for aki_data, fm_data in zip(aki_dataframes, fm_dataframes):
for key, value in aki_format_attributes.items():
aki_data[value] = fm_data[key]
# offset coordinates (lon:lat) used for offsetting the selected focal mechanisms
offset_coordinates = {-117.70:35.78}
fm_aki_format_offset['offset_longitude'] = offset_coordinates.keys()
fm_aki_format_offset['offset_latitude'] = offset_coordinates.values()
# adds label data for plotting above each focal mechanism beachball
fm_aki_format_offset['label_mag'] = fm_label
# path to the "small" EQ focal mechanisms
focal_mechanisms_lesser_output = os.path.join(self.main_dir, 'Data', 'focal_mechanism_data_less_than_M{}.csv').format(magnitude_filter)
# path to the "large" EQ focal mechanisms
focal_mechanisms_greater_output = os.path.join(self.main_dir, 'Data', 'focal_mechanism_data_greater_than_or_equal_to_M{}.csv').format(magnitude_filter)
focal_mechanism_outputs = [focal_mechanisms_lesser_output, focal_mechanisms_greater_output]
# writes the aki format dataframes to csv files, leaving out the header. This is important because if the header is included, fig.meca will plot the header text
for dataframe, path in zip(aki_dataframes, focal_mechanism_outputs):
dataframe.to_csv(path, sep='\t', index=False, header=None)
return fm_depths, fm_mags, magnitude_filter
Creating The Legend
When the focal mechanisms “beach-balls” are eventually plotted on a map, the size differences will need to be given context with a legend. For this kind of legend, the PyGMT legend function needs to be passed a postscript file.
The following demo script function creates a postscript file for the legend:
# Creates a postscript legend file
def Create_Legend(self, fm_mags):
import os
import math
# File path to postscript legend file to be created
focal_mechanism_legend = os.path.join(self.main_dir, 'Data', 'focal_mechanism_legend.txt')
# Creates a postscript file and writes some explainer, legend header text, and the column format
with open(focal_mechanism_legend, 'w') as f:
f.write(
'# G is vertical gap, V is vertical line, N sets # of columns,\n'
'# D draws horizontal line, H is header, L is column header,\n'
'# S is symbol
\n' '\n' 'H 12p,Helvetica Magnitude\n' 'D 0.1i 1p\n' 'G 0.05i\n' 'N 1\n' '\n' ) # sets the index to the minimum magnitude of the unified dataset. This is used to set the minimum size for the circles used to represent magnitude within the legend, as well as label them i = math.floor(fm_mags.min()) # sets the vertical space between each magnitude circle/text line v_spacer = 0.08 # creteas a legend entry for each integer magnitude within the dataset and appends it to the legend file with open(focal_mechanism_legend, 'a') as f: while i <= round(fm_mags.max(), 0):\ # replicates the exponential scaling of the magnitide and linear scaling of the fig.meca function so that the magnitude circles in the legend are the same size as the plotted focal mechanisms size = 0.1*(0.1*(2**i)) size = round(size, 2) # creates a legend of transparent circles that are labeld with the associated magnitude f.write( 'S 0.20i c {} [email protected] 0.5p 0.6i M{}\n' 'G {}i\n' .format(size, i, v_spacer) ) # exponentially scales the vertical spacing between each legend entry so that they are evenly spaced in this particular demo v_spacer = v_spacer**0.5 i += 1
Building The Map
Now that the focal mechanism data and legend have been properly prepared, it’s time to plot everything on a map.
The following demo script function plots focal mechanism “beach-balls” over an automatically downloaded 3 arc second Digital Elevation Model from the Shuttle Radar Topography Mission. The “beach-balls” are scaled exponentially in size by the magnitude and colored by depth, with magnitude, and depth legends plotted along the map axes. An inset map displaying the study region relative to the greater region around it is plotted in the upper right corner of the map:
+
+
/
legend_position = 'JMR+w1.05i+o0.15i/-0.02i' #
+
#legend_box = '+gwhite+p1p' ### INSET MAP ----------------------------------------------------------------------------- # extent of the inset map i_region = [-130, -105, 27, 45] i_projection = 'M1.75i' # sets the border types, line thickness, and color
. In this case 1 = national boundaries, 2 = state boundaries within America i_borders = ['1/0.8p,gray40', '2/0.8p,gray40'] # sets the color of the land i_land = 'gray' # sets the shorline type, line thickness, and color
. In this case, 1 = coastline i_shorelines = '1/0.1p,gray40' i_pen = '1p,black' # sets the water color i_water = 'azure2' # shifts the inset map in the x-direction. Positive numbers shift up, negative numbers shifts down.
i_xshift = '12.5c' # shifts the inset map in the y-direction. Positive numbers shift up, negative numbers shifts down.
i_yshift ='11.0c' # first two list elements of 'rectangle' are the longitude and latitude of the bottom left corner of the rectangle, and the last two elements are the longitude and latitude of the uppper right corner. rectangle = [[region[0], region[2], region[1], region[3]]] # sets style of plot as rectagle (r) with first two list items as the longitude and latitude of the bottom left corner of the rectangle, and the last two list items as the longitude and latitude of the uppper right corner (+s) i_style = 'r+s' # sets an auto-detrmined frame with ticks in intervals of 10 (a10), then displays the frame with ticks on the top and right sides of the map (NE), but doesn't show ticks for the left and bottoms sides (sw) i_frame = ['a10', 'swNE'] ### Save File ----------------------------------------------------------------------------- # map save name save_name = os.path.join(self.main_dir, 'Results', 'Focal_Mechanism_Demo.png') # sets the transparency of the white space around the map save_transparency = False # Establishes a figure to hold map layers fig = pygmt.Figure() # Forces Lat/Lon to display as degree decimal pygmt.config(FORMAT_GEO_MAP = 'D') # Forces Lat/Lon to display with two decimal places
, as opposed the default that limits by significan digits
pygmt.config(FORMAT_FLOAT_OUT = '%.2f') # Forces the map frame annotation to be smaller pygmt.config(FONT_ANNOT_PRIMARY = '10p,Helvetica,black') # Forces the scale font to be smaller pygmt.config(FONT_LABEL = '10p,Helvetica,black') # Basemap layer fig.basemap( region = region, projection = projection, frame = frame ) # topography layer fig.grdimage( # downloads a 3 arc second shuttle radar topography mission digital elevation model of the region extent '@srtm_relief_03s', # applies shading to the grid shading = True, cmap = grid_cmap, ) depth_series = [fm_depths.min(), fm_depths.max()] # makes a color palette table to color the events by depth pygmt.makecpt( cmap = depth_color, series = depth_series ) for focal_mechanism in focal_mechanisms: fig.meca( # requires a file to be provided directly to the spec parameter if offseting the focal mechanisms is desired spec = focal_mechanism, # sets the focal mechanism convention convention = convention, # linear scaling factor with respect to M5 earthquakes scale = fm_scale, # sets the focal mechanisms to be colored by the previously created color palette table C = True, # controls whether to offset the focal mechanims offset = fm_offset, ) # Plots the color bar for depth fig.colorbar( frame = cbar_frame ) # plots the legend fig.legend( spec = legend, position = legend_position, ) # plots the map scale fig.basemap( map_scale = map_scale ) # Plots a mini coast map as an offset inset map fig.coast( region = i_region, projection = i_projection, frame = i_frame, land = i_land, borders = i_borders, shorelines = i_shorelines, water = i_water, xshift = i_xshift, yshift = i_yshift ) # Plots a rectangle of the study area in the inset map fig.plot( data = rectangle, style = i_style, pen = i_pen, projection = i_projection ) # Saves a copy of the generated figure fig.savefig(save_name, transparent = save_transparency) print('Map Saved!') ">
# Plots the map
def Plot_Map(self, fm_depths, magnitude_filter):
import os
import pygmt
print('Building map...')
### BASEMAP -----------------------------------------------------------------------------
# sets the region to be displayed in the map
region = [-118.0, -117.15, 35.5, 36.1]
# map projection (Mercator):
projection = 'M6i' # frame annotations: [,
,
] frame = ['SWne', 'xa', 'ya'] # map scale - displays a fancy (f) distance scale (lScale) in kilometers (u). Format of
/
/
+
+
Result
I hope that this has been helpful for anyone looking to plot and offset focal mechanisms with PyGMT.
Complete Code
< magnitude_filter)] # focal mechanisms with earthquake magnitude greater than or equal to the magnitude filter df_offset_mag_filtered_focal_mecha =df_focal_mechanisms[(df_focal_mechanisms.MAG >= magnitude_filter)] # reads the unmodifed magnitude of the focal mechanisms that will be offset into a list so it can be later used to label them fm_label = [] for i in df_offset_mag_filtered_focal_mecha.MAG: fm_label.append(i) fm_dataframes = [df_mag_filtered_focal_mecha, df_offset_mag_filtered_focal_mecha] # iterates through each focal mechanism dataframe and scales them expoentially so that they will be plotted as exponentially larger sizes, and resets the index in place while not including the old index as a column which is done for tidyness for dataframe in fm_dataframes: dataframe.reset_index(inplace=True, drop=True) dataframe['MAG'] = 0.1 * (2 ** dataframe.MAG) # new dataframes to hold only the columns used for the "aki" focal mechanism format fm_aki_format = pd.DataFrame() fm_aki_format_offset = pd.DataFrame() aki_dataframes = [fm_aki_format, fm_aki_format_offset] # column names of the aki format, with keys being columns from the focal mechanism dataframes and values being the new column names to be used aki_format_attributes = {'LON':'longitude', 'LAT':'latitude', 'DEPTH':'depth', 'STRIKE':'strike', 'DIP':'dip', 'RAKE':'rake', 'MAG':'magnitude'} # add the desired columns for the 'aki' focal mechanism format to the new dataframes in the desired order for aki_data, fm_data in zip(aki_dataframes, fm_dataframes): for key, value in aki_format_attributes.items(): aki_data[value] = fm_data[key] # offset coordinates (lon:lat) used for offsetting the selected focal mechanisms offset_coordinates = {-117.70:35.78} fm_aki_format_offset['offset_longitude'] = offset_coordinates.keys() fm_aki_format_offset['offset_latitude'] = offset_coordinates.values() # adds label data for plotting above each focal mechanism beachball fm_aki_format_offset['label_mag'] = fm_label # path to the "small" EQ focal mechanisms focal_mechanisms_lesser_output = os.path.join(self.main_dir, 'Data', 'focal_mechanism_data_less_than_M{}.csv').format(magnitude_filter) # path to the "large" EQ focal mechanisms focal_mechanisms_greater_output = os.path.join(self.main_dir, 'Data', 'focal_mechanism_data_greater_than_or_equal_to_M{}.csv').format(magnitude_filter) focal_mechanism_outputs = [focal_mechanisms_lesser_output, focal_mechanisms_greater_output] # writes the aki format dataframes to csv files, leaving out the header. This is important because if the header is included, fig.meca will plot the header text for dataframe, path in zip(aki_dataframes, focal_mechanism_outputs): dataframe.to_csv(path, sep='\t', index=False, header=None) return fm_depths, fm_mags, magnitude_filter # Creates a postscript legend file def Create_Legend(self, fm_mags): import os import math # File path to postscript legend file to be created focal_mechanism_legend = os.path.join(self.main_dir, 'Data', 'focal_mechanism_legend.txt') # Creates a postscript file and writes some explainer, legend header text, and the column format with open(focal_mechanism_legend, 'w') as f: f.write( '# G is vertical gap, V is vertical line, N sets # of columns,\n' '# D draws horizontal line, H is header, L is column header,\n' '# S is symbol
\n' '\n' 'H 12p,Helvetica Magnitude\n' 'D 0.1i 1p\n' 'G 0.05i\n' 'N 1\n' '\n' ) # sets the index to the minimum magnitude of the unified dataset. This is used to set the minimum size for the circles used to represent magnitude within the legend, as well as label them i = math.floor(fm_mags.min()) # sets the vertical space between each magnitude circle/text line v_spacer = 0.08 # creteas a legend entry for each integer magnitude within the dataset and appends it to the legend file with open(focal_mechanism_legend, 'a') as f: while i <= round(fm_mags.max(), 0):\ # replicates the exponential scaling of the magnitide and linear scaling of the fig.meca function so that the magnitude circles in the legend are the same size as the plotted focal mechanisms size = 0.1*(0.1*(2**i)) size = round(size, 2) # creates a legend of transparent circles that are labeld with the associated magnitude f.write( 'S 0.20i c {} [email protected] 0.5p 0.6i M{}\n' 'G {}i\n' .format(size, i, v_spacer) ) # exponentially scales the vertical spacing between each legend entry so that they are evenly spaced in this particular demo v_spacer = v_spacer**0.5 i += 1 # Plots the map def Plot_Map(self, fm_depths, magnitude_filter): import os import pygmt print('Building map...') ### BASEMAP ----------------------------------------------------------------------------- # sets the region to be displayed in the map
region = [-118.0, -117.15, 35.5, 36.1] # map projection (Mercator):
projection = 'M6i' # frame annotations: [,
,
] frame = ['SWne', 'xa', 'ya'] # map scale - displays a fancy (f) distance scale (lScale) in kilometers (u). Format of
/
/
+
+
46 Dec 01, 2022
15 Dec 10, 2022
1 Oct 21, 2021
12 Nov 07, 2022
17.1k Dec 31, 2022
2.3k Jan 02, 2023
5 Jan 03, 2022
2.4k Jan 07, 2023
4 May 18, 2022
10 Jun 01, 2022
12 Aug 24, 2022
11 Dec 25, 2022
432 Dec 30, 2022
21 Nov 20, 2022
3 Jun 18, 2022
12 Oct 20, 2022
502 Dec 28, 2022
336 Dec 20, 2022
2 Dec 19, 2022
5 Sep 16, 2022