Turn arias duration utility into a module, added rod velocity module

bbp/comps/arias_duration.py

@@ -21,9 +21,17 @@
 from __future__ import division, print_function
 
 # G2CMSS = 980.665 # Convert g to cm/s/s
+import os
+import sys
 import math
 import scipy.integrate
 
+# Import Broadband modules
+import bband_utils
+import install_cfg
+import bbp_formatter
+from station_list import StationList
+
 # Converting to cm units. Use approximation to g
 G_TO_CMS = 981.0 # %(cm/s)
 
@@ -146,19 +154,31 @@ def ad_from_acc(a_in_peer_file, a_out_ad):
     else:
         ia_norm = [0.0 for i_acc in arias_intensity]
 
-    # Define the time for AI=5%, 75%, 95%
+    # Define the time for AI=5%, 20%, 75%, 80%, 95%
     time_ai5 = 0
     for i in range(pts):
         if ia_norm[i] >= 5:
             break
         time_ai5 = dt * i
 
+    time_ai20 = 0
+    for i in range(pts):
+        if ia_norm[i] >= 20:
+            break
+        time_ai20 = dt * i
+
     time_ai75 = 0
     for i in range(pts):
         if ia_norm[i] >= 75:
             break
         time_ai75 = dt * i
 
+    time_ai80 = 0
+    for i in range(pts):
+        if ia_norm[i] >= 80:
+            break
+        time_ai80 = dt * i
+
     time_ai95 = 0
     for i in range(pts):
         if ia_norm[i] >= 95:
@@ -168,6 +188,7 @@ def ad_from_acc(a_in_peer_file, a_out_ad):
     # Now, calculate the arias intervals 5% to 75% and 5% to 95%
     time_5_75 = time_ai75 - time_ai5
     time_5_95 = time_ai95 - time_ai5
+    time_20_80 = time_ai80 - time_ai20
 
     #print("Arias Intervals: 5 to 75 % 6f (secs), 5 to 95 % 6f (secs) " %
     #      (time_5_75, time_5_95))
@@ -197,8 +218,8 @@ def ad_from_acc(a_in_peer_file, a_out_ad):
     outfile.write("# Arias Intensities from input accel: %s\n" % a_in_peer_file)
     outfile.write("# Peak Arias Intensity (cm/sec): %f (secs): %f\n" %
                   (arias_intensity_max, (arias_intensity_index * dt)))
-    outfile.write("# Arias Intervals: T5-75 %f (s), T5-95 %f (secs)\n" %
-                  (time_5_75, time_5_95))
+    outfile.write("# Arias Intervals: T5-75 %f (s), T5-95 %f (s), T20-80 %f (s)\n" %
+                  (time_5_75, time_5_95, time_20_80))
     outfile.write("# Seconds Accel (g) "
                   "Arias Intensity (cm/s) "
                   "ADNormalized (%)\n")
@@ -207,4 +228,93 @@ def ad_from_acc(a_in_peer_file, a_out_ad):
         outfile.write("% 8f % 8f % 8f % 8f\n" %
                       (dt_vals[i], acc[i], arias_intensity[i], ia_norm[i]))
     outfile.close()
-    return 0
+    return arias_intensity_max, time_5_75, time_5_95, time_20_80;
+
+class AriasDuration(object):
+    """
+    BBP module implementation of arias duration
+    """
+
+    def __init__(self, i_r_stations, sim_id=0):
+        """
+        Initializes class variables
+        """
+        self.sim_id = sim_id
+        self.r_stations = i_r_stations
+
+    def run(self):
+        print("AriasDuration".center(80, '-'))
+        #
+        # convert input bbp acc files to peer format acc files
+        #
+
+        install = install_cfg.InstallCfg.getInstance()
+        sim_id = self.sim_id
+        sta_base = os.path.basename(os.path.splitext(self.r_stations)[0])
+        self.log = os.path.join(install.A_OUT_LOG_DIR, str(sim_id),
+                                "%d.araisduration_%s.log" % (sim_id, sta_base))
+        a_statfile = os.path.join(install.A_IN_DATA_DIR,
+                                  str(sim_id),
+                                  self.r_stations)
+        a_tmpdir = os.path.join(install.A_TMP_DATA_DIR, str(sim_id))
+        a_outdir = os.path.join(install.A_OUT_DATA_DIR, str(sim_id))
+
+        #
+        # Make sure the tmp and out directories exist
+        #
+        bband_utils.mkdirs([a_tmpdir, a_outdir], print_cmd=False)
+
+        slo = StationList(a_statfile)
+        site_list = slo.getStationList()
+
+        for site in site_list:
+            stat = site.scode
+            print("==> Processing station: %s" % (stat))
+            bbpfile = os.path.join(a_outdir,
+                                   "%d.%s.acc.bbp" % (sim_id, stat))
+
+            # Now we need to convert to peer format
+            out_n_acc = os.path.join(a_tmpdir,
+                                     "%d.%s.peer_n.acc" % (sim_id, stat))
+            out_e_acc = os.path.join(a_tmpdir,
+                                     "%d.%s.peer_e.acc" % (sim_id, stat))
+            out_z_acc = os.path.join(a_tmpdir,
+                                     "%d.%s.peer_z.acc" % (sim_id, stat))
+            bbp_formatter.bbp2peer(bbpfile, out_n_acc, out_e_acc, out_z_acc, accel=True)
+
+            # Duration output files
+            out_n_arias = os.path.join(a_outdir,
+                                     "%d.%s_N.arias" % (sim_id, stat))
+            out_e_arias = os.path.join(a_outdir,
+                                     "%d.%s_E.arias" % (sim_id, stat))
+            out_z_arias = os.path.join(a_outdir,
+                                     "%d.%s_Z.arias" % (sim_id, stat))
+
+            # compute each one
+            n_max, n_time_5_75, n_time_5_95, n_time_20_80 = ad_from_acc(out_n_acc, out_n_arias)
+            e_max, e_time_5_75, e_time_5_95, e_time_20_80 = ad_from_acc(out_e_acc, out_e_arias)
+            z_max, z_time_5_75, z_time_5_95, z_time_20_80 = ad_from_acc(out_z_acc, out_z_arias)
+            geo_max = math.sqrt(n_max*e_max)
+            geo_time_5_75 = math.sqrt(n_time_5_75*e_time_5_75)
+            geo_time_5_95 = math.sqrt(n_time_5_95*e_time_5_95)
+            geo_time_20_80 = math.sqrt(n_time_20_80*e_time_20_80)
+
+            # write summary file
+            out_duration = os.path.join(a_outdir,
+                                     "%d.%s.ard" % (sim_id, stat))
+            outfile = open(out_duration, "w")
+            outfile.write("# Component  Peak Arias (cm/s)  T5-75 (s)  T5-95 (s)  T20-80 (s)\n")
+            outfile.write("N %f %f %f %f\n" % (n_max, n_time_5_75, n_time_5_95, n_time_20_80))
+            outfile.write("E %f %f %f %f\n" % (e_max, e_time_5_75, e_time_5_95, e_time_20_80))
+            outfile.write("Z %f %f %f %f\n" % (z_max, z_time_5_75, z_time_5_95, z_time_20_80))
+            outfile.write("GEOM %f %f %f %f\n" % (geo_max, geo_time_5_75, geo_time_5_95, geo_time_20_80))
+            outfile.close()
+
+        # All done!
+        print("AriasDuration Completed".center(80, '-'))
+
+if __name__ == '__main__':
+    print("Testing Module: %s" % (os.path.basename(sys.argv[0])))
+    ME = AriasDuration(sys.argv[1], sim_id=int(sys.argv[2]))
+    ME.run()
+    sys.exit(0)

bbp/comps/bbp_formatter.py

@@ -173,7 +173,7 @@ def peer2bbp(in_peer_n_file, in_peer_e_file, in_peer_z_file, out_bbp_file):
     # Lastly, close the file
     bbp_file.close()
 
-def bbp2peer(in_bbp_file, out_peer_n_file, out_peer_e_file, out_peer_z_file):
+def bbp2peer(in_bbp_file, out_peer_n_file, out_peer_e_file, out_peer_z_file, accel=True):
     """
     Convert bbp file into three peer files for use by RotD50 and
     other programs that input PEER format seismograms
@@ -233,9 +233,14 @@ def bbp2peer(in_bbp_file, out_peer_n_file, out_peer_e_file, out_peer_z_file):
                                                   "Error in conversion.")
             else:
                 dt_vals.append(dt)
-                n_vals.append(float(elems[1]) / bband_utils.G2CMSS)
-                e_vals.append(float(elems[2]) / bband_utils.G2CMSS)
-                z_vals.append(float(elems[3]) / bband_utils.G2CMSS)
+                if accel:
+                    n_vals.append(float(elems[1]) / bband_utils.G2CMSS)
+                    e_vals.append(float(elems[2]) / bband_utils.G2CMSS)
+                    z_vals.append(float(elems[3]) / bband_utils.G2CMSS)
+                else:
+                    n_vals.append(float(elems[1]))
+                    e_vals.append(float(elems[2]))
+                    z_vals.append(float(elems[3]))
 
     # Prepare to write 6 colume format
     n_file = open(out_peer_n_file, "w")
@@ -255,11 +260,16 @@ def bbp2peer(in_bbp_file, out_peer_n_file, out_peer_e_file, out_peer_z_file):
         e_file.write(line)
         z_file.write(line)
 
-    n_file.write("Acceleration in g\n")
+    if accel:
+        n_file.write("Acceleration in g\n")
+        e_file.write("Acceleration in g\n")
+        z_file.write("Acceleration in g\n")
+    else:
+        n_file.write("Velicity in cm/sec\n")
+        e_file.write("Velicity in cm/sec\n")
+        z_file.write("Velicity in cm/sec\n")
     n_file.write("  %d   %1.6f   NPTS, DT\n" % (npts, dt))
-    e_file.write("Acceleration in g\n")
     e_file.write("  %d   %1.6f   NPTS, DT\n" % (npts, dt))
-    z_file.write("Acceleration in g\n")
     z_file.write("  %d   %1.6f   NPTS, DT\n" % (npts, dt))
 
     cur_line = 0

bbp/comps/module.py

@@ -33,6 +33,7 @@
 from uc_site import UCSite
 from wcc_siteamp import WccSiteamp
 from rotd50 import RotD50
+from rotd_vel import RotDVel
 from fas import FAS
 from obs_seismograms import ObsSeismograms
 from copy_seismograms import CopySeismograms
@@ -59,6 +60,7 @@
 from irikura_gen_srf import IrikuraGenSrf
 from irikura_hf import IrikuraHF
 from seismo_soil import SeismoSoil
+from arias_duration import AriasDuration
 
 class Module(object):
     def __init__(self):

bbp/comps/rotd_vel.py

@@ -0,0 +1,164 @@
+#!/usr/bin/env python
+"""
+Copyright 2010-2017 University Of Southern California
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+ http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+"""
+from __future__ import division, print_function
+
+# Import Python modules
+import os
+import sys
+
+# Import Broadband modules
+import bband_utils
+import install_cfg
+import bbp_formatter
+from station_list import StationList
+
+class RotDVel(object):
+    """
+    BBP module implementation of rotd50 for velocity provided by UCB.
+    Rotd50 inputs seismograms and outputs response spectra
+    """
+    @staticmethod
+    def do_rotd_vel(workdir, peer_input_e_file, peer_input_n_file,
+                  peer_input_z_file, output_rotd_file,
+                  logfile):
+        """
+        This function runs the rotd50 command inside workdir, using
+        the inputs and outputs specified
+        """
+        install = install_cfg.InstallCfg.getInstance()
+
+        # Make sure we don't have absolute path names
+        peer_input_e_file = os.path.basename(peer_input_e_file)
+        peer_input_n_file = os.path.basename(peer_input_n_file)
+        peer_input_z_file = os.path.basename(peer_input_z_file)
+        output_rotd_file = os.path.basename(output_rotd_file)
+
+        # Save cwd, change back to it at the end
+        old_cwd = os.getcwd()
+        os.chdir(workdir)
+
+        # Make sure we remove the output files first or Fortran will
+        # complain if they already exist
+        try:
+            os.unlink(output_rotd_file)
+        except OSError:
+            pass
+
+        #
+        # write config file for rotd100 program
+        rd50_conf = open("rotd100_inp.cfg", 'w')
+        # This flag indicates inputs acceleration
+        rd50_conf.write("2 interp flag\n")
+        # This flag indicate we are processing two input files (horizontals)
+        rd50_conf.write("1 Npairs\n")
+        # Number of headers in the file
+        rd50_conf.write("6 Nhead\n")
+        rd50_conf.write("%s\n" % peer_input_e_file)
+        rd50_conf.write("%s\n" % peer_input_n_file)
+        rd50_conf.write("%s\n" % output_rotd_file)
+        # Close file
+        rd50_conf.close()
+
+        progstring = ("%s/rotd100 >> %s 2>&1" % (install.A_UCB_BIN_DIR, logfile))
+        bband_utils.runprog(progstring, abort_on_error=True, print_cmd=False)
+
+        # Restore working directory
+        os.chdir(old_cwd)
+
+    def __init__(self, i_r_stations, sim_id=0):
+        """
+        Initializes class variables
+        """
+        self.sim_id = sim_id
+        self.r_stations = i_r_stations
+
+    def run(self):
+        print("RotDVel".center(80, '-'))
+        #
+        # convert input bbp acc files to peer format acc files
+        #
+
+        install = install_cfg.InstallCfg.getInstance()
+        sim_id = self.sim_id
+        sta_base = os.path.basename(os.path.splitext(self.r_stations)[0])
+        self.log = os.path.join(install.A_OUT_LOG_DIR, str(sim_id),
+                                "%d.rotdvel_%s.log" % (sim_id, sta_base))
+        a_statfile = os.path.join(install.A_IN_DATA_DIR,
+                                  str(sim_id),
+                                  self.r_stations)
+        a_tmpdir = os.path.join(install.A_TMP_DATA_DIR, str(sim_id))
+        a_outdir = os.path.join(install.A_OUT_DATA_DIR, str(sim_id))
+
+        #
+        # Make sure the tmp and out directories exist
+        #
+        bband_utils.mkdirs([a_tmpdir, a_outdir], print_cmd=False)
+
+        slo = StationList(a_statfile)
+        site_list = slo.getStationList()
+
+        for site in site_list:
+            stat = site.scode
+            print("==> Processing station: %s" % (stat))
+
+            # Create path names and check if their sizes are within bounds
+            nsfile = os.path.join(a_tmpdir,
+                                  "%d.%s.000" % (sim_id, stat))
+            ewfile = os.path.join(a_tmpdir,
+                                  "%d.%s.090" % (sim_id, stat))
+            udfile = os.path.join(a_tmpdir,
+                                  "%d.%s.ver" % (sim_id, stat))
+            bbpfile = os.path.join(a_outdir,
+                                   "%d.%s.vel.bbp" % (sim_id, stat))
+
+            bband_utils.check_path_lengths([nsfile, ewfile, udfile],
+                                           bband_utils.GP_MAX_FILENAME)
+
+            cmd = ("%s/wcc2bbp " % (install.A_GP_BIN_DIR) +
+                   "wcc2bbp=0 nsfile=%s ewfile=%s udfile=%s < %s >> %s 2>&1" %
+                   (nsfile, ewfile, udfile, bbpfile, self.log))
+            bband_utils.runprog(cmd, abort_on_error=True, print_cmd=False)
+
+            # Now we need to convert to peer format
+            out_n_vel = os.path.join(a_tmpdir,
+                                     "%d.%s.peer_n.vel" % (sim_id, stat))
+            out_e_vel = os.path.join(a_tmpdir,
+                                     "%d.%s.peer_e.vel" % (sim_id, stat))
+            out_z_vel = os.path.join(a_tmpdir,
+                                     "%d.%s.peer_z.vel" % (sim_id, stat))
+            bbp_formatter.bbp2peer(bbpfile, out_n_vel, out_e_vel, out_z_vel, accel=False)
+
+            # Let's have rotD50 create these output files
+            out_rotd_base = "%d.%s.rdvel" % (sim_id, stat)
+            tmp_rotd = os.path.join(a_tmpdir, out_rotd_base)
+            out_rotd = os.path.join(a_outdir, out_rotd_base)
+
+            # Run the rotD50 program
+            self.do_rotd_vel(a_tmpdir, out_e_vel, out_n_vel, out_z_vel,
+                           out_rotd, self.log)
+
+            cmd = "cp %s %s" % (tmp_rotd, out_rotd)
+            bband_utils.runprog(cmd, abort_on_error=True, print_cmd=False)
+
+        # All done!
+        print("RotDVel Completed".center(80, '-'))
+
+if __name__ == '__main__':
+    print("Testing Module: %s" % (os.path.basename(sys.argv[0])))
+    ME = RotDVel(sys.argv[1], sim_id=int(sys.argv[2]))
+    ME.run()
+    sys.exit(0)