ResponseSpectrum.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Ground motion model (GMM) response spectra subplots
%
% Author: Peter Powers (pmpowers@usgs.gov)
%         Demi Girot (dgirot@usgs.gov)
%
% 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
% 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 reponse 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;

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;

% 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, pga);

subplot(2,1,2)
sigmas = readSpectra(urlbase, gmms_2023, in).response.sigmas.data;
plotSigmas(sigmas, in, epi, pga);

%% Shared Functions

function plotGmms(gmm_data, in, epi, pga)

    plot_handles = [];
    legend_labels = {};

    % loop means response array
    for j = 1 : length(gmm_data) 
    
        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)
        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',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)
                epi_ys = gmm_tree(k).values.sa;
                
                phe = loglog(gmm_xs, epi_ys, ...
                    'LineWidth',1, ...
                    'LineStyle',':', ...
                    'Color',color); 

                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
    
    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.007 10]);
    set(gca,'FontSize',20);
    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, pga)

    plot_handles = [];
    legend_labels = {};

    % loop sigmas response array
    for m = 1 : length(gmm_data) 
    
        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;
        end

        % Plot the total spectrum
        ph = loglog(sig_xs, sig_ys, ...
            'LineWidth',3, ...
            'LineStyle','-', ...
            'Color',in.colors(m));
        
        color = get(ph,'Color');
        hold on;
        
        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)
                epi_ys = sig_tree(n).values.sa;
                
                phe = loglog(sig_xs, epi_ys, ...
                    'LineWidth',1, ...
                    'LineStyle',':', ...
                    'Color',color); 
                
                if n == 1
                    plot_handles = [plot_handles phe];
                    legend_labels{end+1} = [gmm_data(m).label ' ε'];
                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.007 10]);
    set(gca,'FontSize',20);
    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

% 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