notebook.community

Edit and run

Reading SW4 output

In this tutorial we shall go over several examples to show how to read SW4 output.





In [2]:



    


import os
import sys

import pySW4 as sw4
import obspy
import matplotlib.pyplot as plt
import numpy as np

Images





In [3]:



    


f = '../sw4-v1.1/examples/rfile/results/berkeley.cycle=00000.z=0.topo.sw4img'
topo = sw4.read_image(f)
topo.plot(cmap='terrain')



    


















    
Out[3]:








(<matplotlib.figure.Figure at 0x11c4bb890>,
 <matplotlib.axes._subplots.AxesSubplot at 0x11c4bbd50>,
 <matplotlib.colorbar.Colorbar at 0x11fc837d0>)

Plotting with symbols

pySW4 can parse the input file and plot source and station locations automatically.





In [4]:



    


f = '../sw4-v1.1/examples/rfile/results/berkeley.cycle=00000.z=0.p.sw4img'
input_, output_ = sw4.read_metadata('../sw4-v1.1/examples/rfile/berkeley.sw4',
                                    '../sw4-v1.1/examples/rfile/mon.txt')
image = sw4.read_image(f, input_)

print image
for patch in image.patches:
    print patch
    
fig, ax, cb = image.plot()



    


















    






        Image information :
        ----------------- :
                 Filename : berkeley.cycle=00000.z=0.p.sw4img
                Precision : <type 'numpy.float32'>
        Number of patches : 1
                  Time, s : 0.0
                    Plane : z
               Coordinate : 0.0
                     Mode : Vp, m/s
Curvilinear grid included : False
             Image extent : (-10.0, 11010.0, -10.0, 12010.0)
            Creation time : Tue Sep 27 11:55:41 2016

Patch information :
----------------- :
           Number : 0
       Spacing, m : 20.0
          Zmin, m : 0.0
           Extent : (-10.0, 11010.0, -10.0, 12010.0)
               ni : 601
               nj : 551
              Max : 2580.0
              Min : 768.0
              STD : 458.01373291
              RMS : 1578.19506836

Velocity magnitude plots





In [5]:



    


f = '../sw4-v1.1/examples/rfile/results/berkeley.cycle=01402.z=0.magdudt.sw4img'
input_, output_ = sw4.read_metadata('../sw4-v1.1/examples/rfile/berkeley.sw4',
                                    '../sw4-v1.1/examples/rfile/mon.txt')
image = sw4.read_image(f, input_)
    
fig, ax, cb = image.plot(vmax='3', cmap='jet')

Horizontal Peak Ground Velocity (HPGV) plots





In [6]:



    


f = '../sw4-v1.1/examples/rfile/results/berkeley.cycle=14016.z=0.hmaxdudt.sw4img'
input_, output_ = sw4.read_metadata('../sw4-v1.1/examples/rfile/berkeley.sw4',
                                    '../sw4-v1.1/examples/rfile/mon.txt')
image = sw4.read_image(f, input_)
    
fig, ax, cb = image.plot(vmax='3', cmap='hot_r')

Plotting cross-section patches





In [7]:



    


f = '../sw4-v1.1/examples/rfile/results/berkeley.cycle=00000.x=2500.p.sw4img'
image = sw4.read_image(f, input_)

print image
for patch in image.patches:
    print patch

print image.curvilinear_grid_patch
    
image.patches[1].plot()
image.patches[0].plot()



    


















    






        Image information :
        ----------------- :
                 Filename : berkeley.cycle=00000.x=2500.p.sw4img
                Precision : <type 'numpy.float32'>
        Number of patches : 2
                  Time, s : 0.0
                    Plane : x
               Coordinate : 2500.0
                     Mode : Vp, m/s
Curvilinear grid included : True
             Image extent : (-10.0, 11010.0, -330.42556762695312, 5010.0)
            Creation time : Tue Sep 27 11:55:43 2016

Patch information :
----------------- :
           Number : 0
       Spacing, m : 20.0
          Zmin, m : 2000.0
           Extent : (-10.0, 11010.0, 1990.0, 5010.0)
               ni : 551
               nj : 151
              Max : 5551.25
              Min : 3512.5
              STD : 542.043334961
              RMS : 5088.8828125

Patch information :
----------------- :
           Number : 1
       Spacing, m : 20.0
          Zmin, m : -320.425567627
           Extent : (-10.0, 11010.0, -330.42556762695312, 2010.0)
               ni : 551
               nj : 114
              Max : 5045.0
              Min : 768.0
              STD : 959.145141602
              RMS : 3736.68237305

Patch information :
----------------- :
           Number : curvilinear grid
       Spacing, m : 20.0
          Zmin, m : -320.425567627
           Extent : (-10.0, 11010.0, -330.42556762695312, 2010.0)
               ni : 551
               nj : 114
              Max : 2000.0
              Min : -320.425567627
              STD : 617.937072754
              RMS : 1099.14575195











    
Out[7]:








(<matplotlib.figure.Figure at 0x121142ad0>,
 <matplotlib.axes._subplots.AxesSubplot at 0x122c39690>,
 <matplotlib.colorbar.Colorbar at 0x122680cd0>)

Plotting the entire cross-section





In [8]:



    


fig, ax, cb = image.plot(projection_distance=1000)
ax.set_ylim(image.extent[3], image.extent[2]-300)



    


















    
Out[8]:








(5010.0, -630.42556762695312)










    
























In [9]:



    


f = '../sw4-v1.1/examples/rfile/results/berkeley.cycle=01402.x=2500.magdudt.sw4img'
image = sw4.read_image(f, input_)

fig, ax, cb = image.plot(vmax='3', projection_distance=1000)
ax.set_ylim(image.extent[3], image.extent[2]-300)



    


















    
Out[9]:








(5010.0, -630.42556762695312)










    
























In [10]:



    


f = '../sw4-v1.1/examples/rfile/results/berkeley.cycle=00000.y=2500.p.sw4img'
image = sw4.read_image(f, input_)

print image
for patch in image.patches:
    print patch

print image.curvilinear_grid_patch
    
image.patches[1].plot()
image.patches[0].plot()
fig, ax, cb = image.plot(projection_distance=1800)
ax.set_ylim(image.extent[3], image.extent[2]-300)



    


















    






        Image information :
        ----------------- :
                 Filename : berkeley.cycle=00000.y=2500.p.sw4img
                Precision : <type 'numpy.float32'>
        Number of patches : 2
                  Time, s : 0.0
                    Plane : y
               Coordinate : 2500.0
                     Mode : Vp, m/s
Curvilinear grid included : True
             Image extent : (-10.0, 12010.0, -423.19338989257812, 5010.0)
            Creation time : Fri Sep 23 11:49:59 2016

Patch information :
----------------- :
           Number : 0
       Spacing, m : 20.0
          Zmin, m : 2000.0
           Extent : (-10.0, 12010.0, 1990.0, 5010.0)
               ni : 601
               nj : 151
              Max : 5213.4375
              Min : 3532.5
              STD : 445.164398193
              RMS : 4466.79492188

Patch information :
----------------- :
           Number : 1
       Spacing, m : 20.0
          Zmin, m : -413.193389893
           Extent : (-10.0, 12010.0, -423.19338989257812, 2010.0)
               ni : 601
               nj : 114
              Max : 3687.5
              Min : 1800.0
              STD : 501.37600708
              RMS : 2910.00366211

Patch information :
----------------- :
           Number : curvilinear grid
       Spacing, m : 20.0
          Zmin, m : -413.193389893
           Extent : (-10.0, 12010.0, -423.19338989257812, 2010.0)
               ni : 601
               nj : 114
              Max : 2000.0
              Min : -413.193389893
              STD : 666.885375977
              RMS : 1090.12866211











    
Out[10]:








(5010.0, -723.19338989257812)

Seismograms

Station array or line

Reading an array or a line of stations added to the input file by station_array or station_line.





In [11]:



    


name = 'surface'
path = '../sw4-v1.1/examples/rfile/results'
array = sw4.read_stations(name, path, 'velocity')



    











In [12]:



    


for sta in array:
    x, y, z = (sta.stats.coordinates.x,
               sta.stats.coordinates.y,
               sta.stats.coordinates.z)
               
    plt.plot(y, x,'kv', mfc='none')



    


















    
























In [13]:



    


array[10].plot()



    


















    






/Users/shahar/anaconda2/lib/python2.7/site-packages/matplotlib/cbook.py:136: MatplotlibDeprecationWarning: The axisbg attribute was deprecated in version 2.0. Use facecolor instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)

Plot time slices of the Z component through all stations of the array





In [14]:



    


x, y = input_.source[0].x, input_.source[0].y
z_data = array.get_data('*Z')

for time in (3.5, 5.0, 6.2):
    fig, ax = plt.subplots()
    
    step = int(time / array.delta)
    values_t = z_data[:, :, step]
    max_ = np.abs(values_t).max()
    im = ax.imshow(values_t, interpolation='nearest', extent=array.extent,
               cmap='seismic', vmin=-max_, vmax=max_, origin='lower')
    plt.colorbar(im)
    ax.set_title('time: {:.3f}, s'.format(array.times()[step]))
    
    ax.plot(y, x, 'k*', mfc='none', ms=10)
    
    [ax.plot(y_, x_, 'kv', mfc='none') for x_ in array.xi for y_ in array.yi ]

Single station

Reading and plotting with obspy





In [15]:



    


st = obspy.read('../sw4-v1.1/examples/rfile/results/BDSN.BKS*', 'SAC')
st.plot(size=(800, 400))

Reading and plotting with pySW4





In [16]:



    


name = 'BDSN.BKS'
path = '../sw4-v1.1/examples/rfile/results'
st = sw4.read_stations(name, path, 'velocity')

starttime = st.starttime
st.trim(endtime=starttime + 20)
fig, ax = st.plot_traces(lw=1)

Content source: shaharkadmiel/pySW4

Similar notebooks:

notebook.community | gallery | about