Added 2015 Illapel earthquake example

illapel2015_chile/Makefile

@@ -0,0 +1,68 @@
+# Makefile for Clawpack code in this directory.
+# This version only sets the local files and frequently changed
+# options, and then includes the standard makefile pointed to by CLAWMAKE.
+CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common
+
+# See the above file for details and a list of make options, or type
+#   make .help
+# at the unix prompt.
+
+
+# Adjust these variables if desired:
+# ----------------------------------
+
+CLAW_PKG = geoclaw                  # Clawpack package to use
+EXE = xgeoclaw                 # Executable to create
+SETRUN_FILE = setrun.py        # File containing function to make data
+OUTDIR = _output               # Directory for output
+SETPLOT_FILE = setplot.py      # File containing function to set plots
+PLOTDIR = _plots               # Directory for plots
+
+# Environment variable FC should be set to fortran compiler, e.g. gfortran
+
+# Compiler flags can be specified here or set as an environment variable
+FFLAGS ?= 
+
+# ---------------------------------
+# package sources for this program:
+# ---------------------------------
+
+GEOLIB = $(CLAW)/geoclaw/src/2d/shallow
+include $(GEOLIB)/Makefile.geoclaw
+
+# ---------------------------------------
+# package sources specifically to exclude
+# (i.e. if a custom replacement source 
+#  under a different name is provided)
+# ---------------------------------------
+
+EXCLUDE_MODULES = \
+
+EXCLUDE_SOURCES = \
+
+# ----------------------------------------
+# List of custom sources for this program:
+# ----------------------------------------
+
+
+MODULES = \
+
+SOURCES = \
+  $(CLAW)/riemann/src/rpn2_geoclaw.f \
+  $(CLAW)/riemann/src/rpt2_geoclaw.f \
+  $(CLAW)/riemann/src/geoclaw_riemann_utils.f \
+
+#-------------------------------------------------------------------
+# Include Makefile containing standard definitions and make options:
+include $(CLAWMAKE)
+
+# Construct the topography data
+.PHONY: topo all
+topo:
+	python maketopo.py
+
+all: 
+	$(MAKE) topo
+	$(MAKE) .plots
+	$(MAKE) .htmls
+

illapel2015_chile/README.md

@@ -0,0 +1,80 @@
+# 2015 Illapel Earthquake Report
+
+This example contains code for running a simulation of the tsunami generated by the 2015 Illapel earthquake on the Chilean coast.
+
+## The Earthquake
+
+*Sources: USGS (https://earthquake.usgs.gov/earthquakes/eventpage/us20003k7a/executive), NOAA (https://www.ngdc.noaa.gov/hazard/16sep2015.html)*
+
+On September 16, 2015, an earthquake with a magnitude of 8.3 occurred on the boundary between the Nazca and South America tectonic plates, north of the rupture zone from the 2010 Chile earthquake. It produced significant aftershocks, some of which exceeded magnitude 7.0, and it was reported to have caused 15 deaths and 14 injuries.
+
+In this example, we model the tsunami resulting from the earthquake. We validate our model using data from two tide gauges in the region.
+
+## Topography/Bathymetry Data
+
+*Topography data was sourced from the ETOPO 10-arc-minute model (Esri ASCII format) here: http://depts.washington.edu/clawpack/geoclaw/topo/etopo/etopo10min120W60W60S0S.asc*
+
+In the example, [maketopo.py](maketopo.py) downloads the data directly from source into the scratch directory within GeoClaw. 
+
+```python
+# maketopo.py, lines 26-27
+url = 'http://depts.washington.edu/clawpack/geoclaw/topo/etopo/' + topo_fname
+clawpack.clawutil.data.get_remote_file(url, output_dir=scratch_dir, file_name=topo_fname, verbose=True)
+```
+
+If you want to run a simulation with more refined topography, simply download the file into the scratch directory (or edit the download URL) and edit the appropriate filenames within [maketopo.py](maketopo.py) and [setrun.py](setrun.py). If the geographical bounds or topography file format are different, note that these parameters must be changed in [setrun.py](setrun.py).
+
+## Moving Topography Data
+
+*Moving topography (dtopo) data was sourced from the USGS finite fault model for the earthquake here: https://earthquake.usgs.gov/product/finite-fault/us20003k7a/us/1539809967421/basic_inversion.param*
+
+As with topography data, moving topography (dtopo) data is directly downloaded. 
+
+```python
+# maketopo.py, lines 57-59
+param_fname = 'basic_inversion.param'
+csv_fname = 'chile2015.csv'
+get_csv('https://earthquake.usgs.gov/product/finite-fault/us20003k7a/us/1539809967421/' + param_fname, param_fname, csv_fname)
+```
+In order to produce a CSV file of the subfaults readable by the `clawpack.geoclaw.dtopotools` module, the data file is reformatted in [get_dtopo_csv.py](get_dtopo_csv.py). Therefore, to use different finite fault data, edit this file alongside the download URL in [maketopo.py](maketopo.py) (and potentially the file format in [setrun.py](setrun.py)).
+
+## GeoClaw Results
+
+The simulation was set to run 5 hours from the origin time of the earthquake, with a maximum of 6 adaptive mesh refinement layers. Areas of high refinement were manually implemented around tide gauges. Surface plots and tide gauge plots are generated by [setplot.py](setplot.py), and they show predicted tide residuals ranging from -0.075 to +0.125. Overall, these values match NOAA tide data relatively well.
+
+## Validation Results
+
+*Tide data was downloaded directly from the NOAA database for the earthquake here: https://www.ngdc.noaa.gov/hazard/dart/2015chile.html*
+
+We test our simulation against tide residual values at two DART stations: gauge no. 32412 (SW of Lima, Peru) and gauge no. 32402 (W of Caldera, Chile). The simulation matched the residual patterns at the closer gauge (32402) relatively well but those at farther gauge (32412) relatively poorly.
+
+For additional testing, one can add more gauges for plotting. To do so, the download must be specified in [maketopo.py](maketopo.py) and the plotting in [setplot.py](setplot.py).
+
+## Conclusion
+
+Topography and moving topography data for the 2015 Illapel earthquake were obtained from the internet and used by GeoClaw to run a simulation of the regional sea surface. Tide residuals were compared with observed NOAA data and demonstrated good local predictiveness but poor faraway predictiveness. Better topography data and adaptive mesh refinement customization may yield improved results. More tide gauges may provide a better picture of the simulation's predictive scope.
+
+## Installation
+
+Make sure to have the following dependencies downloaded in order to run this example:
+- [Clawpack](https://www.clawpack.org/installing.html)
+- [pandas](https://pandas.pydata.org/)
+
+## Usage
+
+Make sure you have GNU Make installed, and run the following terminal commands inside of the `illapel2015_chile` directory:
+
+```
+make clobber
+make .data
+make xgeoclaw
+make .output
+make .plots
+```
+
+To optimize simulation speed, consider [setting OpenMP compiler flags for multiprocessing](https://www.clawpack.org/openmp.html).
+
+---
+
+For questions, please contact Ashwin Padaki at [email protected].
+

illapel2015_chile/format_gauge.py

@@ -0,0 +1,14 @@
+"""Script for formatting NOAA gauge data for the Chile 2015 earthquake"""
+def format_gauges_chile2015(gaugeinfo):
+    import pandas as pd
+    for gaugeno, offset in gaugeinfo:
+        df = pd.read_csv('dart{}.txt'.format(gaugeno), delim_whitespace=True)
+        df.columns = ['time', 'year', 'month', 'day', 'hours', 'minutes', 'seconds', 'raw_obs', 'fitted_tidal', 'residual', 'blank']
+        # only store time and residual data
+        df = df[['time', 'residual']]
+        # earthquake origin time (from USGS): 2015-09-16 22:54:32.860 (UTC)
+        origin_time = 259 + (22 * 60**2 + 54*60 + 32.860)/(24 * 60**2)
+        df.time = df.time.map(lambda time : 24*60**2 * (time-origin_time))
+        # option to include manual offsets to correct for equilibria that are not at sea level
+        df.residual = df.residual.map(lambda res : res - offset)
+        df.to_csv('{}_notide.txt'.format(gaugeno), index=False, header=False, sep=' ')

illapel2015_chile/get_dtopo_csv.py

@@ -0,0 +1,22 @@
+"""Script for formatting USGS finite fault data into a CSV file readable by GeoClaw's dtopotools module"""
+def get_csv(url, param_fname, csv_fname):
+    import clawpack.clawutil.data
+    import os    
+    try:
+        CLAW = os.environ['CLAW']
+    except:
+        raise Exception("*** Must first set CLAW enviornment variable")
+    scratch_dir = os.path.join(CLAW, 'geoclaw', 'scratch')
+    clawpack.clawutil.data.get_remote_file(url, output_dir=scratch_dir, file_name=param_fname, verbose=True)
+
+    import pandas as pd
+    df = pd.read_csv(os.path.join(scratch_dir,param_fname), skiprows=9, delim_whitespace=True)
+    df.columns = ['latitude', 'longitude', 'depth', 'slip', 'rake', 'strike', 'dip', 'rupture_time', 'rise_time', 'fall_time', 'mo']
+    # length/width can be automated, but for convenience, they're manually written here
+    df['length'] = 12.00
+    df['width'] = 8.80
+    # rigidity = (moment seismicity) / (slip * area)
+    # cm * km * km = 10^4 m^3
+    df['mu'] = df['mo'] / (10**4 * df['length'] * df['width'] * df['slip']) 
+    df.drop(['fall_time', 'mo'], axis=1, inplace=True)
+    df.to_csv(csv_fname, index=False)

illapel2015_chile/maketopo.py

@@ -0,0 +1,153 @@
+"""
+Script for downloading gauge and topography data, and generating differential topography data
+Call functions with makeplots==True to create plots of topo, slip, and dtopo.
+"""
+
+from __future__ import absolute_import
+from __future__ import print_function
+import os
+
+import clawpack.clawutil.data
+
+try:
+    CLAW = os.environ['CLAW']
+except:
+    raise Exception("*** Must first set CLAW enviornment variable")
+
+scratch_dir = os.path.join(CLAW, 'geoclaw', 'scratch')
+
+def get_topo(makeplots=False):
+    """
+    Retrieve the topo file from the GeoClaw repository.
+    """
+    from clawpack.geoclaw import topotools
+    # download 10 arc-minute topography data in region around Chile
+    topo_fname = 'etopo10min120W60W60S0S.asc'
+    url = 'http://depts.washington.edu/clawpack/geoclaw/topo/etopo/' + topo_fname
+    clawpack.clawutil.data.get_remote_file(url, output_dir=scratch_dir, file_name=topo_fname, verbose=True)
+
+    if makeplots:
+        from matplotlib import pyplot as plt
+        topo = topotools.Topography(os.path.join(scratch_dir,topo_fname), topo_type=2)
+        topo.plot()
+        fname = os.path.splitext(topo_fname)[0] + '.png'
+        plt.savefig(fname)
+        print("Created ",fname)
+
+# generate side-by-side fault slip and seafloor deformation contour plots
+def plot_subfaults_dZ(t, fig, fault, dtopo, xlower, xupper, ylower, yupper, xylim, dz_max):
+    fig.clf()
+    ax1 = fig.add_subplot(121)
+    ax2 = fig.add_subplot(122)
+    fault.plot_subfaults(axes=ax1, slip_color=True,
+                        slip_time=t, xylim=xylim)
+    dtopo.plot_dZ_colors(axes=ax2, t=t, cmax_dZ=dz_max)
+    ax2.set_xlim(xlower,xupper)
+    ax2.set_ylim(ylower,yupper)
+    return fig
+
+# create dtopo file from finite fault data    
+def make_dtopo(makeplots=False):
+    from clawpack.geoclaw import dtopotools
+    from matplotlib import pyplot as plt
+    import numpy
+    from get_dtopo_csv import get_csv
+
+    # specify the URL for the USGS finite fault data file and properly format into a CSV file
+    param_fname = 'basic_inversion.param'
+    csv_fname = 'chile2015.csv'
+    get_csv('https://earthquake.usgs.gov/product/finite-fault/us20003k7a/us/1539809967421/' + param_fname, param_fname, csv_fname)
+
+    dtopo_fname = os.path.join(scratch_dir, "dtopo_usgs150916.tt3")
+
+    # manually specify the bounds for the fault model (leave unchanged)
+    xlower = -73
+    xupper = -70.5
+    ylower = -33
+    yupper = -29.5
+    xylim = [xlower,xupper,ylower,yupper]
+
+    input_units = {
+        'depth' : 'km',
+        'slip' : 'cm', 
+        'length' : 'km',
+        'width' : 'km',
+        'mu' : 'Pa'
+    }
+
+    fault = dtopotools.CSVFault()
+    fault.read(csv_fname,input_units=input_units, coordinate_specification='noaa sift')
+    fault.rupture_type = 'dynamic'
+
+    print ("%s subfaults read in " % len(fault.subfaults))
+
+    if os.path.exists(dtopo_fname):
+        print("*** Not regenerating dtopo file (already exists): %s" \
+                    % dtopo_fname)
+    else:
+        print("Using Okada model to create dtopo file")
+        points_per_degree = 60
+        mx = int((xupper - xlower)*points_per_degree + 1)
+        my = int((yupper - ylower)*points_per_degree + 1)
+        x = numpy.linspace(xlower, xupper, mx)
+        y = numpy.linspace(ylower, yupper, my)
+
+        # list of the times (in seconds) when the dtopo data is captured
+        # here, `times` consists of the interval [0,200] divided evenly into 20 markers (0 and 200 included)
+        # note: ensure that the maximum value in `times` exceeds the maximum rupture time of any subfault
+        times = numpy.linspace(0,200,20)
+
+        # create and write dtopo file
+        fault.create_dtopography(x,y,times)
+        dtopo = fault.dtopo
+        dtopo.write(dtopo_fname, dtopo_type=3)
+
+        import clawpack.visclaw.JSAnimation.JSAnimation_frametools as J
+        plotdir = '_fault_slip_plots'
+        J.make_plotdir(plotdir, clobber=True)
+        fig = plt.figure(figsize=(12,5))
+        dzmax = abs(dtopo.dZ).max()
+
+        # generate fault slip and deformation plots
+        for k,t in enumerate(times):
+            plot_subfaults_dZ(t,fig, fault, dtopo, xlower, xupper, ylower, yupper, xylim, dzmax)
+            J.save_frame(k, plotdir = plotdir, verbose=True)
+
+    if makeplots:
+        if fault.dtopo is None:
+            # read in the pre-existing file:
+            print("Reading in dtopo file...")
+            dtopo = dtopotools.DTopography()
+            dtopo.read(dtopo_fname, dtopo_type=3)
+            x = dtopo.x
+            y = dtopo.y
+        plt.figure(figsize=(12,7))
+        ax1 = plt.subplot(121)
+        ax2 = plt.subplot(122)
+        fault.plot_subfaults(axes=ax1,slip_color=True)
+        ax1.set_xlim(x.min(),x.max())
+        ax1.set_ylim(y.min(),y.max())
+        dtopo.plot_dZ_colors(1.,axes=ax2)
+        fname = os.path.splitext(os.path.split(dtopo_fname)[-1])[0] + '.png'
+        plt.savefig(fname)
+        print("Created ",fname)
+
+# grabs gauge data from NOAA and formats it
+def get_gauge():
+    # find more gauge numbers here: https://www.ngdc.noaa.gov/hazard/dart/2015chile.html
+    gaugeinfo = []
+    # gaugeinfo.append([gauge_no, offset]) where `offset` is the nonzero equilibrium in the data
+    gaugeinfo.append([32412, -0.03])
+    gaugeinfo.append([32402, 0.00])
+
+    for gaugeno, offset in gaugeinfo:
+        url = 'https://www.ngdc.noaa.gov/hazard/data/DART/20150916_chile/dart{}_20150915to20150921_meter.txt'.format(gaugeno)
+        clawpack.clawutil.data.get_remote_file(url, output_dir='.', file_name='dart{}.txt'.format(gaugeno), verbose=True)
+
+    from format_gauge import format_gauges_chile2015
+    format_gauges_chile2015(gaugeinfo)
+
+if __name__=='__main__':
+    get_topo(False)
+    make_dtopo(False)
+    get_gauge()