nhlib demo

nhlib-demo.py

import nhlib, nhlib.site, nhlib.calc


fault_mfd = nhlib.mfd.EvenlyDiscretizedMFD(min_mag=5.5, bin_width=1,
                                           occurrence_rates=[0.1, 0.01, 0.001])
fault_trace = nhlib.geo.Line([nhlib.geo.Point(*coords)
                              for coords in [(10, 7), (10.4, 7), (10.6, 6.8)]])
fault = nhlib.source.SimpleFaultSource(
    'f1', 'fault 1', tectonic_region_type=nhlib.const.TRT.STABLE_CONTINENTAL,
    mfd=fault_mfd, rupture_mesh_spacing=5, rupture_aspect_ratio=1.5,
    magnitude_scaling_relationship=nhlib.scalerel.WC1994(),
    upper_seismogenic_depth=1.0, lower_seismogenic_depth=20.0,
    fault_trace=fault_trace, dip=60, rake=90)


from nhlib.gsim.chiou_youngs_2008 import ChiouYoungs2008
gsims = {nhlib.const.TRT.STABLE_CONTINENTAL: ChiouYoungs2008()}


import numpy
lons = numpy.mgrid[9.5:11.1:0.02]
lats = numpy.mgrid[6.4:7.3:0.02]
sites = [nhlib.site.Site(nhlib.geo.Point(lon, lat), vs30=760,
                         vs30measured=False, z1pt0=100, z2pt5=5)
         for lat in lats for lon in lons]
sitecol = nhlib.site.SiteCollection(sites)


imt = nhlib.imt.PGA()
imls = 10 ** numpy.linspace(-4, 1)
imts = {imt: imls}
time_span = 5
curves = nhlib.calc.hazard_curves_poissonian([fault], sitecol, imts, time_span,
                                             gsims=gsims, truncation_level=None)
curves_pga = curves[imt]


from matplotlib import pyplot as plt
iml = 30  # index in imls array
plt.imshow(curves_pga[:, iml].reshape((len(lats), len(lons))),
           aspect='equal', origin='lower',
           extent=(lons.min(), lons.max(), lats.min(), lats.max()))
plt.xlabel('longitude')
plt.ylabel('latitude')
plt.colorbar().set_label('Probability of PGA to exceed %.3f g' % imls[iml])
plt.show()