Comparisons of Ground Motion Predictions for EEW
This code package is written to compare ground motion predictions made from different algorithms. The original intent is to compare the source models coming out of geodetic algorithms for ShakeAlert, the different seismic algorithms, and the ground truth (real PGA/PGV values recorded at seismic stations).
The usage of this code package is fully dictated by the .props file. If the name of the properties file is changed, only the following line towards the top of GMCOMP.py needs to be modified:
props = Properties('gmcomp.props')
As such, maintaining the integrity of the .props file is key. The variables should have no spaces on a line such that:
variable=value
Since there are lots of variables and paths, my suggestion is to copy your original .props file and create a new one for each earthquake run you want to perform such that:
earthquake1.props, earthquake2.props, etc.
Output file names will be modified according to specific properties you use such as earthquake name, gmpe used, gmice used, definition of peak motion, and algorithm name, so you don't need to make a new properties file for small changes such as gmpe used since the output files will not overwrite others for the same earthquake.
If everything is setup correctly, you should be able to run the code with:
python3 GMCOMP.py
Prerequisites
In order to run GMCOMP, you will need several Python3 packages. On OSX, use pip3 for all of these, on other operating systems it will be dependent on your package manager, but these are fairly standard packages.
obspy
shapely
numpy
For plotting, I have set up several GMT scripts. The most annoying thing with GMT is ensuring your path is correct. I have added the GMT path as a variable in the props file; this should be the same as in your .bashrc file. If your version of GMT is finnicky, simply turn off plotting in the props file.
Ground Motion Models (GMMs):
There are two types of GMMs within GMCOMP, Ground Motion Prediction Equations (GMPEs) and Ground Motion Intensity Conversion Equations (GMICEs). The two Python packages are aptly called GMPE.py and GMICE.py.
GMPE.py - There are currently 4 GMPEs within the main package that compute PGA or PGV. The name of the GMPE is stored in the properties dictionary (read from the .props file), and is called within GMPE.py through props.getgmpe(). To perform ground motion predictions, the following is called:
Y = GMPE.Y(M, Rrup, Rjb, Rx, Ztor, Zhyp, VS30, Hw, W, dip, rake, mode, props)
and the return is either the values of PGA or PGV at the locations specified. Here, the locations of interest are embedded within the different distance measures and site conditions.
The variables are:
M - magnitude (single value)
Rrup - Rupture distance (km, a numpy array, n by 1, where n is the number of sites/locations you are interested in)
Rjb - Joyner Boore distance (km, array)
Rx - Distance to the updip projection of the fault plane (km, array)
Ztor - Depth to the top of the rupture plane (km, array)
Zhyp - hypocentral depth (km, single value)
VS30 - Shear wave velocity in upper 30 meters at site (m/s, array)
Hw - Hanging wall flag (1 for on hanging wall or 0 for not, array)
W - width of fault plane (km, single value)
dip - dip of the fault plane (degrees, single value)
rake - rake of the earthwauke (degrees, single value)
mode - 0 for PGA, 1 for PGV. PGV is not fully operational
props - the dictionary of properties.
The 4 different GMPEs are:
cy08 - Chiou and Youngs (2008)
cb14 - Campbell and Bozorgnia (2014)
ba08 - Boore and Atkinson (2008)
bradley13 - Bradley (2013)
GMICE.py - There are two GMICEs embedded within GMICE.py. They both can take in either the PGA or PGV values and convert to instrumental MMI. The particular GMICE equation used is embedded in the properties dictionary, accessible through props.getgmice(). The usage is the following:
MMI = GMICE.MMI_Y(Y, mode, props)
where Y is either the PGA or PGV values, mode is 0 for PGA, 1 for PGV, and props is the dictionary of properties. The two GMICE equations are as follows:
wald99 - Wald et al. (1999)
wgrw12 - Worden et al. (2012)
MSG_xmlreader.py
The xmlreader section contains 5 different xml readers for source information. The most basic is the cireader, which takes in a point source model. An example coreinfo message is below:
<event_message alg_vers="3.1.0-2018-04-09" category="live" instance="E2@eew-bk-dev1" message_type="new" orig_sys="elarms" ref_id="0" ref_src="" timestamp="2018-04-25T20:53:27Z" version="0">
<core_info id="1497">
<mag units="Mw">5.6094</mag>
<mag_uncer units="Mw">0.3366</mag_uncer>
<lat units="deg">38.2250</lat>
<lat_uncer units="deg">0.1000</lat_uncer>
<lon units="deg">-122.2866</lon>
<lon_uncer units="deg">0.1000</lon_uncer>
<depth units="km">8.0000</depth>
<depth_uncer units="km">5.0000</depth_uncer>
<orig_time units="UTC">2018-04-25T20:53:24.268Z</orig_time>
<orig_time_uncer units="sec">1.6751</orig_time_uncer>
<likelihood>0.9256</likelihood>
<num_stations>5</num_stations>
</core_info>
</event_message>
The output of the cireader are two lists, the first the mess_overview which contains the information from the event_message tag and second the mess_details which contains the source parameters such as location, magnitude, and origin time.
The next two subroutines both read finite fault information except the second one only reads a line source finite fault message (ffreader and ffreader_line).
GMCOMP_coord_tools.py
The coord_tools package performs standard coordinate transformations. All subroutines are array capable. Not all of them are used, but the following subroutines are available:
GMCOMP_coord_tools.lla2ecef(lat,lon,alt) - lon,lat,alt to ITRF cartesian
GMCOMP_coord_tools.ecef2lla(x,y,z) - ITRF Cartesian to lon,lat,alt
GMCOMP_coord_tools.dxyz2dneu(dx,dy,dz,lat,lon) - displacement in xyz to neu
GMCOMP_coord_tools.ll2utm(lon,lat,lon0,lat0) - lon/lat to UTM
GMCOMP_coord_tools.utm2ll(UTMeasting,UTMnorthing,lon0,lat0) - UTM to lon/lat
Plotting
In the .props file, there are three variables related to plotting. plotson is either yes or no to turn on plotting, gmtpath controls the path to your GMT bin folder, and gmtversion gives the overarching version of GMT (right now only 5 is supported but I will write GMT4 scripts eventually).
On Macs, the easiest way to get GMT5 is through homebrew such that:
brew install gmt@5
This will place your GMT path in:
/usr/local/Cellar/gmt/"version"/bin
The main plotting script is GMCOMP_gmtplots.py and this will be called at the end of GMCOMP.py. A number of GMT scripts will be created in your output directory (outputdir in .props), and you can modify these to change formatting of figures according to your needs. You can simply go to your output directory, modify, and rerun the following from the command line in the main GMCOMP directory (since there are path references that wouldn't work within the output directory):
./plot_mmibias
This creates MMI Bias PDFs as a function of warning time. If the warning threshold is not exceeded, I use theoretical S wave traveltimes derived from the Obspy taup package.
./plot_cdf_fp
This creates a pseudo-cdf and false-positive chart. What I mean by pseudo-cdf is the plot is the cumulative number of warned stations as a function of warning time, however, that number can decrease as the model results are updated. This would assume optimal cancellations of alerts, which is not the current paradigm.
./plot_quad (First alert, video capable)
This creates a plot of observed versus predicted MMI. The reason I call it a quad plot is the four quadrants can be segregated into true positive (TP), true negative (TN), false positive (FP), and false negative (FN). The dots are color-coded by warning time, fully saturated black dots represent no or negative warning time. The vertical and horizontal lines are modified by the mmiwarnthreshold defined in the .props file.
./plot_warntime
This plots the observed MMI versus the warning time. The color coding is by MMI bias (predicted minus observed). The horizontal line is modified by the mmiwarnthreshold defined in the .props file.
./plot_costratio (First alert, video capable)
This plots the cost savings performance metric, Q as a function of the cost ratio, r.
./plot_mmibias_map (First alert, video capable)
This plots the MMI bias of the first alert for an algorithm.
./plot_warntime_map
This plots the warning time map for a given algorithm based on the first alert times.
./plot_mmiwarntimedensity
This plots the probability of warning times based upon a given mmi band. I have set this to 0.5 MMI units between the warning threshold and MMI=10.
