diff --git a/notebooks/geoclaw/Okada.ipynb b/notebooks/geoclaw/Okada.ipynb new file mode 100644 index 00000000..ba3dd32f --- /dev/null +++ b/notebooks/geoclaw/Okada.ipynb @@ -0,0 +1,365 @@ +{ + "metadata": { + "name": "", + "signature": "sha256:3ed22d5adb293714950468a32da406b8bea77be84677ba08703fab8a9606d0e4" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# The Okada model for computing seafloor deformation\n", + "\n", + "This notebook is part of the GeoClaw documentation, see \n", + " - For general Clawpack documentation,\n", + " - for more documentation on the use of the Okada model in GeoClaw. \n", + " - for information on obtaining the *apps* repository. The source for this notebook is maintained in `$CLAW/apps/notebooks/geoclaw/Okada.ipynb`.\n", + " - : In GeoClaw the Okada model is implemented in `dtopotools.SubFault.okada`, \n", + " \n", + "## Contents\n", + "\n", + "- Examples\n", + "- Typical subduction zone thrust event\n", + "- Varying the strike\n", + "- Varying the dip\n", + "- Varying the rake" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Introduction\n", + "\n", + "The \"Okada model\" takes as input the slip on a rectangular patch of a fault deep in the earth and produces the resulting deformation of the earth's surface. For tsunami modeling, we are particularly interested in the vertical motion of the seafloor above the fault. \n", + "\n", + "The Okada model makes many assumptions that only approximate reality, in particular:\n", + " - The earth is modeled as a perfectly elastic half space. The slip is uniform on a rectangular patch and a Greens' function solution to the half space problem is integrated to obtain the vertical deformation at the surface of the half plane. The seafloor topography is ignored in computing the Okada deformation (which is then generally added to the real topography to get the modified topography during the earthquake).\n", + " - The half space is assumed to be isotropic with uniform elastic moduli. The Poisson ratio is generally taken to be 0.25 and the shear modulus (or rigidity) is constant, often taken to be $4\\times 10^{11}$ dyne/cm$^2 = 40$ GPa.\n", + " \n", + "Complex earthquake fault surfaces are generally approximated by a subdividing into a set of rectangular subfaults, each of which might have different orientation and slip properties. The Okada model is applied to each and then the resulting surface deformations are summed. The Okada model is linear, so that if a single planar fault is split into subfaults the resulting surface deformation should be independent of the number of pieces its split into.\n", + "\n", + "The GeoClaw *dtopotools.Fault* class provides tools for specifying a fault as a collection of *dtopotools.SubFault* objects, and an object of this class has a method *create_dtopography* that applies the Okada model and returns a *dtopotools.DTopography* object expressing the surface deformation. In this notebook we illustrate the Okada model by working with a fault that consists of a single subfault.\n", + "\n", + "A subfault is specified via the following standard parameters:\n", + " - *length* and *width* of the fault plane (specified in meters below),\n", + " - *latitude* and *longitude* of some point on the fault plane, typically\n", + " either the centroid or the center of the top (shallowest edge),\n", + " - *depth* of the specified point below the sea floor (in meters below),\n", + " - *strike*, the orientation of the top edge, measured in degrees\n", + " clockwise from North, between 0 and 360. The fault plane dips downward\n", + " to the right when moving along the top edge in the strike direction.\n", + " - *dip*, angle at which the plane dips downward from the top edge, a\n", + " positive angle between 0 and 90 degrees.\n", + " - *rake*, the angle in the fault plane in which the slip occurs,\n", + " measured in degrees counterclockwise from the strike direction.\n", + " Between -180 and 180.\n", + " - *slip > 0*, the distance (measured in meters below) the hanging-wall block moves\n", + " relative to the footwall block, in the direction specified by the rake.\n", + " The \"hanging-wall block\" is the one above the dipping fault plane (or to the\n", + " right if you move along the fault in the strike direction).\n", + " \n", + "\n", + "For other descriptions and illustrations, see e.g.\n", + " - \n", + " - \n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Import modules and define some utility functions used below..." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "%pylab inline" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from clawpack.geoclaw import dtopotools\n", + "from clawpack.visclaw.JSAnimation import IPython_display\n", + "import clawpack.visclaw.JSAnimation.JSAnimation_frametools as J" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def set_fault(strike, dip, rake, depth):\n", + " \"\"\"\n", + " Set the subfault parameters.\n", + " Most are fixed for the examples below, \n", + " and only the strike, dip, and rake will be varied.\n", + " \"\"\"\n", + " subfault = dtopotools.SubFault()\n", + " subfault.strike = strike\n", + " subfault.dip = dip\n", + " subfault.rake = rake\n", + " subfault.length = 100.e3\n", + " subfault.width = 50.e3\n", + " subfault.depth = depth\n", + " subfault.slip = 1.\n", + " subfault.longitude = 0.\n", + " subfault.latitude = 0.\n", + " subfault.coordinate_specification = \"top center\"\n", + "\n", + " fault = dtopotools.Fault()\n", + " fault.subfaults = [subfault]\n", + " return fault, subfault\n", + "\n", + "# Create a sample fault and print out some information about it...\n", + "fault, subfault = set_fault(0,0,0,5e3)\n", + "print \"This sample fault has %s meter of slip over a %s by %s km patch\" \\\n", + " % (subfault.slip,subfault.length/1e3,subfault.width/1e3)\n", + "print \"With shear modulus %4.1e Pa the seismic moment is %4.1e\" % (subfault.mu, subfault.Mo())\n", + "print \" corresponding to an earthquake with moment magnitude %s\" % fault.Mw()\n", + "print \"The depth at the top edge of the fault plane is %s km\" % (subfault.depth/1e3)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def plot_okada(strike, dip, rake, depth, verbose=False):\n", + " \"\"\"\n", + " Make 3 plots to illustrate the Okada solution.\n", + " \"\"\"\n", + "\n", + " fault,subfault = set_fault(strike, dip, rake, depth)\n", + " ax1 = subplot(2,2,1)\n", + " ax2 = subplot(2,2,2)\n", + " ax3 = subplot(2,2,3)\n", + " ax4 = subplot(2,2,4)\n", + "\n", + " # Subfault projection on surface on ax1:\n", + " ax = fault.plot_subfaults(axes=ax1, plot_rake=True, xylim=[-.5,1.5, -1,1])\n", + " text(0.6,0.8,\"Strike = %5.1f\" % strike, fontsize=12)\n", + " text(0.6,0.6,\"Dip = %5.1f\" % dip, fontsize=12)\n", + " text(0.6,0.4,\"Rake = %5.1f\" % rake, fontsize=12)\n", + " text(0.6,0.2,\"Depth = %5.1f km\" % (depth/1e3), fontsize=12)\n", + " ax1.set_ylabel('latitude (degrees)')\n", + "\n", + " # Depth profile on ax3:\n", + " z_top = -subfault.fault_plane_centers[0][2] / 1e3 # convert to km\n", + " z_bottom = -subfault.fault_plane_centers[2][2] / 1e3 # convert to km\n", + " ax3.plot([0,cos(subfault.dip*pi/180.)*subfault.width/1.e3], [z_top, z_bottom])\n", + " ax3.set_xlim(-50,150)\n", + " ax3.set_ylim(-55,0)\n", + " ax3.set_xlabel('distance orthogonal to strike')\n", + " ax3.set_ylabel('depth (km)')\n", + " ax3.set_title('Depth profile')\n", + " \n", + " \n", + " # Grid to use for evaluating and plotting dz\n", + " x = numpy.linspace(-0.5, 1., 101)\n", + " y = numpy.linspace(-1., 1., 101)\n", + " times = [1.]\n", + "\n", + " # color map of deformation dz on ax2:\n", + " fault.create_dtopography(x,y,times,verbose=verbose)\n", + " dtopo = fault.dtopo\n", + " dtopo.plot_dz_colors(t=1., axes=ax2)\n", + " \n", + " # transect of dz on ax4:\n", + " dz = dtopo.DZ[-1,50,:]\n", + " ax4.plot(x,dz)\n", + " ax4.set_ylim(-0.5,0.5)\n", + " ax4.set_title('Transect of dz along y=0')\n", + " ax4.set_xlabel('Longitude (degrees)')\n", + " ax4.set_ylabel('Seafloor deformation (m)')\n", + " " + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "## Examples\n", + "\n", + "
\n", + "### Typical subduction zone thrust event\n", + "\n", + "For a subduction zone earthquake, the rake is generally close to 90 degrees. In the first example below we take it to be 80 degrees -- note that the green line on the fault plane extending from the centroid in the rake direction is 80 degrees counterclockwise from north (the strike direction, since *strike = 0*). \n", + "\n", + "Note that the fault plane plot shows the projection of the fault plane on the flat surface. It dips down to the right as seen in the depth profile. The *dip* is 10 degrees and the fault extends from a depth of 5 km at the up-dip edge to about 14 km at the downdip edge. Note that above the width of the subfault is set to 50 km and that the \"Fault planes\" plot axes are in degrees. This fault is located at the equator where 1 degree is approximate 111 km in both directions.\n", + "\n", + "The rock above the fault plane (the hanging block) is moving to the left relative to the lower block (as shown by the green line) and hence the rock above the fault will be compressed to the left and bulge upwards. The rock to the right is under tension and the surface dips downwards as a result. This can be seen in the \"seafloor deformation\" plot, where the colors indicate meters of vertical motion. \n", + "\n", + "The *slip* on this subfault was set to 1 meter above. Since the Okada model is linear, changing the slip to some other value would simply multiply the deformation by the same factor everywhere." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "fig=figure(figsize=(10,10))\n", + "plot_okada(strike=0, dip=10, rake=80, depth=5e3, verbose=True)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "## Varying the strike\n", + "\n", + "Changing the strike simply rotates the fault and the resulting deformation. The animation below illustrates this...\n" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "plotdir = 'okada_strike_plots'\n", + "J.make_plotdir(plotdir, clobber=True)\n", + "fig=figure(figsize=(10,10))\n", + "for k,strike in enumerate(linspace(0,340,18)):\n", + " clf()\n", + " plot_okada(strike, 10., 80., 5e3)\n", + " J.save_frame(k, plotdir=plotdir, verbose=False)\n", + " \n", + "clf()\n", + "anim = J.make_anim(plotdir)\n", + "anim" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "## Varying the dip\n", + "\n", + "When *dip=0* the fault plane is horizontal. If the *rake* is near 90 degrees we expect compression and uplift to the left and subsidence to the right, as illustrated above. \n", + "\n", + "When the *dip=90* the fault plane is vertical with the hanging-wall block to the right moving upward and the footwall block to the left moving downward. We then expect to see uplift on the right and subsidence on the left, opposite to what is seen when *dip=0*. \n", + "\n", + "Note how this transition takes place as the dip is changed..." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "plotdir = 'okada_dip_plots'\n", + "J.make_plotdir(plotdir, clobber=True)\n", + "fig=figure(figsize=(10,10))\n", + "for k,dip in enumerate(linspace(0,90,10)):\n", + " clf()\n", + " plot_okada(0., dip, 90., 5e3)\n", + " J.save_frame(k, plotdir=plotdir, verbose=False)\n", + " \n", + "clf()\n", + "anim = J.make_anim(plotdir)\n", + "anim" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "## Varying the rake\n", + "\n", + "If we fix the *dip* at 10 degrees and vary the direction of slip on the fault plane (the *rake*), we get the following patterns:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "plotdir = 'okada_rake_plots'\n", + "J.make_plotdir(plotdir, clobber=True)\n", + "fig=figure(figsize=(9,9))\n", + "for k,rake in enumerate(linspace(-90,90,19)):\n", + " clf()\n", + " plot_okada(0., 10., rake, 5e3)\n", + " J.save_frame(k, plotdir=plotdir, verbose=False)\n", + " \n", + "clf()\n", + "anim = J.make_anim(plotdir)\n", + "anim" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "## Varying the depth\n", + "\n", + "If the fault surface is near the surface the surface deformation will be more concentrated near the fault plane than if the fault is deeper, as the next animation illustrates.\n", + "\n", + "Note that the grid used in this example for evaluting the seafloor deformation *dz* does not extend out far enough to capture all of the deformation, particularly for deeper faults!" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "plotdir = 'okada_depth_plots'\n", + "J.make_plotdir(plotdir, clobber=True)\n", + "fig=figure(figsize=(9,9))\n", + "for k,depth in enumerate(arange(0,40,2)*1e3):\n", + " clf()\n", + " plot_okada(0., 10., 80, depth)\n", + " J.save_frame(k, plotdir=plotdir, verbose=False)\n", + " \n", + "clf()\n", + "anim = J.make_anim(plotdir)\n", + "anim" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [], + "language": "python", + "metadata": {}, + "outputs": [] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file diff --git a/notebooks/geoclaw/data/alaska1964.csv b/notebooks/geoclaw/data/alaska1964.csv new file mode 100644 index 00000000..96396161 --- /dev/null +++ b/notebooks/geoclaw/data/alaska1964.csv @@ -0,0 +1,13 @@ +Longitude,Latitude,Depth(km),Length,Width,Strike,Dip,Rake,Slip,UnitSrc +204.89500,55.97000,17.94,100.0,50.0,236.000,15.000,90.0,19.5000,acsza31 +205.34000,55.59800,5.00,100.0,50.0,236.000,15.000,90.0,19.5000,acszb31 +206.20800,56.47300,17.94,100.0,50.0,236.000,15.000,90.0,24.5000,acsza32 +206.65800,56.10000,5.00,100.0,50.0,236.000,15.000,90.0,24.5000,acszb32 +207.53700,56.97500,17.94,100.0,50.0,236.000,15.000,90.0,29.5000,acsza33 +207.99300,56.60300,5.00,100.0,50.0,236.000,15.000,90.0,29.5000,acszb33 +208.93710,57.51240,17.94,100.0,50.0,236.000,15.000,90.0,34.5000,acsza34 +209.40000,57.14000,5.00,100.0,50.0,236.000,15.000,90.0,34.5000,acszb34 +210.25970,58.04410,17.94,100.0,50.0,230.000,15.000,90.0,34.5000,acsza35 +210.80000,57.70000,5.00,100.0,50.0,230.000,15.000,90.0,34.5000,acszb35 +211.32490,58.65650,17.94,100.0,50.0,218.000,15.000,90.0,34.5000,acsza36 +212.00000,58.38000,5.00,100.0,50.0,218.000,15.000,90.0,34.5000,acszb36 diff --git a/notebooks/geoclaw/data/timedep.csv b/notebooks/geoclaw/data/timedep.csv new file mode 100644 index 00000000..c6ee2d70 --- /dev/null +++ b/notebooks/geoclaw/data/timedep.csv @@ -0,0 +1,13 @@ +Longitude,Latitude,Depth(km),Length,Width,Strike,Dip,Rake,Slip,rupture_time,rise_time +204.89500,55.97000,17.94,100.0,50.0,236.000,15.000,90.0,19.5000,0,1 +205.34000,55.59800,5.00,100.0,50.0,236.000,15.000,90.0,19.5000,1,1 +206.20800,56.47300,17.94,100.0,50.0,236.000,15.000,90.0,24.5000,2,1 +206.65800,56.10000,5.00,100.0,50.0,236.000,15.000,90.0,24.5000,3,1 +207.53700,56.97500,17.94,100.0,50.0,236.000,15.000,90.0,29.5000,4,1 +207.99300,56.60300,5.00,100.0,50.0,236.000,15.000,90.0,29.5000,5,1 +208.93710,57.51240,17.94,100.0,50.0,236.000,15.000,90.0,34.5000,6,1 +209.40000,57.14000,5.00,100.0,50.0,236.000,15.000,90.0,34.5000,7,1 +210.25970,58.04410,17.94,100.0,50.0,230.000,15.000,90.0,34.5000,8,1 +210.80000,57.70000,5.00,100.0,50.0,230.000,15.000,90.0,34.5000,9,1 +211.32490,58.65650,17.94,100.0,50.0,218.000,15.000,90.0,34.5000,10,1 +212.00000,58.38000,5.00,100.0,50.0,218.000,15.000,90.0,34.5000,11,1 diff --git a/notebooks/geoclaw/dtopotools_examples.ipynb b/notebooks/geoclaw/dtopotools_examples.ipynb new file mode 100644 index 00000000..4bf04fd0 --- /dev/null +++ b/notebooks/geoclaw/dtopotools_examples.ipynb @@ -0,0 +1,567 @@ +{ + "metadata": { + "name": "", + "signature": "sha256:4127453e482935e99ae339e2cb581b846947bdc3026ea299194089edd8ee64fb" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# GeoClaw dtopotools examples\n", + "\n", + "This notebook contains some examples of working with the *clawpack.geoclaw.dtopotools* module. These tools facilitate creating and manipulating the *dtopo* files that are required as GeoClaw input to specify moving topography, in particular the seafloor motion due to a subduction zone earthquake. All the examples in this notebook are of this nature: slip on a fault surface (or in geneneral a collection of rectangular planar subfaults making up the fault) is converted to seafloor deformation by applying the Okada model. \n", + "\n", + "See for general documentation of these tools and their use in the context of GeoClaw, and for more information on the Clawpack suite of software." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Contents\n", + "\n", + " - Reading a csv file specifying the subfaults\n", + " - Plot the subfault geometry and slip\n", + " - Plot the seafloor displacement\n", + " - Create the dtopo file needed by GeoClaw\n", + " - Read in an existing dtopo file\n", + " - Create a fault using the NOAA SIFT database directly\n", + " - Time-dependent fault rupture\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Setup notebook and environment:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "%pylab inline" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from clawpack.geoclaw import dtopotools, topotools\n", + "import os\n", + "CLAW = os.environ['CLAW']\n", + "datadir = os.path.join(CLAW,'geoclaw','scratch') # directory for some sample data files" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "\n", + "## Reading a csv file specifying the subfaults\n", + "\n", + "Read in a csv file where the first line labels the columns and the remaining lines give the parameters for each subfault. This is a model of the 1964 Alaska earthquake using 12 subfaults from the NOAA SIFT database of unit sources. \n", + "\n", + "See and\n", + "[gica2937.pdf](http://www.pmel.noaa.gov/pubs/PDF/gica2937/gica2937.pdf).\n", + "\n", + "The full database is available at and a later example in this notebook illustrates using the *dtopotools.SiftFault* class to access this directly. But this example shows how a csv file can be specified. The first line labels the columns with standard attributes needed to specify each subfault. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here is the raw file:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "subfault_fname = 'data/alaska1964.csv'\n", + "print open(subfault_fname).read()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that the last column has the name of the unit source from the NOAA SIFT database and will be ignored. The other columns have labels that are expected. Specify the units used for each and the coordinate system for the (longitude,latitude) and depth when reading the file in:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "input_units = {\"length\":\"km\", \"width\":\"km\", \"depth\":\"km\", \"slip\":\"m\"}\n", + "fault = dtopotools.CSVFault()\n", + "fault.read(subfault_fname, input_units=input_units, coordinate_specification=\"noaa sift\")" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You could skip directly to Create the dtopo file needed by GeoClaw, but first we'll explore this fault model a bit..." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Seismic moment `Mo` and Moment Magnitude `Mw` can be computed using the `Fault` class methods:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "print \"The seismic moment is %g N-m\" % fault.Mo()\n", + "print \"The Moment magnitude is %g\" % fault.Mw()\n", + "print \" (Assuming the rigidity mu of all subfaults is the default value %g Pa)\"\\\n", + " % fault.subfaults[0].mu" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "### Plot the subfault geometry and slip\n", + "\n", + "First with a vector showing the rake (direction of slip):" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "fault.plot_subfaults(plot_rake=True)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot again showing the magnitude of slip on each subfault:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "fault.plot_subfaults(slip_color=True)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the coastline by downloading a file containing the x-y coordinates of shorelines based on etopo 4 minute data:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from clawpack.geoclaw.util import get_remote_file\n", + "filename = 'pacific_shorelines_east_4min.npy'\n", + "url = 'http://www.geoclaw.org/topo/' + filename\n", + "get_remote_file(url=url, output_dir=datadir, force=True, verbose=True)\n", + "shorelines_file = os.path.join(datadir, filename)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "fault.plot_subfaults()\n", + "shore = load(shorelines_file)\n", + "plot(shore[:,0], shore[:,1], 'g')\n", + "axis([170,240,40,65])\n", + "gca().set_aspect(1./cos(55*pi/180.))" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "### Create the dtopo file needed by GeoClaw:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Generate a dtopo grid by applying the Okada model to each subfault. First we need to specify the region over which the seafloor deformation should be specified and the resolution of the grid:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "xlower = 203.5\n", + "xupper = 214. # approximate - adjusted below\n", + "ylower = 54.5\n", + "yupper = 60. # approximate - adjusted below\n", + "\n", + "# dtopo parameters:\n", + "points_per_degree = 60 # 1 minute resolution\n", + "dx = 1./points_per_degree\n", + "mx = int((xupper - xlower)/dx + 1)\n", + "xupper = xlower + (mx-1)*dx\n", + "my = int((yupper - ylower)/dx + 1)\n", + "yupper = ylower + (my-1)*dx\n", + "\n", + "x = linspace(xlower,xupper,mx)\n", + "y = linspace(ylower,yupper,my)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now apply Okada to create the static deformation at a single time $t = 1$ second:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "dtopo = fault.create_dtopography(x,y,times=[1.], verbose=True)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can save this deformation as a `dtopo` file with various `dtopo_type`s recognized by GeoClaw. The most compact is `dtopo_type==3`, which specifies a header followed by all the `dz` data:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "dtopo_fname = os.path.join(datadir, 'alaska1964.tt3')\n", + "dtopo.write(dtopo_fname, dtopo_type=3)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Read the file in and print just the header. (The remaining lines contain all the data.)" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "lines = open(dtopo_fname).readlines()\n", + "for k in range(9):\n", + " print lines[k][:-1]" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "### Plot the seafloor displacement" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "dtopo.plot_dZ_colors(t=1)\n", + "\n", + "# add shoreline and Anchorage for orientation:\n", + "plot(shore[:,0], shore[:,1], 'g')\n", + "plot([210.1],[61.2],'ro')\n", + "text(210.3,61.5,'Anchorage')\n", + "axis([xlower,xupper,ylower,62])" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "\n", + "## Read in an existing dtopo file:\n", + "\n", + "An existing dtopo file can be read in for plotting purposes or to further manipulate it. To illustrate, we read in the file we just created..." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "dtopo2 = dtopotools.DTopography()\n", + "dtopo2.read(dtopo_fname, dtopo_type=3)\n", + "\n", + "# Check that this data looks right:\n", + "assert len(dtopo2.x) == 631, \"*** length of x is wrong\"\n", + "assert len(dtopo2.y) == 331, \"*** length of y is wrong\"\n", + "dz_max = abs(dtopo2.dZ).max()\n", + "assert abs(dz_max - 15.368266081250006) < 1e-5, \"*** dz_max is wrong: %g\" % dz_max\n", + "print \"Looks ok\"" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "\n", + "## Create a fault using the NOAA SIFT database directly\n", + "\n", + "The example above showed how to read a csv file with arbitrary columns. Since this particular fault is actually specified in terms of the NOAA SIFT unit sources, another option to create the same fault is to use the *dtopotools.SiftFault* class, which takes as an argument a dictionary *sift_slip* specifying the unit sources to be used and the slip on each. \n", + "\n", + "To illustrate, here we specify only the northern two and southern two unit sources from the example above, but with a dictionary of 12 unit sources we could recreate the full fault." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "sift_slips = {'acsza31':19.5, 'acszb31':19.5, 'acsza36':34.5, 'acszb36':34.5}\n", + "f2 = dtopotools.SiftFault(sift_slips)\n", + "f2.plot_subfaults(slip_color=True)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "\n", + "## Time-dependent fault rupture\n", + "\n", + "This is synthetic data -- the same set of faults as before but with rupture occuring from south to north. Note that new columns `rupture_time` and `rise_time` have been added." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "subfault_fname = 'data/timedep.csv'\n", + "print open(subfault_fname).read()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Read in this csv file and set the parameters for the desired dtopo file:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "input_units = {\"length\":\"km\", \"width\":\"km\", \"depth\":\"km\", \n", + " \"slip\":\"m\", \"mu\":\"dyne/cm^2\"}\n", + "\n", + "fault = dtopotools.CSVFault()\n", + "fault.read(subfault_fname, input_units=input_units, \n", + " coordinate_specification=\"noaa sift\")\n", + "fault.rupture_type = 'dynamic'\n", + "\n", + "print \"%s subfaults read in \" % len(fault.subfaults)\n", + "\n", + "xlower = 203.5\n", + "xupper = 214.\n", + "ylower = 54.5\n", + "yupper = 60.\n", + "xylim = [xlower,xupper,ylower,yupper]\n", + "\n", + "# dtopo parameters:\n", + "points_per_degree = 60 # 1 minute resolution\n", + "mx = int((xupper - xlower)*points_per_degree + 1)\n", + "my = int((yupper - ylower)*points_per_degree + 1)\n", + "x = linspace(xlower,xupper,mx)\n", + "y = linspace(ylower,yupper,my)\n", + "print \"dZ arrays will have shape %s by %s\" % (len(y),len(x))\n" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Specify the desired output times, in this case 6 times between $t=0$ and $t=13$, when the rupture motion is complete. Then apply the Okada model at each time (to the set of all subfaults with slips computed at each time based on the `rupture_time` and `rise_time` for that subfault). Then `dtopo_dz_list` will be a list of 6 deformation arrays:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "times = linspace(0,13,6)\n", + "dtopo = fault.create_dtopography(x,y,times=times)\n", + "print \"Created array of dtopo deformations with shape\", dtopo.dZ.shape" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Create an animation of the fault rupture\n", + "\n", + "First define a function that will show time-dependent slip on the fault planes and resulting time-dependent seafloor motion in two side-by-side axes in the same figure:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# for setting color scale:\n", + "print \"maximum abs(dz) over the full rupture time:\", \\\n", + " abs(dtopo.dZ).max()\n", + "\n", + "# Read in shoreline segments to give context in plots:\n", + "#shoreline_fname = os.path.join(datadir,'tohoku_shoreline_1min.npy')\n", + "#shoreline_xy = load(shoreline_fname)\n", + "shoreline_xy = shore\n", + "\n", + "# Incorporate this function in dtopotools to replace animate_dz_colors?\n", + "def plot_subfaults_dZ(t, fig):\n", + " fig.clf()\n", + " ax1 = fig.add_subplot(121)\n", + " ax2 = fig.add_subplot(122)\n", + " fault.plot_subfaults(axes=ax1, slip_color=True, \n", + " slip_time=t, xylim=xylim)\n", + " dtopo.plot_dZ_colors(axes=ax2, t=t, cmax_dZ=dz_max)\n", + " ax1.plot(shoreline_xy[:,0],shoreline_xy[:,1],'g')\n", + " ax2.plot(shoreline_xy[:,0],shoreline_xy[:,1],'g')\n", + " ax2.set_xlim(xlower,xupper)\n", + " ax2.set_ylim(ylower,yupper)\n", + " return fig\n" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The loop below makes each plot and then saves it as a png file:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from clawpack.visclaw.JSAnimation import IPython_display\n", + "import clawpack.visclaw.JSAnimation.JSAnimation_frametools as J\n", + "plotdir = '_plots'\n", + "J.make_plotdir(plotdir, clobber=True)\n", + "fig = plt.figure(figsize=(12,5))\n", + "\n", + "for k,t in enumerate(times):\n", + " plot_subfaults_dZ(t,fig)\n", + " J.save_frame(k, verbose=True)\n", + " \n", + "print \"Final frame will be displayed below:\"" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Combine the png files into an animation and display it:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "anim = J.make_anim(plotdir)\n", + "anim" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [], + "language": "python", + "metadata": {}, + "outputs": [] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file diff --git a/notebooks/geoclaw/topotools_examples.ipynb b/notebooks/geoclaw/topotools_examples.ipynb new file mode 100644 index 00000000..b01a7f59 --- /dev/null +++ b/notebooks/geoclaw/topotools_examples.ipynb @@ -0,0 +1,293 @@ +{ + "metadata": { + "name": "", + "signature": "sha256:9d29b7123140c2485c9a142657795e0267923d24de47f7f36b580c70e38c584e" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# GeoClaw topotools examples\n", + "\n", + "This notebook contains some examples of working with the *clawpack.geoclaw.topotools* module. These tools facilitate creating and manipulating the *topo* files that are required as GeoClaw input to specify topography and bathymetry (underwater topography).\n", + "\n", + "See for general documentation of these tools and their use in the context of GeoClaw, and for more information on the Clawpack suite of software." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Contents\n", + "\n", + " - Fetching a topo file from the web\n", + " - Reading a topo file\n", + " - Plotting topography data\n", + " - Cropping topography" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Setup notebook and environment:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "%pylab inline" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from clawpack.geoclaw import topotools\n", + "import os\n", + "CLAW = os.environ['CLAW']\n", + "datadir = os.path.join(CLAW,'geoclaw','scratch') # directory for some sample data files" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "## Fetching a topo file from the web\n", + "\n", + "Many GeoClaw examples are set up to use topography files that have already been created and archived on the web, e.g. the example found in `$CLAW/geoclaw/examples/tsunami/chile2010` uses a topo file that can be obtained by these commands (which are also found in the Python script `maketopo.py` found in that directory):" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from clawpack.geoclaw import util\n", + "filename = 'etopo10min120W60W60S0S.asc'\n", + "url = 'http://www.geoclaw.org/topo/etopo/' + filename\n", + "util.get_remote_file(url=url, output_dir=datadir, force=True, verbose=True)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If all you want to do is use this topo file in a GeoClaw run, you do not need to use any further Python tools." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "\n", + "## Reading a topo file\n", + "\n", + "In order to plot or manipulate the topo data, we first read the file we just downloaded into a `topotools.Topography` object. To do so, we must know how the data is formatted. This file is in the format `topo_type==2` as described at :" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "topo_path = os.path.join(datadir, filename)\n", + "topo = topotools.Topography()\n", + "topo.read(topo_path, topo_type=2)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can now do various things with the data. First let's print out some basic attributes:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "print \"The extent of the data in longitude and latitude: \"\n", + "print topo.extent" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "print \"The shapes of 1-dimensional arrays of longitude x and latitude y values:\", topo.x.shape, topo.y.shape\n", + "print \"The shapes of 2-dimensional arrays X,Y and the topography data Z:\", topo.Z.shape" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "From the filename you might guess this is 10-arcminute data, we can check that it is:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "print \"topo.delta = \",topo.delta\n", + "print \"10 arcminutes is 1/6 degree = %8.6f degree\" % (1./6.)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "## Plotting topography data\n", + "\n", + "A simple plot with various defaults used can be obtained simply by:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "topo.plot()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Various arguments can be supplied for more control, e.g.:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "topo.plot(axes=None, region_extent=(-90,-60,-60,-40), coastlines=False)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Of course you can plot the data any way you want using the X, Y, and Z attributes, e.g." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "contourf(topo.X, topo.Y, topo.Z, [-7000, -5000, -4000, -3000, -2000, 0])\n", + "colorbar()\n", + "# rescale aspect ratio based on mean latitude so distances more correct:\n", + "gca().set_aspect(1.0 / cos(pi / 180.0 * topo.y.mean())) " + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "## Cropping topography\n", + "\n", + "Sometimes it is useful to crop a large topo file to create a smaller one that contains a sub-region, e.g." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "filter_region = (-90, -60, -60, -40)\n", + "topo2 = topo.crop(filter_region)\n", + "topo2.Z.shape" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "topo2.plot()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The cropped topography can be written out to a new file, which will be smaller than the original:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "topo2_path = os.path.join(datadir,'tip_of_south_america.tt3')\n", + "topo2.write(topo2_path, topo_type=3)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "print \"Size of the new file: %s bytes\" % os.stat(topo2_path).st_size\n", + "print \"Size of the original file: %s bytes\" % os.stat(topo_path).st_size" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To obtain topography for other regions, see the links at " + ] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file