Skip to content
Snippets Groups Projects
Commit c1e2a50d authored by Powers, Peter M.'s avatar Powers, Peter M.
Browse files

Merge branch 'matlab_scripts_update' into 'main' 

Added updated RS GMvD GMvMw scripts

See merge request !229
parents 1d56f74d 81dc2a27
No related branches found
No related tags found
1 merge request!229Added updated RS GMvD GMvMw scripts
Pipeline #445819 passed
Stage: build
Stage: trigger
Hide whitespace changes
Inline Side-by-side
etc/matlab/GroundMotionvsDistance.m
+ 46
− 25
View file @ c1e2a50d
......@@ -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);
......
etc/matlab/GroundMotionvsMagnitude.m
+ 45
− 22
View file @ c1e2a50d
......@@ -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);
......
etc/matlab/ResponseSpectrum.m
+ 116
− 34
View file @ c1e2a50d
......@@ -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) + ...
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment