Newer
Older
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Ground motion model (GMM) ground motion vs distance subplots
%
% Author: Peter Powers (pmpowers@usgs.gov)
% Demi Girot (dgirot@usgs.gov)
%
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
%
% This script demonstrates how to call web services for the ground motion
% models implemented in USGS National Sesmic Hazard Models and plot the
% results. A JSON web service response is automatically converted to a
% struct by the Matlab webread() function. This script shows how to plot
% the GMMs requested and optionally the underlying models of epistemic
% uncertainty.
%
% Users are STRONGLY encouraged to run the web services locally so as to
% not overburden USGS servers. See the repository README for details on how
% to do so: https://code.usgs.gov/ghsc/nshmp/nshmp-ws
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot a figure of GMM ground motion vs distance, sigmas, and epistemic
% uncertainty
clearvars
% The root url for the response spectra web service
urlbase = "https://earthquake.usgs.gov/ws/nshmp/data/gmm/distance?";
% If you are running nshmp-ws locally, use this URL:
% urlbase = "http://localhost:8080/gmm/distance?";
% Struct of ground motion model parameters
in.Mw = 7.5;
in.dip = 90;
in.rake = 0;
in.width = 14;
in.rJB = 10;
in.rMax = 300;
in.rMin = 0.1;
in.rRup = 10.3;
in.rX = 10;
in.vs30 = 760;
in.vsInf = true;
in.zHyp = 7.5;
in.zTor = 0.5;
% Uncomment to include any of these additional parameters in the struct
% in.z2p5 = 1;
% in.z1p0 = 1;
% in.zSed = 1;
% Input imt value from PGA to 10.0s spectral acceleration
in.imt = 'PGA';
in.colors = ["#D95319" "#0072BD" "#7E2F8E" "#77AC30" "#EDB120"];
gmms_2023 = ["ASK_14" "BSSA_14" "CB_14" "CY_14"];
%%
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
figure
plot_handles = [];
% Set to true to show epistemic branches
epi = false;
subplot(2,1,1)
means = readGMvD(urlbase, gmms_2023, in).response.means.data;
plotGmms(means, in, epi);
subplot(2,1,2)
sigmas = readGMvD(urlbase, gmms_2023, in).response.sigmas.data;
plotSigmas(sigmas, in, epi);
%% Shared Functions
function plotGmms(gmm_data, in, epi)
plot_handles = [];
legend_labels = {};
% loop means response array
for j = 1 : length(gmm_data)
gmm_xs = gmm_data(j).data.xs;
gmm_ys = gmm_data(j).data.ys;
gmm_tree = gmm_data(j).tree;
% if 0.01 absent use PGA (0.001 is used for PGA)
if gmm_xs(1) == 0.001 && gmm_xs(2) ~= 0.01
gmm_xs(1) = 0.01;
end
% Plot the total spectrum
ph = loglog(gmm_xs, gmm_ys, ...
'LineWidth',3, ...
'LineStyle','-');
color = get(ph,'Color');
hold on;
plot_handles = [plot_handles ph];
legend_labels{end+1} = gmm_data(j).label;
% Plot epistemic spectra, if present
if epi && ~isempty(gmm_tree)
for k = 1 : length(gmm_tree)
epi_ys = gmm_tree(k).values;
phe = loglog(gmm_xs, epi_ys, ...
'LineWidth',1, ...
'LineStyle',':', ...
if k == 1
plot_handles = [plot_handles phe];
legend_labels{end+1} = [gmm_data(j).label ' ε'];
end
end
end
end
grid on
xlabel('Distance (km)','FontSize',16)
ylabel('Median Ground Motion (g)','FontSize',16)
title(['M',num2str(in.Mw),', ',num2str(in.rRup),'km'],'FontSize',20)
xlim([0.1 300]);
set(gca,'FontSize',20);
set(gca,'XTick',[0.1 1 10 100]);
set(gca,'XTickLabel',{'0.1','1','10','100'});
legend( ...
plot_handles,legend_labels, ...
'Location','best', ...
'NumColumns',2, ...
'FontSize',12, ...
'Interpreter', 'none');
end
function plotSigmas(gmm_data, in, epi)
plot_handles = [];
legend_labels = {};
% loop sigmas response array
sig_xs = gmm_data(l).data.xs;
sig_ys = gmm_data(l).data.ys;
sig_tree = gmm_data(l).tree;
% if 0.01 absent use PGA (0.001 is used for PGA)
if sig_xs(1) == 0.001 && sig_xs(2) ~= 0.01
sig_xs(1) = 0.01;
end
% Plot the total spectrum
ph = loglog(sig_xs, sig_ys, ...
'LineWidth',3, ...
'LineStyle','-');
color = get(ph,'Color');
hold on;
plot_handles = [plot_handles ph];
legend_labels{end+1} = gmm_data(l).label;
% Plot epistemic spectra, if present
if epi && ~isempty(sig_tree)
for m = 1 : length(sig_tree)
epi_ys = sig_tree(m).values;
phe = loglog(sig_xs, epi_ys, ...
'LineWidth',1, ...
'LineStyle',':', ...
'Color',color);
if m == 1
plot_handles = [plot_handles phe];
legend_labels{end+1} = [gmm_data(l).label ' ε'];
end
end
end
end
grid on
xlabel('Distance (km)','FontSize',16)
ylabel('Standard Deviation','FontSize',16)
title(['M',num2str(in.Mw),', ',num2str(in.rRup),'km'],'FontSize',20)
xlim([0.1 300]);
set(gca,'FontSize',20);
set(gca,'XTick',[0.1 1 10 100]);
set(gca,'XTickLabel',{'0.1','1','10','100'})
legend( ...
plot_handles,legend_labels, ...
'Location','best', ...
'NumColumns',2, ...
'FontSize',12, ...
'Interpreter', 'none');
end
% Read reponse spectra from a NSHM web service for the supplied
% ground motion models (GMMs) and GMM input parameters
function response = readGMvD(urlbase, gmms, in)
url = urlbase;
for i = 1 : size(gmms, 2)
if i == 1
url = url + "gmm=" + gmms(i);
else
url = url + "&gmm=" + gmms(i);
end
end
url = url + ...
"&Mw=" + num2str(in.Mw) + ...
"&dip=" + num2str(in.dip) + ...
"&rake=" + num2str(in.rake) + ...
"&width=" + num2str(in.width) + ...
"&rJB=" + num2str(in.rJB) + ...
"&rMax=" + num2str(in.rMax) + ...
"&rMin=" + num2str(in.rMin) + ...
"&rRup=" + num2str(in.rRup) + ...
"&rX=" + num2str(in.rX) + ...
"&vs30=" + num2str(in.vs30) + ...
"&vsInf=" + string(in.vsInf) + ...
"&zHyp=" + num2str(in.zHyp) + ...
"&zTop=" + num2str(in.zTor) + ...
"&imt=" + string(in.imt);
if isfield(in, 'z2p5')
url = url + num2str(in.z2p5);
end
if isfield(in, 'z1p0')
url = url + num2str(in.z1p0);
end
if isfield(in, 'zSed')
url = url + num2str(in.zSed);
end
response = webread(url);
end