From 5cd98c559352c72e8275151819ec7757cdac0332 Mon Sep 17 00:00:00 2001 From: Kevin Milner Date: Fri, 5 Nov 2021 16:32:23 -0700 Subject: [PATCH 1/4] added additional periods to RotD100 calculation --- bbp/src/ucb/rotd50/rotd100.f | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/bbp/src/ucb/rotd50/rotd100.f b/bbp/src/ucb/rotd50/rotd100.f index 25f60486..e712afca 100644 --- a/bbp/src/ucb/rotd50/rotd100.f +++ b/bbp/src/ucb/rotd50/rotd100.f @@ -30,22 +30,22 @@ program Calc_RotD50 real dt10 real x0(MAXPTS), y0(MAXPTS), u(MAXPTS), y2(MAXPTS) real rspTH1(MAXPTS), rspTH2(MAXPTS), rsp1(MAXPTS), rsp2(MAXPTS) - real x(MAXPTS), y(MAXPTS), rsp_Period(63), famp15(3) + real x(MAXPTS), y(MAXPTS), rsp_Period(65), famp15(3) real rotangle, w(200), sa(1000), workArray(1000), rotD50(3,200), rotD100(3,200), psa5E(200), psa5N(200) integer rD100ang(3,200), rD50ang(3,200) real saUnsort(1000) real damping complex cu1(MAXPTS) - data RSP_Period / 0.010, 0.011, 0.012, 0.013, 0.015, 0.017, 0.020, 0.022, 0.025, 0.029, - 1 0.032, 0.035, 0.040, 0.045, 0.050, 0.055, 0.060, 0.065, 0.075, 0.085, - 2 0.100, 0.110, 0.120, 0.130, 0.150, 0.170, 0.200, 0.220, 0.240, 0.260, - 3 0.280, 0.300, 0.350, 0.400, 0.450, 0.500, 0.550, 0.600, 0.650, 0.750, - 4 0.850, 1.000, 1.100, 1.200, 1.300, 1.500, 1.700, 2.000, 2.200, 2.400, - 5 2.600, 2.800, 3.000, 3.500, 4.000, 4.400, 5.000, 5.500, 6.000, 6.500, - 6 7.500, 8.500, 10.000 / + data RSP_Period / 0.0001, 0.001, 0.010, 0.011, 0.012, 0.013, 0.015, 0.017, 0.020, 0.022, + 1 0.025, 0.029, 0.032, 0.035, 0.040, 0.045, 0.050, 0.055, 0.060, 0.065, + 2 0.075, 0.085, 0.100, 0.110, 0.120, 0.130, 0.150, 0.170, 0.200, 0.220, + 3 0.240, 0.260, 0.280, 0.300, 0.350, 0.400, 0.450, 0.500, 0.550, 0.600, + 4 0.650, 0.750, 0.850, 1.000, 1.100, 1.200, 1.300, 1.500, 1.700, 2.000, + 5 2.200, 2.400, 2.600, 2.800, 3.000, 3.500, 4.000, 4.400, 5.000, 5.500, + 6 6.000, 6.500, 7.500, 8.500, 10.000 / - nFreq = 63 + nFreq = 65 damping = 0.05 dt_max = 0.001 From ed9347c7d6c86e985af91103d2902c74cc1e53e1 Mon Sep 17 00:00:00 2001 From: Kevin Milner Date: Fri, 5 Nov 2021 16:33:11 -0700 Subject: [PATCH 2/4] new rotd velocity module --- bbp/comps/bbp_formatter.py | 24 ++++-- bbp/comps/rotd_vel.py | 164 +++++++++++++++++++++++++++++++++++++ 2 files changed, 181 insertions(+), 7 deletions(-) create mode 100644 bbp/comps/rotd_vel.py diff --git a/bbp/comps/bbp_formatter.py b/bbp/comps/bbp_formatter.py index 217f6578..0279e0a0 100644 --- a/bbp/comps/bbp_formatter.py +++ b/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 diff --git a/bbp/comps/rotd_vel.py b/bbp/comps/rotd_vel.py new file mode 100644 index 00000000..0650d987 --- /dev/null +++ b/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) From 881ff581b9e00bc558ea860f405d0799e95d5b1d Mon Sep 17 00:00:00 2001 From: Kevin Milner Date: Fri, 5 Nov 2021 16:33:54 -0700 Subject: [PATCH 3/4] turned arias duration tool into a module --- bbp/comps/arias_duration.py | 118 ++++++++++++++++++++++++++++++++++-- 1 file changed, 114 insertions(+), 4 deletions(-) diff --git a/bbp/comps/arias_duration.py b/bbp/comps/arias_duration.py index d16e85ab..0bf8da5e 100644 --- a/bbp/comps/arias_duration.py +++ b/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) From 940cf6b5130fb361ae8747742dcdcd311ceb97d7 Mon Sep 17 00:00:00 2001 From: Kevin Milner Date: Fri, 5 Nov 2021 16:34:21 -0700 Subject: [PATCH 4/4] imports for new rotd vel and arias duration modules --- bbp/comps/module.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/bbp/comps/module.py b/bbp/comps/module.py index bd792406..89dc25a9 100644 --- a/bbp/comps/module.py +++ b/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):