Merge pull request #10148 from gem/gmm_basin

Refactoring to homogenise GMM basin terms

openquake/hazardlib/gsim/abrahamson_2014.py

@@ -240,7 +240,7 @@ def _get_site_response_term(C, imt, vs30, sa1180):
     return site_resp_term
 
 
-def _get_soil_depth_term(region, C, z1pt0, vs30):
+def _get_basin_term(region, C, z1pt0, vs30):
     """
     Compute and return soil depth term.  See page 1042.
     """
@@ -392,7 +392,7 @@ def _get_sa_at_1180(region, C, imt, ctx):
             _get_site_response_term(C, imt, vs30_1180, ref_iml) +
             _get_hanging_wall_term(C, ctx) +
             _get_top_of_rupture_depth_term(C, imt, ctx) +
-            _get_soil_depth_term(region, C, fake_z1pt0, vs30_1180) +
+            _get_basin_term(region, C, fake_z1pt0, vs30_1180) +
             _get_regional_term(region, C, imt, vs30_1180, ctx.rrup))
 
 def get_epistemic_sigma(ctx):
@@ -487,7 +487,7 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
             # f5 = _get_site_response_term(C, imt, ctx.vs30, sa1180)
             # f6 = _get_top_of_rupture_depth_term(C, imt, ctx)
             # f7 = _get_faulting_style_term(C, ctx)
-            # f10 =_get_soil_depth_term(self.region, C, ctx.z1pt0, ctx.vs30)
+            # f10 = _get_basin_term(self.region, C, ctx.z1pt0, ctx.vs30)
             # fre = _get_regional_term(self.region, C, imt, ctx.vs30, ctx.rrup)
 
             # get the mean value
@@ -496,8 +496,7 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
                        _get_site_response_term(C, imt, ctx.vs30, sa1180) +
                        _get_top_of_rupture_depth_term(C, imt, ctx) +
                        _get_faulting_style_term(C, ctx) +
-                       _get_soil_depth_term(self.region, C, ctx.z1pt0,
-                                            ctx.vs30))
+                       _get_basin_term(self.region, C, ctx.z1pt0, ctx.vs30))
 
             mean[m] += _get_regional_term(
                 self.region, C, imt, ctx.vs30, ctx.rrup)

openquake/hazardlib/gsim/abrahamson_gulerce_2020.py

@@ -288,7 +288,7 @@ def get_reference_basin_depth(region, vs30):
     return np.exp(ln_zref)
 
 
-def get_basin_depth_scaling(C, region, vs30, z25):
+def _get_basin_term(C, region, vs30, z25):
     """
     Returns the basin depth scaling term, applicable only for the Cascadia
     and Japan regions, defined in equations 3.9 - 3.11 and corrected in the
@@ -353,7 +353,7 @@ def get_mean_acceleration(C, trt, region, ctx, pga1000, apply_adjustment):
             get_rupture_depth_scaling_term(C, trt, ctx) +
             get_inslab_scaling_term(C, trt, region, ctx.mag, ctx.rrup) +
             get_site_amplification_term(C, region, ctx.vs30, pga1000) +
-            get_basin_depth_scaling(C, region, ctx.vs30, z25))
+            _get_basin_term(C, region, ctx.vs30, z25))
 
 
 def _get_f2(t1, t2, t3, t4, alpha, period):

openquake/hazardlib/gsim/abrahamson_silva_2008.py

@@ -157,7 +157,7 @@ def _compute_large_distance_term(C, ctx):
     return large_distance_term
 
 
-def _compute_soil_depth_term(C, imt, z1pt0, vs30):
+def _get_basin_term(C, imt, z1pt0, vs30):
     """
     Compute and return soil depth model term, that is the 9-th term in
     equation 1, page 74. The calculation of this term is explained in
@@ -189,7 +189,7 @@ def _compute_imt1100(C_PGA, ctx):
             _compute_hanging_wall_term(C_PGA, ctx) +
             _compute_top_of_rupture_depth_term(C_PGA, ctx) +
             _compute_large_distance_term(C_PGA, ctx) +
-            _compute_soil_depth_term(C_PGA, imt, ctx.z1pt0, vs30_1100) +
+            _get_basin_term(C_PGA, imt, ctx.z1pt0, vs30_1100) +
             # this is the site response term in case of vs30=1100
             ((C_PGA['a10'] + C_PGA['b'] * CONSTS['n']) *
              np.log(vs30_star / C_PGA['VLIN'])))
@@ -474,7 +474,7 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
                        _compute_hanging_wall_term(C, ctx) +
                        _compute_top_of_rupture_depth_term(C, ctx) +
                        _compute_large_distance_term(C, ctx) +
-                       _compute_soil_depth_term(C, imt, ctx.z1pt0, ctx.vs30))
+                       _get_basin_term(C, imt, ctx.z1pt0, ctx.vs30))
 
             sig[m], tau[m], phi[m] = _get_stddevs(C, C_PGA, pga1100, ctx)
 

openquake/hazardlib/gsim/afshari_stewart_2016.py

@@ -132,15 +132,21 @@ def get_magnitude_term(C, ctx):
     return term
 
 
+def _get_basin_term(C, ctx, region):
+    """
+    Return the basin term (equation 9)
+    """
+    dz1 = ctx.z1pt0 - np.exp(_get_lnmu_z1(region, ctx.vs30))
+    fb = C['c5'] * dz1
+    fb[dz1 > CONSTANTS["dz1ref"]] = (C["c5"] * CONSTANTS["dz1ref"])
+    return fb
+
+
 def get_site_amplification(region, C, ctx):
     """
     Returns the site amplification term
     """
-    # Gets delta normalised z1
-    dz1 = ctx.z1pt0 - np.exp(_get_lnmu_z1(region, ctx.vs30))
-    f_s = C["c5"] * dz1
-    # Calculates site amplification term
-    f_s[dz1 > CONSTANTS["dz1ref"]] = (C["c5"] * CONSTANTS["dz1ref"])
+    f_s = _get_basin_term(C, ctx, region) # First get basin term
     idx = ctx.vs30 > CONSTANTS["v1"]
     f_s[idx] += (C["c4"] * np.log(CONSTANTS["v1"] / C["vref"]))
     idx = np.logical_not(idx)

openquake/hazardlib/gsim/aristeidou_2023.py

@@ -142,7 +142,7 @@ def _get_site_amplification_term(C, ctx):
     return f_s
 
 
-def _get_basin_response_term(C, ctx):
+def _get_basin_term(C, ctx):
     """
     Returns the basin response term defined in equation (9), p. 1611
     """
@@ -258,7 +258,7 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
             f_s = _get_site_amplification_term(C, ctx)
 
             # Basin-effects correction function:
-            f_basin = _get_basin_response_term(C, ctx)
+            f_basin = _get_basin_term(C, ctx)
 
             # log of the functional form considered
             mean[m] = np.squeeze(C["a"] + f_m + f_d + f_sof + f_s + f_basin)

openquake/hazardlib/gsim/bahrampouri_2021_duration.py

@@ -61,18 +61,25 @@ def _get_path_term(C, ctx):
     return fpath
 
 
-def _get_site_term(C, ctx):
+def _get_basin_term(C, ctx):
     """
-    Implementing Eqs. 5, 6 and 12
+    Get the basin term
     """
     mean_z1pt0 = (np.exp(((-5.23 / 2.) * np.log((ctx.vs30 ** 2. +
                   412.39 ** 2.) / (1360 ** 2. + 412.39 ** 2.)))-0.9))
     delta_z1pt0 = np.round(ctx.z1pt0 - mean_z1pt0, 4)
-    fsite = []
-    for i, value in enumerate(delta_z1pt0):
-        s = (np.round(C['s1'] * np.log(min(ctx.vs30[i], 600.) / 600.) +
-             C['s2']*min(delta_z1pt0[i], 250.0) + C['s3'], 4))
-        fsite.append(s)
+    return C['s2'] * np.minimum(delta_z1pt0, 250.0)
+
+
+def _get_site_term(C, ctx):
+    """
+    Implementing Eqs. 5, 6 and 12
+    """
+    fbasin = _get_basin_term(C, ctx)
+
+    fsite = np.round(C['s1'] * np.log(np.clip(ctx.vs30, None, 600.0) / 600.0) +
+                     fbasin + C['s3'], 4)
+
     return fsite
 
 

openquake/hazardlib/gsim/bayless_abrahamson_2018.py

@@ -86,7 +86,7 @@ def _linear_site_response(C, ctx):
     return fsl
 
 
-def _soil_depth_scaling(C, ctx):
+def _get_basin_term(C, ctx):
     """ Compute the soil depth scaling term """
     # Set the c11 coefficient - See eq.13b at page 2093
     c11 = np.ones_like(ctx.vs30) * C['c11a']
@@ -190,7 +190,7 @@ def _get_mean_linear(C, ctx):
             _linear_site_response(C, ctx) +
             _ztor_scaling(C, ctx) +
             _normal_fault_effect(C, ctx) +
-            _soil_depth_scaling(C, ctx))
+            _get_basin_term(C, ctx))
     return mean
 
 

openquake/hazardlib/gsim/boore_2014.py

@@ -40,7 +40,7 @@
     "v2": 300.0}
 
 
-def _get_basin_depth_term(region, C, ctx, period):
+def _get_basin_term(region, C, ctx, period):
     """
     In the case of the base model the basin depth term is switched off.
     Therefore we return an array of zeros.
@@ -204,7 +204,7 @@ def _get_site_scaling(kind, region, C, pga_rock, ctx, period, rjb):
     if kind == 'stewart':
         fbd = 0.  # there is no basin term in the Stewart models
     else:
-        fbd = _get_basin_depth_term(region, C, ctx, period)
+        fbd = _get_basin_term(region, C, ctx, period)
     return flin + fnl + fbd
 
 

openquake/hazardlib/gsim/boore_2020.py

@@ -101,7 +101,7 @@ def _get_pga_on_rock(C, rup, dists):
 def _get_site_scaling(C, pga_rock, sites):
     """
     Returns the site-scaling term (equation 5), broken down into a
-    linear scaling, a nonlinear scaling and a basin scaling term
+    linear scaling and a nonlinear scaling
     """
     flin = _get_linear_site_term(C, sites.vs30)
     fnl = _get_nonlinear_site_term(C, sites.vs30, pga_rock)

openquake/hazardlib/gsim/bozorgnia_campbell_2016.py

@@ -45,7 +45,7 @@ def _get_anelastic_attenuation_term(sgn, C, rrup):
     return f_atn
 
 
-def _get_basin_response_term(SJ, C, z2pt5):
+def _get_basin_term(SJ, C, z2pt5):
     """
     Returns the basin response term, f_sed, defined in equation 20
 
@@ -150,7 +150,7 @@ def get_mean_values(SJ, sgn, C, ctx):
             _get_style_of_faulting_term(C, ctx) +
             _get_hanging_wall_term(C, ctx) +
             _get_shallow_site_response_term(SJ, C, ctx.vs30) +
-            _get_basin_response_term(SJ, C, temp_z2pt5) +
+            _get_basin_term(SJ, C, temp_z2pt5) +
             _get_hypocentral_depth_term(C, ctx) +
             _get_fault_dip_term(C, ctx) +
             _get_anelastic_attenuation_term(sgn, C, ctx.rrup))

openquake/hazardlib/gsim/bradley_2013.py

@@ -46,6 +46,13 @@
      (172.6220, -43.5233)])
 
 
+def _get_basin_term(C, z1pt0):
+    fb1 = C['phi5'] * (1.0 - 1.0 / np.cosh(
+        C['phi6'] * (z1pt0 - C['phi7']).clip(0, np.inf)))
+    fb2 = C['phi8'] / np.cosh(0.15 * (z1pt0 - 15).clip(0, np.inf))
+    return fb1 + fb2
+
+
 def _adjust_mean_model(region, in_cshm, in_cbd, imt_per, b13_mean):
     dL2L = dS2S = np.array(np.zeros(np.shape(b13_mean)))
     # If the site is in the CBD polygon, get dL2L and dS2S terms
@@ -357,11 +364,7 @@ def _get_mean(ctx, C, ln_y_ref, exp1, exp2, v1):
 
     Implements eq. 5
     """
-    # we do not support estimating of basin depth and instead
-    # rely on it being available (since we require it).
-    z1pt0 = ctx.z1pt0
-
-    # we consider random variables being zero since we want
+    # We consider random variables being zero since we want
     # to find the exact mean value.
     eta = epsilon = 0
 
@@ -373,10 +376,7 @@ def _get_mean(ctx, C, ln_y_ref, exp1, exp2, v1):
         + C['phi2'] * (exp1 - exp2)
         * np.log((np.exp(ln_y_ref) + C['phi4']) / C['phi4'])
         # third line
-        + C['phi5']
-        * (1.0 - 1.0 / np.cosh(
-            C['phi6'] * (z1pt0 - C['phi7']).clip(0, np.inf)))
-        + C['phi8'] / np.cosh(0.15 * (z1pt0 - 15).clip(0, np.inf))
+        + _get_basin_term(C, ctx.z1pt0)
         # fourth line
         + eta + epsilon)
 

openquake/hazardlib/gsim/campbell_bozorgnia_2008.py

@@ -27,7 +27,7 @@
 from openquake.hazardlib.imt import PGA, PGV, PGD, CAV, SA
 
 
-def _compute_basin_response_term(C, z2pt5):
+def _get_basin_term(C, z2pt5):
     """
     Returns the basin response term (equation 12, page 146)
     """
@@ -76,7 +76,7 @@ def _compute_imt1100(C, ctx, get_pga_site=False):
                      _compute_distance_term(C, ctx) +
                      _compute_style_of_faulting_term(C, ctx) +
                      _compute_hanging_wall_term(C, ctx) +
-                     _compute_basin_response_term(C, ctx.z2pt5) +
+                     _get_basin_term(C, ctx.z2pt5) +
                      fsite)
     # If PGA at the site is needed then remove factor for rock and
     # re-calculate on correct site condition
@@ -321,7 +321,7 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
                        _compute_style_of_faulting_term(C, ctx) +
                        _compute_hanging_wall_term(C, ctx) +
                        _compute_shallow_site_response(C, ctx, pga1100) +
-                       _compute_basin_response_term(C, ctx.z2pt5))
+                       _get_basin_term(C, ctx.z2pt5))
 
             # If it is necessary to ensure that Sa(T) >= PGA
             # (see previous comment)

openquake/hazardlib/gsim/campbell_bozorgnia_2014.py

@@ -59,7 +59,7 @@ def _get_anelastic_attenuation_term(C, rrup):
     return f_atn
 
 
-def _get_basin_response_term(SJ, C, z2pt5):
+def _get_basin_term(SJ, C, z2pt5):
     """
     Returns the basin response term defined in equation 20
     """
@@ -309,7 +309,7 @@ def get_mean_values(SJ, C, ctx, a1100=None):
             _get_style_of_faulting_term(C, ctx) +
             _get_hanging_wall_term(C, ctx) +
             _get_shallow_site_response_term(SJ, C, temp_vs30, a1100) +
-            _get_basin_response_term(SJ, C, temp_z2pt5) +
+            _get_basin_term(SJ, C, temp_z2pt5) +
             _get_hypocentral_depth_term(C, ctx) +
             _get_fault_dip_term(C, ctx) +
             _get_anelastic_attenuation_term(C, ctx.rrup))

openquake/hazardlib/gsim/can20/can_shm6_active_crust.py

@@ -86,7 +86,7 @@ def _get_site_scaling_ba14(kind, region, C, pga_rock, sites, period, rjb):
         BSSA14_1100[0], BSSA14_2000[0], sites.vs30[sites.vs30 > 1100], imt)
 
     fnl = BA14._get_nonlinear_site_term(C, sites.vs30, pga_rock)
-    fbd = BA14._get_basin_depth_term(region, C, sites, period)  # returns 0
+    fbd = BA14._get_basin_term(region, C, sites, period)  # returns 0
     fbd = 0.0
 
     return flin + fnl + fbd
@@ -205,7 +205,7 @@ class CanadaSHM6_ActiveCrust_ChiouYoungs2014(ChiouYoungs2014):
     """
     #: Required site parameters are Vs30, Vs30 measured flag
     #: and Z1.0.
-    REQUIRES_SITES_PARAMETERS = {'vs30', 'vs30measured'}
+    REQUIRES_SITES_PARAMETERS = {'vs30', 'vs30measured', 'z1pt0'}
 
     def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
         """

openquake/hazardlib/gsim/can20/can_shm6_interface.py

@@ -260,7 +260,7 @@ class CanadaSHM6_Interface_AtkinsonMacias2009(AtkinsonMacias2009):
     """
 
     REQUIRES_DISTANCES = {'rrup', 'rjb'}
-    REQUIRES_SITES_PARAMETERS = {'vs30', 'z1pt0'}
+    REQUIRES_SITES_PARAMETERS = {'vs30'}
     DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, PGV, SA}
 
     def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):

openquake/hazardlib/gsim/chao_2020.py

@@ -123,7 +123,7 @@ def _fvs30(geology, C, ctx):
                     C['c27'] if geology else C['c28'])
 
 
-def _fz1pt0(C, ctx):
+def _get_basin_term(C, ctx):
     """
     z1pt0 factor.
     """
@@ -217,7 +217,7 @@ def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
             sa1180 = np.exp(med + math.log(1180/s['vs30_ref']) * C['c24'])
             med += _fc(C, imt, ctx.vs30, sa1180)
             med += np.log(ctx.vs30 / s['vs30_ref']) * C['c24']
-            med += _fz1pt0(C, ctx)
+            med += _get_basin_term(C, ctx)
 
             sig[m], tau[m], phi[m] = get_stddevs(self.SBCR, C, ctx.mag)
 

openquake/hazardlib/gsim/chiou_youngs_2008.py

@@ -26,6 +26,19 @@
 from openquake.hazardlib.imt import PGA, PGV, SA
 
 
+def _get_basin_term(C, z1pt0):
+    """
+    Return the basin term describing effects of deep sediment sites and shallow
+    sediment sites through z1pt0 
+    """
+    # Equation 3.10
+    deep_soil = C['phi5'] * (1.0 - 1.0 / np.cosh(C['phi6'] * (
+        z1pt0 - C['phi7']).clip(0, np.inf)))
+    # Equation 3.11
+    shallow_soil = C['phi8'] / np.cosh(0.15 * (z1pt0 - 15).clip(0, np.inf))
+    return deep_soil + shallow_soil
+
+
 def _get_ln_y_ref(ctx, C):
     """
     Get an intensity on a reference soil.
@@ -98,10 +111,7 @@ def _get_mean(ctx, C, ln_y_ref, exp1, exp2):
         + C['phi2'] * (exp1 - exp2)
         * np.log((np.exp(ln_y_ref) + C['phi4']) / C['phi4'])
         # third line
-        + C['phi5']
-        * (1.0 - 1.0 / np.cosh(
-            C['phi6'] * (z1pt0 - C['phi7']).clip(0, np.inf)))
-        + C['phi8'] / np.cosh(0.15 * (z1pt0 - 15).clip(0, np.inf))
+        + _get_basin_term(C, z1pt0)
         # fourth line
         + eta + epsilon
     )