# -*- coding: utf-8 -*-
"""
Python module to handle SW4 images of Maps or Cross-Sections.
.. module:: image
:author:
Shahar Shani-Kadmiel (s.shanikadmiel@tudelft.nl)
Omry Volk (omryv@post.bgu.ac.il)
Tobias Megies (megies@geophysik.uni-muenchen.de)
:copyright:
Shahar Shani-Kadmiel (s.shanikadmiel@tudelft.nl)
Omry Volk (omryv@post.bgu.ac.il)
Tobias Megies (megies@geophysik.uni-muenchen.de)
:license:
This code is distributed under the terms of the
GNU Lesser General Public License, Version 3
(https://www.gnu.org/copyleft/lesser.html)
"""
from __future__ import absolute_import, print_function, division
from warnings import warn
import inspect
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import patheffects as path_effects
from mpl_toolkits.axes_grid1 import make_axes_locatable
from obspy.core.util import AttribDict
try:
from obspy.imaging.cm import obspy_divergent as cmap_divergent
from obspy.imaging.cm import obspy_divergent_r as cmap_divergent_r
from obspy.imaging.cm import obspy_sequential as cmap_sequential
from obspy.imaging.cm import obspy_sequential_r as cmap_sequential_r
except ImportError:
cmap_divergent = None
cmap_divergent_r = None
cmap_sequential = None
cmap_sequential_r = None
from ..sw4_metadata import Inputfile
from ..headers import (
IMAGE_HEADER_DTYPE, PATCH_HEADER_DTYPE, IMAGE_PLANE,
IMAGE_MODE_DISPLACEMENT, IMAGE_MODE_VELOCITY, IMAGE_PRECISION,
STF)
import copy
[docs]class Image():
"""
A class to hold SW4 image files.
Initialize an empty Image object, preferentially specifying the
input file used to run the simulation.
Parameters
----------
input_file : str or :class:`~obspy.core.util.attribdict.AttribDict`
Input file (already parsed or filename) used to compute the
image output.
stf : {'displacement', 'velocity'}
Only needed if no metadata from original input_file is used.
"""
CMAP = {"divergent" : cmap_divergent,
"divergent_r" : cmap_divergent_r,
"sequential" : cmap_sequential,
"sequential_r" : cmap_sequential_r}
MPL_SCATTER_PROPERTIES = {
"source" : {"s" : 200,
"marker" : "*",
"edgecolors" : "k",
"facecolors" : "",
"alpha" : 1,
"linewidths" : 1.5},
"rec" : {"s" : 200,
"marker" : "v",
"edgecolors" : "k",
"facecolors" : "",
"alpha" : 1,
"linewidths" : 1.5}
}
MPL_SCATTER_PATH_EFFECTS = [
path_effects.Stroke(linewidth=1.5 + 1, foreground='w'),
path_effects.Normal()]
def __init__(self, input_file=None, stf=None):
self.patches = []
if type(input_file) in (str, Inputfile):
self.input_file = Inputfile(input_file)
if not self.input_file:
warn('File {} was not parsed and might not '
'be an SW4 inputfile'.format(input_file))
self.input_file = None
else:
self.input_file = None
try:
stf_ = STF[self.input_file.source[0].type].type
except AttributeError:
stf_ = None
if (stf and stf_ and
stf != stf_):
msg = ('Overriding user specified source time function '
'type ({}) with the one found in input file '
'({}).')
warn(msg.format(stf, stf_))
self.stf = stf_
else:
self.stf = stf
# set mode code mapping, depending on the type of
# source time function
if self.stf == 'displacement':
self._mode_dict = IMAGE_MODE_DISPLACEMENT
elif self.stf == 'velocity':
self._mode_dict = IMAGE_MODE_VELOCITY
else:
msg = ("Unrecognized 'stf': '{}'")
raise ValueError(msg.format(self.stf))
self.filename = None
self._precision = None
self.number_of_patches = None
self.time = None
self._plane = None
self.coordinate = None
self._mode = None
self.gridinfo = None
self.extent = None
self.creation_time = None
[docs] def _read_patches(self, f):
"""
Read SW4 patch data and store it in a list of :class:`~.Patch`
objects under :obj:`~.Image.patches`.
Parameters
----------
f : file
Open file handle of SW4 image file (at correct position).
"""
patch_info = np.fromfile(
f, PATCH_HEADER_DTYPE, self.number_of_patches)
for i, header in enumerate(patch_info):
patch = Patch(number=i, image=self)
patch._set_header(header)
data = np.fromfile(f, self.precision, patch.ni * patch.nj)
data = data.reshape(patch.nj, patch.ni)
patch._set_data(data)
if (self.gridinfo and
patch.number == self.number_of_patches - 1):
patch.is_curvilinear = True
else:
patch.is_curvilinear = False
self.patches.append(patch)
if self.gridinfo == 1:
self._read_curvilinear_grid(f)
[docs] def _read_curvilinear_grid(self, f):
"""
Read the last bit of the SW4 image file in case a curvilinear
grid is found.
Parameters
----------
f : file
Open file handle of SW4 image file (at correct position).
"""
grid = Patch(number='curvilinear grid', image=self)
# get data from the corresponding image patch
patch = self.patches[-1]
grid.h = patch.h
grid.ib = patch.ib
grid.ni = patch.ni
grid.jb = patch.jb
grid.nj = patch.nj
# read data from file
data = np.fromfile(f, self.precision, grid.ni * grid.nj)
data = data.reshape(grid.nj, grid.ni)
zmin = data.min()
grid.zmin = zmin
grid._set_data(data)
# update curvilinear patch and grid extent
extent = list(grid.extent)
extent[2] = zmin - (grid.h / 2.0)
extent[3] = grid._max + (grid.h / 2.0)
extent = tuple(extent)
grid.extent = extent
self.patches[-1].zmin = zmin
self.patches[-1].extent = extent
self.curvilinear_grid_patch = grid
[docs] def _calc_global_min_max(self):
"""
Calculate ``min``, ``max``, ``rms``, and ``ptp``.
"""
_max = []
_min = []
_rms = []
for patch in self.patches:
_max += [patch._max]
_min += [patch._min]
_rms += [patch._rms]
self._max = max(_max)
self._min = min(_min)
self._rms = max(_rms)
self._ptp = self._max - self._min
if self.is_cross_section:
x1, x2 = self.patches[-1].extent[:2]
y1 = min([patch.extent[2] for patch in self.patches])
y2 = max([patch.extent[3] for patch in self.patches])
self.extent = (x1, x2, y1, y2)
else:
self.extent = self.patches[-1].extent
[docs] def plot(self, patches=None, ax=None, vmin='min', vmax='max',
colorbar=True, cmap=None, interpolation='nearest',
origin='lower', projection_distance=np.inf, **kwargs):
"""
Plot all (or specific) patches in :class:`~.Image`.
Parameters
----------
patches : list of int
Patches to plot
Other Parameters
----------------
ax : :class:`~matplotlib.axes.Axes`
Use existing axes.
vmin : float
Manually set minimum of color scale.
vmax : str or float
Used to clip the coloring of the data at the set value.
Default is 'max' which shows all the data. If ``float``,
the colorscale saturates at the given value. Finally, if a
string is passed (other than 'max'), it is casted to float
and used as an ``rms`` multiplier. For instance, if
``vmax='3'``, clipping is done at 3.0\*rms of the data.
colorbar : bool or str
If ``colorbar`` is a string, that string is used to override
the automatic label chosen based on the image header. To
Supress plotting of the colorbar set to ``False``.
cmap : :class:`~matplotlib.colors.Colormap`
Colormap for the plot
interpolation : str
Acceptable values are 'none', 'nearest', 'bilinear',
'bicubic', 'spline16', 'spline36', 'hanning', 'hamming',
'hermite', 'kaiser', 'quadric', 'catrom', 'gaussian',
'bessel', 'mitchell', 'sinc', 'lanczos'.
origin : str
Places the origin at the 'lower' (default)
or 'upper' left corner of the plot.
projection_distance : float
Threshold distance from the plane coordinate to include
symbols of stations and sources. These are orthogonally
*projected* onto the plotted 2D plane. By default everything
is included but this can cause too many symbols to be
plotted, obscuring the image.
Examples
--------
>>> my_image.plot() # plots all patches
>>> my_image.plot(patches=[0, 2]) # plots first and third patch
"""
if ax is None:
fig, ax = plt.subplots()
else:
fig = None
# set vmin, vmax
if vmax is None or vmax is 'max':
clip = self._max
elif type(vmax) in [float, int]:
clip = vmax
else:
clip = float(vmax) * self._rms
if vmin is None or vmin is 'min':
vmin = self._min
elif type(vmin) in [int, float]:
pass
elif vmin in [None, False]:
vmin = -clip
vmax = clip
# plot patches
if patches is None:
for patch in self.patches:
im = patch.plot(
ax=ax, vmin=vmin, vmax=vmax, colorbar=False,
cmap=cmap, interpolation=interpolation, origin=origin,
**kwargs)
else:
for i in patches:
im = self.patches[i].plot(
ax=ax, vmin=vmin, vmax=vmax, colorbar=False,
cmap=cmap, interpolation=interpolation, origin=origin,
**kwargs)
# plot the colorbar
if colorbar:
# set colorbar extend mode
if vmin > self._min and vmax < self._max:
extend = 'both'
elif vmin > self._min:
extend = 'min'
elif vmax < self._max:
extend = 'max'
else:
extend = 'neither'
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="3%", pad=0.1)
if type(colorbar) is str:
colorbar_label = colorbar
else:
colorbar_label = "{name} [{unit}]".format(
name=self.quantity_name,
unit=self.quantity_unit)
cb = plt.colorbar(im, cax=cax, extend=extend,
label=colorbar_label)
# invert Z axis for cross-section plots and certain
# quantities that usually increase with depths
if (self.is_cross_section and
self._mode in (4, 7, 8)):
cax.invert_yaxis()
else:
cb = None
# labels
if self._plane == 0:
xlabel = "Y"
ylabel = "Z"
elif self._plane == 1:
xlabel = "X"
ylabel = "Z"
elif self._plane == 2:
xlabel = "Y"
ylabel = "X"
ax.set_xlabel(xlabel + ' [m]')
ax.set_ylabel(ylabel + ' [m]')
title = self.filename.rsplit('.', 1)[0]
if len(title) > 40:
title = '...' + title[-40:]
ax.set_title("{}\n{}={} t={:.2f} seconds".format(
title, self.plane, self.coordinate,
self.time), y=1.03, fontsize=12)
# plot receivers, source etc.
if self.input_file:
try: # get surface elevation to correct depth coordinate
xi = self.curvilinear_grid_patch.x
elev = self.curvilinear_grid_patch.data[0]
except AttributeError: # no curvilinear grid... no correction
xi = None
elev = None
for key, kwargs in self.MPL_SCATTER_PROPERTIES.items():
coordinates = self.input_file.get_coordinates(
key, xi, elev,
self._plane, self.coordinate, projection_distance)
if not coordinates:
continue
x, y = coordinates
collection = ax.scatter(
x, y, **kwargs)
collection.set_path_effects(
self.MPL_SCATTER_PATH_EFFECTS)
# reset axes limits to data, they sometimes get changed by the
# scatter plot
ax.axis(self.extent)
# invert Z axis if not a map view
if self.is_cross_section:
ax.invert_yaxis()
return fig, ax, cb
def get_source_coordinates(self):
try:
return self._get_plot_coordinates_from_input("source")
except KeyError:
return None
@property
def is_cross_section(self):
if self._plane in (0, 1):
return True
elif self._plane == 2:
return False
@property
def is_divergent(self):
if self._cmap_type in ("divergent", "divergent_r"):
return True
return False
@property
def _cmap_type(self):
try:
return self._mode_dict[self._mode]['cmap_type']
except KeyError:
return None
@property
def precision(self):
try:
return IMAGE_PRECISION[self._precision]
except KeyError:
return None
@property
def plane(self):
try:
return IMAGE_PLANE[self._plane]
except KeyError:
return None
@property
def type(self):
if self.is_cross_section:
return 'cross-section'
else:
return 'map'
@property
def quantity_name(self):
try:
return self._mode_dict[self._mode]['name']
except KeyError:
return None
@property
def quantity_altname(self):
try:
return self._mode_dict[self._mode]['altname']
except KeyError:
return None
@property
def quantity_symbol(self):
try:
return self._mode_dict[self._mode]['symbol']
except KeyError:
return None
@property
def quantity_altsymbol(self):
try:
return self._mode_dict[self._mode]['altsymbol']
except KeyError:
return None
@property
def quantity_unit(self):
try:
return self._mode_dict[self._mode]['unit']
except KeyError:
return None
def __str__(self):
string = (
' Image information :\n'
' ----------------- :\n'
' Filename : {}\n'
' Precision : {}\n'
' Number of patches : {}\n'
' Time, s : {}\n'
' Plane : {}\n'
' Coordinate : {}\n'
' Mode : {}, {}\n'
'Curvilinear grid included : {}\n'
' Image extent : {}\n'
' Creation time : {}\n'
).format(
str(self.filename).rsplit('/', 1)[-1], self.precision,
self.number_of_patches, self.time,
self.plane, self.coordinate, self.quantity_symbol,
self.quantity_unit, True if self.gridinfo else False,
self.extent, self.creation_time
)
return string
[docs] def copy(self):
"""
Return a deepcopy of `self`.
"""
return copy.deepcopy(self)
[docs]class Patch():
"""
A class to hold SW4 patch data.
Initialize an empty Patch object, preferentially specifying the
parent :class:`.Image`.
Parameters
----------
image : :class:`.Image`
Parent Image object.
number : int
Patch index in parent image (starts at `0`).
"""
def __init__(self, image=None, number=None):
self.number = number
self._image = image # link back to the image this patch belongs to
self.h = None
self.zmin = None
self.ib = None
self.ni = None
self.jb = None
self.nj = None
self.data = None
self.extent = None
def _set_data(self, data):
self.data = data
if self._image.is_cross_section:
self.extent = (
0 - (self.h / 2.0),
(self.ni - 1) * self.h + (self.h / 2.0),
self.zmin - (self.h / 2.0),
self.zmin + (self.nj - 1) * self.h + (self.h / 2.0))
else:
self.data = self.data.T
self.extent = (
0 - (self.h / 2.0),
(self.nj - 1) * self.h + (self.h / 2.0),
0 - (self.h / 2.0),
(self.ni - 1) * self.h + (self.h / 2.0))
self._min = data.min()
self._max = data.max()
self._std = data.std()
self._rms = np.sqrt(np.mean(data**2))
self.shape = data.shape
[docs] def plot(self, ax=None, vmin=None, vmax=None, colorbar=True, cmap=None,
interpolation='nearest', origin='lower', **kwargs):
"""
Plot patch.
.. note:: Should not really be used directly by the user but
rather called by the :meth:`~.Image.plot` method.
Parameters
----------
ax : :class:`~matplotlib.axes.Axes`
Use existing axes.
vmin : float
Manually set minimum of color scale.
vmax : float
Manually set maximum of color scale.
colorbar : bool or str
If ``colorbar`` is a string, that string is used to override
the automatic label chosen based on the image header. To
Supress plotting of the colorbar set to ``False``.
cmap : :class:`~matplotlib.colors.Colormap`
Colormap for the plot
interpolation : str
Acceptable values are 'none', 'nearest', 'bilinear',
'bicubic', 'spline16', 'spline36', 'hanning', 'hamming',
'hermite', 'kaiser', 'quadric', 'catrom', 'gaussian',
'bessel', 'mitchell', 'sinc', 'lanczos'.
origin : str
Places the origin at the 'lower' (default)
or 'upper' left corner of the plot.
Other Parameters
----------------
kwargs : :func:`~matplotlib.pyplot.imshow`
"""
caller = inspect.stack()[1][3] # find out who made the call
if not ax:
fig, ax = plt.subplots()
else:
fig = None
if not cmap:
cmap = self._image.CMAP[self._image._cmap_type]
ax.set_aspect(1)
if self.is_curvilinear:
y = self._image.curvilinear_grid_patch.data
x = (np.ones_like(y) * self.x)
im = ax.pcolormesh(x, y, self.data,
shading='flat', edgecolors='None',
vmin=vmin, vmax=vmax, cmap=cmap)
ax.axis(self.extent)
elif not self.is_curvilinear:
im = ax.imshow(self.data, extent=self.extent,
vmin=vmin, vmax=vmax, origin=origin,
interpolation=interpolation, cmap=cmap, **kwargs)
# invert Z axis if not a map view
if self._image.is_cross_section:
ax.invert_yaxis()
if caller is not 'plot': # do this only if the user calls
if colorbar:
if vmin > self._min and vmax < self._max:
extend = 'both'
elif vmin > self._min:
extend = 'min'
elif vmax < self._max:
extend = 'max'
else:
extend = 'neither'
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="3%", pad=0.1)
if type(colorbar) is str:
colorbar_label = colorbar
else:
colorbar_label = "{name} [{unit}]".format(
name=self._image.quantity_name,
unit=self._image.quantity_unit)
cb = plt.colorbar(im, cax=cax, extend=extend,
label=colorbar_label)
# invert Z axis for cross-section plots and certain
# quantities that usually increase with depths
if self._image.is_cross_section and \
self._image._mode in (4, 7, 8):
cax.invert_yaxis()
else:
cb = None
if self._image._plane == 0:
xlabel = "Y"
ylabel = "Z"
elif self._image._plane == 1:
xlabel = "X"
ylabel = "Z"
elif self._image._plane == 2:
xlabel = "Y"
ylabel = "X"
ax.set_xlabel(xlabel + ' [m]')
ax.set_ylabel(ylabel + ' [m]')
# setup the title so that it is not too long
# I guess for some purposes you would want the
# entire path to show so we might want to just
# brake it into lines?
# For now I am truncating it at -40: chars and
# getting rid of the .sw4img extention... That's
# kind of obvious.
title = self._image.filename.rsplit('.', 1)[0]
if len(title) > 40:
title = ('...' + title[-40:]
+ ' : patch no. {}'.format(self.number))
ax.set_title("{}\n{}={} t={:.2f} seconds".format(
title, self._image.plane, self._image.coordinate,
self._image.time), y=1.03, fontsize=12)
return fig, ax, cb
elif caller is 'plot': # do this only if Image.plot() calls
return im
def __str__(self):
string = (
'Patch information :\n'
'----------------- :\n'
' Number : {}\n'
' Spacing, m : {}\n'
' Zmin, m : {}\n'
' Extent : {}\n'
' ni : {}\n'
' nj : {}\n'
' Max : {}\n'
' Min : {}\n'
' STD : {}\n'
' RMS : {}\n'
).format(
self.number, self.h, self.zmin, self.extent, self.ni, self.nj,
self._max, self._min, self._std, self._rms
)
return string
@property
def x(self):
return np.linspace(*self.extent[:2], num=self.ni)
[docs] def copy(self):
"""
Return a deepcopy of `self`.
"""
return copy.deepcopy(self)
[docs]def read_image(filename='random', input_file=None,
stf="displacement", verbose=False):
"""
Read image data, cross-section or map into a
:class:`.Image` object.
Parameters
----------
filename : str
If no filename is passed, by default, a random image is
generated. if filename is ``None``, an empty :class:`.Image`
object is returned.
input_file : str or AttribDict
Input file (already parsed or filename) used to compute the
image output.
stf : str
`'displacement'` or `'velocity'`. Only needed if no metadata
from original input_file is used.
verbose : bool
If set to ``True``, print some information while reading the
file.
Returns
-------
:class:`.Image`
An :class:`~.Image` object with a list of :class:`~.Patch`
objects.
"""
image = Image(stf=stf,
input_file=input_file)
image.filename = filename
if filename is 'random': # generate random data, populate objects
image = _create_random_image(
stf=stf)
elif filename is None:
pass
else:
if not filename.endswith('.sw4img'):
msg = ("Using 'read_image()' on file with uncommon file "
"extension: '{}'.").format(filename)
warn(msg)
with open(image.filename, 'rb') as f:
image._read_header(f)
image._read_patches(f)
image._calc_global_min_max()
return image
def _create_random_image(stf="displacement"):
image = Image(stf=stf)
image.filename = None
image.number_of_patches = 1
image._precision = 4
image.cycle = 0
image.time = 0
image._min = 0
image._max = 0
image._rms = 0
image._ptp = 0
image.patches = [_create_random_patch()]
return image
def _create_random_patch():
patch = Patch()
patch.ni = 100
patch.nj = 200
patch.data = 2 * (np.random.rand(patch.ni, patch.nj) - 0.5)
patch.number = 0
patch.h = 100.0
patch.zmin = 0
patch.extent = (0, patch.nj * patch.h, patch.zmin + patch.ni * patch.h,
patch.zmin)
patch._min = patch.data.min()
patch._max = patch.data.max()
patch._rms = np.sqrt(np.mean(patch.data**2))
patch.ib = 1
patch.jb = 1
return patch