Cascadia_Rising_2016.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Earthquake parameters (2014 exercise)\n",
    "\n",
    "Location: N 45.73, W 125.12\n",
    "\n",
    "Depth: 2.0 km\n",
    "\n",
    "Prelim Magnitude: M 9.0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#import globalETAS as gep\n",
    "import numpy\n",
    "import scipy\n",
    "import math\n",
    "import os\n",
    "import sys\n",
    "import pylab as plt\n",
    "import datetime as dtm\n",
    "import pytz\n",
    "import multiprocessing as mpp\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "#\n",
    "import matplotlib.dates as mpd\n",
    "from mpl_toolkits.basemap import Basemap as Basemap\n",
    "#from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas\n",
    "#from geographiclib.geodesic import Geodesic as ggp\n",
    "#\n",
    "#\n",
    "import ANSStools as atp\n",
    "import contours2kml\n",
    "import globalETAS as gep\n",
    "import globalETAS_scripts as gesp\n",
    "import eq_params as eqp\n",
    "tzutc=pytz.timezone('UTC')\n",
    "#\n",
    "# define earthquake parameters.\n",
    "eq_lat = 45.74\n",
    "eq_lon = -125.12\n",
    "eq_z = 2.0   #km\n",
    "eq_m = 9.0\n",
    "#\n",
    "#ETAS resolution, etc. parameters:\n",
    "n_processes = mpp.cpu_count()\n",
    "n_contours=25\n",
    "d_lat,d_lon = .1, .1\n",
    "etas_range_factor=25\n",
    "#etas_range_factor=5\n",
    "etas_range_padding=2.25\n",
    "#\n",
    "mainshock = {'lat':eq_lat, 'lon':eq_lon, 'depth': eq_z, 'mag':eq_m}\n",
    "#\n",
    "# and let's have a go at replicating the rupture from the shakemap.\n",
    "rupture_poly = [[-127.,49.5], [-125.5,46.], [-125., 40.], [-124.,40.], [-123., 47.], [-125.0, 49.75]]\n",
    "if rupture_poly[-1]!=rupture_poly[0]: rupture_poly+=[rupture_poly[0]]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## make a quick map that shows our rupture area, and maybe more stuff too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# lons=[-132.5, -115.0], lats=[35.0, 53.]\n",
    "lat_min = 35\n",
    "lat_max=53\n",
    "lon_min=-132.\n",
    "lon_max = -115.\n",
    "#\n",
    "time_mainshock = dtm.datetime(2016,6,7, 8,30,0, tzinfo=pytz.timezone('US/Pacific'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def angle_between_2d(v_seg1, v_seg2, deg_rad='rad'):\n",
    "    #v0 = np.array(p0) - np.array(p1)\n",
    "    #v1 = np.array(p2) - np.array(p1)\n",
    "    #\n",
    "    v_seg_1=numpy.array(v_seg1)\n",
    "    if isinstance(v_seg_1[0], int) or isinstance(v_seg_1[0], float): v_seg_1 = [numpy.zeros(len(v_seg_1)),v_seg_1]\n",
    "    v_seg_2=numpy.array(v_seg2)\n",
    "    if isinstance(v_seg_2[0], int) or isinstance(v_seg_2[0], float): v_seg_2 = [numpy.zeros(len(v_seg_2)),v_seg_2]\n",
    "    #\n",
    "    v0 = v_seg_1[1]-v_seg_1[0]\n",
    "    v1 = v_seg_2[1]-v_seg_2[0]\n",
    "    #\n",
    "    angle = numpy.math.atan2(numpy.linalg.det([v0,v1]), numpy.dot(v0,v1))\n",
    "    #\n",
    "    if deg_rad=='rad':\n",
    "        return angle\n",
    "    else:\n",
    "        return numpy.degrees(angle)\n",
    "#\n",
    "def get_rotation_matrix(v_seg1, v_seg2):\n",
    "    theta = angle_between_2d(v_seg1, v_seg2)\n",
    "    #\n",
    "    return numpy.array([[numpy.cos(theta), -numpy.sin(theta)], \n",
    "                         [numpy.sin(theta),  numpy.cos(theta)]])\n",
    "#\n",
    "#def move_and_rotate_coordinates(coords, )\n",
    "#def vector_to_vector_rotation(vc1 = [1.,0.], vc2 = [0,1], ab=[1,1]):\n",
    "#    # after doing some algebra and then translating to linear albebra, we get this matrix to transform one kn\n",
    "\n",
    "\n",
    "#\n",
    "def rotate_points(XY, theta=0., center=[0.,0.], deg=True):\n",
    "    if deg: theta = numpy.radians(theta)\n",
    "    center = numpy.array(center)\n",
    "    #XY\n",
    "    X,Y = (numpy.array(x) for x in zip(*(XY-center)))    \n",
    "    #\n",
    "    R = numpy.array([[numpy.cos(theta), -numpy.sin(theta)], \n",
    "                         [numpy.sin(theta),  numpy.cos(theta)]])\n",
    "    #\n",
    "    #XY_prime = numpy.dot(R, list(zip(*(XY-center))))\n",
    "    X_prime, Y_prime = numpy.dot(R, (X,Y))\n",
    "    #\n",
    "    #XY_prime.shape=(int(XY_prime.size/2),2)\n",
    "    X_prime += center[0]\n",
    "    Y_prime += center[1]\n",
    "    #\n",
    "    return numpy.array(list(zip(X_prime, Y_prime)))\n",
    "#\n",
    "def draw_basic_map():\n",
    "    plt.figure(0, figsize=(10,10))\n",
    "\n",
    "    bm = Basemap(projection='cyl',llcrnrlat=lat_min,urcrnrlat=lat_max,\\\n",
    "                llcrnrlon=lon_min,urcrnrlon=lon_max, resolution='i')\n",
    "    bm.drawcoastlines()\n",
    "    bm.fillcontinents(color='coral',lake_color='aqua')\n",
    "    # draw parallels and meridians.\n",
    "    bm.drawparallels(numpy.arange(int(lat_min), int(lat_max), 2.))\n",
    "    bm.drawmeridians(numpy.arange(int(lon_min), int(lon_max), 2.))\n",
    "    bm.drawmapboundary(fill_color='aqua')\n",
    "    plt.title(\"Cylindrical Projection\")\n",
    "    #\n",
    "    #print('poly: ', rupture_poly)\n",
    "    bm.plot(*zip(*rupture_poly), latlon=True, marker='', ls='-', lw=2.5)\n",
    "    bm.plot([eq_lon], [eq_lat], marker='*',ls='', ms=18, color='k', zorder=10)\n",
    "    bm.plot([eq_lon], [eq_lat], marker='*',ls='', ms=16, color='r', zorder=11)\n",
    "    #\n",
    "    return bm\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## For the mainshock, use the Tohoku event. Copy the Tohoku sequence, rotate and translate to the Cascadia region. if we really want to be smart about it, we should probably also do a mirror transformation along the rupture axis, since the subduction is in the opposite direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#print('angle: ', angle_between_2d([[1.,1.], [2,1]], [[1,0.], [1,4]]))\n",
    "#print(eqp.tohoku_ETAS_prams)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "tpp = eqp.tohoku_ETAS_prams\n",
    "# get mainshock:\n",
    "tpp['lats']=[30.0, 43.]\n",
    "mc=tpp['mc']\n",
    "mc=4.0\n",
    "tohoku_cat = atp.catfromANSS(lon=tpp['lons'], lat=tpp['lats'], minMag=mc, \n",
    "                             dates0=[dtm.datetime(2005,3,1, tzinfo=tzutc), dtm.datetime(2012,1,1,tzinfo=tzutc)],\n",
    "                             Nmax=None, fout=None, rec_array=True)\n",
    "# we ran the original data sets on approximately 16 May 2016; if we want to add supplementary ETAS, we should\n",
    "# populate the local catalog to that date.\n",
    "local_cat_date = dtm.datetime.now(tzutc)\n",
    "local_cat_date = dtm.datetime(2016,5,16, tzinfo=tzutc)\n",
    "local_cat = atp.catfromANSS(lon=[lon_min, lon_max], lat=[lat_min, lat_max], minMag=mc,\n",
    "                            dates0=[local_cat_date-dtm.timedelta(days=5*365), local_cat_date],\n",
    "                            Nmax=None, fout=None, rec_array=True)\n",
    "#                            \n",
    "max_m = max(tohoku_cat['mag'])\n",
    "j_tohoku, tohoku_ms_rw = [[j,rw] for j,rw in enumerate(tohoku_cat) if rw['mag']==max_m][0]\n",
    "print('tohoku_ms: ', j_tohoku, tohoku_ms_rw)\n",
    "#\n",
    "# translation vector from source to target location\n",
    "dxdy = [mainshock['lon']-tohoku_ms_rw['lon'], mainshock['lat']-tohoku_ms_rw['lat']]\n",
    "print('dxdy: ', dxdy)\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tohoku_cat_prime = tohoku_cat.copy()\n",
    "local_cat_prime = local_cat.copy()\n",
    "#\n",
    "theta = 50\n",
    "#\n",
    "XY_prime = rotate_points(numpy.array(list(zip( (tohoku_cat_prime['lon']), (tohoku_cat_prime['lat'])))), theta, center=[tohoku_ms_rw['lon'], tohoku_ms_rw['lat']])\n",
    "XY_prime += dxdy\n",
    "#\n",
    "# and a symmetry operation:\n",
    "do_symm_opp=False\n",
    "if do_symm_opp:\n",
    "    ms_center = [mainshock['lon'], mainshock['lat']]\n",
    "    XY_prime-=ms_center\n",
    "    XY_prime = rotate_points(XY_prime, -10, center=[0,0])\n",
    "    #print('xy: ', XY_prime[0:5])\n",
    "    Ms = [[-1., 0], [0.,1.]]\n",
    "    XY_prime = numpy.array([numpy.dot(Ms, rw) for rw in XY_prime])\n",
    "    XY_prime = rotate_points(XY_prime, 10, center=(0,0)) + ms_center\n",
    "    #XY_prime = numpy.dot(Ms,)+ms_center\n",
    "#\n",
    "tohoku_cat_prime['lon'], tohoku_cat_prime['lat'] = numpy.array(list(zip(*XY_prime)))\n",
    "# now, get the tohoku datetime. move up all the tohoku events to (around) the 7 june cascadia event time.\n",
    "# let's also pull some events from today around the region and construct a composite catalog.\n",
    "t_tohoku = dtm.datetime(*(x for x in tohoku_ms_rw['event_date'].tolist().timetuple()[0:6]), tzinfo=tzutc)\n",
    "tohoku_time_delta = time_mainshock-t_tohoku\n",
    "print('tohoku time_delta: ', tohoku_time_delta)\n",
    "for j,dt in enumerate(tohoku_cat_prime['event_date']):\n",
    "    tohoku_cat_prime['event_date'][j] = numpy.datetime64(dt.tolist() + tohoku_time_delta)\n",
    "    tohoku_cat_prime['event_date_float'][j] = mpd.date2num(dt.tolist() + tohoku_time_delta)\n",
    "#\n",
    "# and bump up the local catalog to the event date.\n",
    "#local_time_delta = time_mainshock-dtm.datetime.now(tzutc)\n",
    "local_time_delta = time_mainshock-local_cat_date\n",
    "#print('ltd: ', local_time_delta.days)\n",
    "for j,dt in enumerate(local_cat_prime['event_date']):\n",
    "    local_cat_prime['event_date'][j] = numpy.datetime64(dt.tolist() + local_time_delta)\n",
    "    local_cat_prime['event_date_float'][j] = mpd.date2num(dt.tolist() + local_time_delta)\n",
    "#\n",
    "bm = draw_basic_map()\n",
    "bm.plot(tohoku_cat_prime['lon'], tohoku_cat_prime['lat'], marker='.', color='b', ls='', alpha=.5,zorder=6)\n",
    "bm.plot(local_cat_prime['lon'], local_cat_prime['lat'], marker='.', color='g', ls='', alpha=.5, zorder=7)\n",
    "#\n",
    "# Take a look at what we've got so far:\n",
    "#print('dtype', tohoku_cat_prime.dtype.names)\n",
    "plt.figure()\n",
    "plt.plot(tohoku_cat_prime['event_date_float'], tohoku_cat_prime['mag'])\n",
    "plt.plot(local_cat_prime['event_date_float'], local_cat_prime['mag'])\n",
    "plt.plot([mpd.date2num(dtm.datetime.now(tzutc))], [6], marker='*', color='r', ls='', ms=15, zorder=6)\n",
    "plt.plot([mpd.date2num(time_mainshock)], [6], marker='*', color='m', ls='', ms=15, zorder=6)\n",
    "#plt.plot(tohoku_cat_prime['event_date'], tohoku_cat_prime['mag'])\n",
    "#plt.plot(local_cat_prime['event_date'], local_cat_prime['mag'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# now, make an ETAS-ready catalog:\n",
    "# TODO: set lats=[], lons=[] to get better definition of map boundaries.\n",
    "# now, combine the Tohoku and Local catalogs and process for ETAS.\n",
    "combined_catalog = numpy.core.records.fromarrays(zip(*tohoku_cat_prime.tolist() + local_cat_prime.tolist()), dtype=local_cat_prime.dtype)\n",
    "combined_catalog.sort(order='event_date_float')\n",
    "etas_cat = gep.make_ETAS_catalog_mpp(incat=combined_catalog)    # and otherwise, the defaults should be ok.\n",
    "#\n",
    "print('ETAS catalog processed.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(tohoku_ms_rw, tohoku_ms_rw.dtype.names)\n",
    "print(type(tohoku_ms_rw['event_date'].tolist()))\n",
    "print(time_mainshock)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bm = draw_basic_map()\n",
    "bm.plot(etas_cat['lon'], etas_cat['lat'], ls='', marker='.', ms=2, color='r')\n",
    "    #\n",
    "#draw_map()\n",
    "    #\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, load up an ETAS object. maybe we can develop some sample seed catalogs from which to base earthquake orientation, and maybe create some examples with Nepal-like not co-located clustering. maybe grab the Nepal sequence and rotate/scale it to the Cascadia fault?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# now, generate ETAS from our new catalog\n",
    "# resolution and other parameters are set in opening cell.\n",
    "#\n",
    "t_now = time_mainshock+dtm.timedelta(hours=1)\n",
    "#\n",
    "#etas = gep.ETAS_mpp_handler_xyz(n_processes = n_processes, *args, **kwargs)\n",
    "# catalog=None, lats=None, lons=None, mc=2.5, mc_etas=None, d_lon=.1, d_lat=.1, bin_lon0=0., bin_lat0=0., \n",
    "#       etas_range_factor=10.0, etas_range_padding=.25, etas_fit_factor=1.0, t_0=dtm.datetime(1990,1,1, tzinfo=tz_utc),\n",
    "#       t_now=dtm.datetime.now(tzutc), transform_type='equal_area', transform_ratio_max=2.5, cat_len=5.*365., \n",
    "#      calc_etas=True, n_contours=15, etas_cat_range=None, etas_xyz_range=None, p_cat=1.1, q_cat=1.5, p_etas=None,**kwargs):\n",
    "#\n",
    "# there may be a problem with pickling the catalog objects when we provide a catalog.\n",
    "# TODO: set t_now to define the forcast date. we'll want to make a list of t_now and create aftershock export products\n",
    "# for each t_now.\n",
    "etas = gep.ETAS_mpp_handler_xyz(n_processes = n_processes, t_now=t_now, catalog=etas_cat, lons=[lon_min, lon_max], lats=[lat_min, lat_max], mc=None,\n",
    "                               d_lat=d_lat, d_lon=d_lon, etas_range_factor=etas_range_factor,\n",
    "                                etas_range_padding=etas_range_padding, etas_fit_factor=1.0,\n",
    "                                n_contours=n_contours)\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fignum=0\n",
    "plt.figure()\n",
    "figsize=(10,10)\n",
    "color_map='hot'\n",
    "bm2=etas.make_etas_contour_map(n_contours=n_contours, fignum=fignum, fig_size=figsize, kml_contours_bottom=0., \n",
    "                           kml_contours_top=1.0, alpha=.6, map_resolution='i', map_projection='cyl', map_cmap=color_map)\n",
    "#                           ,\n",
    "#                           contour_fig_file=png_file, contour_kml_file=kml_file, alpha_kml=.5, refresh_etas=False)\n",
    "plt.title('CR2016: Simulated ETAS, %s' % t_now)\n",
    "etas.cm.plot(*zip(*rupture_poly), latlon=True, marker='', ls='-', lw=2.5) \n",
    "plt.savefig('/home/myoder/Dropbox/Research/CascadiaRising2016/ETAS_CR2016_%s_hot.png' % str(t_now))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### ... and we should be able to use the globalETAS scripts to produce a set of ETAS, including export files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# etas_outputs(n_processes = None, output_path='etas_outputs/', kml_file='etas_kml.kml', png_file=None,\n",
    "# xyz_file=None, fignum=0, color_map='jet', *args, **kwargs)\n",
    "time_deltas = [-dtm.timedelta(days=10), -dtm.timedelta(days=1), dtm.timedelta(hours=2), dtm.timedelta(hours=8), dtm.timedelta(hours=16),\n",
    "               dtm.timedelta(days=1), dtm.timedelta(hours=36), dtm.timedelta(days=2), dtm.timedelta(days=2,hours=12),\n",
    "               dtm.timedelta(days=3), dtm.timedelta(days=3,hours=12), dtm.timedelta(days=4)]\n",
    "#\n",
    "# some near-field supplementary runs:\n",
    "time_deltas = [dtm.timedelta(minutes=x) for x in (15, 30, 60, 90)]\n",
    "\n",
    "output_path='%s/Dropbox/Research/CascadiaRising2016/outputs_hires' % os.getenv('HOME')\n",
    "#etas = gep.ETAS_mpp_handler_xyz(n_processes = n_processes, t_now=t_now, catalog=etas_cat, lons=[lon_min, lon_max], lats=[lat_min, lat_max], mc=None,\n",
    "#                               d_lat=d_lat, d_lon=d_lon, etas_range_factor=etas_range_factor,\n",
    "#                                etas_range_padding=etas_range_padding, etas_fit_factor=1.0,\n",
    "#                                n_contours=n_contours)\n",
    "kwds = {'n_processes':n_processes, 't_now':t_now, 'catalog':etas_cat, 'lons':[lon_min, lon_max], 'lats':[lat_min, lat_max],\n",
    "        'mc':None, 'd_lat':d_lat, 'd_lon':d_lon, 'etas_range_factor':etas_range_factor, \n",
    "        'etas_range_padding':etas_range_padding, 'etas_fit_factor':1.0, 'n_contours':n_contours}\n",
    "#\n",
    "for td in time_deltas:\n",
    "    t_now = time_mainshock+td\n",
    "    #\n",
    "    kml_file = 'Cascadia_Rising_Exercise_2016_etas_%s_kml.kml' % str(t_now)\n",
    "    kwds.update({'t_now':t_now})\n",
    "    #\n",
    "    etas_output = gesp.etas_outputs(output_path=output_path, kml_file=kml_file, png_file=None,\n",
    "                                   xyz_file=None, fignum=0, color_map='jet', **kwds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}