From aafa096e7688d2fded6edc434be25465fcb737bd Mon Sep 17 00:00:00 2001
From: Peter Powers <pmpowers@usgs.gov>
Date: Thu, 19 Jan 2023 14:23:22 -0700
Subject: [PATCH] readme updates; matlab example

---
 README.md                     |  28 ++++--
 etc/matlab/README.md          |   7 ++
 etc/matlab/ResponseSpectrum.m | 156 ++++++++++++++++++++++++++++++++++
 3 files changed, 185 insertions(+), 6 deletions(-)
 create mode 100644 etc/matlab/README.md
 create mode 100644 etc/matlab/ResponseSpectrum.m

diff --git a/README.md b/README.md
index 66edd7e..e68146a 100644
--- a/README.md
+++ b/README.md
@@ -27,13 +27,29 @@ After startup, web services and documentation are available at <http://localhost
 
 #### Customizing Code
 
-It is not uncommon for users to want to customize the behavior of the code. Whereas *nshmp-ws*
-contains code to run web services, ground motion model (GMM) calculations are handled in the
-dependent library [*nshmp-lib*](https://code.usgs.gov/ghsc/nshmp/nshmp-lib). To use a local,
-modified version of *nshmp-lib*, one need only to set an environment variable
+Whereas *nshmp-ws* contains code to run web services, ground motion model (GMM) calculations
+are handled in the dependent library [*nshmp-lib*](https://code.usgs.gov/ghsc/nshmp/nshmp-lib).
+To use a local, modified version of *nshmp-lib*, set an environment variable
 `NSHMP_LIB_LOCAL=true` and *nshmp-ws* will look for *nshmp-lib* in a directory adjacent to
-*nshmp-haz*. If *nshmp-lib* is located somewhere else, modify the path specified in
-[`gradle/dependencies.gradle`](../../gradle/dependencies.gradle).
+*nshmp-ws*. If *nshmp-lib* is located somewhere else, modify the path specified in
+[`gradle/dependencies.gradle`](../../gradle/dependencies.gradle). When using a local version
+of *nshmp-lib*, first build the *nshmp-lib* project using `./gradlew fatJar` so that
+required dependencies are included.
+
+Summary of steps to check out and start running code for local development:
+
+bash
+```
+# --> Set a NSHMP_LIB_LOCAL=true environment variable on your system
+cd yourProjectDirectory
+git clone https://code.usgs.gov/ghsc/nshmp/nshmp-lib.git
+cd nshmp-lib
+./gradlew fatJar
+cd ..
+git clone https://code.usgs.gov/ghsc/nshmp/nshmp-ws.git
+cd nshmp-ws
+./gradlew run
+```
 
 ### Docker
 
diff --git a/etc/matlab/README.md b/etc/matlab/README.md
new file mode 100644
index 0000000..d0f7324
--- /dev/null
+++ b/etc/matlab/README.md
@@ -0,0 +1,7 @@
+# Using Ground Motion Models (GMMs) from Matlab
+
+The scripts in this directory provide examples of how to access web services for NSHM GMMs.
+
+Users are *strongly* encouraged to run the web services locally so as to not overburden USGS
+servers. See the [README](../../README.md) at the root of this repository for instructions on
+how to do so.
diff --git a/etc/matlab/ResponseSpectrum.m b/etc/matlab/ResponseSpectrum.m
new file mode 100644
index 0000000..a3d0b7b
--- /dev/null
+++ b/etc/matlab/ResponseSpectrum.m
@@ -0,0 +1,156 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Ground motion model (GMM) response spectra plots
+%
+% Author: Peter Powers (pmpowers@usgs.gov)
+%
+% Updated: 01/19/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 GMM reponse spectra 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-haz 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"];
+means = readSpectra(urlbase, gmms_2023, in).response.means.data;
+plotGmms(means, in, epi);
+
+%% Shared Functions
+
+function plotGmms(gmm_data, in, epi)
+
+    plot_handles = [];
+    legend_labels = {};
+
+    % loop means_new 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','southwest', ...
+        '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
-- 
GitLab