Repository navigation
Site_Factors() rejects Vs30 in (400, 450] m/s with z1 > 750 m, although the 450 m/s, 900 m data ships and loads #104
Description
Found while reviewing #74. Pre-existing on main (b764057); #74 does not change the behavior.
Problem
Site_Factors._range_check() marks Vs30 > 400 and z1 > 750 as an invalid Vs30-z1 combination:
PySeismoSoil/PySeismoSoil/class_site_factors.py
Lines 775 to 788 in b764057
This may be intended; if so, the reason for 750 m is not written down: the code comment only says that a stiff site cannot have a deep basin. The facts:
- The other seven rules match the shipped data. For each rule "Vs30 > V and z1 > Z", the next Vs30 on the grid has data down to exactly Z. The first rule is the exception: Vs30 = 450 m/s has data down to 900 m.
- 750 m is not on the z1 grid (
[8, 16, 24, 36, 75, 150, 300, 450, 600, 900]); all the other thresholds are. - The
450_900amplification and phase files ship and load without NaN.Site_Factors(450, 900, 0.3)still raises, and so doesSite_Factors(430, 800, 0.3), whose four corner grids (400 and 450 m/s, 600 and 900 m) all have data.
test_range_check() asserts the current behavior (same name in #74):
PySeismoSoil/tests/test_class_site_factors.py
Lines 20 to 21 in b764057
from pathlib import Path
import numpy as np
from PySeismoSoil.class_site_factors import Site_Factors as SF
sf = SF(400, 900, 0.3)
amp_dir, phase_dir = Path(sf.dir_amplif), Path(sf.dir_phase)
def deepest_z1(vs30): # deepest z1 on the grid that has data for this Vs30
return max(z1 for z1 in SF.z1_array if (amp_dir / f'{vs30}_{z1:03d}_freq.csv').exists())
# Each rule "Vs30 > V and z1 > Z" vs. the data of the next Vs30 on the grid
rules = [(400, 750), (450, 600), (550, 450), (600, 300), (650, 150), (750, 75), (800, 36), (850, 16)]
for v, z in rules:
v_next = SF.Vs30_array[SF.Vs30_array.index(v) + 1]
print(f'rule z1 > {z} for Vs30 > {v}; deepest data for Vs30 = {v_next}: {deepest_z1(v_next)}')
# rule z1 > 750 for Vs30 > 400; deepest data for Vs30 = 450: 900
# rule z1 > 600 for Vs30 > 450; deepest data for Vs30 = 500: 600
# rule z1 > 450 for Vs30 > 550; deepest data for Vs30 = 600: 450
# rule z1 > 300 for Vs30 > 600; deepest data for Vs30 = 650: 300
# rule z1 > 150 for Vs30 > 650; deepest data for Vs30 = 700: 150
# rule z1 > 75 for Vs30 > 750; deepest data for Vs30 = 800: 75
# rule z1 > 36 for Vs30 > 800; deepest data for Vs30 = 850: 36
# rule z1 > 16 for Vs30 > 850; deepest data for Vs30 = 900: 16
print(750 in SF.z1_array) # False
# The 450_900 files ship (10 amplification files, 3 phase files) and load without NaN
print(len(list(amp_dir.glob('450_900_*'))), len(list(phase_dir.glob('450_900_*')))) # 10 3
_, amp = SF._query('amplif', 450, 900, 0.3, data_dir=sf.dir_amplif)
_, phase = SF._query('phase', 450, 900, 0.3, method='eq_hh', data_dir=sf.dir_phase)
print(amp[:3].round(3), bool(np.isnan(amp).any() or np.isnan(phase).any())) # [1.064 1.061 1.061] False
# Site_Factors() rejects that grid point, and points whose 4 corners (400/450 m/s, 600/900 m) all have data
for vs30, z1 in [(450, 900), (430, 800), (450, 750)]:
try:
SF(vs30, z1, 0.3)
print(vs30, z1, 'accepted')
except ValueError as e:
print(vs30, z1, type(e).__name__, e)
# 450 900 ValueError Vs30 and z1 combination not valid. (The `lenient` option does not apply to this type of issue.)
# 430 800 ValueError Vs30 and z1 combination not valid. (The `lenient` option does not apply to this type of issue.)
# 450 750 acceptedFix
- If the 750 m limit is intended, add a comment next to the rule that says where it comes from. Add a sentence to the
Site_Factorsdocstring saying that Vs30 in (400, 450] m/s is limited to z1 <= 750 m, even though data down to 900 m exists. - If it is not intended, remove the
(Vs30_in_mps > 400 and z1_in_m > 750)rule. Intest_range_check(), change the assertion for(450, 751)to expect[], and add(450, 900)and(451, 601)cases. With the rule removed,get_amplification(),get_amplification(Fourier=False)andget_phase_shift()return results without NaN for (450, 900), (430, 800) and (401, 751).
Activity