ResponseSpectrum.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Ground motion model (GMM) response spectra subplots
%
% Author: Peter Powers (pmpowers@usgs.gov)
%         Demi Girot (dgirot@usgs.gov)
%
% Updated: 05/16/2023
%
% 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 subplots for GMM response spectra, sigmas, and 
% epistemic uncertainty

clearvars

% The root url for the response spectra web service
urlbase = "https://earthquake.usgs.gov/ws/nshmp/data/gmm/spectra?";

% If you are running nshmp-ws locally, use this URL:
% urlbase = "http://localhost:8080/gmm/spectra?";

% Struct of ground motion model parameters
in.Mw = 7.5;
in.dip = 90;
in.rake = 0;
in.width = 15;
in.rJB = 50;
in.rRup = 50;
in.rX = 50;
in.vs30 = 760;
in.zHyp = 7.5;
in.zTor = 10;

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 = readSpectra(urlbase, gmms_2023, in).response.means.data;
plotGmms(means, in, epi);

subplot(2,1,2)
sigmas = readSpectra(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)
                id = gmm_tree(k).id;
                epi_ys = gmm_tree(k).values;
                
                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];
            end
        end
    end
    
    grid on

    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]);

    set(gca,'FontSize',20);
    set(gca,'XTick',[0.01 0.1 1 10]);
    set(gca,'XTickLabel',{'0.01','0.1','1','10'});
    
    legend( ...
        plot_handles,legend_labels, ...
        'Location','best', ...
        'Interpreter', 'none');
end

function plotSigmas(gmm_data, in, epi)

    plot_handles = [];
    legend_labels = {};

    % 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_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;
        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(m).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;
                
                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];
            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]);

    set(gca,'FontSize',20);
    set(gca,'XTick',[0.01 0.1 1 10]);
    set(gca,'XTickLabel',{'0.01','0.1','1','10'});
    
    legend( ...
        plot_handles,legend_labels, ...
        'Location','best', ...
        'Interpreter', 'none');
end

% Read reponse spectra from a NSHM web service for the supplied
% ground motion models (GMMs) and GMM input parameters
function response = readSpectra(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) + ...
    "&rRup=" + num2str(in.rRup) + ...
    "&rX=" + num2str(in.rX) + ...
    "&vs30=" + num2str(in.vs30) + ...
    "&zHyp=" + num2str(in.zHyp) + ...
    "&zTop=" + num2str(in.zTor);

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