Using GMPEs with OpenQuake.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preamble"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this exercise we will look at how to run the GMPEs in OpenQuake's extensive GMPE library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "# Import our usual Numpy and Matplotlib tools\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Which GMPEs are available to me?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If at any time you want to see which GMPEs are available on your version of OpenQuake you can just run this command.\n",
    "\n",
    "There tell you the names you need to use when calling the GMPEs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the function to tell me which GMPEs are available\n",
    "from openquake.hazardlib.gsim import get_available_gsims\n",
    "\n",
    "# Print out the list\n",
    "for gmpe in get_available_gsims():\n",
    "    print(gmpe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <i>Everything in Context</i>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make every GMPE run in the same way, OpenQuake needs you to organise your inputs into <b>source</b> parameters, <b>path (distance)</b> paramters and <b>site</b> parameters\n",
    "\n",
    "<b>source</b> parameters refer to the properties of the rupture, e.g. the magnitude, the rake (or style-of faulting),\n",
    "the top of rupture depth, hypocentral depth. These should always be scalar values (i.e. a single number)\n",
    "\n",
    "<b>distances</b> refer to the source-to-site distance metrics (e.g. epicentral distsnce, hypocentral distance, Joyner-Boore distance etc. These must always be vector values (i.e. an array of numbers)\n",
    "\n",
    "<b>site</b> refers to the site properties (e.g. $V_{S30}$, basin depth}, and should always be an array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openquake.hazardlib.gsim.base import DistancesContext, SitesContext, RuptureContext"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, let's look at a simple example. We will take a simple GMPE (Bindi et al., 2014) and calculate the peak ground acceleration and 1.0 s spectral acceleration at a site 20 km Joyner-Boore distance from a $M_W$ 6.0 normal-faulting earthquakes, for a $V_{S30}$ of 500 m/s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openquake.hazardlib.gsim.bindi_2014 import BindiEtAl2014Rjb\n",
    "from openquake.hazardlib.imt import PGA, SA\n",
    "from openquake.hazardlib import const"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now set up the scenario"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is quite common to encounter a style-of-faulting term, e.g. normal, strike-slip, reverse.\n",
    "\n",
    "OpenQuake requires that style of faulting is described by the rake angle of the fault, which follows this convention (Aki & Richards, 2002):\n",
    "\n",
    "Normal = $-90^{\\circ}$\n",
    "\n",
    "Strike-slip = $0^{\\circ}$ or $180^{\\circ}$\n",
    "\n",
    "Reverse = $90^{\\circ}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rctx = RuptureContext()\n",
    "rctx.mag = 6.0\n",
    "rctx.rake = -90.0\n",
    "\n",
    "dctx = DistancesContext()\n",
    "dctx.rjb = np.array([20.0])\n",
    "\n",
    "sctx = SitesContext()\n",
    "sctx.vs30 = np.array([500.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gmpe = BindiEtAl2014Rjb()\n",
    "imt = PGA()\n",
    "median_pga, [stddev_pga] = gmpe.get_mean_and_stddevs(sctx, rctx, dctx, imt, [const.StdDev.TOTAL])\n",
    "print(\"Median PGA = %.5f, Total Std. Dev = %.5f\" % (np.exp(median_pga), stddev_pga))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imt = SA(1.0)\n",
    "median_sa1, [stddev_sa1] = gmpe.get_mean_and_stddevs(sctx, rctx, dctx, imt, [const.StdDev.TOTAL])\n",
    "print(\"Median Sa(1.0) = %.5f, Total Std. Dev = %.5f\" % (np.exp(median_sa1), stddev_sa1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How does the PGA attenuate with distance?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dctx = DistancesContext()\n",
    "dctx.rjb = np.arange(0., 201., 1.)\n",
    "\n",
    "sctx = SitesContext()\n",
    "sctx.vs30 = 500. * np.ones_like(dctx.rjb) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imt = PGA()\n",
    "median_pga, [stddev_pga] = gmpe.get_mean_and_stddevs(sctx, rctx, dctx, imt, [const.StdDev.TOTAL])\n",
    "plt.figure()\n",
    "plt.plot(dctx.rjb, np.exp(median_pga), \"k-\", lw=2)\n",
    "plt.plot(dctx.rjb, np.exp(median_pga + stddev_pga), \"r--\")\n",
    "plt.plot(dctx.rjb, np.exp(median_pga - stddev_pga), \"r--\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How does the GMPE vary with Spectral Period (for 1st scenario)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dctx = DistancesContext()\n",
    "dctx.rjb = np.array([20.0])\n",
    "\n",
    "sctx = SitesContext()\n",
    "sctx.vs30 = np.array([500.])\n",
    "\n",
    "periods = np.logspace(np.log10(0.05), np.log10(3.0), 100)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "median_pga = np.empty(len(periods))\n",
    "stddev_pga = np.empty(len(periods))\n",
    "\n",
    "# Loop over each period\n",
    "for i, period in enumerate(periods):\n",
    "    median, [stddev] = gmpe.get_mean_and_stddevs(sctx, rctx, dctx, SA(period), [const.StdDev.TOTAL])\n",
    "    median_pga[i] = median[0]\n",
    "    stddev_pga[i] = stddev[0]\n",
    "\n",
    "# Plot the results\n",
    "plt.figure()\n",
    "plt.semilogx(periods, np.exp(median_pga), \"k-\", lw=2)\n",
    "plt.semilogx(periods, np.exp(median_pga + stddev_pga), \"r--\")\n",
    "plt.semilogx(periods, np.exp(median_pga - stddev_pga), \"r--\")\n",
    "plt.grid(True)\n",
    "plt.xlim(0.05, 3.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now it's your turn ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you do the following?\n",
    "\n",
    "1. Show how the PGA changes with magnitude for the same site (20 km from the earthquake source with Vs30 500 m/s)\n",
    "\n",
    "2. Show how the PGA changes with Vs30 for the Mw 6.0 earthquake at a site 20 km from the normal faulting source"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What's the damage?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You are building a new building at the site that is 20 km from the normal fault that can produce a Mw 6.5 earthquake, with the Vs30 = 300 m/s. The engineer tells you that the building can withstand 0.15 g without sustaining any damage, but at 0.5 g it will collapse!\n",
    "\n",
    "If the earthquake were to occur, what would be the probability that the building:\n",
    "1. sustains no damage\n",
    "2. collapses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import a Scipy tool for the normal distribution\n",
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rctx = RuptureContext()\n",
    "rctx.mag = 6.5\n",
    "rctx.rake = -90.0\n",
    "\n",
    "dctx = DistancesContext()\n",
    "dctx.rjb = np.array([20.0])\n",
    "\n",
    "sctx = SitesContext()\n",
    "sctx.vs30 = np.array([300.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imt = PGA()\n",
    "gmpe = BindiEtAl2014Rjb()\n",
    "median_pga, [stddev] = gmpe.get_mean_and_stddevs(sctx, rctx, dctx, imt, [const.StdDev.TOTAL])\n",
    "print(\"Median PGA = %.5f and Std. Deviation = %.4f\" % (np.exp(median_pga), stddev))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the \"no damage\" case - how many standard deviations above or below the median is this 0.15 g threshold?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accel = np.logspace(-3., 1., 1000.)\n",
    "cdf = norm.cdf(np.log(accel), loc=median_pga[0], scale=stddev[0])\n",
    "plt.figure(figsize=(8,8))\n",
    "plt.semilogx(accel, cdf, \"r-\")\n",
    "plt.grid()\n",
    "plt.xlabel(\"g\", fontsize=16)\n",
    "plt.ylabel(\"Probability\", fontsize=16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of standard deviations above/below mean\n",
    "num_stddevs = (np.log(0.15) - median_pga) / stddev\n",
    "print(\"0.15 g is %.3f standard deviations above the median\" % num_stddevs)\n",
    "prob_no_damage = norm.cdf(num_stddevs, loc=0., scale=1.)\n",
    "print(\"The building has a %.4f probability of observing no damage\" % prob_no_damage)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What is the probability of collapse?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_stddevs = (np.log(0.5) - median_pga) / stddev\n",
    "print(\"0.5 g is %.3f standard deviations above the median\" % num_stddevs)\n",
    "# Note the use of the \"survival function [sf]\": 1.0 - cdf\n",
    "prob_collapse = norm.sf(num_stddevs, loc=0., scale=1.)\n",
    "print(\"The building has a %.4f probability of collapsing\" % prob_collapse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}