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)

Bases: object

A class to hold rfile block header and data section

_read_block_data_section(f, precision)

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()

Return values of Qp for block

qs()

Return values of Qs for block

rho()

Return values of density for block

vp()

Return values of vp for block

vs()

Return values of vs for block

x()

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

y()

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

z()

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)

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)

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)

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)

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')

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()

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')

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)

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)

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()

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)

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)

Read rfile into Model object.

See Model for documentation.

read_block_hdr(f)

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)

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)

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)

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)

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)

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