hw_gr.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Applied Seismology, GEOS 626, University of Alaska Fairbanks\n",
    "\n",
    "- Gutenberg-Richter frequency-magnitude relationship\n",
    "- explores the relationship between cumulative and regular distributions using the GCMT catalog as an example\n",
    "- see Pacheco, Scholz, Sykes (1992 Nature) for details on the two segments of the GCMT distribution"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.io\n",
    "import warnings\n",
    "\n",
    "from lib_seis import seis2GR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# script settings\n",
    "\n",
    "warnings.filterwarnings('ignore')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem 1-1\n",
    "\n",
    "# Load GCMT catalog - Note: currently just loads Matlab database, may want to change later\n",
    "# M: 6 x n set of moment tensors\n",
    "# M0: 1 x n set of scalar seismic moments (derived from moment tensors)\n",
    "# Mw: n x 1 set of moment magnitudes (derived from moment tensors)\n",
    "# dep: n x 1 set of depths\n",
    "# eid: n x 1 set of event IDs\n",
    "# lat: n x 1 set of latitudes\n",
    "# lon: n x 1 set of longitudes\n",
    "# otime: n x 1 set of origin times, in seconds (relative to Matlab t=0) Note: This may be weird in python\n",
    "\n",
    "data = scipy.io.loadmat('./data/cmtall_sub.mat') # Make sure that cmtall_sub.mat is in your working directory \n",
    "print(data.keys()) # This allows you to see what variables are in your dictionary\n",
    "\n",
    "Mw = data['Mw']\n",
    "dep = data['dep']\n",
    "lat = data['lat']\n",
    "lon = data['lon']\n",
    "otime = data['otime']\n",
    "# To see the shape of any variable -> Mw.shape\n",
    "\n",
    "tran_yr = (max(otime) - min(otime)) / 365.25 #tran_yr = \"time range in years\"\n",
    "tran_yr = np.round(tran_yr, decimals = 3)\n",
    "\n",
    "# To see what variables and modules at any time -> %whos\n",
    "print('GCMT catalog duration is %.3f years' %tran_yr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the catalog, colored by depth\n",
    "isort = np.argsort(dep, axis=0)\n",
    "\n",
    "lon[lon<0] = lon[lon<0] + 360 # center map on 180 deg\n",
    "title = 'GCMT catalog (' + str(len(dep)) + ' events), colored by depth (km)'\n",
    "\n",
    "plt.figure(figsize=(14,7),dpi=80)\n",
    "plt.scatter(lon[isort],lat[isort],s=0.5,c=dep[isort],cmap='viridis',vmin=0,vmax=600)\n",
    "plt.xlim(min(lon),max(lon))\n",
    "plt.ylim(min(lat),max(lat))\n",
    "plt.xlabel('Longitude')\n",
    "plt.ylabel('Latitude')\n",
    "plt.title(title)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem 1-2\n",
    "\n",
    "dmag = 0.1 # Magnitude bin width\n",
    "\n",
    "# Note: seis2GR is a homebrewed function. You must have it in your working directory\n",
    "# Note: Python is 0-indexed, i.e. 'bin 0' is the first bin in the array\n",
    "\n",
    "N, Ninc, Medges = seis2GR(Mw,dmag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# try the histogram plotting option that is available\n",
    "idisplay=0\n",
    "ifigure=1\n",
    "N, Ninc, Medges = seis2GR(Mw,dmag,idisplay,ifigure)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem 1-3\n",
    "# Start your work here\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example of plotting multiple items on log-scaled axes\n",
    "# 1) Generate fake data and best-fitting line\n",
    "# 2) Plot lines\n",
    "\n",
    "n = 50\n",
    "x = np.linspace(4,10,num=n)\n",
    "b = 1\n",
    "a = 9\n",
    "e = 0.1\n",
    "err = e*(-1+2*np.random.randn(n))\n",
    "Nd = 10**(a-b*x + err)\n",
    "N  = 10**(a-b*x)\n",
    "N2 = 10**(a-1-b*x)\n",
    "\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.semilogy(x,Nd,color='blue',marker='v',linestyle='none',ms=8,mfc='white',label='fake data') # points\n",
    "# After this, you do not need to specify semilogy\n",
    "plt.plot(x,N,color='k',label='cumulative fit') # a line\n",
    "plt.plot(x,N2,'r',label='incremental fit') # a different line\n",
    "plt.legend(loc='upper right')\n",
    "plt.title('An example of plotting multiple things with a log-y axis')\n",
    "plt.xticks([4, 5, 6, 7, 8, 9, 10])\n",
    "plt.yticks(10.0**np.arange(-1,6),['0.1','1','10','100','1000','10000','100000']) #numpy.arange does include the 'stop' value\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "name": "python3",
   "language": "python",
   "display_name": "Python 3 (ipykernel)"
  },
  "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}