From 81dc2a272535e56bd9b0fb6f530646896e74b460 Mon Sep 17 00:00:00 2001
From: Girot <dgirot@usgs.gov>
Date: Fri, 21 Jun 2024 16:46:27 -0600
Subject: [PATCH] Added updated RS GMvD GMvMw scripts
---
etc/matlab/GroundMotionvsDistance.m | 71 ++++++++-----
etc/matlab/GroundMotionvsMagnitude.m | 67 ++++++++----
etc/matlab/ResponseSpectrum.m | 150 +++++++++++++++++++++------
3 files changed, 207 insertions(+), 81 deletions(-)
diff --git a/etc/matlab/GroundMotionvsDistance.m b/etc/matlab/GroundMotionvsDistance.m
index a309373..20754b4 100644
--- a/etc/matlab/GroundMotionvsDistance.m
+++ b/etc/matlab/GroundMotionvsDistance.m
@@ -5,7 +5,7 @@
% Author: Peter Powers (pmpowers@usgs.gov)
% Demi Girot (dgirot@usgs.gov)
%
-% Updated: 05/23/2023
+% Updated: 05/13/2024
%
% This script demonstrates how to call web services for the ground motion
% models implemented in USGS National Sesmic Hazard Models and plot the
@@ -51,14 +51,20 @@ in.zTor = 0.5;
% 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"];
+
+%%
figure
plot_handles = [];
% Set to true to show epistemic branches
epi = false;
-gmms_2023 = ["ASK_14" "BSSA_14" "CB_14" "CY_14"];
-
subplot(2,1,1)
means = readGMvD(urlbase, gmms_2023, in).response.means.data;
plotGmms(means, in, epi);
@@ -88,6 +94,7 @@ function plotGmms(gmm_data, in, epi)
% Plot the total spectrum
ph = loglog(gmm_xs, gmm_ys, ...
+ 'Color',in.colors(j), ...
'LineWidth',3, ...
'LineStyle','-');
@@ -96,20 +103,21 @@ function plotGmms(gmm_data, in, epi)
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)
- id = gmm_tree(k).id;
epi_ys = gmm_tree(k).values;
phe = loglog(gmm_xs, epi_ys, ...
'LineWidth',1, ...
'LineStyle',':', ...
- 'Color',color);
+ 'Color',color);
- plot_handles = [plot_handles phe];
- legend_labels{end+1} = [gmm_data(j).label ' ' id];
+ if k == 1
+ plot_handles = [plot_handles phe];
+ legend_labels{end+1} = [gmm_data(j).label ' ε'];
+ end
end
end
end
@@ -120,12 +128,16 @@ function plotGmms(gmm_data, in, epi)
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');
+ legend( ...
+ plot_handles,legend_labels, ...
+ 'Location','best', ...
+ 'NumColumns',2, ...
+ 'FontSize',12, ...
+ 'Interpreter', 'none');
end
function plotSigmas(gmm_data, in, epi)
@@ -134,11 +146,11 @@ function plotSigmas(gmm_data, in, epi)
legend_labels = {};
% loop sigmas response array
- for m = 1 : length(gmm_data)
+ for l = 1 : length(gmm_data)
- sig_xs = gmm_data(m).data.xs;
- sig_ys = gmm_data(m).data.ys;
- sig_tree = gmm_data(m).tree;
+ 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
@@ -147,6 +159,7 @@ function plotSigmas(gmm_data, in, epi)
% Plot the total spectrum
ph = loglog(sig_xs, sig_ys, ...
+ 'Color',in.colors(l), ...
'LineWidth',3, ...
'LineStyle','-');
@@ -154,21 +167,22 @@ function plotSigmas(gmm_data, in, epi)
hold on;
plot_handles = [plot_handles ph];
- legend_labels{end+1} = gmm_data(m).label;
-
+ legend_labels{end+1} = gmm_data(l).label;
+
% Plot epistemic spectra, if present
if epi && ~isempty(sig_tree)
- for n = 1 : length(sig_tree)
- id = sig_tree(n).id;
- epi_ys = sig_tree(n).values;
+ for m = 1 : length(sig_tree)
+ epi_ys = sig_tree(m).values;
phe = loglog(sig_xs, epi_ys, ...
'LineWidth',1, ...
'LineStyle',':', ...
'Color',color);
-
- plot_handles = [plot_handles phe];
- legend_labels{end+1} = [gmm_data(m).label ' ' id];
+
+ if m == 1
+ plot_handles = [plot_handles phe];
+ legend_labels{end+1} = [gmm_data(l).label ' ε'];
+ end
end
end
end
@@ -179,12 +193,16 @@ function plotSigmas(gmm_data, in, epi)
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');
+ 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
@@ -192,6 +210,7 @@ end
function response = readGMvD(urlbase, gmms, in)
url = urlbase;
+
for i = 1 : size(gmms, 2)
if i == 1
url = url + "gmm=" + gmms(i);
@@ -199,6 +218,7 @@ for i = 1 : size(gmms, 2)
url = url + "&gmm=" + gmms(i);
end
end
+
url = url + ...
"&Mw=" + num2str(in.Mw) + ...
"&dip=" + num2str(in.dip) + ...
@@ -212,7 +232,8 @@ url = url + ...
"&vs30=" + num2str(in.vs30) + ...
"&vsInf=" + string(in.vsInf) + ...
"&zHyp=" + num2str(in.zHyp) + ...
- "&zTop=" + num2str(in.zTor);
+ "&zTop=" + num2str(in.zTor) + ...
+ "&imt=" + string(in.imt);
if isfield(in, 'z2p5')
url = url + num2str(in.z2p5);
diff --git a/etc/matlab/GroundMotionvsMagnitude.m b/etc/matlab/GroundMotionvsMagnitude.m
index bbe0b25..595914f 100644
--- a/etc/matlab/GroundMotionvsMagnitude.m
+++ b/etc/matlab/GroundMotionvsMagnitude.m
@@ -5,7 +5,7 @@
% Author: Peter Powers (pmpowers@usgs.gov)
% Demi Girot (dgirot@usgs.gov)
%
-% Updated: 05/23/2023
+% Updated: 05/13/2024
%
% This script demonstrates how to call web services for the ground motion
% models implemented in USGS National Sesmic Hazard Models and plot the
@@ -51,14 +51,20 @@ in.zTor = 0.5;
% 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"];
+
+%%
figure
plot_handles = [];
% Set to true to show epistemic branches
epi = false;
-gmms_2023 = ["ASK_14" "BSSA_14" "CB_14" "CY_14"];
-
subplot(2,1,1)
means = readGMvMw(urlbase, gmms_2023, in).response.means.data;
plotGmms(means, in, epi);
@@ -89,7 +95,8 @@ function plotGmms(gmm_data, in, epi)
% Plot the total spectrum
ph = loglog(gmm_xs, gmm_ys, ...
'LineWidth',3, ...
- 'LineStyle','-');
+ 'LineStyle','-', ...
+ 'Color',in.colors(j));
color = get(ph,'Color');
hold on;
@@ -100,7 +107,6 @@ function plotGmms(gmm_data, in, epi)
% Plot epistemic spectra, if present
if epi && ~isempty(gmm_tree)
for k = 1 : length(gmm_tree)
- id = gmm_tree(k).id;
epi_ys = gmm_tree(k).values;
phe = loglog(gmm_xs, epi_ys, ...
@@ -108,8 +114,10 @@ function plotGmms(gmm_data, in, epi)
'LineStyle',':', ...
'Color',color);
- plot_handles = [plot_handles phe];
- legend_labels{end+1} = [gmm_data(j).label ' ' id];
+ if k == 1
+ plot_handles = [plot_handles phe];
+ legend_labels{end+1} = [gmm_data(j).label ' ε'];
+ end
end
end
end
@@ -121,7 +129,12 @@ function plotGmms(gmm_data, in, epi)
title(['M',num2str(in.Mw),', ',num2str(in.rRup),'km'],'FontSize',20)
set(gca,'FontSize',20);
- legend(plot_handles,legend_labels,'Location','best');
+ legend( ...
+ plot_handles,legend_labels, ...
+ 'Location','best', ...
+ 'NumColumns',2, ...
+ 'FontSize',12, ...
+ 'Interpreter', 'none');
end
function plotSigmas(gmm_data, in, epi)
@@ -130,11 +143,11 @@ function plotSigmas(gmm_data, in, epi)
legend_labels = {};
% loop sigmas response array
- for m = 1 : length(gmm_data)
+ for l = 1 : length(gmm_data)
- sig_xs = gmm_data(m).data.xs;
- sig_ys = gmm_data(m).data.ys;
- sig_tree = gmm_data(m).tree;
+ 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
@@ -144,27 +157,29 @@ function plotSigmas(gmm_data, in, epi)
% Plot the total spectrum
ph = loglog(sig_xs, sig_ys, ...
'LineWidth',3, ...
- 'LineStyle','-');
+ 'LineStyle','-', ...
+ 'Color',in.colors(l));
color = get(ph,'Color');
hold on;
plot_handles = [plot_handles ph];
- legend_labels{end+1} = gmm_data(m).label;
+ legend_labels{end+1} = gmm_data(l).label;
% Plot epistemic spectra, if present
if epi && ~isempty(sig_tree)
- for n = 1 : length(sig_tree)
- id = sig_tree(n).id;
- epi_ys = sig_tree(n).values;
+ for m = 1 : length(sig_tree)
+ epi_ys = sig_tree(m).values;
phe = loglog(sig_xs, epi_ys, ...
'LineWidth',1, ...
'LineStyle',':', ...
'Color',color);
-
- plot_handles = [plot_handles phe];
- legend_labels{end+1} = [gmm_data(m).label ' ' id];
+
+ if m == 1
+ plot_handles = [plot_handles phe];
+ legend_labels{end+1} = [gmm_data(l).label ' ε'];
+ end
end
end
end
@@ -176,7 +191,12 @@ function plotSigmas(gmm_data, in, epi)
title(['M',num2str(in.Mw),', ',num2str(in.rRup),'km'],'FontSize',20)
set(gca,'FontSize',20);
- legend(plot_handles,legend_labels,'Location','best');
+ 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
@@ -184,6 +204,7 @@ end
function response = readGMvMw(urlbase, gmms, in)
url = urlbase;
+
for i = 1 : size(gmms, 2)
if i == 1
url = url + "gmm=" + gmms(i);
@@ -191,6 +212,7 @@ for i = 1 : size(gmms, 2)
url = url + "&gmm=" + gmms(i);
end
end
+
url = url + ...
"&Mw=" + num2str(in.Mw) + ...
"&dip=" + num2str(in.dip) + ...
@@ -204,7 +226,8 @@ url = url + ...
"&vs30=" + num2str(in.vs30) + ...
"&vsInf=" + string(in.vsInf) + ...
"&zHyp=" + num2str(in.zHyp) + ...
- "&zTop=" + num2str(in.zTor);
+ "&zTop=" + num2str(in.zTor) + ...
+ "&imt=" + string(in.imt);
if isfield(in, 'z2p5')
url = url + num2str(in.z2p5);
diff --git a/etc/matlab/ResponseSpectrum.m b/etc/matlab/ResponseSpectrum.m
index 3ceec53..8b9d84a 100644
--- a/etc/matlab/ResponseSpectrum.m
+++ b/etc/matlab/ResponseSpectrum.m
@@ -5,7 +5,7 @@
% Author: Peter Powers (pmpowers@usgs.gov)
% Demi Girot (dgirot@usgs.gov)
%
-% Updated: 05/16/2023
+% Updated: 05/13/2024
%
% This script demonstrates how to call web services for the ground motion
% models implemented in USGS National Sesmic Hazard Models and plot the
@@ -20,8 +20,7 @@
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Plot a figure of subplots for GMM response spectra, sigmas, and
-% epistemic uncertainty
+%% Plot a figure of GMM reponse spectra, sigmas, and epistemic uncertainty
clearvars
@@ -43,25 +42,31 @@ in.vs30 = 760;
in.zHyp = 7.5;
in.zTor = 10;
+in.colors = ["#D95319" "#0072BD" "#7E2F8E" "#77AC30" "#EDB120"];
+
+gmms_2023 = ["ASK_14" "BSSA_14" "CB_14" "CY_14"];
+
+%%
figure
plot_handles = [];
% Set to true to show epistemic branches
epi = false;
-gmms_2023 = ["ASK_14" "BSSA_14" "CB_14" "CY_14"];
+% Set to true to show PGA
+pga = false;
subplot(2,1,1)
means = readSpectra(urlbase, gmms_2023, in).response.means.data;
-plotGmms(means, in, epi);
+plotGmms(means, in, epi, pga);
subplot(2,1,2)
sigmas = readSpectra(urlbase, gmms_2023, in).response.sigmas.data;
-plotSigmas(sigmas, in, epi);
+plotSigmas(sigmas, in, epi, pga);
%% Shared Functions
-function plotGmms(gmm_data, in, epi)
+function plotGmms(gmm_data, in, epi, pga)
plot_handles = [];
legend_labels = {};
@@ -69,8 +74,9 @@ function plotGmms(gmm_data, in, epi)
% 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_xs = gmm_data(j).data.sa.xs;
+ gmm_ys = gmm_data(j).data.sa.ys;
+ gmm_pga = gmm_data(j).data.pga;
gmm_tree = gmm_data(j).tree;
% if 0.01 absent use PGA (0.001 is used for PGA)
@@ -81,27 +87,63 @@ function plotGmms(gmm_data, in, epi)
% Plot the total spectrum
ph = loglog(gmm_xs, gmm_ys, ...
'LineWidth',3, ...
- 'LineStyle','-');
+ 'LineStyle','-', ...
+ 'Color',in.colors(j));
color = get(ph,'Color');
hold on;
-
+
plot_handles = [plot_handles ph];
legend_labels{end+1} = gmm_data(j).label;
+
+ % Plot PGA, if present
+ if pga && ~isempty(gmm_pga)
+
+ phpga = loglog(0.0085,gmm_pga, ...
+ 'LineStyle','none', ...
+ 'Marker','square', ...
+ 'MarkerFaceColor',color, ...
+ 'MarkerSize',12);
+
+ % Uncomment to show in legend
+ % plot_handles = [plot_handles phpga];
+ % legend_labels{end+1} = [gmm_data(j).label ' PGA'];
+ end
% Plot epistemic spectra, if present
if epi && ~isempty(gmm_tree)
for k = 1 : length(gmm_tree)
- id = gmm_tree(k).id;
- epi_ys = gmm_tree(k).values;
+ epi_ys = gmm_tree(k).values.sa;
phe = loglog(gmm_xs, epi_ys, ...
'LineWidth',1, ...
'LineStyle',':', ...
'Color',color);
- plot_handles = [plot_handles phe];
- legend_labels{end+1} = [gmm_data(j).label ' ' id];
+ if k == 1
+ plot_handles = [plot_handles phe];
+ legend_labels{end+1} = [gmm_data(j).label ' ε'];
+ end
+ end
+ end
+
+ % Plot epistemic PGA, if present
+ if pga && ~isempty(gmm_tree)
+ for l = 1 : length(gmm_tree)
+ gmm_tree_pga = gmm_tree(l).values.pga;
+
+ phpga = loglog(0.0085,gmm_tree_pga, ...
+ 'LineWidth',2, ...
+ 'LineStyle','none', ...
+ 'Color',color, ...
+ 'Marker','square', ...
+ 'MarkerSize',7); hold on
+
+ % Uncomment to show in legend
+ % if l == 1
+ % plot_handles = [plot_handles phpga];
+ % legend_labels{end+1} = [gmm_data(j).label ' ε PGA'];
+ % end
end
end
end
@@ -111,19 +153,20 @@ function plotGmms(gmm_data, in, epi)
xlabel('Spectral Period (s)','FontSize',16)
ylabel('Median Ground Motion (g)','FontSize',16)
title(['M',num2str(in.Mw),', ',num2str(in.rRup),'km'],'FontSize',20)
- xlim([0.01 10]);
-
+ xlim([0.007 10]);
set(gca,'FontSize',20);
- set(gca,'XTick',[0.01 0.1 1 10]);
- set(gca,'XTickLabel',{'0.01','0.1','1','10'});
+ set(gca,'XTick',[0.001 0.01 0.1 1 10]);
+ set(gca,'XTickLabel',{'0.001','0.01','0.1','1','10'});
legend( ...
plot_handles,legend_labels, ...
'Location','best', ...
+ 'NumColumns',2, ...
+ 'FontSize',12, ...
'Interpreter', 'none');
end
-function plotSigmas(gmm_data, in, epi)
+function plotSigmas(gmm_data, in, epi, pga)
plot_handles = [];
legend_labels = {};
@@ -131,10 +174,11 @@ function plotSigmas(gmm_data, in, epi)
% loop sigmas response array
for m = 1 : length(gmm_data)
- sig_xs = gmm_data(m).data.xs;
- sig_ys = gmm_data(m).data.ys;
+ sig_xs = gmm_data(m).data.sa.xs;
+ sig_ys = gmm_data(m).data.sa.ys;
+ sig_pga = gmm_data(m).data.pga;
sig_tree = gmm_data(m).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;
@@ -143,7 +187,8 @@ function plotSigmas(gmm_data, in, epi)
% Plot the total spectrum
ph = loglog(sig_xs, sig_ys, ...
'LineWidth',3, ...
- 'LineStyle','-');
+ 'LineStyle','-', ...
+ 'Color',in.colors(m));
color = get(ph,'Color');
hold on;
@@ -151,37 +196,72 @@ function plotSigmas(gmm_data, in, epi)
plot_handles = [plot_handles ph];
legend_labels{end+1} = gmm_data(m).label;
+ % Plot PGA, if present
+ if pga && ~isempty(sig_pga)
+
+ phpga = loglog(0.0085,sig_pga, ...
+ 'LineStyle','none', ...
+ 'Marker','square', ...
+ 'MarkerFaceColor',color, ...
+ 'MarkerSize',12);
+
+ % Uncomment to show in legend
+ % plot_handles = [plot_handles phpga];
+ % legend_labels{end+1} = [gmm_data(m).label ' PGA'];
+ end
+
% Plot epistemic spectra, if present
if epi && ~isempty(sig_tree)
for n = 1 : length(sig_tree)
- id = sig_tree(n).id;
- epi_ys = sig_tree(n).values;
+ epi_ys = sig_tree(n).values.sa;
phe = loglog(sig_xs, epi_ys, ...
'LineWidth',1, ...
'LineStyle',':', ...
'Color',color);
-
- plot_handles = [plot_handles phe];
- legend_labels{end+1} = [gmm_data(m).label ' ' id];
+
+ if n == 1
+ plot_handles = [plot_handles phe];
+ legend_labels{end+1} = [gmm_data(m).label ' ε'];
+ end
end
end
- end
+ % Plot epistemic PGA, if present
+ if pga && ~isempty(sig_tree)
+ for o = 1 : length(sig_tree)
+ sig_tree_pga = sig_tree(o).values.pga;
+
+ phpga = loglog(0.0085,sig_tree_pga, ...
+ 'LineWidth',2, ...
+ 'LineStyle','none', ...
+ 'Color',color, ...
+ 'Marker','square', ...
+ 'MarkerSize',7); hold on
+
+ % Uncomment to show in legend
+ % if o == 1
+ % plot_handles = [plot_handles phpga];
+ % legend_labels{end+1} = [gmm_data(m).label ' ε PGA'];
+ % end
+ end
+ end
+ end
+
grid on
xlabel('Spectral Period (s)','FontSize',16)
ylabel('Standard Deviation','FontSize',16)
title(['M',num2str(in.Mw),', ',num2str(in.rRup),'km'],'FontSize',20)
- xlim([0.01 10]);
-
+ xlim([0.007 10]);
set(gca,'FontSize',20);
- set(gca,'XTick',[0.01 0.1 1 10]);
- set(gca,'XTickLabel',{'0.01','0.1','1','10'});
+ set(gca,'XTick',[0.001 0.01 0.1 1 10]);
+ set(gca,'XTickLabel',{'0.001','0.01','0.1','1','10'});
legend( ...
plot_handles,legend_labels, ...
'Location','best', ...
+ 'FontSize',12, ...
'Interpreter', 'none');
end
@@ -190,6 +270,7 @@ end
function response = readSpectra(urlbase, gmms, in)
url = urlbase;
+
for i = 1 : size(gmms, 2)
if i == 1
url = url + "gmm=" + gmms(i);
@@ -197,6 +278,7 @@ for i = 1 : size(gmms, 2)
url = url + "&gmm=" + gmms(i);
end
end
+
url = url + ...
"&Mw=" + num2str(in.Mw) + ...
"&dip=" + num2str(in.dip) + ...
--
GitLab