Python integration

New in version 5.3.0!

Python is a rapidly glowing ecosystem for data analysis and visualization in seismology. Followings are some tips & tricks for those who would like to handle the input/output of the OpenSWPC in python. The following example uses

• ObsPy
• numpy
• xarray
• PyGMT

under the Jupyter Notebook environment.

Input parameter files

A python module OpenSWPC/src/tools/swpc.py contains several utility functions to handle input parameter to the OpenSWPC. To use it, first add this directory to the module path:

import sys
sys.path.append('/path/to/OpenSWPC/src/tools/')

Then, the module can be loaded as follows:

import swpc

A function prm_new() create the default parameter set. These parameters are equivalent to that defined in example/input.inf.

prm0 = swpc.prm_new()

In python, parameter data are stored in a simple dictionary as shown below. Thus, the user can easily change the parameter values if necessary.

print(prm0)

{
 'title': 'swpc', 'odir': './out', 'ntdec_r': 50, 'strict_mode': False, 
 'nproc_x': 2, 'nproc_y': 2, 'nx': 384, 'ny': 384, 'nz': 384, 'nt': 1000, 
 'dx': 0.5, 'dy': 0.5, 'dz': 0.5, 'dt': 0.02, 'vcut': 1.5, 
 'xbeg': -96.0, 'ybeg': -96.0, 'zbeg': -10.0, 'tbeg': 0.0, 'clon': 139.7604, 'clat': 35.7182, 
 'phi': 0.0, 'fq_min': 0.02, 'fq_max': 2.0, 'fq_ref': 1.0, 'snp_format': 'netcdf', 
 'xy_ps%sw': False, 'xz_ps%sw': True, 'yz_ps%sw': False, 'fs_ps%sw': False, 
 'ob_ps%sw': True, 'xy_v%sw': False, 'xz_v%sw': True, 'yz_v%sw': False, 
 'fs_v%sw': False, 'ob_v%sw': True, 'xy_u%sw': False, 'xz_u%sw': True, 
 'yz_u%sw': False, 'fs_u%sw': False, 'ob_u%sw': True, 
 'z0_xy': 7.0, 'x0_yz': 0.0, 'y0_xz': 0.0, 'ntdec_s': 5, 
 'idec': 2, 'jdec': 2, 'kdec': 2, 'sw_wav_v': True, 'sw_wav_u': False, 
 'sw_wav_stress': False, 'sw_wav_strain': False, 'ntdec_w': 5, 
 'st_format': 'xy', 'fn_stloc': './example/stloc.xy', 'wav_format': 'sac', 
 'wav_calc_dist': False, 'stf_format': 'xym0ij', 'stftype': 'kupper', 
 'fn_stf': './example/source.dat', 'sdep_fit': 'asis', 'bf_mode': False, 
 'pw_mode': False, 'pw_ztop': 100.0, 'pw_zlen': 30.0, 'pw_ps': 'p', 
 'pw_strike': 0.0, 'pw_dip': 0.0, 'pw_rake': 0.0, 'abc_type': 'pml', 
 'na': 20, 'stabilize_pml': False, 'vmodel_type': 'lhm', 'is_ocean': True, 
 'topo_flatten': False, 'munk_profile': True, 'earth_flattening': False, 
 'vp0': 5.0, 'vs0': 3.0, 'rho0': 2.7, 'qp0': 200, 'qs0': 200, 'topo0': 0, 
 'dir_grd': '${DATASET}/vmodel/ejivsm/', 'fn_grdlst': './example/grd.lst', 
 'node_grd': 0, 'fn_lhm': 'example/lhm.dat', 'dir_rmed': './in/', 
 'fn_grdlst_rmed': './example/grd.lst', 'rhomin': 1.0, 'fn_rmed0': 'dummy.ns', 
 'is_ckp': False, 'ckpdir': './out/ckp', 'ckp_interval': 1000000, 
 'ckp_time': 1000000.0, 'ckp_seq': True, 'green_mode': False, 'green_stnm': 'st01', 
 'green_cmp': 'z', 'green_trise': 1.0, 'green_bforce': False, 'green_maxdist': 550.0,
 'green_fmt': 'llz', 'fn_glst': 'example/green.lst', 'stopwatch_mode': False, 
 'benchmark_mode': False, 'ipad': 0, 'jpad': 0, 'kpad': 0
}

Also, one can read the parameter file by using prm_read() function. The following

prm1 = swpc.prm_read('/path/to/OpenSWPC/example/input.inf')

Apply some modification to the prm1 ...

prm1['title'] = 'new-simulation'
prm1['nx'] = 1024
prm1['ny'] = 1024
prm1['nz'] = 1024

To save the parameter to use in OpenSWPC, use prm_print() function. In default, this function exports all parameters to standard output. The followings are an example to save the parameter to file:

with open('swpc.prm', 'w') as fp:
    swpc.prm_print(prm1, io=fp)

Confirm the modified parameter is stored in the specified file:

! head swpc.prm

title              =   'new-simulation'
odir               =   './out'
ntdec_r            =   50
strict_mode        =   .false.
nproc_x            =   2
nproc_y            =   2
nx                 =   1024
ny                 =   1024
nz                 =   1024
nt                 =   1000

Output waveforms

ObsPy is a de-facto standard on reading seismograms and applying basic manipulation of seismic data. Since the ObsPy can directly read the SAC-formatted seismogram, one can easily import the output waveform data generated by the OpenSWPC without any additional processings.

Note that some user reported that the ObsPy occasionaly fails to specify the format when it reads SAC files produced by OpenSWPC properly. It is recommended to add format='sac' option to obspy.read() function to explicitly tell this is SAC-formatted file.

ls /path/to/OpenSWPC/output/wav/swpc_N.N.AAKH.Vz.sac

/path/to/OpenSWPC/output/wav/swpc_N.N.AAKH.Vz.sac

import obspy

tr = obspy.read('/path/to/OpenSWPC/output/wav/*.Vz.sac', format='sac')

# set distance to stats
# original distance header in stats.sac.dist is not read by record section plot function in obspy
for t in tr:
    t.stats.distance = t.stats.sac.dist * 1000

tr.plot(type='section', scale=5, offset_min=0, offset_max=500*1000, recordlength=200, reduce = 7).show()

tr[10].spectrogram(wlen=60, log=True, title='Spectrogram example')

Output snapshots

2D snapshot data produced by OpenSWPC follows the CF metadata conventions of NetCDF format. There seems several libraries to read NetCDF files, and here is an exaple using xarray which has an affinity to the PyGMT.

fn_ob = '/path/to/OpenSWPC/output/swpc_N.ob.v.nc'

import xarray as xr

xr_ob = xr.open_dataset(fn_ob)

xarray dataset contains multiple list (array) data, as the NetCDF data format:

xr_ob

<xarray.Dataset>
Dimensions:  (x: 1067, y: 533, t: 600)
Coordinates:
  * x        (x) float32 -799.2 -797.8 -796.2 -794.8 ... 795.2 796.8 798.2 799.8
  * y        (y) float32 -399.2 -397.8 -396.2 -394.8 ... 394.2 395.8 397.2 398.8
  * t        (t) float32 0.025 1.025 2.025 3.025 ... 596.0 597.0 598.0 599.0
    lon      (y, x) float32 ...
    lat      (y, x) float32 ...
Data variables:
    rho      (y, x) float32 ...
    lambda   (y, x) float32 ...
    mu       (y, x) float32 ...
    topo     (y, x) float32 ...
    Vx       (t, y, x) float32 ...
    Vy       (t, y, x) float32 ...
    Vz       (t, y, x) float32 ...
    max-V    (y, x) float32 ...
    max-H    (y, x) float32 ...
    max-A    (y, x) float32 ...
Attributes: (12/26)
    generated_by:  SWPC
    codetype:      SWPC_3D
    hdrver:        6
    title:         swpc_N
    exedate:       1664430044
    ns1:           1067
    ...            ...
    evdp:          43.18
    evx:           -320.6624
    evy:           -66.1461
    clon:          142.0
    clat:          39.0
    phi:           0.0

xarray.Dataset

• Dimensions:
  • x: 1067
  • y: 533
  • t: 600
• Coordinates: (5)
  • x
    (x)
    float32
    -799.2 -797.8 ... 798.2 799.8

    long_name :

    x

    standard_name :

    projection_x_coordinate

    units :

    km

    actual_range :

    [-799.25 799.75]

    array([-799.25, -797.75, -796.25, ...,  796.75,  798.25,  799.75],
          dtype=float32)

  • y
    (y)
    float32
    -399.2 -397.8 ... 397.2 398.8

    long_name :

    y

    standard_name :

    projection_y_coordinate

    units :

    km

    actual_range :

    [-399.25 398.75]

    array([-399.25, -397.75, -396.25, ...,  395.75,  397.25,  398.75],
          dtype=float32)

  • t
    (t)
    float32
    0.025 1.025 2.025 ... 598.0 599.0

    long_name :

    t

    units :

    s

    array([2.50000e-02, 1.02500e+00, 2.02500e+00, ..., 5.97025e+02, 5.98025e+02,
           5.99025e+02], dtype=float32)

  • lon
    (y, x)
    float32
    ...

    long_name :

    Longitude

    units :

    degrees_east

    actual_range :

    [136.8374 147.15616]

    [568711 values with dtype=float32]

  • lat
    (y, x)
    float32
    ...

    long_name :

    Latitude

    units :

    degrees_north

    actual_range :

    [31.725807 46.200157]

    [568711 values with dtype=float32]

• Data variables: (10)
  • rho
    (y, x)
    float32
    ...

    long_name :

    rho

    units :

    10^3 kg/cm^3

    actual_range :

    [2.25 3.2 ]

    [568711 values with dtype=float32]

  • lambda
    (y, x)
    float32
    ...

    long_name :

    lambda

    units :

    10^9 Pa

    actual_range :

    [10.002981 50.560623]

    [568711 values with dtype=float32]

  • mu
    (y, x)
    float32
    ...

    long_name :

    mu

    units :

    10^9 Pa

    actual_range :

    [ 5.004911 64.35583 ]

    [568711 values with dtype=float32]

  • topo
    (y, x)
    float32
    ...

    long_name :

    Topography

    units :

    m

    actual_range :

    [-9315.041 3645.4475]

    [568711 values with dtype=float32]

  • Vx
    (t, y, x)
    float32
    ...

    long_name :

    Vx

    units :

    m/s

    actual_range :

    [-5.082945 4.5914364]

    [341226600 values with dtype=float32]

  • Vy
    (t, y, x)
    float32
    ...

    long_name :

    Vy

    units :

    m/s

    actual_range :

    [-5.6842103 6.2938366]

    [341226600 values with dtype=float32]

  • Vz
    (t, y, x)
    float32
    ...

    long_name :

    Vz

    units :

    m/s

    actual_range :

    [-2.0089784 2.3452911]

    [341226600 values with dtype=float32]

  • max-V
    (y, x)
    float32
    ...

    long_name :

    Maximum amplitude of the vertical component

    units :

    m/s

    actual_range :

    [2.6916168e-05 2.3838887e+00]

    [568711 values with dtype=float32]

  • max-H
    (y, x)
    float32
    ...

    long_name :

    Maximum amplitude of the horizontal components

    units :

    m/s

    actual_range :

    [7.022181e-05 6.767687e+00]

    [568711 values with dtype=float32]

  • max-A
    (y, x)
    float32
    ...

    long_name :

    Maximum amplitude of the vector motion

    units :

    m/s

    actual_range :

    [7.5203607e-05 7.0090966e+00]

    [568711 values with dtype=float32]

• Attributes: (26)

  generated_by :

  SWPC

  codetype :

  SWPC_3D

  hdrver :

  6

  title :

  swpc_N

  exedate :

  1664430044

  ns1 :

  1067

  ns2 :

  533

  beg1 :

  -799.25

  beg2 :

  -399.25

  na1 :

  6

  na2 :

  6

  ds1 :

  1.5

  ds2 :

  1.5

  nmed :

  6

  nsnp :

  3

  coordinate :

  ob

  datatype :

  v3

  dt :

  1.0

  evlo :

  141.2653

  evla :

  36.1083

  evdp :

  43.18

  evx :

  -320.6624

  evy :

  -66.1461

  clon :

  142.0

  clat :

  39.0

  phi :

  0.0

Each data can be obtained by appending data varialbe name after . of the object:

xr_ob.topo.data

array([[-3902.103 , -3902.103 , -3902.103 , ...,   841.1037,   841.1037,
          841.1037],
       [-3902.103 , -3902.103 , -3902.103 , ...,   841.1037,   841.1037,
          841.1037],
       [-3902.103 , -3902.103 , -3902.103 , ...,   841.1037,   841.1037,
          841.1037],
       ...,
       [-6059.085 , -6059.085 , -6059.085 , ..., -3449.3516, -3449.3516,
        -3449.3516],
       [-6059.085 , -6059.085 , -6059.085 , ..., -3449.3516, -3449.3516,
        -3449.3516],
       [-6059.085 , -6059.085 , -6059.085 , ..., -3449.3516, -3449.3516,
        -3449.3516]], dtype=float32)

xr_ob.lon.data

array([[137.7886 , 137.788  , 137.78738, ..., 136.83992, 136.83867,
        136.8374 ],
       [137.80438, 137.80377, 137.80316, ..., 136.85922, 136.85797,
        136.85672],
       [137.82018, 137.81956, 137.81895, ..., 136.87854, 136.87729,
        136.87604],
       ...,
       [146.17456, 146.17517, 146.17578, ..., 147.11504, 147.11627,
        147.11752],
       [146.19035, 146.19096, 146.19157, ..., 147.13434, 147.13559,
        147.13684],
       [146.20613, 146.20674, 146.20737, ..., 147.15364, 147.1549 ,
        147.15616]], dtype=float32)

xr_ob.lat.data

array([[31.725807, 31.7393  , 31.752792, ..., 46.056538, 46.06998 ,
        46.08342 ],
       [31.72633 , 31.739822, 31.753315, ..., 46.05741 , 46.070854,
        46.084293],
       [31.726849, 31.740343, 31.753836, ..., 46.05828 , 46.07172 ,
        46.085163],
       ...,
       [31.727022, 31.740515, 31.754007, ..., 46.058567, 46.07201 ,
        46.085453],
       [31.726503, 31.739996, 31.753489, ..., 46.0577  , 46.071144,
        46.084583],
       [31.725983, 31.739475, 31.752968, ..., 46.056828, 46.07027 ,
        46.083714]], dtype=float32)

... and time slice of the wavefield can be reffered by the sequential time step number like

xr_ob.Vz[50].data

array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]], dtype=float32)

Note that the OpenSWPC performs simulations in the Cartesian coordinate. The result is mapped to the geographical cooridnate upon output, but the data is not evenly-spaced in longitude-latitude system. Thus, interpolation and/or resampling may be necessary to plot the wavefield on a map.

import pygmt
import numpy as np

# make interpolated grid data by nearneighbor function

region = [136.5, 147.5, 31.5, 46.5]
spacing = 0.01
radius = 0.05

# topography grid
grd_topo = pygmt.nearneighbor(x = xr_ob.lon.data.ravel(), 
                              y = xr_ob.lat.data.ravel(),
                              z = xr_ob.topo.data.ravel(), 
                              region = region, spacing = spacing, search_radius = radius) 

# wavefield at 50th time step
grd_wave = pygmt.nearneighbor(x = xr_ob.lon.data.ravel(), 
                              y = xr_ob.lat.data.ravel(),
                              z = xr_ob.Vz[50].data.ravel(), 
                              region = region, spacing = spacing, search_radius = radius) 

fig = pygmt.Figure()

with pygmt.config(MAP_FRAME_TYPE  = 'plain'):

    # background topography
    pygmt.makecpt(cmap='gray', series=[-2000, 2000, 500], continuous=True)
    fig.grdimage(grid=grd_topo, projection='N142/10c', region=region, frame=['WSen', 'xaf', 'yaf'])

    # coastline
    fig.coast(resolution = 'h', area_thresh = 100, shorelines='thin,black')

    # wavefield with tranparency
    pygmt.makecpt(cmap='polar',series = [-0.1, 0.1, 0.02], continuous=True, background=True)
    fig.grdimage(grid=grd_wave, transparency=40)

fig.show()