Merge pull request #396 from GEMScienceTools/update_sub_plotting

Update subduction plotting tools

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)
 

openquake/mbi/sub/geojson_from_profiles.py

@@ -11,18 +11,21 @@ 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 = []
-    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])])
-        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)
+    print(f'profiles written to {fname_output}')
 
 
 descr = 'The pattern to the text files with the profiles'

openquake/mbi/sub/make_cs_coords.py

@@ -0,0 +1,30 @@
+#!/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
+
+    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)
+
+
+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)

openquake/mbi/sub/plot_cross_sections_map.py

@@ -0,0 +1,25 @@
+#!/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):
+    """
+    Plots map of cross sections with earthquake data
+
+    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)
+
+
+main.config_fname = 'config file to datasets'
+main.cs_file = 'existing cross sections details'
+
+if __name__ == "__main__":
+    sap.run(main)

openquake/mbi/sub/plot_multiple_cross_sections.py

@@ -0,0 +1,26 @@
+#!/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):
+    """
+    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
+
+    Example:
+
+        oqm sub plot_multiple_cross_sections cs_profiles.cs pdf 
+    """
+
+    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)

openquake/mbi/sub/srcmod_to_json.py

@@ -0,0 +1,44 @@
+import pandas as pd
+import geopandas as gpd
+from shapely.geometry import Point
+
+from openquake.baselib import sap
+
+
+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'
+
+    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")
+    print(f'Written to {outfi}')
+
+main.fippath = 'path to .fsp file'
+
+if __name__ == "__main__":
+    sap.run(main)

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),

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,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])
@@ -121,12 +118,16 @@ 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]
             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 +136,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 +181,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)

openquake/sub/cross_sections.py

@@ -244,13 +244,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,
@@ -315,16 +315,19 @@ 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.):
+    def set_gcmt(self, filename, gcmt_mag=0.0, bffer=75.):
         """
         :parameter cmt_cat:
             Name of a file in the .ndk format
         """
         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']
 
@@ -342,7 +345,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):
@@ -378,7 +381,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.):
         """
@@ -413,9 +416,8 @@ 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)
 
 
 class Trench: