[SciPy-Dev] 2D histogram: request for plotting variable bin size
Sturla Molden
sturla@molden...
Thu Jan 31 09:40:03 CST 2013
Hi Frank
I use Python to plot wavelet spectrograms with logarithmically scaled
frequencies. I e.g. interpolate the spectrum to the frequencies I like
to plot, and pass that to pcolor or contourf. Usually I plot with a
linear y-axis first, and then use set_yscale('log'). Usually the y-axis
must be flipped too.
Then I e.g. end up with something like this, plotting with 250
logartihmically scaled frequencies from 10 Hz to 10 kHz:
cwt,scale= wavelet(signal, sampling_rate)
time = np.arange(signal.shape[0])/sampling_rate
power = np.abs(cwt)**2/scale # bias corrected power
periods = wavelet_to_fourier(scale)
spectrum_freqs = np.logspace(np.log10(10), np.log10(10000), 250)
spectrum_periods = spectrum_freqs**-1
spectrum = interpolate_spectrum(spectrum_periods, power, periods)
X,Y = np.meshgrid(time, spectrum_freqs)
cs = plt.contourf(X,Y,spectrum, cmap=cm.jet, origin=origin)
ax = plt.gca()
ax.set_yscale('log')
ylim = ax.get_ylim()
ax.set_ylim(ylim[::-1])
import numpy as np
from bisect import bisect_left
def interpolate_spectrum(spectrum_periods, wavelet_power, wavelet_periods):
""" linear interpolation for plotting spectrum at given frequencies
"""
n0 = spectrum_periods.shape[0]
n1 = wavelet_power.shape[1]
m = wavelet_power.shape[0]
spectrum_power = np.zeros((n0,n1))
for k,sp in enumerate(spectrum_periods):
i = bisect_left(wavelet_periods,sp)
j = i-1
i = 0 if i<0 else i
i = m-1 if i > m-1 else i
j = 0 if j<0 else j
j = m-1 if j > m-1 else j
eps = 1E-5 # avoid log of zero
di = 1./(np.abs(sp - wavelet_periods[i]) + eps)
dj = 1./(np.abs(sp - wavelet_periods[j]) + eps)
spectrum_power[k,:] = (di*wavelet_power[i,:] +
dj*wavelet_power[j,:])/(di+dj)
return spectrum_power
Not sure if this helps or not, but it might.
Sturla
On 31.01.2013 16:00, Frank Breitling wrote:
> Hi Skipper,
>
> I would like to use Python for plotting dynamic radio spectra e.g. as
> found at
>
> http://science.nasa.gov/media/medialibrary/1998/12/20/ast22dec98_1_resources/leonid_dynspectrum.jpg
>
> (Image is also attached).
> Given the high resolution in time and frequency I believe broken_barh is
> not suited.
>
> But thanks for your suggestion
>
> Frank
>
>
> On 31.01.2013 15:17, Skipper Seabold wrote:
>> On Thu, Jan 31, 2013 at 9:01 AM, Frank Breitling <fbreitling@aip.de
>> <mailto:fbreitling@aip.de>> wrote:
>>
>> Hello,
>>
>> http://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram2d.html
>>
>> explains the usage of 2D histograms.
>> Although it is possible to construct 2D histograms with a variable bin
>> size, there is no simple way to plot them.
>> Therefore I would like to request an implementation to imshow that
>> allows the visualization of 2D histograms with variable bin size. I
>> imagine a syntax like
>>
>> imshow(historgram2d, xedges=xedges, yedges=yedges)
>>
>> where xedges and yedges are the bin edges along the x- and y-axis.
>>
>> I have attached a short program and its output which constructs a
>> numpy.historgram2d with variable bin width (xedges=[0,1,3]) and shows
>> its bin content with imshow.
>> I would highly appreciate if imshow would be extended to represent the
>> histogram correctly.
>> If there is an alternative solution, I would be interested to learn
>> about it as well.
>>
>>
>> This may be a discussion for matplotlib-user, but have you had a look
>> at trying to use broken_barh for your needs?
>>
>> http://matplotlib.org/examples/pylab_examples/broken_barh.html
>>
>> Skipper
>>
>>
>> _______________________________________________
>> SciPy-Dev mailing list
>> SciPy-Dev@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
>
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
More information about the SciPy-Dev
mailing list