You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
705 lines
24 KiB
705 lines
24 KiB
from __future__ import print_function |
|
from builtins import zip |
|
from builtins import range |
|
|
|
import numpy as num |
|
import time |
|
from matplotlib.animation import FuncAnimation |
|
|
|
from pyrocko.gui.snuffling import Param, Snuffling, Choice, Switch |
|
from scipy.signal import fftconvolve, lfilter, hilbert |
|
from scipy.interpolate import UnivariateSpline |
|
from pyrocko import orthodrome as ortho |
|
from pyrocko import parstack |
|
from pyrocko import util |
|
from pyrocko import trace |
|
from pyrocko import model |
|
import logging |
|
|
|
|
|
logger = logging.getLogger('pyrocko.gui.snufflings.fk_parstack.py') |
|
d2r = num.pi/180. |
|
km = 1000. |
|
|
|
|
|
def search_max_block(n_maxsearch, data): |
|
''' |
|
Find indices of maxima of *data* in groups of *n_maxsearch* |
|
samples. |
|
|
|
returns an array of indices. |
|
|
|
If length of *data* is not a multiple of *n_maxsearch*, *data* will be |
|
padded. |
|
''' |
|
n = len(data) |
|
n_dim2 = (int(n/n_maxsearch)+1) * n_maxsearch |
|
n_missing = n_dim2 - n |
|
if n % n_maxsearch != 0: |
|
a = num.pad(data, [0, n_missing], mode='minimum') |
|
else: |
|
a = data |
|
|
|
return num.argmax(a.reshape((-1, n_maxsearch)), axis=1) +\ |
|
num.arange(0, (n+n_missing)/n_maxsearch) * n_maxsearch |
|
|
|
|
|
def instantaneous_phase(signal): |
|
analytic_signal = hilbert(signal) |
|
return num.unwrap(num.angle(analytic_signal)) |
|
|
|
|
|
def get_instantaneous_frequency(signal, fs): |
|
inst_phase = instantaneous_phase(signal) |
|
return (num.diff(inst_phase) / (2.0*num.pi) * fs) |
|
|
|
|
|
def get_center_station(stations, select_closest=False): |
|
''' gravitations center of *stations* list.''' |
|
n = len(stations) |
|
slats = num.empty(n) |
|
slons = num.empty(n) |
|
for i, s in enumerate(stations): |
|
slats[i] = s.lat |
|
slons[i] = s.lon |
|
|
|
center_lat, center_lon = ortho.geographic_midpoint(slats, slons) |
|
|
|
if select_closest: |
|
center_lats = num.ones(n)*center_lat |
|
center_lons = num.ones(n)*center_lon |
|
dists = ortho.distance_accurate50m_numpy( |
|
center_lats, center_lons, slats, slons) |
|
|
|
return stations[num.argmin(dists)] |
|
else: |
|
return model.Station(center_lat, center_lon) |
|
|
|
|
|
def get_theoretical_backazimuth(event, stations, center_station): |
|
return (ortho.azimuth_numpy( |
|
event.lat, event.lon, center_station.lat, center_station.lon) |
|
+ 180.) % 360. |
|
|
|
|
|
def get_shifts(stations, center_station, bazis, slownesses): |
|
''' shape = (len(bazi), len(slow))''' |
|
lats = num.array([s.lat for s in stations]) |
|
lons = num.array([s.lon for s in stations]) |
|
lat0 = num.array([center_station.lat] * len(stations)) |
|
lon0 = num.array([center_station.lon] * len(stations)) |
|
ns, es = ortho.latlon_to_ne_numpy(lat0, lon0, lats, lons) |
|
station_vector = num.array((ns, es)).T |
|
bazis = bazis * d2r |
|
shifts = num.zeros((len(bazis)*len(slownesses), len(stations))) |
|
ishift = 0 |
|
for ibazi, bazi in enumerate(bazis): |
|
s_vector = num.array((num.cos(bazi), num.sin(bazi))) |
|
for islowness, slowness in enumerate(slownesses): |
|
shifts[ishift] = station_vector.dot(s_vector) * slowness |
|
ishift += 1 |
|
|
|
return shifts |
|
|
|
|
|
def to_db(d): |
|
return 10*num.log10(d/num.max(d)) |
|
|
|
|
|
def lowpass_array(ydata_array, deltat, order, corner, demean=True, axis=1): |
|
''' |
|
Apply butterworth highpass to the trace. |
|
|
|
:param order: order of the filter |
|
:param corner: corner frequency of the filter |
|
|
|
Mean is removed before filtering. |
|
''' |
|
|
|
(b, a) = trace._get_cached_filter_coefs( |
|
order, [corner*2.0*deltat], btype='low') |
|
data = ydata_array.astype(num.float64) |
|
if len(a) != order+1 or len(b) != order+1: |
|
logger.warn( |
|
'Erroneous filter coefficients returned by ' |
|
'scipy.signal.butter(). You may need to downsample the ' |
|
'signal before filtering.') |
|
if demean: |
|
data -= num.mean(data, axis=1)[None].T |
|
|
|
return lfilter(b, a, data) |
|
|
|
|
|
def highpass_array(ydata_array, deltat, order, corner, demean=True, axis=1): |
|
''' |
|
Apply butterworth highpass to the trace. |
|
|
|
:param order: order of the filter |
|
:param corner: corner frequency of the filter |
|
|
|
Mean is removed before filtering. |
|
''' |
|
(b, a) = trace._get_cached_filter_coefs( |
|
order, [corner*2.0*deltat], btype='high') |
|
data = ydata_array.astype(num.float64) |
|
if len(a) != order+1 or len(b) != order+1: |
|
logger.warn( |
|
'Erroneous filter coefficients returned by ' |
|
'scipy.signal.butter(). You may need to downsample the ' |
|
'signal before filtering.') |
|
if demean: |
|
data -= num.mean(data, axis=1)[None].T |
|
|
|
return lfilter(b, a, data) |
|
|
|
|
|
def value_to_index(value, range_min, range_max, range_delta, clip=True): |
|
''' map a(n array of) *values* to its' index in a continuous data range |
|
defined by *range_min*, *range_max* and *range_delta*. |
|
''' |
|
indices = num.round((value-range_min)/range_delta) |
|
if clip: |
|
indices = num.clip(indices, 0, (range_max-range_min)/range_delta) |
|
|
|
return num.asarray(indices, dtype=num.int) |
|
|
|
|
|
class FK(Snuffling): |
|
''' |
|
<html> |
|
<head> |
|
<style type="text/css"> |
|
body { margin-left:10px }; |
|
</style> |
|
</head> |
|
<body> |
|
<h1 align='center'>FK ANALYSIS</h1> |
|
<p> |
|
Performs delay and sum in the time domain. |
|
<u>Usage</u><br> |
|
|
|
- Load station information at startup <br> |
|
|
|
- Zoom into the data until you see only data you desire to analyse or |
|
use extended markers to selected time regions for analysis<br> |
|
|
|
- Press the 'Run' button <br> |
|
|
|
</p> |
|
The slowness is in units s/km. |
|
<p> |
|
if <b>Show</n> is activated, three images will be genereated for each |
|
processing block: A polar plot which shows the maximum coherence found |
|
along the time axis for all slownesses and back-azimuths. The blue dot |
|
within that figure indicates the position of the maximum. Based on that |
|
slowness/back-azimuth the other two coherence maps are generated which |
|
show the coherence in the slowness and back-azimuth domain for that |
|
specific maximum of that processing block.<br> |
|
|
|
Picinbono, et. al, 1997, On Instantaneous Amplitude and Phase of Signals, |
|
552 IEEE TRANSACTIONS ON SIGNAL PROCESSING, 45, 3, March 1997 |
|
</body> |
|
</html> |
|
''' |
|
|
|
def setup(self): |
|
self.set_name('FK (parstack)') |
|
self.add_parameter(Param( |
|
'max slowness [s/km]', 'slowness_max', 0.2, 0., 1.)) |
|
self.add_parameter(Param( |
|
'min slowness [s/km]', 'slowness_min', 0.01, 0., 1.)) |
|
self.add_parameter(Param( |
|
'delta slowness [s/km]', 'slowness_delta', 0.002, 0., 0.2)) |
|
self.add_parameter(Param( |
|
'delta backazimut', 'delta_bazi', 2, 1, 20)) |
|
self.add_parameter(Param( |
|
'Increment [s]', 'tinc', 60., 0.5, 60., |
|
high_is_none=True)) |
|
self.add_parameter(Param( |
|
'Smoothing length [N]', 'ntaper', 0, 0, 30, low_is_none=True)) |
|
self.add_parameter(Param( |
|
'Maximum search factor', 'search_factor', 1, 0, 3)) |
|
self.add_parameter(Choice( |
|
'Use channels', 'want_channel', '*', |
|
['*', '*Z', '*E', '*N', 'SHZ', 'BHZ', 'p0'])) |
|
self.add_parameter( |
|
Choice('method', 'method', 'stack', ['stack', 'correlate']) |
|
) |
|
self.add_parameter(Switch('Show', 'want_all', True)) |
|
self.add_parameter(Switch('Phase weighted stack', 'want_pws', False)) |
|
self.set_live_update(False) |
|
self.irun = 0 |
|
self.figs2draw = [] |
|
|
|
def new_figure(self, title=''): |
|
'''Return a new Figure instance''' |
|
fig_frame = self.pylab(name='FK: %s (%i)' % |
|
(title, self.irun), get='figure_frame') |
|
self.figs2draw.append(fig_frame.gcf()) |
|
return self.figs2draw[-1] |
|
|
|
def draw_figures(self): |
|
''' Draw all new figures and clear list.''' |
|
for fig in self.figs2draw: |
|
fig.canvas.draw() |
|
|
|
self.figs2draw = [] |
|
|
|
def call(self): |
|
|
|
self.cleanup() |
|
figs = [] |
|
azi_theo = None |
|
method = {'stack': 0, |
|
'correlate': 2}[self.method] |
|
|
|
bazis = num.arange(0., 360.+self.delta_bazi, self.delta_bazi) |
|
slownesses = num.arange(self.slowness_min/km, |
|
self.slowness_max/km, |
|
self.slowness_delta/km) |
|
n_bazis = len(bazis) |
|
n_slow = len(slownesses) |
|
|
|
viewer = self.get_viewer() |
|
event = viewer.get_active_event() |
|
|
|
stations = self.get_stations() |
|
stations_dict = dict(zip([viewer.station_key(s) for s in stations], |
|
stations)) |
|
|
|
traces_pile = self.get_pile() |
|
deltats = traces_pile.deltats.keys() |
|
if len(deltats) > 1: |
|
self.fail('sampling rates differ in dataset') |
|
else: |
|
deltat_cf = deltats[0] |
|
|
|
tinc_use = self.get_tinc_use(precision=deltat_cf) |
|
|
|
if self.ntaper: |
|
taper = num.hanning(int(self.ntaper)) |
|
else: |
|
taper = None |
|
|
|
frames = None |
|
t1 = time.time() |
|
|
|
# make sure that only visible stations are used |
|
use_stations = stations |
|
center_station = get_center_station(use_stations, select_closest=True) |
|
print('Center station: ', center_station) |
|
shift_table = get_shifts( |
|
stations=use_stations, |
|
center_station=center_station, |
|
bazis=bazis, |
|
slownesses=slownesses) |
|
|
|
shifts = num.round(shift_table / deltat_cf).astype(num.int32) |
|
|
|
# padding from maximum shift of traces: |
|
npad = num.max(num.abs(shifts)) |
|
tpad = npad * deltat_cf |
|
|
|
# additional padding for cross over fading |
|
npad_fade = 0 |
|
tpad_fade = npad_fade * deltat_cf |
|
|
|
npad += npad_fade |
|
tpad += tpad_fade |
|
|
|
frames = None |
|
tinc_add = tinc_use or 0 |
|
|
|
def trace_selector(x): |
|
return util.match_nslc('*.*.*.%s' % self.want_channel, x.nslc_id) |
|
|
|
for traces in self.chopper_selected_traces( |
|
tinc=tinc_use, tpad=tpad, fallback=True, |
|
want_incomplete=False, trace_selector=trace_selector): |
|
|
|
if len(traces) == 0: |
|
self.fail('No traces matched') |
|
continue |
|
|
|
# should be correct |
|
t_min = traces[0].tmin |
|
t_max = traces[0].tmax |
|
|
|
use_stations = [] |
|
for tr in traces: |
|
try: |
|
use_stations.append(stations_dict[viewer.station_key(tr)]) |
|
except KeyError: |
|
self.fail('no trace %s' % ('.'.join(tr.nslc_id))) |
|
|
|
shift_table = get_shifts( |
|
stations=use_stations, |
|
center_station=center_station, |
|
bazis=bazis, |
|
slownesses=slownesses) |
|
|
|
shifts = num.round(shift_table / deltat_cf).astype(num.int32) |
|
|
|
wmin = traces[0].tmin |
|
wmax = wmin + tinc_add |
|
|
|
iwmin = int(round((wmin-wmin) / deltat_cf)) |
|
iwmax = int(round((wmax-wmin) / deltat_cf)) |
|
lengthout = iwmax - iwmin |
|
arrays = num.zeros((len(traces), lengthout + npad*2)) |
|
|
|
for itr, tr in enumerate(traces): |
|
tr = tr.copy() |
|
if viewer.highpass: |
|
tr.highpass(4, viewer.highpass, demean=True) |
|
else: |
|
tr.ydata = num.asarray( |
|
tr.ydata, dtype=num.float) - num.mean(tr.ydata) |
|
if viewer.lowpass: |
|
tr.lowpass(4, viewer.lowpass) |
|
|
|
arrays[itr] = tr.get_ydata() |
|
|
|
# if viewer.highpass: |
|
# arrays = highpass_array( |
|
# arrays, deltat_cf, 4, viewer.highpass) |
|
# if viewer.lowpass: |
|
# arrays = lowpass_array(arrays, deltat_cf, 4, viewer.lowpass) |
|
|
|
_arrays = [] |
|
for itr, tr in enumerate(traces): |
|
if taper is not None: |
|
ydata = fftconvolve(arrays[itr], taper, mode='same') |
|
else: |
|
ydata = arrays[itr] |
|
_arrays.append(num.asarray(ydata, dtype=num.float64)) |
|
arrays = _arrays |
|
|
|
offsets = num.array( |
|
[int(round((tr.tmin-wmin) / deltat_cf)) for tr in traces], |
|
dtype=num.int32) |
|
|
|
ngridpoints = len(bazis)*len(slownesses) |
|
weights = num.ones((ngridpoints, len(traces))) |
|
|
|
frames, ioff = parstack.parstack( |
|
arrays, offsets, shifts, weights, method, |
|
offsetout=iwmin, |
|
lengthout=lengthout, |
|
result=frames, |
|
impl='openmp') |
|
|
|
# theoretical bazi |
|
if event is not None: |
|
azi_theo = get_theoretical_backazimuth( |
|
event, use_stations, center_station) |
|
print('theoretical azimuth %s degrees' % (azi_theo)) |
|
|
|
print('processing time: %s seconds' % (time.time()-t1)) |
|
|
|
if frames is None: |
|
self.fail('Could not process data!') |
|
return |
|
|
|
frames_reshaped = frames.reshape((n_bazis, n_slow, lengthout)) |
|
times = num.linspace(t_min-tpad_fade, t_max+tpad_fade, lengthout) |
|
max_powers = num.max(frames, axis=0) |
|
|
|
# power maxima in blocks |
|
i_max_blocked = search_max_block( |
|
n_maxsearch=int(npad*self.search_factor), data=max_powers) |
|
|
|
max_powers += (num.min(max_powers)*-1) |
|
max_powers /= num.max(max_powers) |
|
max_powers *= max_powers |
|
weights = max_powers[i_max_blocked] |
|
block_max_times = times[i_max_blocked] |
|
|
|
_argmax = num.argmax(frames, axis=0) |
|
imax_bazi_all, imax_slow_all = num.unravel_index( |
|
_argmax, dims=(n_bazis, n_slow)) |
|
|
|
local_max_bazi = bazis[imax_bazi_all][i_max_blocked] |
|
local_max_slow = slownesses[imax_slow_all][i_max_blocked]*km |
|
|
|
k_north = num.sin(local_max_bazi * d2r) * local_max_slow |
|
k_east = num.cos(local_max_bazi * d2r) * local_max_slow |
|
|
|
smooth = 4e7 |
|
|
|
spline_north = UnivariateSpline( |
|
block_max_times, k_north, w=weights, |
|
s=smooth |
|
) |
|
|
|
spline_east = UnivariateSpline( |
|
block_max_times, k_east, w=weights, |
|
s=smooth, |
|
) |
|
|
|
k_north_fit = spline_north(times) |
|
k_east_fit = spline_east(times) |
|
|
|
bazi_fitted = num.arctan2(k_east_fit, k_north_fit) / d2r |
|
bazi_fitted -= 90. |
|
bazi_fitted *= -1. |
|
bazi_fitted[num.where(bazi_fitted<0.)] += 360. |
|
|
|
spline_slow = UnivariateSpline( |
|
block_max_times, |
|
local_max_slow, |
|
w=weights, |
|
) |
|
|
|
slow_fitted = spline_slow(times) |
|
i_bazi_fitted = value_to_index( |
|
bazi_fitted, 0., 360., self.delta_bazi) |
|
|
|
i_slow_fitted = value_to_index( |
|
slow_fitted, self.slowness_min, self.slowness_max, |
|
self.slowness_delta) |
|
|
|
i_shift = num.ravel_multi_index( |
|
num.vstack((i_bazi_fitted, i_slow_fitted)), |
|
(n_bazis, n_slow), |
|
) |
|
|
|
stack_trace = num.zeros(lengthout) |
|
i_base = num.arange(lengthout, dtype=num.int) + npad |
|
for itr, tr in enumerate(traces): |
|
isorting = num.clip( |
|
i_base-shifts[i_shift, itr], npad, lengthout+npad) |
|
stack_trace += tr.ydata[isorting] |
|
|
|
beam_tr = trace.Trace( |
|
tmin=t_min+tpad, ydata=stack_trace, deltat=deltat_cf) |
|
|
|
self.add_trace(beam_tr) |
|
|
|
if self.want_all: |
|
|
|
# --------------------------------------------------------- |
|
# maxima search |
|
# --------------------------------------------------------- |
|
fig1 = self.new_figure('Max Power') |
|
nsubplots = 1 |
|
ax = fig1.add_subplot(nsubplots, 1, 1) |
|
ax.plot(num.max(frames, axis=0)) |
|
# -------------------------------------------------------------- |
|
# coherence maps |
|
# -------------------------------------------------------------- |
|
|
|
max_time = num.amax(frames, axis=0) |
|
imax_time = num.argmax(max_time) |
|
best_frame = num.amax(frames, axis=1) |
|
imax_bazi_slow = num.argmax(best_frame) |
|
imax_bazi, imax_slow = num.unravel_index( |
|
num.argmax(best_frame), |
|
dims=(n_bazis, n_slow)) |
|
|
|
fig2 = self.new_figure('Slowness') |
|
data = frames_reshaped[imax_bazi, :, :] |
|
data_max = num.amax(frames_reshaped, axis=0) |
|
|
|
ax = fig2.add_subplot(211) |
|
ax.set_title('Global maximum slize') |
|
ax.set_ylabel('slowness [s/km]') |
|
ax.plot(times[imax_time], slownesses[imax_slow]*km, 'b.') |
|
ax.pcolormesh(times, slownesses*km, data) |
|
|
|
ax = fig2.add_subplot(212, sharex=ax, sharey=ax) |
|
ax.set_ylabel('slowness [s/km]') |
|
ax.pcolormesh(times, slownesses*km, data_max) |
|
ax.set_title('Maximum') |
|
|
|
# highlight block maxima |
|
ax.plot(block_max_times, local_max_slow, 'wo') |
|
|
|
ax.plot(times, num.clip( |
|
slow_fitted, self.slowness_min, self.slowness_max) |
|
) |
|
|
|
fig3 = self.new_figure('Back-Azimuth') |
|
data = frames_reshaped[:, imax_slow, :] |
|
data_max = num.amax(frames_reshaped, axis=1) |
|
|
|
ax = fig3.add_subplot(211, sharex=ax) |
|
ax.set_title('Global maximum slize') |
|
ax.set_ylabel('back-azimuth') |
|
ax.pcolormesh(times, bazis, data) |
|
ax.plot(times[imax_time], bazis[imax_bazi], 'b.') |
|
|
|
ax = fig3.add_subplot(212, sharex=ax, sharey=ax) |
|
ax.set_ylabel('back-azimuth') |
|
ax.set_title('Maximum') |
|
ax.pcolormesh(times, bazis, data_max) |
|
|
|
# highlight block maxima |
|
ax.plot(block_max_times, local_max_bazi, 'wo') |
|
ax.plot(times, num.clip(bazi_fitted, 0, 360.)) |
|
|
|
# xfmt = md.DateFormatter('%Y-%m-%d %H:%M:%S') |
|
# ax.xaxis.set_major_formatter(xfmt) |
|
# fig.autofmt_xdate() |
|
# fig.subplots_adjust(hspace=0) |
|
|
|
semblance = best_frame.reshape((n_bazis, n_slow)) |
|
|
|
fig4 = self.new_figure('Max') |
|
theta, r = num.meshgrid(bazis, slownesses) |
|
theta *= (num.pi/180.) |
|
|
|
ax = fig4.add_subplot(111, projection='polar') |
|
m = ax.pcolormesh(theta.T, r.T*km, to_db(semblance)) |
|
|
|
ax.plot(bazis[imax_bazi]*d2r, slownesses[imax_slow]*km, 'o') |
|
|
|
bazi_max = bazis[imax_bazi]*d2r |
|
slow_max = slownesses[imax_slow]*km |
|
ax.plot(bazi_max, slow_max, 'b.') |
|
ax.text(0.5, 0.01, 'Maximum at %s degrees, %s s/km' % |
|
(num.round(bazi_max, 1), slow_max), |
|
transform=fig4.transFigure, |
|
horizontalalignment='center', |
|
verticalalignment='bottom') |
|
|
|
if azi_theo: |
|
ax.arrow(azi_theo/180.*num.pi, num.min(slownesses), 0, |
|
num.max(slownesses), alpha=0.5, width=0.015, |
|
edgecolor='black', facecolor='green', lw=2, |
|
zorder=5) |
|
|
|
self.adjust_polar_axis(ax) |
|
fig4.colorbar(m) |
|
|
|
# --------------------------------------------------------- |
|
# CF and beam forming |
|
# --------------------------------------------------------- |
|
fig5 = self.new_figure('Beam') |
|
nsubplots = 4 |
|
nsubplots += self.want_pws |
|
|
|
ax_raw = fig5.add_subplot(nsubplots, 1, 1) |
|
ax_shifted = fig5.add_subplot(nsubplots, 1, 2) |
|
ax_beam = fig5.add_subplot(nsubplots, 1, 3) |
|
ax_beam_new = fig5.add_subplot(nsubplots, 1, 4) |
|
|
|
axkwargs = dict(alpha=0.3, linewidth=0.3, color='grey') |
|
|
|
ybeam = num.zeros(lengthout) |
|
ybeam_weighted = num.zeros(lengthout) |
|
for i, (shift, array) in enumerate(zip(shifts.T, arrays)): |
|
ax_raw.plot(times, array[npad: -npad], **axkwargs) |
|
ishift = shift[imax_bazi_slow] |
|
ax_shifted.plot( |
|
times, array[npad-ishift: -npad-ishift], **axkwargs) |
|
|
|
ydata = traces[i].get_ydata()[npad-ishift: -npad-ishift] |
|
ybeam += ydata |
|
|
|
# calculate phase weighting |
|
if self.want_pws: |
|
ph_inst = instantaneous_phase(ydata) |
|
ybeam_weighted += num.abs(num.exp(ph_inst))**4 |
|
|
|
ax_beam_new.plot(stack_trace) |
|
ax_beam_new.set_title('continuous mode') |
|
# ax_beam.plot(times, ybeam, color='black') |
|
ax_beam.plot(ybeam, color='black') |
|
ax_raw.set_title('Characteristic Function') |
|
ax_shifted.set_title('Shifted CF') |
|
ax_beam.set_title('Linear Stack') |
|
|
|
if self.want_pws: |
|
ax_playground = fig5.add_subplot(nsubplots, 1, 5) |
|
ax_playground.plot(ybeam*ybeam_weighted/len(arrays)) |
|
ax_playground.set_title('Phase Weighted Stack') |
|
|
|
# ----------------------------------------------------------- |
|
# polar movie: |
|
# ----------------------------------------------------------- |
|
fig6 = self.new_figure('Beam') |
|
self.polar_movie( |
|
fig=fig6, |
|
frames=frames, |
|
times=times, |
|
theta=theta.T, |
|
r=r.T*km, |
|
nth_frame=2, |
|
n_bazis=n_bazis, |
|
n_slow=n_slow, |
|
) |
|
|
|
self.draw_figures() |
|
|
|
self.irun += 1 |
|
|
|
def polar_movie(self, fig, frames, times, theta, r, nth_frame, |
|
n_bazis, n_slow): |
|
frame_artists = [] |
|
# progress_artists = [] |
|
ax = fig.add_subplot(111, projection='polar') |
|
iframe_min = 0 |
|
iframe_max = len(times)-1 |
|
vmin = num.min(frames) |
|
vmax = num.max(frames) |
|
self.adjust_polar_axis(ax) |
|
|
|
def update(iframe): |
|
# if iframe is not None: |
|
if False: |
|
frame = frames[:, iframe] |
|
''' |
|
if not progress_artists: |
|
progress_artists[:] = [axes2.axvline( |
|
tmin_frames - t0 + deltat_cf * iframe, |
|
color=scolor('scarletred3'), |
|
alpha=0.5, |
|
lw=2.)] |
|
|
|
else: |
|
progress_artists[0].set_xdata( |
|
tmin_frames - t0 + deltat_cf * iframe) |
|
''' |
|
else: |
|
frame = frames[:, iframe].reshape((n_bazis, n_slow)) |
|
|
|
frame_artists[:] = [ax.pcolormesh(theta, r, frame, vmin=vmin, |
|
vmax=vmax)] |
|
|
|
# return frame_artists + progress_artists + static_artists |
|
return frame_artists |
|
|
|
axf = FuncAnimation( # noqa |
|
fig, update, |
|
frames=list( |
|
range(iframe_min, iframe_max+1))[::nth_frame] + [None], |
|
interval=20., |
|
repeat=False, |
|
blit=True) |
|
|
|
fig.canvas.draw() |
|
|
|
def adjust_polar_axis(self, ax): |
|
ax.set_theta_zero_location('N') |
|
ax.set_theta_direction(-1) |
|
ax.set_xticks([0, num.pi/2., num.pi, 3*num.pi/2]) |
|
ax.set_xticklabels(['N', 'E', 'S', 'W']) |
|
|
|
def get_tinc_use(self, precision=1.): |
|
''' |
|
Set increment time for data processing. |
|
''' |
|
if self.tinc is not None: |
|
tinc = self.tinc |
|
else: |
|
tmin, tmax = self.get_selected_time_range(fallback=True) |
|
if tmin == tmax: |
|
# selected event marker |
|
tmin, tmax = self.get_viewer().get_time_range() |
|
tinc = tmax-tmin |
|
return num.floor(tinc/precision) * precision |
|
|
|
|
|
def __snufflings__(): |
|
return [FK()]
|
|
|