Time series analysis to compute surface response spectrum and site amplification functions using velocity profiles from geopsy.
In [1]:
import re
import matplotlib.pyplot as plt
import numpy as np
import pysra
%matplotlib inline
In [2]:
# Increased figure sizes
plt.rcParams['figure.dpi'] = 120
In [3]:
def iter_geopsy_profiles(fname):
"""Read a Geopsy formatted text file created by gpdcreport."""
with open(fname) as fp:
next(fp)
while True:
try:
line = next(fp)
except StopIteration:
break
m = re.search(r'Layered model (\d+): value=([0-9.]+)', line)
id, score = m.groups()
count = int(next(fp))
d = {
'id': id,
'score': score,
'layers': [],
}
cols = ['thickness', 'vel_comp', 'vel_shear', 'density']
for _ in range(count):
values = [float(p) for p in next(fp).split()]
d['layers'].append(dict(zip(cols, values)))
yield d
In [4]:
fname = 'data/NIS090.AT2'
ts = pysra.motion.TimeSeriesMotion.load_at2_file(fname)
ts.accels
Out[4]:
In [5]:
calc = pysra.propagation.LinearElasticCalculator()
In [6]:
freqs = np.logspace(-1, 2, num=500)
outputs = pysra.output.OutputCollection([
pysra.output.ResponseSpectrumOutput(
# Frequency
freqs,
# Location of the output
pysra.output.OutputLocation('outcrop', index=0),
# Damping
0.05),
pysra.output.ResponseSpectrumRatioOutput(
# Frequency
freqs,
# Location in (denominator),
pysra.output.OutputLocation('outcrop', index=-1),
# Location out (numerator)
pysra.output.OutputLocation('outcrop', index=0),
# Damping
0.05),
])
In [7]:
fname = 'data/best100_GM_linux.txt'
for geopsy_profile in iter_geopsy_profiles(fname):
profile = pysra.site.Profile([
pysra.site.Layer(
pysra.site.SoilType(
'soil-%d' % i, l['density'] / pysra.site.GRAVITY,
damping=0.05), l['thickness'], l['vel_shear'])
for i, l in enumerate(geopsy_profile['layers'])
])
# Use 1% damping for the half-space
profile[-1].soil_type.damping = 0.01
# Compute the waves from the last layer
calc(ts, profile, profile.location('outcrop', index=-1))
# Compute the site amplification
outputs(calc)
In [8]:
for o in outputs:
fig, ax = plt.subplots()
ax.plot(o.refs, o.values, color='C0', alpha=0.6, linewidth=0.2)
ax.set(xlabel=o.xlabel, xscale='log', ylabel=o.ylabel)
fig.tight_layout();
In [ ]: