From a8a213f0b5c3b128206ec94a3f6e4e38af537d72 Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Fri, 2 Feb 2024 17:07:49 +0100
Subject: [PATCH 01/19] adding function to create cs coords ready for plotting

---
 openquake/sub/cross_sections.py              | 10 ++--
 openquake/sub/make_cs_coords.py              | 61 ++++++++++++++++++++
 openquake/sub/plotting/plot_cross_section.py |  6 +-
 3 files changed, 69 insertions(+), 8 deletions(-)
 create mode 100755 openquake/sub/make_cs_coords.py

diff --git a/openquake/sub/cross_sections.py b/openquake/sub/cross_sections.py
index 33b90a155..448cc8dee 100644
--- a/openquake/sub/cross_sections.py
+++ b/openquake/sub/cross_sections.py
@@ -217,13 +217,13 @@ def set_slab1pt0(self, filename, bffer=2.0):
             slab1pt0[idx[0], 0] = slab1pt0[idx[0], 0] - 360.
         if qual == 0:
             minlo, maxlo, minla, maxla, qual = self.csec.get_mm(2.0)
-            idxslb = self.csec.get_grd_nodes_within_buffer(slab1pt0[:, 0],
+            idxslb, dst = self.csec.get_grd_nodes_within_buffer(slab1pt0[:, 0],
                                                            slab1pt0[:, 1],
                                                            bffer,
                                                            minlo, maxlo,
                                                            minla, maxla)
         if qual == 1:
-            idxslb = self.csec.get_grd_nodes_within_buffer_idl(slab1pt0[:, 0],
+            idxslb, dst = self.csec.get_grd_nodes_within_buffer_idl(slab1pt0[:, 0],
                                                                slab1pt0[:, 1],
                                                                bffer,
                                                                minlo, maxlo,
@@ -288,7 +288,7 @@ def set_litho_moho_depth(self, filename, bffer=100.):
         if idxl is not None and len(idxl):
             boo = numpy.zeros_like(dataa[:, 0], dtype=int)
             boo[idxl[0]] = 1
-            self.litho = numpy.squeeze(dataa[idxl, :])
+            self.litho = numpy.squeeze(dataa[idxl[0], :])
 
     def set_gcmt(self, filename, bffer=75.):
         """
@@ -351,7 +351,7 @@ def set_topo(self, filename, bffer=0.25):
         if idxb is not None and len(idxb):
             boo = numpy.zeros_like(datab[:, 0], dtype=int)
             boo[idxb[0]] = 1
-            self.topo = numpy.squeeze(datab[idxb, :])
+            self.topo = numpy.squeeze(datab[idxb[0], :])
 
     def set_volcano(self, filename, bffer=75.):
         """
@@ -386,7 +386,7 @@ def set_volcano(self, filename, bffer=75.):
         if idxv is not None and len(idxv):
             voo = numpy.zeros_like(vulc[:, 0], dtype=int)
             voo[idxv[0]] = 1
-            self.volc = numpy.squeeze(vulc[idxv, :])
+            self.volc = numpy.squeeze(vulc[idxv[0], :])
         fin.close()
         print(self.volc)
 
diff --git a/openquake/sub/make_cs_coords.py b/openquake/sub/make_cs_coords.py
new file mode 100755
index 000000000..c94e52500
--- /dev/null
+++ b/openquake/sub/make_cs_coords.py
@@ -0,0 +1,61 @@
+#!/usr/bin/env python
+# ------------------- The OpenQuake Model Building Toolkit --------------------
+# Copyright (C) 2022 GEM Foundation
+#           _______  _______        __   __  _______  _______  ___   _
+#          |       ||       |      |  |_|  ||  _    ||       ||   | | |
+#          |   _   ||   _   | ____ |       || |_|   ||_     _||   |_| |
+#          |  | |  ||  | |  ||____||       ||       |  |   |  |      _|
+#          |  |_|  ||  |_|  |      |       ||  _   |   |   |  |     |_
+#          |       ||      |       | ||_|| || |_|   |  |   |  |    _  |
+#          |_______||____||_|      |_|   |_||_______|  |___|  |___| |_|
+#
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Affero General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public License for more
+# details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+# -----------------------------------------------------------------------------
+# vim: tabstop=4 shiftwidth=4 softtabstop=4
+# coding: utf-8
+
+import os
+import glob
+import pandas as pd
+
+from openquake.baselib import sap
+from openquake.hazardlib.geo.geodetic import azimuth
+
+def make_cs_coords(cs_dir, outfi):
+    cs_files = glob.glob(f'{cs_dir}/*csv')
+    lines = []
+    for fi in cs_files:
+        sz = os.path.getsize(fi)
+        if sz == 0:
+            continue
+
+        df = pd.read_csv(fi, sep=' ', names=["lon", "lat", "depth"])
+        az = azimuth(df.lon[0], df.lat[0], df.lon.values[-1], df.lat.values[-1])
+
+        csid = fi.split('/')[-1][3:].replace('.csv','')
+        line = f'{df.lon[0]} {df.lat[0]} 700.0 300.0 {az} {csid} cs.ini \n'
+        lines.append(line) 
+    
+    os.remove(outfi) if os.path.exists(outfi) else None
+
+    with open(outfi, 'w') as f:
+        for line in lines:
+            f.write(line)
+    print(f'Written to {outfi}')
+
+make_cs_coords.cs_dir = 'directory with cross section coordinates'
+make_cs_coords.outfi = 'output filename'
+
+if __name__ == "__main__":
+    sap.run(make_cs_coords)
diff --git a/openquake/sub/plotting/plot_cross_section.py b/openquake/sub/plotting/plot_cross_section.py
index a6a8b0672..f23ccaa87 100755
--- a/openquake/sub/plotting/plot_cross_section.py
+++ b/openquake/sub/plotting/plot_cross_section.py
@@ -345,7 +345,7 @@ def _plot_volc(axes, csda):
     ola = csda.csec.ola
     patches = []
 
-    if (len(csda.volc)-1) >= 1:
+    if (len(csda.volc)-1) >= 2:
         vuls = geodetic_distance(olo, ola,
                                  csda.volc[:, 0],
                                  csda.volc[:, 1])
@@ -564,7 +564,7 @@ def plt_cs(olo, ola, depp, lnght, strike, ids, ini_filename):
 
     config = configparser.ConfigParser()
     config.read(ini_filename)
-    fname_trench = config['data']['trench_axis_filename']
+#    fname_trench = config['data']['trench_axis_filename']
     fname_eqk_cat = config['data']['catalogue_pickle_filename']
     fname_slab = config['data']['slab1pt0_filename']
     fname_crust = config['data']['crust1pt0_filename']
@@ -574,7 +574,7 @@ def plt_cs(olo, ola, depp, lnght, strike, ids, ini_filename):
     fname_litho = config['data']['litho_filename']
     fname_volc = config['data']['volc_filename']
 
-    csda.set_trench_axis(fname_trench)
+#    csda.set_trench_axis(fname_trench)
     cat = pickle.load(open(fname_eqk_cat, 'rb'))
     csda.set_catalogue(cat)
     if re.search('[a-z]', fname_slab):

From 54fe050babb49a8078926da777d78fa07aa5dd29 Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Fri, 2 Feb 2024 18:20:21 +0100
Subject: [PATCH 02/19] volcanos

---
 openquake/sub/plotting/plot_cross_section.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/openquake/sub/plotting/plot_cross_section.py b/openquake/sub/plotting/plot_cross_section.py
index f23ccaa87..76bba90c4 100755
--- a/openquake/sub/plotting/plot_cross_section.py
+++ b/openquake/sub/plotting/plot_cross_section.py
@@ -345,7 +345,11 @@ def _plot_volc(axes, csda):
     ola = csda.csec.ola
     patches = []
 
-    if (len(csda.volc)-1) >= 2:
+    
+
+    if (len(csda.volc)-1) >= 1:
+        if str(csda.volc.shape) == '(2,)':
+            csda.volc = csda.volc.reshape(1,2)
         vuls = geodetic_distance(olo, ola,
                                  csda.volc[:, 0],
                                  csda.volc[:, 1])
@@ -354,6 +358,7 @@ def _plot_volc(axes, csda):
             patches.append(square)
 
     else:
+        breakpoint()
         vuls = geodetic_distance(olo, ola,
                                  csda.volc[0],
                                  csda.volc[1])

From 644114784a97c2da37ffb18848e1c0a8453e48cd Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Fri, 16 Feb 2024 15:28:21 +0100
Subject: [PATCH 03/19] updating plotting tools for subduction sources

---
 openquake/mbi/sub/make_cs_coords.py           | 25 +++++++
 openquake/mbi/sub/plot_cross_sections_map.py  | 23 +++++++
 .../sub/create_multiple_cross_sections.py     | 29 +++++----
 openquake/sub/make_cs_coords.py               | 14 ++--
 openquake/sub/plotting/plot_cross_section.py  |  2 +-
 .../plotting/plot_multiple_cross_sections.py  |  3 +-
 .../plot_multiple_cross_sections_map.py       | 65 ++++++++++++-------
 openquake/sub/tests/cross_section_test.py     | 17 ++++-
 openquake/sub/tests/expected/cs_file.cs       |  2 +
 requirements_linux.txt                        |  1 +
 requirements_macos_arm64.txt                  |  1 +
 requirements_macos_x86_64.txt                 |  1 +
 requirements_win64.txt                        |  1 +
 13 files changed, 139 insertions(+), 45 deletions(-)
 create mode 100755 openquake/mbi/sub/make_cs_coords.py
 create mode 100755 openquake/mbi/sub/plot_cross_sections_map.py
 create mode 100644 openquake/sub/tests/expected/cs_file.cs

diff --git a/openquake/mbi/sub/make_cs_coords.py b/openquake/mbi/sub/make_cs_coords.py
new file mode 100755
index 000000000..d29ac211c
--- /dev/null
+++ b/openquake/mbi/sub/make_cs_coords.py
@@ -0,0 +1,25 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+from openquake.baselib import sap
+from openquake.sub.make_cs_coords import make_cs_coords
+
+def main(cs_dir, outfi, ini_fname, cs_length=300., cs_depth=300.):
+    """
+    Creates file with parameters needed to plot cross sections.
+    Output file is a list for each cross section with the format:
+
+        lon lat length depth azimuth id <config>.ini
+    """
+
+    make_cs_coords(cs_dir, outfi, ini_fname, cs_length, cs_depth)
+
+
+make_cs_coords.cs_dir = 'directory with cross section coordinates'
+make_cs_coords.outfi = 'output filename'
+make_cs_coords.ini_fname = 'name of ini file specifying data paths'
+make_cs_coords.cs_length = 'length of cross sections (default 300)'
+make_cs_coords.cs_depth = 'depth extent of cross sections (default 300 km)'
+
+if __name__ == '__main__':
+    sap.run(main)
diff --git a/openquake/mbi/sub/plot_cross_sections_map.py b/openquake/mbi/sub/plot_cross_sections_map.py
new file mode 100755
index 000000000..e03d72cef
--- /dev/null
+++ b/openquake/mbi/sub/plot_cross_sections_map.py
@@ -0,0 +1,23 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+from openquake.baselib import sap
+from openquake.sub.plotting.plot_multiple_cross_sections_map import plot
+
+
+def main(config_fname, cs_file=None):
+    """
+    Creates file with parameters needed to plot cross sections.
+    Output file is a list for each cross section with the format:
+
+        lon lat length depth azimuth id <config>.ini
+    """
+
+    plot(config_fname, cs_file)
+
+
+main.config_fname = 'config file to datasets'
+main.cs_file = 'existing cross sections details'
+
+if __name__ == "__main__":
+    sap.run(main)
diff --git a/openquake/sub/create_multiple_cross_sections.py b/openquake/sub/create_multiple_cross_sections.py
index 4c4340625..11f700743 100644
--- a/openquake/sub/create_multiple_cross_sections.py
+++ b/openquake/sub/create_multiple_cross_sections.py
@@ -12,7 +12,7 @@
 from matplotlib.collections import PatchCollection
 
 from openquake.sub.cross_sections import Trench
-
+from openquake.baselib import sap
 
 def get_cs(trench, ini_filename, cs_len, cs_depth, interdistance, qual,
            fname_out_cs='cs_traces.cs'):
@@ -81,8 +81,6 @@ def plot(trench, cat, cs_dict, interdistance):
     midlo = (minlo+maxlo)/2
     midla = (minla+maxla)/2
 
-    # Plot the traces of cross-sections
-    ts = trench.resample(interdistance)
 
     _ = plt.figure(figsize=(12, 9))
 
@@ -95,7 +93,7 @@ def plot(trench, cat, cs_dict, interdistance):
     #
     # Draw paralleles and meridians with labels
     # labels = [left,right,top,bottom]
-    m.drawcoastlines()
+#    m.drawcoastlines()
     m.drawmeridians(numpy.arange(numpy.floor(minlo/10.)*10,
                                  numpy.ceil(maxlo/10.)*10, 5.),
                     labels=[False, False, False, True])
@@ -121,12 +119,17 @@ def plot(trench, cat, cs_dict, interdistance):
 
     x, y = m(trench.axis[:, 0], trench.axis[:, 1])
     plt.plot(x, y, '-g', linewidth=2, zorder=10)
-    x, y = m(ts.axis[:, 0], ts.axis[:, 1])
     plt.plot(x, y, '--y', linewidth=4, zorder=20)
 
-    for key in cs_dict:
-        cs = cs_dict[key]
-        if cs is not None:
+    # Plot the traces of cross-sections
+    if interdistance != 0:
+        ts = trench.resample(interdistance)
+        x, y = m(ts.axis[:, 0], ts.axis[:, 1])
+
+    else:
+        for key in cs_dict:
+            cs = cs_dict[key]
+        #if cs is not None:
             x, y = m(cs.plo, cs.pla)
             plt.plot(x, y, ':r', linewidth=2, zorder=20)
             text = plt.text(x[0], y[0], '%s' % key, ha='center',
@@ -135,17 +138,18 @@ def plot(trench, cat, cs_dict, interdistance):
                             va='center', size=10, zorder=30)
             text.set_path_effects([PathEffects.withStroke(linewidth=3,
                                                           foreground="w")])
+    m.drawcoastlines(zorder=26)
     plt.show()
 
 
-def main(argv):
+def main(config_fname):
     """
-    argv[0] is the .ini file
+    config_fname is the .ini file
     """
 
     # Parse .ini file
     config = configparser.ConfigParser()
-    config.read(argv[0])
+    config.read(config_fname)
     fname_trench = config['data']['trench_axis_filename']
     fname_eqk_cat = config['data']['catalogue_pickle_filename']
     cs_length = float(config['section']['lenght'])
@@ -179,6 +183,7 @@ def main(argv):
     if False:
         plot(trench, cat, cs_dict, interdistance)
 
+main.config = 'config file for creating cross sections from trench axis'
 
 if __name__ == "__main__":
-    main(sys.argv[1:])
+    sap.run(main)
diff --git a/openquake/sub/make_cs_coords.py b/openquake/sub/make_cs_coords.py
index c94e52500..f6c8c052f 100755
--- a/openquake/sub/make_cs_coords.py
+++ b/openquake/sub/make_cs_coords.py
@@ -32,8 +32,8 @@
 from openquake.baselib import sap
 from openquake.hazardlib.geo.geodetic import azimuth
 
-def make_cs_coords(cs_dir, outfi):
-    cs_files = glob.glob(f'{cs_dir}/*csv')
+def make_cs_coords(cs_dir, outfi, ini_fname, cs_length=300., cs_depth=300.):
+    cs_files = sorted(glob.glob(f'{cs_dir}/*csv'))
     lines = []
     for fi in cs_files:
         sz = os.path.getsize(fi)
@@ -44,9 +44,10 @@ def make_cs_coords(cs_dir, outfi):
         az = azimuth(df.lon[0], df.lat[0], df.lon.values[-1], df.lat.values[-1])
 
         csid = fi.split('/')[-1][3:].replace('.csv','')
-        line = f'{df.lon[0]} {df.lat[0]} 700.0 300.0 {az} {csid} cs.ini \n'
-        lines.append(line) 
-    
+        line = f'{df.lon[0]} {df.lat[0]} {cs_length} {cs_depth} '
+        line += f'{az:.4} {csid} {ini_fname} \n'
+        lines.append(line)
+
     os.remove(outfi) if os.path.exists(outfi) else None
 
     with open(outfi, 'w') as f:
@@ -56,6 +57,9 @@ def make_cs_coords(cs_dir, outfi):
 
 make_cs_coords.cs_dir = 'directory with cross section coordinates'
 make_cs_coords.outfi = 'output filename'
+make_cs_coords.ini_fname = 'name of ini file specifying data paths'
+make_cs_coords.cs_length = 'length of cross sections (default 300)'
+make_cs_coords.cs_depth = 'depth extent of cross sections (default 300 km)'
 
 if __name__ == "__main__":
     sap.run(make_cs_coords)
diff --git a/openquake/sub/plotting/plot_cross_section.py b/openquake/sub/plotting/plot_cross_section.py
index 76bba90c4..c78f7e03b 100755
--- a/openquake/sub/plotting/plot_cross_section.py
+++ b/openquake/sub/plotting/plot_cross_section.py
@@ -182,7 +182,7 @@ def _plot_slab1pt0(axes, csda):
     iii = numpy.argsort(slb_dst)
     if len(iii) > 2:
         plt.plot(slb_dst[iii], -1*slb_dep[iii], '-b', linewidth=3, zorder=30)
-        plt.text(slb_dst[iii[-1]], -1*slb_dep[iii[-1]], 'Slab1.0', fontsize=8)
+        plt.text(slb_dst[iii[-1]], -1*slb_dep[iii[-1]], 'Slab2.0', fontsize=8)
 
 
 def _plot_np_intersection(axes, csda):
diff --git a/openquake/sub/plotting/plot_multiple_cross_sections.py b/openquake/sub/plotting/plot_multiple_cross_sections.py
index 8aa56268b..01e7ca9f1 100755
--- a/openquake/sub/plotting/plot_multiple_cross_sections.py
+++ b/openquake/sub/plotting/plot_multiple_cross_sections.py
@@ -5,9 +5,10 @@
 import sys
 from plot_cross_section import plt_cs
 import matplotlib.pyplot as plt
+from openquake.baselib import sap
 
 
-def main(argv):
+def pcs(argv):
     if len(argv) < 2:
         folder = './'
     else:
diff --git a/openquake/sub/plotting/plot_multiple_cross_sections_map.py b/openquake/sub/plotting/plot_multiple_cross_sections_map.py
index 6809c9c97..35589295d 100755
--- a/openquake/sub/plotting/plot_multiple_cross_sections_map.py
+++ b/openquake/sub/plotting/plot_multiple_cross_sections_map.py
@@ -5,6 +5,8 @@
 import numpy
 import configparser
 import pickle
+
+from openquake.baselib import sap
 from openquake.sub.cross_sections import CrossSection, Trench
 from openquake.sub.create_multiple_cross_sections import plot as pcs
 
@@ -26,43 +28,56 @@ def read_cs(filename):
     return cs_dict
 
 
-def main(argv):
+def plot(config_file, cs_file=None):
     """
     """
 
     # Parse .ini file
     config = configparser.ConfigParser()
-    config.read(argv[0])
-    fname_trench = config['data']['trench_axis_filename']
-    fname_eqk_cat = config['data']['catalogue_pickle_filename']
-    interdistance = float(config['section']['interdistance'])
+    config.read(config_file)
+    if cs_file == None:
+        fname_trench = config['data']['trench_axis_filename']
+        intdist = float(config['section']['interdistance'])
+
+
+        # Load trench axis
+        fin = open(fname_trench, 'r')
+        lotmp = []
+        latmp = []
+        for line in fin:
+            aa = re.split('\\s+', re.sub('^\\s+', '', line))
+            lotmp.append(float(aa[0]))
+            latmp.append(float(aa[1]))
+        fin.close()
+        qual = 0
+        if (min(lotmp)/max(lotmp) < 0.) & (min(lotmp) < -150.):
+            # MN: 'qual' assigned but never used
+            qual = 1
+            lotmp = (x+360. if x < 0. else x for x in lotmp)
+        trench = list(zip(lotmp, latmp))
+        trench = Trench(numpy.array(trench))
+
+    else:
+        # Read cs file
+        cs_dict = read_cs(cs_file)
+        lotmp = [cs_dict[c].olo for c in cs_dict]
+        latmp = [cs_dict[c].ola for c in cs_dict]
+        trench = list(zip(lotmp, latmp))
+        trench = Trench(numpy.array(trench))
+        intdist = 0
 
-    # Load trench axis
-    fin = open(fname_trench, 'r')
-    lotmp = []
-    latmp = []
-    for line in fin:
-        aa = re.split('\\s+', re.sub('^\\s+', '', line))
-        lotmp.append(float(aa[0]))
-        latmp.append(float(aa[1]))
-    fin.close()
-    qual = 0
-    if (min(lotmp)/max(lotmp) < 0.) & (min(lotmp) < -150.):
-        # MN: 'qual' assigned but never used
-        qual = 1
-        lotmp = (x+360. if x < 0. else x for x in lotmp)
-    trench = list(zip(lotmp, latmp))
-    trench = Trench(numpy.array(trench))
 
     # Load catalogue
+    fname_eqk_cat = config['data']['catalogue_pickle_filename']
     cat = pickle.load(open(fname_eqk_cat, 'rb'))
 
-    # Read cs file
-    cs_dict = read_cs(argv[1])
 
     # Plotting
-    pcs(trench, cat, cs_dict, interdistance)
+    pcs(trench, cat, cs_dict, interdistance=intdist)
 
+plot.config_file = 'config file to datasets'
+plot.cs_file = 'existing cross sections details'
 
 if __name__ == "__main__":
-    main(sys.argv[1:])
+    #main(sys.argv[1:])
+    sap.run(plot)
diff --git a/openquake/sub/tests/cross_section_test.py b/openquake/sub/tests/cross_section_test.py
index 014f038f6..ef27ebec7 100644
--- a/openquake/sub/tests/cross_section_test.py
+++ b/openquake/sub/tests/cross_section_test.py
@@ -11,7 +11,8 @@
 
 from openquake.hazardlib.geo.point import Point
 from openquake.hazardlib.geo.line import Line
-
+from openquake.sub.tests import __file__ as tests__init__
+from openquake.sub.make_cs_coords import make_cs_coords
 
 BASE_PATH = os.path.dirname(__file__)
 tmp = 'data/crust/crust_except.xyz'
@@ -136,3 +137,17 @@ def test01(self):
         pnts = numpy.array([[10.3, 45.0], [9.0, 45.0]])
         line = Line([Point(lo, la) for lo, la in zip(lons, lats)])
         get_min_distance(line, pnts)
+
+
+class edge2profileTest(unittest.TestCase):
+
+    def test_cs_generator(self):
+        """
+        Tests cross section files generated from exisiting profiles
+        """
+
+        reference_file = os.path.join(BASE_PATH, 'expected', 'cs_file.cs')
+        cs_dir = os.path.join(BASE_PATH, 'data', 'cs')
+        outfi = os.path.join('/tmp/cs_file.cs')
+        make_cs_coords(cs_dir, outfi, 'cs.ini')
+        assert filecmp.cmp(outfi, reference_file) == True
diff --git a/openquake/sub/tests/expected/cs_file.cs b/openquake/sub/tests/expected/cs_file.cs
new file mode 100644
index 000000000..5982d08f2
--- /dev/null
+++ b/openquake/sub/tests/expected/cs_file.cs
@@ -0,0 +1,2 @@
+10.0 45.0 300.0 300.0 34.85 3 cs.ini 
+10.1 45.0 300.0 300.0 34.85 4 cs.ini 
diff --git a/requirements_linux.txt b/requirements_linux.txt
index 879d23487..c015d3e22 100644
--- a/requirements_linux.txt
+++ b/requirements_linux.txt
@@ -31,3 +31,4 @@ reportlab
 netCDF4
 igraph
 pyvista
+basemap
diff --git a/requirements_macos_arm64.txt b/requirements_macos_arm64.txt
index 68e3a3094..df9fcfab6 100644
--- a/requirements_macos_arm64.txt
+++ b/requirements_macos_arm64.txt
@@ -32,3 +32,4 @@ reportlab
 netCDF4
 igraph
 pyvista
+basemap
diff --git a/requirements_macos_x86_64.txt b/requirements_macos_x86_64.txt
index 205f55047..d1ae32779 100644
--- a/requirements_macos_x86_64.txt
+++ b/requirements_macos_x86_64.txt
@@ -32,3 +32,4 @@ reportlab
 netCDF4
 igraph
 pyvista
+basemap
diff --git a/requirements_win64.txt b/requirements_win64.txt
index f20000c77..06a6615f7 100644
--- a/requirements_win64.txt
+++ b/requirements_win64.txt
@@ -31,3 +31,4 @@ reportlab
 netCDF4
 igraph
 pyvista
+basemap

From c99052327ba56c0117d031772033d84989930b38 Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Mon, 19 Feb 2024 10:57:11 +0100
Subject: [PATCH 04/19] small improvements

---
 openquake/mbi/sub/geojson_from_profiles.py     |  3 ++-
 openquake/sub/cross_sections.py                |  3 +--
 openquake/sub/plotting/plot_cross_section.py   | 15 ++++++++-------
 .../plotting/plot_multiple_cross_sections.py   | 18 +++++++++---------
 4 files changed, 20 insertions(+), 19 deletions(-)

diff --git a/openquake/mbi/sub/geojson_from_profiles.py b/openquake/mbi/sub/geojson_from_profiles.py
index 0626d8b75..22b893529 100755
--- a/openquake/mbi/sub/geojson_from_profiles.py
+++ b/openquake/mbi/sub/geojson_from_profiles.py
@@ -19,7 +19,8 @@ def main(pattern: str, fname_output: str = "profiles.geojson"):
         print(fname)
         dat = np.loadtxt(fname)
         tmp = LineString([(x, y) for x, y in zip(dat[:, 0], dat[:, 1])])
-        features.append(Feature(geometry=tmp))
+        prop = {'csid': fname.split('_')[1].replace('.csv','')}
+        features.append(Feature(geometry=tmp, properties=prop))
     feature_collection = FeatureCollection(features)
     with open(fname_output, 'w') as f:
         dump(feature_collection, f)
diff --git a/openquake/sub/cross_sections.py b/openquake/sub/cross_sections.py
index 448cc8dee..146febb83 100644
--- a/openquake/sub/cross_sections.py
+++ b/openquake/sub/cross_sections.py
@@ -315,7 +315,7 @@ def set_gcmt(self, filename, bffer=75.):
                                                              minlo, maxlo,
                                                              minla, maxla)
         if idxs is not None:
-            cmt_cat.select_catalogue_events(idxs)
+            cmt_cat.select_catalogue_events(idxs[0])
             self.gcmt = cmt_cat
 
     def set_topo(self, filename, bffer=0.25):
@@ -388,7 +388,6 @@ def set_volcano(self, filename, bffer=75.):
             voo[idxv[0]] = 1
             self.volc = numpy.squeeze(vulc[idxv[0], :])
         fin.close()
-        print(self.volc)
 
 
 class Trench:
diff --git a/openquake/sub/plotting/plot_cross_section.py b/openquake/sub/plotting/plot_cross_section.py
index c78f7e03b..390826fa9 100755
--- a/openquake/sub/plotting/plot_cross_section.py
+++ b/openquake/sub/plotting/plot_cross_section.py
@@ -345,7 +345,7 @@ def _plot_volc(axes, csda):
     ola = csda.csec.ola
     patches = []
 
-    
+
 
     if (len(csda.volc)-1) >= 1:
         if str(csda.volc.shape) == '(2,)':
@@ -358,12 +358,13 @@ def _plot_volc(axes, csda):
             patches.append(square)
 
     else:
-        breakpoint()
-        vuls = geodetic_distance(olo, ola,
-                                 csda.volc[0],
-                                 csda.volc[1])
-        square = Rectangle((vuls, -10.0), 7, 12)
-        patches.append(square)
+        print(csda.volc)
+    # breakpoint()
+       # vuls = geodetic_distance(olo, ola,
+       #                          csda.volc[0],
+       #                          csda.volc[1])
+       # square = Rectangle((vuls, -10.0), 7, 12)
+       # patches.append(square)
 
     vv = PatchCollection(patches, zorder=6, color='red', edgecolors='red')
     vv.set_alpha(0.85)
diff --git a/openquake/sub/plotting/plot_multiple_cross_sections.py b/openquake/sub/plotting/plot_multiple_cross_sections.py
index 01e7ca9f1..e33df6cef 100755
--- a/openquake/sub/plotting/plot_multiple_cross_sections.py
+++ b/openquake/sub/plotting/plot_multiple_cross_sections.py
@@ -3,18 +3,16 @@
 import re
 import os
 import sys
-from plot_cross_section import plt_cs
 import matplotlib.pyplot as plt
 from openquake.baselib import sap
+from openquake.sub.plotting.plot_cross_section import plt_cs
 
 
-def pcs(argv):
-    if len(argv) < 2:
-        folder = './'
-    else:
-        folder = argv[1]
+def pcs(cs_fname, out_folder=None):
+    if out_folder == None:
+        out_folder = './'
 
-    fin = open(argv[0], 'r')
+    fin = open(cs_fname, 'r')
     for line in fin:
 
         aa = re.split('\\s+', line.rstrip())
@@ -28,11 +26,13 @@ def pcs(argv):
 
         fig = plt_cs(olo, ola, depth, lnght, strike, ids, ini_fle)
         name = 'section_%s.pdf' % (ids)
-        path = os.path.join(folder, name)
+        path = os.path.join(out_folder, name)
         fig.savefig(path, bbox_inches='tight')
         print('Created %s' % (path))
         plt.close()
 
+pcs.cs_fname = 'file with cross sections details'
+pcs.out_folder = 'place to store the pdfs'
 
 if __name__ == "__main__":
-    main(sys.argv[1:])
+    sap.run(pcs)

From 572659fc3a331c5c5443973e7789d0ecdc5ce878 Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Mon, 19 Feb 2024 10:57:53 +0100
Subject: [PATCH 05/19] adding mbi part

---
 .../mbi/sub/plot_multiple_cross_sections.py   | 23 +++++++++++++++++++
 1 file changed, 23 insertions(+)
 create mode 100755 openquake/mbi/sub/plot_multiple_cross_sections.py

diff --git a/openquake/mbi/sub/plot_multiple_cross_sections.py b/openquake/mbi/sub/plot_multiple_cross_sections.py
new file mode 100755
index 000000000..ce2fecb60
--- /dev/null
+++ b/openquake/mbi/sub/plot_multiple_cross_sections.py
@@ -0,0 +1,23 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+from openquake.baselib import sap
+from openquake.sub.plotting.plot_multiple_cross_sections import pcs
+
+
+def main(cs_file, output_folder=None):
+    """
+    Creates file with parameters needed to plot cross sections.
+    Output file is a list for each cross section with the format:
+
+        lon lat length depth azimuth id <config>.ini
+    """
+
+    pcs(cs_file, output_folder)
+
+
+main.cs_file = 'file with cross sections details'
+main.output_folder = 'place to store the pdfs'
+
+if __name__ == "__main__":
+    sap.run(main)

From bac899521f76be647eb9d784e0cfca2c41032b9f Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MBP-001.lan>
Date: Tue, 20 Feb 2024 17:20:21 +0100
Subject: [PATCH 06/19] adding script to convert ruptures from srcmod to
 geojson

---
 openquake/mbi/sub/srcmod_to_json.py | 31 +++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)
 create mode 100644 openquake/mbi/sub/srcmod_to_json.py

diff --git a/openquake/mbi/sub/srcmod_to_json.py b/openquake/mbi/sub/srcmod_to_json.py
new file mode 100644
index 000000000..c076c874e
--- /dev/null
+++ b/openquake/mbi/sub/srcmod_to_json.py
@@ -0,0 +1,31 @@
+import pandas as pd
+import geopandas as gpd
+from shapely.geometry import Point
+
+from openquake.baselib import sap
+
+
+def main(datapath, eventid):
+    fin = f'{datapath}/seismicity/srcmod/{eventid}.fsp'
+    root = fin.split('/')[-1].replace('.fsp','')
+    outfi = f'{root}.geojson'
+
+    file1 = open(fin, 'r')
+    Lines = file1.readlines()
+
+    lons, lats, depths, slips = [],[],[],[]
+    for line in Lines:
+        if (line[0] != '%') & (line != ' \n'):
+            parts = line.strip().split()
+            lons.append(parts[1]); lats.append(parts[0])
+            depths.append(float(parts[4])); slips.append(float(parts[5]))
+    df = pd.DataFrame({'lon': lons, 'lat': lats, 'depth': depths, 'slip': slips})
+
+    gdf = gpd.GeoDataFrame(df, geometry=[Point(xy) for xy in zip(df['lon'], df['lat'])], crs="epsg:4326")
+    gdf.to_file(outfi, driver="GeoJSON")
+
+main.datapath = 'path to gem_hazard_data'
+main.eventid = 'name of event to convert'
+
+if __name__ == "__main__":
+    sap.run(main)

From 2988b93b8bc6fd0bd735cf918d4c268a0351b326 Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Wed, 21 Feb 2024 11:29:46 +0100
Subject: [PATCH 07/19] adding option to filter gcmt by magnitude

---
 openquake/sub/cross_sections.py                        | 5 ++++-
 openquake/sub/plotting/plot_cross_section.py           | 8 +++++++-
 openquake/sub/plotting/plot_multiple_cross_sections.py | 2 ++
 3 files changed, 13 insertions(+), 2 deletions(-)

diff --git a/openquake/sub/cross_sections.py b/openquake/sub/cross_sections.py
index 146febb83..75f2800e4 100644
--- a/openquake/sub/cross_sections.py
+++ b/openquake/sub/cross_sections.py
@@ -290,7 +290,7 @@ def set_litho_moho_depth(self, filename, bffer=100.):
             boo[idxl[0]] = 1
             self.litho = numpy.squeeze(dataa[idxl[0], :])
 
-    def set_gcmt(self, filename, bffer=75.):
+    def set_gcmt(self, filename, gcmt_mag=0.0, bffer=75.):
         """
         :parameter cmt_cat:
             Name of a file in the .ndk format
@@ -298,6 +298,9 @@ def set_gcmt(self, filename, bffer=75.):
         print('setting gcmt')
         parser = ParseNDKtoGCMT(filename)
         cmt_cat = parser.read_file()
+        # prune to magnitude range
+        mags = cmt_cat.data['magnitude']
+        cmt_cat.select_catalogue_events(mags > gcmt_mag)
         loc = cmt_cat.data['longitude']
         lac = cmt_cat.data['latitude']
 
diff --git a/openquake/sub/plotting/plot_cross_section.py b/openquake/sub/plotting/plot_cross_section.py
index 390826fa9..97d2a12ed 100755
--- a/openquake/sub/plotting/plot_cross_section.py
+++ b/openquake/sub/plotting/plot_cross_section.py
@@ -586,7 +586,13 @@ def plt_cs(olo, ola, depp, lnght, strike, ids, ini_filename):
     if re.search('[a-z]', fname_slab):
         csda.set_slab1pt0(fname_slab)
     csda.set_crust1pt0_moho_depth(fname_crust)
-    csda.set_gcmt(fname_gcmt)
+
+    try:
+        gcmt_mag = float(config['params']['gcmt_mag_lower_limit'])
+    except:
+        gcmt_mag = 0.0
+
+    csda.set_gcmt(fname_gcmt, gcmt_mag=gcmt_mag)
     csda.set_topo(fname_topo)
     csda.set_litho_moho_depth(fname_litho)
     csda.set_volcano(fname_volc)
diff --git a/openquake/sub/plotting/plot_multiple_cross_sections.py b/openquake/sub/plotting/plot_multiple_cross_sections.py
index e33df6cef..13206100a 100755
--- a/openquake/sub/plotting/plot_multiple_cross_sections.py
+++ b/openquake/sub/plotting/plot_multiple_cross_sections.py
@@ -11,6 +11,8 @@
 def pcs(cs_fname, out_folder=None):
     if out_folder == None:
         out_folder = './'
+    if not os.path.exists(out_folder):
+        os.makedirs(out_folder)
 
     fin = open(cs_fname, 'r')
     for line in fin:

From 2607d83d13934293483131593d6cddb993cc754e Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Fri, 23 Feb 2024 09:57:10 +0100
Subject: [PATCH 08/19] handling directories as to not refill existing
 profile/edge directory

---
 openquake/mbi/sub/build_complex_surface.py | 12 ++++++++++++
 openquake/sub/build_complex_surface.py     |  2 ++
 2 files changed, 14 insertions(+)

diff --git a/openquake/mbi/sub/build_complex_surface.py b/openquake/mbi/sub/build_complex_surface.py
index f9822dbf9..0ed8c2382 100755
--- a/openquake/mbi/sub/build_complex_surface.py
+++ b/openquake/mbi/sub/build_complex_surface.py
@@ -1,6 +1,8 @@
 #!/usr/bin/env python
 # coding: utf-8
 
+import os
+
 from openquake.baselib import sap
 from openquake.sub.build_complex_surface import build_complex_surface
 
@@ -11,6 +13,16 @@ def main(in_path, max_sampl_dist, out_path, upper_depth=0,
     Builds edges that can be used to generate a complex fault surface
     starting from a set of profiles
     """
+    # check if output folder is empty
+    if os.path.exists(out_path):
+        tmps = f'\nError: {out_path} already exists and contains profiles and '
+        tmps += '\n edges from a previous run! Specify empty output directory.'
+        print(tmps)
+        sys.exit()
+    # otherwise create directory
+    else:
+        os.makedirs(out_path)
+
     build_complex_surface(in_path, max_sampl_dist, out_path, upper_depth,
                           lower_depth, from_id, to_id)
 
diff --git a/openquake/sub/build_complex_surface.py b/openquake/sub/build_complex_surface.py
index b21f386a7..8732f3ddb 100644
--- a/openquake/sub/build_complex_surface.py
+++ b/openquake/sub/build_complex_surface.py
@@ -71,6 +71,8 @@ def build_complex_surface(in_path, max_sampl_dist, out_path, upper_depth=0,
         tmps += '    output: {0:s}\n'.format(out_path)
         sys.exit()
 
+
+
     # Read the profiles
     sps, dmin, dmax = read_profiles_csv(in_path,
                                         float(upper_depth),

From 0d7c2638a7f40e0ff07ed22c41437afc0ad496be Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Thu, 29 Feb 2024 15:05:51 +0100
Subject: [PATCH 09/19] make print statements more meaningful

---
 openquake/mbi/sub/geojson_from_profiles.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/openquake/mbi/sub/geojson_from_profiles.py b/openquake/mbi/sub/geojson_from_profiles.py
index 22b893529..58943119d 100755
--- a/openquake/mbi/sub/geojson_from_profiles.py
+++ b/openquake/mbi/sub/geojson_from_profiles.py
@@ -14,9 +14,7 @@ def main(pattern: str, fname_output: str = "profiles.geojson"):
     """
 
     features = []
-    print(pattern)
     for fname in glob.glob(pattern):
-        print(fname)
         dat = np.loadtxt(fname)
         tmp = LineString([(x, y) for x, y in zip(dat[:, 0], dat[:, 1])])
         prop = {'csid': fname.split('_')[1].replace('.csv','')}
@@ -24,6 +22,7 @@ def main(pattern: str, fname_output: str = "profiles.geojson"):
     feature_collection = FeatureCollection(features)
     with open(fname_output, 'w') as f:
         dump(feature_collection, f)
+    print(f'profiles written to {fname_output}')
 
 
 descr = 'The pattern to the text files with the profiles'

From fa193189a0e5fb46674c0fb789ff189b67804178 Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Thu, 29 Feb 2024 15:08:59 +0100
Subject: [PATCH 10/19] adding line to explain use

---
 openquake/mbi/sub/geojson_from_profiles.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/openquake/mbi/sub/geojson_from_profiles.py b/openquake/mbi/sub/geojson_from_profiles.py
index 58943119d..94b3fc6b9 100755
--- a/openquake/mbi/sub/geojson_from_profiles.py
+++ b/openquake/mbi/sub/geojson_from_profiles.py
@@ -11,6 +11,9 @@ def main(pattern: str, fname_output: str = "profiles.geojson"):
     """
     Creates a geojson file with all the sections included in the text files
     matching the pattern
+
+    Example use: 
+    oqm sub geojson_from_profiles 'openquake/sub/tests/data/cs_cam/cs*csv'
     """
 
     features = []

From a4999f5f782e3b42e2f03140985d617adc6d79a5 Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Thu, 29 Feb 2024 15:13:02 +0100
Subject: [PATCH 11/19] add exmaple for use

---
 openquake/mbi/sub/make_cs_coords.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/openquake/mbi/sub/make_cs_coords.py b/openquake/mbi/sub/make_cs_coords.py
index d29ac211c..75ca55edf 100755
--- a/openquake/mbi/sub/make_cs_coords.py
+++ b/openquake/mbi/sub/make_cs_coords.py
@@ -10,6 +10,11 @@ def main(cs_dir, outfi, ini_fname, cs_length=300., cs_depth=300.):
     Output file is a list for each cross section with the format:
 
         lon lat length depth azimuth id <config>.ini
+
+    Example use: 
+
+    oqm sub make_cs_coords openquake/sub/tests/data/cs_cam cs_profiles.csv 
+                example.ini 300 300
     """
 
     make_cs_coords(cs_dir, outfi, ini_fname, cs_length, cs_depth)

From 7f19851e34d9f0a3cae5b37d618f73b9815070bf Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Thu, 29 Feb 2024 15:16:59 +0100
Subject: [PATCH 12/19] adding example for making map

---
 openquake/mbi/sub/plot_cross_sections_map.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/openquake/mbi/sub/plot_cross_sections_map.py b/openquake/mbi/sub/plot_cross_sections_map.py
index e03d72cef..a16e0778d 100755
--- a/openquake/mbi/sub/plot_cross_sections_map.py
+++ b/openquake/mbi/sub/plot_cross_sections_map.py
@@ -7,10 +7,12 @@
 
 def main(config_fname, cs_file=None):
     """
-    Creates file with parameters needed to plot cross sections.
-    Output file is a list for each cross section with the format:
+    Plots map of cross sections with earthquake data
 
-        lon lat length depth azimuth id <config>.ini
+    Example:
+        oqm sub plot_cross_sections_map config.ini cs_profiles.cs
+
+    Note: paths in config.ini are relative to cwd
     """
 
     plot(config_fname, cs_file)

From b4123c6af440e4207aa7f3973fa202acdad2df26 Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Thu, 29 Feb 2024 15:21:37 +0100
Subject: [PATCH 13/19] adding correct description to
 plot_multiple_cross_sections.py

---
 openquake/mbi/sub/plot_multiple_cross_sections.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/openquake/mbi/sub/plot_multiple_cross_sections.py b/openquake/mbi/sub/plot_multiple_cross_sections.py
index ce2fecb60..8bccd5452 100755
--- a/openquake/mbi/sub/plot_multiple_cross_sections.py
+++ b/openquake/mbi/sub/plot_multiple_cross_sections.py
@@ -7,10 +7,13 @@
 
 def main(cs_file, output_folder=None):
     """
-    Creates file with parameters needed to plot cross sections.
-    Output file is a list for each cross section with the format:
+    Plots cross section of data for each cross section in
+    cs_file including the data in the ini file referenced
+    in the cs_file lines. Saves pdfs to output_folder
 
-        lon lat length depth azimuth id <config>.ini
+    Example:
+
+        oqm sub plot_multiple_cross_sections cs_profiles.cs pdf 
     """
 
     pcs(cs_file, output_folder)

From 1b14be6cafdeacaf54daf995daf6e1c92b7ec27c Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Thu, 29 Feb 2024 15:29:16 +0100
Subject: [PATCH 14/19] generalize srcmod converter; add example running syntax

---
 openquake/mbi/sub/srcmod_to_json.py | 25 +++++++++++++++++++------
 1 file changed, 19 insertions(+), 6 deletions(-)

diff --git a/openquake/mbi/sub/srcmod_to_json.py b/openquake/mbi/sub/srcmod_to_json.py
index c076c874e..63679b70e 100644
--- a/openquake/mbi/sub/srcmod_to_json.py
+++ b/openquake/mbi/sub/srcmod_to_json.py
@@ -5,8 +5,19 @@
 from openquake.baselib import sap
 
 
-def main(datapath, eventid):
-    fin = f'{datapath}/seismicity/srcmod/{eventid}.fsp'
+def main(fippath):#, eventid):
+    """
+    Converts an event from the sourcemod database of the format
+    .fsp to a geojson file
+
+    Note: first download events from http://equake-rc.info/srcmod/ 
+
+    Example:
+
+        oqm sub srcmod_to_json srcmod_events/s2013SCOTIA01HAYE.fsp
+
+    """
+    fin = f'{fippath}'
     root = fin.split('/')[-1].replace('.fsp','')
     outfi = f'{root}.geojson'
 
@@ -19,13 +30,15 @@ def main(datapath, eventid):
             parts = line.strip().split()
             lons.append(parts[1]); lats.append(parts[0])
             depths.append(float(parts[4])); slips.append(float(parts[5]))
-    df = pd.DataFrame({'lon': lons, 'lat': lats, 'depth': depths, 'slip': slips})
+    df = pd.DataFrame({'lon': lons, 'lat': lats,
+                       'depth': depths, 'slip': slips})
 
-    gdf = gpd.GeoDataFrame(df, geometry=[Point(xy) for xy in zip(df['lon'], df['lat'])], crs="epsg:4326")
+    gdf = gpd.GeoDataFrame(df, geometry=[Point(xy) for xy in \
+            zip(df['lon'], df['lat'])], crs="epsg:4326")
     gdf.to_file(outfi, driver="GeoJSON")
+    print(f'Written to {outfi}')
 
-main.datapath = 'path to gem_hazard_data'
-main.eventid = 'name of event to convert'
+main.fippath = 'path to .fsp file'
 
 if __name__ == "__main__":
     sap.run(main)

From 1b2c97ed767534697921c1606b08730b19553727 Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Thu, 29 Feb 2024 15:32:26 +0100
Subject: [PATCH 15/19] updates to use sap, generalises so can either ask to
 create the profiels or provide them

---
 openquake/sub/create_multiple_cross_sections.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/openquake/sub/create_multiple_cross_sections.py b/openquake/sub/create_multiple_cross_sections.py
index 11f700743..a9120c499 100644
--- a/openquake/sub/create_multiple_cross_sections.py
+++ b/openquake/sub/create_multiple_cross_sections.py
@@ -93,7 +93,6 @@ def plot(trench, cat, cs_dict, interdistance):
     #
     # Draw paralleles and meridians with labels
     # labels = [left,right,top,bottom]
-#    m.drawcoastlines()
     m.drawmeridians(numpy.arange(numpy.floor(minlo/10.)*10,
                                  numpy.ceil(maxlo/10.)*10, 5.),
                     labels=[False, False, False, True])
@@ -129,7 +128,6 @@ def plot(trench, cat, cs_dict, interdistance):
     else:
         for key in cs_dict:
             cs = cs_dict[key]
-        #if cs is not None:
             x, y = m(cs.plo, cs.pla)
             plt.plot(x, y, ':r', linewidth=2, zorder=20)
             text = plt.text(x[0], y[0], '%s' % key, ha='center',

From 5606242431d3088a67d6274afc972fa5a0eac282 Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Thu, 29 Feb 2024 15:34:43 +0100
Subject: [PATCH 16/19] include use note

---
 openquake/sub/make_cs_coords.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/openquake/sub/make_cs_coords.py b/openquake/sub/make_cs_coords.py
index f6c8c052f..ee1f048d4 100755
--- a/openquake/sub/make_cs_coords.py
+++ b/openquake/sub/make_cs_coords.py
@@ -33,6 +33,12 @@
 from openquake.hazardlib.geo.geodetic import azimuth
 
 def make_cs_coords(cs_dir, outfi, ini_fname, cs_length=300., cs_depth=300.):
+    """
+    Creates cs_coords file of format to be used by plotting script, in the
+    case that the profiles have been generated in any other way than by
+    the create_multiple_cross_sections.py
+    """
+
     cs_files = sorted(glob.glob(f'{cs_dir}/*csv'))
     lines = []
     for fi in cs_files:

From 54e61a984a4a5c747a78c646c3176d7c25ab1617 Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Thu, 29 Feb 2024 15:38:53 +0100
Subject: [PATCH 17/19] cleaning

---
 openquake/sub/plotting/plot_cross_section.py | 13 +++++--------
 1 file changed, 5 insertions(+), 8 deletions(-)

diff --git a/openquake/sub/plotting/plot_cross_section.py b/openquake/sub/plotting/plot_cross_section.py
index 97d2a12ed..634663ee1 100755
--- a/openquake/sub/plotting/plot_cross_section.py
+++ b/openquake/sub/plotting/plot_cross_section.py
@@ -359,12 +359,6 @@ def _plot_volc(axes, csda):
 
     else:
         print(csda.volc)
-    # breakpoint()
-       # vuls = geodetic_distance(olo, ola,
-       #                          csda.volc[0],
-       #                          csda.volc[1])
-       # square = Rectangle((vuls, -10.0), 7, 12)
-       # patches.append(square)
 
     vv = PatchCollection(patches, zorder=6, color='red', edgecolors='red')
     vv.set_alpha(0.85)
@@ -570,7 +564,11 @@ def plt_cs(olo, ola, depp, lnght, strike, ids, ini_filename):
 
     config = configparser.ConfigParser()
     config.read(ini_filename)
-#    fname_trench = config['data']['trench_axis_filename']
+    try:
+        fname_trench = config['data']['trench_axis_filename']
+        csda.set_trench_axis(fname_trench)
+    except: 
+        fname_trench = None
     fname_eqk_cat = config['data']['catalogue_pickle_filename']
     fname_slab = config['data']['slab1pt0_filename']
     fname_crust = config['data']['crust1pt0_filename']
@@ -580,7 +578,6 @@ def plt_cs(olo, ola, depp, lnght, strike, ids, ini_filename):
     fname_litho = config['data']['litho_filename']
     fname_volc = config['data']['volc_filename']
 
-#    csda.set_trench_axis(fname_trench)
     cat = pickle.load(open(fname_eqk_cat, 'rb'))
     csda.set_catalogue(cat)
     if re.search('[a-z]', fname_slab):

From 695b4aaa0da044774e2bb9598566134587dd9aa2 Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Thu, 29 Feb 2024 15:40:34 +0100
Subject: [PATCH 18/19] add use case note

---
 openquake/sub/plotting/plot_multiple_cross_sections.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/openquake/sub/plotting/plot_multiple_cross_sections.py b/openquake/sub/plotting/plot_multiple_cross_sections.py
index 13206100a..20a048c0d 100755
--- a/openquake/sub/plotting/plot_multiple_cross_sections.py
+++ b/openquake/sub/plotting/plot_multiple_cross_sections.py
@@ -9,6 +9,12 @@
 
 
 def pcs(cs_fname, out_folder=None):
+    """
+    plots cross sections based on lines in cs_fname of format
+
+        lon lat length depth azimuth id <config>.ini
+
+    """
     if out_folder == None:
         out_folder = './'
     if not os.path.exists(out_folder):

From c69f74f582e2f226b950d910fa43e3f08c90880e Mon Sep 17 00:00:00 2001
From: Kendra Johnson <kjohnson@Kendras-MacBook-Pro.local>
Date: Thu, 29 Feb 2024 15:41:43 +0100
Subject: [PATCH 19/19] making general to plot even if cross sections exist

---
 openquake/sub/plotting/plot_multiple_cross_sections_map.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/openquake/sub/plotting/plot_multiple_cross_sections_map.py b/openquake/sub/plotting/plot_multiple_cross_sections_map.py
index 35589295d..e3e1e4232 100755
--- a/openquake/sub/plotting/plot_multiple_cross_sections_map.py
+++ b/openquake/sub/plotting/plot_multiple_cross_sections_map.py
@@ -79,5 +79,4 @@ def plot(config_file, cs_file=None):
 plot.cs_file = 'existing cross sections details'
 
 if __name__ == "__main__":
-    #main(sys.argv[1:])
     sap.run(plot)