pySW4.prep.rfileIO module

Python module to read and write rfiles.

author:Shahar Shani-Kadmiel (s.shanikadmiel@tudelft.nl)
copyright:Shahar Shani-Kadmiel
license:This code is distributed under the terms of the GNU Lesser General Public License, Version 3 (https://www.gnu.org/copyleft/lesser.html)
class Block(f, b)[source]

Bases: object

A class to hold rfile block header and data section

_read_block_data_section(f, precision)[source]

Private method. Reads the data section of cued file handle. Should not be used by the user. See read_block_data_section() method in the Model class.

qp()[source]

Return values of Qp for block

qs()[source]

Return values of Qs for block

rho()[source]

Return values of density for block

vp()[source]

Return values of vp for block

vs()[source]

Return values of vs for block

x()[source]

Return a vector with x coordinates of the block based on hh. Coordinates are distance in km.

y()[source]

Return a vector with y coordinates of the block based on hh. Coordinates are distance in km.

z()[source]

Return a vector with z coordinates of the block based on hv. Coordinates are elevation below sea level in km (positive is down).

class CrossSection(model, x1=None, x2=None, y1=None, y2=None)[source]

Bases: object

Class to generate and hold a cross-section of the material properties along a line in an existing model (Model).

Parameters:

model : Model

A populated Model object from which to extract the data

x1 : float

Distance in km along the x-axis of start-point.

x2 : float

Distance in km along the x-axis of end-point.

y1 : float

Distance in km along the y-axis of start-point.

y2 : float

Distance in km along the y-axis of end-point.

Returns:

CrossSection

A CrossSection object with data and a plotting method.

plot(property='vp', ax=None, draw_separator=False, **kwargs)[source]

Plot cross-section.

Parameters:

property : {‘vp’, ‘vs’, ‘qp’, ‘qs’}

Which property to plot.

ax : Axes

Use existing axes. If ax=None function returns: a Figure instance, a Axes instance, and a Colorbar instance. Otherwise, only a Colorbar instance is returned.

draw_separator : bool

If set to True, a separator is plotted between blocks.

Other Parameters:
 

vmin : int or float

Used to limit the scale of the data. By default the minimum and maximum of the data is used.

vmax : int or float

Used to limit the scale of the data. By default the minimum and maximum of the data is used.

ylim : tuple

Set the ylim (vertical) extent of the plot.

aspect : str or int or float

Changes the axes aspect ratio.

interpolation : str

Acceptable values are ‘none’, ‘nearest’, ‘bilinear’, ‘bicubic’, ‘spline16’, ‘spline36’, ‘hanning’, ‘hamming’, ‘hermite’, ‘kaiser’, ‘quadric’, ‘catrom’, ‘gaussian’, ‘bessel’, ‘mitchell’, ‘sinc’, ‘lanczos’.

pad : int or float

Horizontal space between the axes and the colorbar.

size : str or float

Width of the colorbar. If float, width is set to size. If string colorbar width is calculated as a fraction of the axes width.

For example:
size='3%' will result in a colorbar who’s width is

3% of the width of the axes.

Similarly:
size='0.3' will result in a colorbar who’s width is

0.3 of the width of the axes.

label : str

Label used for the colorbar. By default the label is set by the property plotted.

kwargs : imshow() propeties.

class Model(filename, block_number=None, verbose=False)[source]

Bases: object

A class to hold rfile header and blocks.

Note

Setting block_number='all' for large models may take a while and take up a lot of memory.

Block number 1 holds the elevation of the topography/bathymetry.

Parameters:

filename : file

Path to rfile

block_number : int, str, None

By default (block_number=None) only the rfile header and the block headers are read into a Model object. No block data sections are read. If block_number='all', all block data sections are read into memory. Therefor, you could read only one block by specifying block_number=? replacing ? with the block number you are interested in.

verbose : bool

If True, information about the read process is printed to stdout.

Returns:

Model

A Model object holding rfile header and block headers and perhaps also block data sections, depending on the block_number value.

get_cross_section(x1=None, x2=None, y1=None, y2=None)[source]

Extract a cross-section from a Model object.

Generate a cross-section of the material properties along a line in an existing model (Model).

Parameters:

model : Model

A populated Model object from which to extract the data

x1, x2 : float

Distance in km along the x-axis of start- and end-points.

y1, y2 : float

Distance in km along the y-axis of start- and end-points.

Returns:

CrossSection

A CrossSection object with data and a plotting method.

get_properties_at_point(x, y, z, property='all')[source]

Get the material properties at a point.

Parameters:

x, y, z : float

Distance in km along the x-, y-, and z-axis of the model. (z - positive is down)

property : {‘all’, ‘rho’, ‘vp’, ‘vs’, ‘qp’, ‘qs’}

‘all’ returns all 5 material properties. Otherwise, only ‘rho’, ‘vp’, ‘vs’, ‘qp’, or ‘qs’ are returned.

Returns:

ndarray or float

One or all of the material properties at point (x,y,z) (rho, vp, vs, qp, qs)

get_topography()[source]

Returns a 2d ndarray with elevation data retrieved from the topography/bathymetry block (the first block).

Note

The origin of the data is the bottom-left corner so if plotted with imshow() set the origin keyword to ‘lower’ or simply use plot_topography() for plotting.

get_z_profile(x, y, property='all')[source]

Get the material properties along a z-profile at location x,y.

Parameters:

x, y : float

Distance in km along the x-, and y-axis of the model.

property : {‘all’, ‘rho’, ‘vp’, ‘vs’, ‘qp’, ‘qs’}

‘all’ returns all 5 material properties. Otherwise, only ‘rho’, ‘vp’, ‘vs’, ‘qp’, or ‘qs’ are returned.

Returns:

2 ndarray ‘s

z : ndarray

Elevation values allong the profile, shape (n,)

properties : ndarray

Material properties (one or all):

(‘rho’, ‘vp’, ‘vs’, ‘qp’, ‘qs’), shape (n,) or (n, 5)

plot_topography(ax=None, **kwargs)[source]

Plot the top surface (topography and bathymetry) of the model.

Parameters:

ax : Axes

Use existing axes. If ax=None function returns a Figure instance, a Axes instance, and a Colorbar instance. Otherwise, only a Colorbar instance is returned.

Other Parameters:
 

vmin, vmax : int or float

Used to limit the scale of the data. By default the minimum and maximum of the data is used.

interpolation : str

Acceptable values are ‘none’, ‘nearest’, ‘bilinear’, ‘bicubic’, ‘spline16’, ‘spline36’, ‘hanning’, ‘hamming’, ‘hermite’, ‘kaiser’, ‘quadric’, ‘catrom’, ‘gaussian’, ‘bessel’, ‘mitchell’, ‘sinc’, ‘lanczos’.

pad : int or float

Horizontal space between the axes and the colorbar.

size : str or float

Width of the colorbar. If float, width is set to size. If string colorbar width is calculated as a fraction of the axes width.

For example:
size='3%' will result in a colorbar who’s width is

3% of the width of the axes.

Similarly:
size='0.3' will result in a colorbar who’s width is

0.3 of the width of the axes.

label : str

Label used for the colorbar. By default the label is set by the property plotted.

kwargs : imshow() propeties.

read_block_data_section(block_number)[source]

Read the data section of the specified block into the current Model object to be stored in memory.

Note

For large models this may take a while and take up a lot of memory.

Parameters:

block_number : int or str

The number of the block you are interested in. Block number 1 holds the elevation of the topography/bathymetry. Blocks 2 and on (2<=block_number<=self.nb) hold the material properties of the model. If block_number='all', all block data sections are read into memory.

z()[source]

Return a vector (ndarray) with z coordinates of the model based on the hv of each block. Coordinates are distance in km along the z-axis of the model which are elevation below sea level (positive is down).

line_func(h1, h2, v1, v2, spacing=1)[source]

Return all coordinates along a straight line from point (h1,`v1`) to point (h2,`v2`) with spacing.

Parameters:

h1, h2 : int or float

Horizontal coordinate of the start-point and end-point.

v1, v1 : int or float

Vertical coordinate of the start-point and end-point.

spacing : int or float

Spacing along the straight line between coordinates. Controls the data-type returned: If spacing type is int, the values in the returned arrays are int. If spacing type is float, the values in the returned arrays are float.

Returns:

2 class:~numpy.ndarray

2 arrays with horizontal coordinates and vertical coordinates.

read(filename, block_number=None, verbose=False)[source]

Read rfile into Model object.

See Model for documentation.

read_block_hdr(f)[source]

Read rfile block header.

Parameters:

f : file

Open file handle in 'rb' mode of the rfile to be read.

Returns:

rfile block header data (7 elements)

hh, hv, z0, nc, ni, nj, nk

hh, hv : numpy.float64

Grid size in the horizontal (x and y) and vertical (z) directions in meters.

z0 : numpy.float64

The base z-level of the block. Not used for the first block which holds the elevation of the topography/ bathymetry.

nc : int

The number of components:

The first block holds the elevation of the topography/ bathymetry, so nc=1. The following blocks have either 3 if only rho, vp, and vs are present (attenuation=0) or 5 if qp and qs are pressent (attenuation=1).

ni, nj, nk : int

Number of grid points in the i, j, k directions.

Because the topography/bathymetry is (only) a function of the horizontal coordinates, the first block must have nk=1.

read_hdr(f)[source]

Read rfile header.

Parameters:

f : file

Open file handle in 'rb' mode of the rfile to be read.

Returns:

rfile header data (9 elements):

magic, precision, attenuation, az, lon0, lat0, mlen, proj_str, nb

magic : int

Determine byte ordering in file.

precision : int

The number of bytes per entry in the data section.

4 - single precision

8 - double precision

attenuation : int

Indicates whether the visco-elastic attenuation parameters QP and QS are included in the data section.

0 - no visco-elastic attenuation parameters included

1 - visco-elastic attenuation parameters included

az : numpy.float64

Angle in degrees between North and the positive x-axis.

lon0, lat0 : numpy.float64

Longitude and Latitude of the origin of the data.

mlen : int

The number of characters in the string proj_str.

proj_str : str

Projection string which is read by the Proj4 library if SW4 was built with Proj4 support.

nb : int

The number of blocks in the data section.

See also

write_hdr

Notes

See the SW4 User Guide for more details about these header parameters.

write_block_hdr(f, hh, hv, z0, nc, ni, nj, nk)[source]

Write rfile block header

Block headers are appended after the rfile header has been written. All block headers are written one after the other.

Parameters:

f : file

Open file handle in 'wa' mode.

hh, hv : numpy.float64

Grid size in the horizontal (x and y) and vertical (z) directions in meters.

z0 : numpy.float64

The base z-level of the block. Not used for the first block which holds the elevation of the topography/ bathymetry.

nc : int

The number of components:

The first block holds the elevation of the topography/ bathymetry, so nc=1. The following blocks must have either 3 if only rho, vp, and vs are present (attenuation=0) or 5 if qp and qs are pressent (attenuation=1).

ni, nj, nk : int

Number of grid points in the i, j, k directions.

Because the topography/bathymetry is (only) a function of the horizontal coordinates, the first block must have nk=1.

write_hdr(f, magic=1, precision=4, attenuation=1, az=0.0, lon0=33.5, lat0=28.0, proj_str=b'+proj=utm +zone=36 +datum=WGS84 +units=m +no_defs', nb=1)[source]

Write rfile header.

Parameters:

f : file

Open file handle in 'wb' mode

magic : int

Determine byte ordering in file. Defaults to 1.

precision : int

The number of bytes per entry in the data section.

4 - single precision (default)

8 - double precision

attenuation : int

Indicates whether the visco-elastic attenuation parameters QP and QS are included in the data section. 0 - no visco-elastic attenuation parameters included 1 - visco-elastic attenuation parameters included (default)

az : float

Angle in degrees between North and the positive x-axis. Defaults to 0. See the SW4 User Guide.

lon0, lat0 : float

Longitude and Latitude of the origin of the data. Defaults to 33.5, 28.0 .

proj_str : str

Projection string which is read by the Proj4 library if SW4 was built with Proj4 support. See the SW4 User Guide and the Proj4 documentation. Defaults to ‘+proj=utm +zone=36 +datum=WGS84 +units=m +no_defs’.

nb : int

The number of blocks in the data section. Must be > 0. Defaults to 1.

write_properties(f, vp, nc, vs=None, rho=None, qp=None, qs=None)[source]

Write material properties at a point i, j in block b.

This is a convenient function to use while looping over i, j in a specific block b for writing out material properties at each index k. At the very least vp should be provided. If only vp is provided, the other properties are calculated using the material_model module.

Parameters:

f : file

Open file handle in 'wa' mode for appending data to the end of a file in construction ot in 'r+b' mode for overwriting existing data.

When overwriting existing data, the user must take care to place the cursor in the right place in the file.

vp : array-like

P wave velocity at indices of k in m/s.

nc : int

Number of components to write out. Either 3 (rho, vp, and vs if attenuation=0) or 5 (also qp and qs if attenuation=1).

vs : array-like, optional

S wave velocity at indices of k in m/s. If not given, vs is calculated from get_vs().

rho : array-like, optional

Density at indices of k in kg/m^3. If not given, rho is calculated from get_rho().

qp : array-like, optional

P quality factor at indices of k. If not given, qp is calculated from get_qp().

qs : array-like, optional

S quality factor at indices of k. If not given, qs is calculated from get_qs().

write_topo_block(f, data, precision=4)[source]
Parameters:

f : file or str

Open file handle in 'wa' mode or path to an rfile for appending data to the end.

data : ndarray

Elevation data in meters below sea level (z positive down).

precision : {4 (default), 8}

The number of bytes per entry in the data section.

4 - single precision

8 - double precision