该博客文章详细介绍了如何使用pycwt库进行小波分析，包括数据预处理、小波变换、逆变换、功率谱计算以及显著性测试。文中提供了官方示例代码，并展示了如何修改代码以绘制小波图和方差图，适用于时间序列分析和信号处理。
安装模块
pycwt这个库在官方文档中没有提供conda的安装方法，官网推荐使用pip进行安装pip install pycwt。
但是，conda也可以用如下方法安装，推荐conda安装
conda install -c conda-forge/label/gcc7 pycwt
官网例子
https://github.com/regeirk/pycwt
"""
In this example we will load the NINO3 sea surface temperature anomaly dataset
between 1871 and 1996. This and other sample data files are kindly provided by
C. Torrence and G. Compo at
<http://paos.colorado.edu/research/wavelets/software.html>.
"""
# We begin by importing the relevant libraries. Please make sure that PyCWT is
# properly installed in your system.
from __future__ import division
import numpy
from matplotlib import pyplot
import pycwt as wavelet
from pycwt.helpers import find
# Then, we load the dataset and define some data related parameters. In this
# case, the first 19 lines of the data file contain meta-data, that we ignore,
# since we set them manually (*i.e.* title, units).
url = 'http://paos.colorado.edu/research/wavelets/wave_idl/nino3sst.txt'
dat = numpy.genfromtxt(url, skip_header=19)
title = 'NINO3 Sea Surface Temperature'
label = 'NINO3 SST'
units = 'degC'
t0 = 1871.0
dt = 0.25 # In years
# We also create a time array in years.
N = dat.size
t = numpy.arange(0, N) * dt + t0
# We write the following code to detrend and normalize the input data by its
# standard deviation. Sometimes detrending is not necessary and simply
# removing the mean value is good enough. However, if your dataset has a well
# defined trend, such as the Mauna Loa CO\ :sub:`2` dataset available in the
# above mentioned website, it is strongly advised to perform detrending.
# Here, we fit a one-degree polynomial function and then subtract it from the
# original data.
p = numpy.polyfit(t - t0, dat, 1)
dat_notrend = dat - numpy.polyval(p, t - t0)
std = dat_notrend.std() # Standard deviation
var = std ** 2 # Variance
dat_norm = dat_notrend / std # Normalized dataset
# The next step is to define some parameters of our wavelet analysis. We
# select the mother wavelet, in this case the Morlet wavelet with
# :math:`\omega_0=6`.
mother = wavelet.Morlet(6)
s0 = 2 * dt # Starting scale, in this case 2 * 0.25 years = 6 months
dj = 1 / 12 # Twelve sub-octaves per octaves
J = 7 / dj # Seven powers of two with dj sub-octaves
alpha, _, _ = wavelet.ar1(dat) # Lag-1 autocorrelation for red noise
# The following routines perform the wavelet transform and inverse wavelet
# transform using the parameters defined above. Since we have normalized our
# input time-series, we multiply the inverse transform by the standard
# deviation.
wave, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(dat_norm, dt, dj, s0, J,
mother)
iwave = wavelet.icwt(wave, scales, dt, dj, mother) * std
# We calculate the normalized wavelet and Fourier power spectra, as well as
# the Fourier equivalent periods for each wavelet scale.
power = (numpy.abs(wave)) ** 2
fft_power = numpy.abs(fft) ** 2
period = 1 / freqs
# We could stop at this point and plot our results. However we are also
# interested in the power spectra significance test. The power is significant
# where the ratio ``power / sig95 > 1``.
signif, fft_theor = wavelet.significance(1.0, dt, scales, 0, alpha,
significance_level=0.95,
wavelet=mother)
sig95 = numpy.ones([1, N]) * signif[:, None]
sig95 = power / sig95
# Then, we calculate the global wavelet spectrum and determine its
# significance level.
glbl_power = power.mean(axis=1)
dof = N - scales # Correction for padding at edges
glbl_signif, tmp = wavelet.significance(var, dt, scales, 1, alpha,
significance_level=0.95, dof=dof,
wavelet=mother)
# We also calculate the scale average between 2 years and 8 years, and its
# significance level.
sel = find((period >= 2) & (period < 8))
Cdelta = mother.cdelta
scale_avg = (scales * numpy.ones((N, 1))).transpose()
scale_avg = power / scale_avg # As in Torrence and Compo (1998) equation 24
scale_avg = var * dj * dt / Cdelta * scale_avg[sel, :].sum(axis=0)
scale_avg_signif, tmp = wavelet.significance(var, dt, scales, 2, alpha,
significance_level=0.95,
dof=[scales[sel[0]],
scales[sel[-1]]],
wavelet=mother)
# Finally, we plot our results in four different subplots containing the
# (i) original series anomaly and the inverse wavelet transform; (ii) the
# wavelet power spectrum (iii) the global wavelet and Fourier spectra ; and
# (iv) the range averaged wavelet spectrum. In all sub-plots the significance
# levels are either included as dotted lines or as filled contour lines.
# Prepare the figure
pyplot.close('all')
pyplot.ioff()
figprops = dict(figsize=(11, 8), dpi=72)
fig = pyplot.figure(**figprops)
# First sub-plot, the original time series anomaly and inverse wavelet
# transform.
ax = pyplot.axes([0.1, 0.75, 0.65, 0.2])
ax.plot(t, iwave, '-', linewidth=1, color=[0.5, 0.5, 0.5])
ax.plot(t, dat, 'k', linewidth=1.5)
ax.set_title('a) {}'.format(title))
ax.set_ylabel(r'{} [{}]'.format(label, units))
# Second sub-plot, the normalized wavelet power spectrum and significance
# level contour lines and cone of influece hatched area. Note that period
# scale is logarithmic.
bx = pyplot.axes([0.1, 0.37, 0.65, 0.28], sharex=ax)
levels = [0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16]
bx.contourf(t, numpy.log2(period), numpy.log2(power), numpy.log2(levels),
extend='both', cmap=pyplot.cm.viridis)
extent = [t.min(), t.max(), 0, max(period)]
bx.contour(t, numpy.log2(period), sig95, [-99, 1], colors='k', linewidths=2,
extent=extent)
bx.fill(numpy.concatenate([t, t[-1:] + dt, t[-1:] + dt,
t[:1] - dt, t[:1] - dt]),
numpy.concatenate([numpy.log2(coi), [1e-9], numpy.log2(period[-1:]),
numpy.log2(period[-1:]), [1e-9]]),
'k', alpha=0.3, hatch='x')
bx.set_title('b) {} Wavelet Power Spectrum ({})'.format(label, mother.name))
bx.set_ylabel('Period (years)')
#
Yticks = 2 ** numpy.arange(numpy.ceil(numpy.log2(period.min())),
numpy.ceil(numpy.log2(period.max())))
bx.set_yticks(numpy.log2(Yticks))
bx.set_yticklabels(Yticks)
# Third sub-plot, the global wavelet and Fourier power spectra and theoretical
# noise spectra. Note that period scale is logarithmic.
cx = pyplot.axes([0.77, 0.37, 0.2, 0.28], sharey=bx)
cx.plot(glbl_signif, numpy.log2(period), 'k--')
cx.plot(var * fft_theor, numpy.log2(period), '--', color='#cccccc')
cx.plot(var * fft_power, numpy.log2(1./fftfreqs), '-', color='#cccccc',
linewidth=1.)
cx.plot(var * glbl_power, numpy.log2(period), 'k-', linewidth=1.5)
cx.set_title('c) Global Wavelet Spectrum')
cx.set_xlabel(r'Power [({})^2]'.format(units))
cx.set_xlim([0, glbl_power.max() + var])
cx.set_ylim(numpy.log2([period.min(), period.max()]))
cx.set_yticks(numpy.log2(Yticks))
cx.set_yticklabels(Yticks)
pyplot.setp(cx.get_yticklabels(), visible=False)
# Fourth sub-plot, the scale averaged wavelet spectrum.
dx = pyplot.axes([0.1, 0.07, 0.65, 0.2], sharex=ax)
dx.axhline(scale_avg_signif, color='k', linestyle='--', linewidth=1.)
dx.plot(t, scale_avg, 'k-', linewidth=1.5)
dx.set_title('d) {}--{} year scale-averaged power'.format(2, 8))
dx.set_xlabel('Time (year)')
dx.set_ylabel(r'Average variance [{}]'.format(units))
ax.set_xlim([t.min(), t.max()])
pyplot.show()
修改后的小波程序
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import pycwt as wavelet
from pycwt.helpers import find
from matplotlib import pyplot
def drwa_cwt(t, period, power, glbl_signif, glbl_power):
# 小波图
ax = plt.axes([0.1, 0.37, 0.65, 0.28])
levels = [0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8]
cs = ax.contourf(t, np.log2(period), np.log2(power), np.log2(levels),
extend='both', cmap=plt.cm.viridis)
plt.colorbar(cs)
extent = [t.min(), t.max(), 0, max(period)]
ax.contour(t, np.log2(period), sig95, [-99, 1], colors='k', linewidths=1,
extent=extent)
ax.fill(np.concatenate([t, t[-1:] + dt, t[-1:] + dt,
t[:1] - dt, t[:1] - dt]),
np.concatenate([np.log2(coi), [1e-9], np.log2(period[-1:]),
np.log2(period[-1:]), [1e-9]]),
'k', alpha=0.3, hatch='x')
ax.set_ylabel('Period (years)')
ax.set_xlabel('Time')
Yticks = 2 ** np.arange(np.ceil(np.log2(period.min())),
np.ceil(np.log2(period.max())))
ax.set_yticks(np.log2(Yticks))
ax.set_yticklabels(Yticks)
ax.set_xlim([t.min(), t.max()])
ax.set_xticks(range(1961, 2021, 10))
# 方差图
bx = pyplot.axes([0.77, 0.37, 0.2, 0.28], sharey=ax)
bx.plot(glbl_signif, np.log2(period), 'k--')
bx.plot(var * glbl_power, np.log2(period), 'k-', linewidth=1.5)
bx.set_xlabel('Power (mm^2)')
bx.set_ylim(np.log2([period.min(), period.max()]))
bx.set_yticks(np.log2(Yticks))
bx.set_yticklabels(Yticks)
pyplot.setp(bx.get_yticklabels(), visible=False)
pyplot.show()
if __name__ == '__main__':
x = np.linspace(1961, 2020, 60)
dat = np.cos(np.pi / 2 * x)
plt.plot(x, dat)
plt.show()
dat = np.array(dat)
t0 = 1961.0 # 开始的时间，以年为单位
dt = 1 # 采样间隔，以年为单位
N = dat.size # 时间序列的长度
t = np.arange(0, N) * dt + t0 # 构造时间序列数组
p = np.polyfit(t - t0, dat, 1) # 线性拟合
dat_notrend = dat - np.polyval(p, t - t0) # 去趋势
std = dat_notrend.std() # 标准差
var = std ** 2 # 方差
dat_norm = dat_notrend / std # 标准化
mother = wavelet.Morlet(6) # Monther Wavelet: Morlet
s0 = -1 # Starting scale, in this case 2 * 0.25 years = 6 months
dj = 1 / 12 # Twelve sub-octaves per octaves
J = -1 # Seven powers of two with dj sub-octaves
alpha, _, _ = wavelet.ar1(dat) # Lag-1 autocorrelation for red noise
wave, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(dat_norm, dt, dj=1 / 12, s0=-1,
J=-1, wavelet='morlet', freqs=None)
iwave = wavelet.icwt(wave, scales, dt, dj, mother) * std
power = (np.abs(wave)) ** 2
fft_power = np.abs(fft) ** 2
period = 1 / freqs
# We could stop at this point and plot our results. However we are also
# interested in the power spectra significance test. The power is significant
# where the ratio ``power / sig95 > 1``.
signif, fft_theor = wavelet.significance(1.0, dt, scales, 0, alpha,
significance_level=0.95,
wavelet=mother)
sig95 = np.ones([1, N]) * signif[:, None]
sig95 = power / sig95
glbl_power = power.mean(axis=1)
dof = N - scales # Correction for padding at edges
glbl_signif, tmp = wavelet.significance(var, dt, scales, 1, alpha,
significance_level=0.95, dof=dof,
wavelet=mother)
sel = find((period >= 2) & (period < 8))
Cdelta = mother.cdelta
scale_avg = (scales * np.ones((N, 1))).transpose()
scale_avg = power / scale_avg # As in Torrence and Compo (1998) equation 24
scale_avg = var * dj * dt / Cdelta * scale_avg[sel, :].sum(axis=0)
scale_avg_signif, tmp = wavelet.significance(var, dt, scales, 2, alpha,
significance_level=0.95,
dof=[scales[sel[0]],
scales[sel[-1]]],
wavelet=mother)
drwa_cwt(t, period, power, glbl_signif, glbl_power)
按y=sin(x)画的结果：
进一步修改
# -*- coding: utf-8 -*-
"""
@Features: 画小波和方差图
@Author:
@Date：2022-05-20
"""
import glob
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import pycwt as wavelet
from pycwt.helpers import find
from matplotlib import pyplot
def cwt_calculate(dat):
t0 = 1961.0 # 开始的时间，以年为单位
dt = 1 # 采样间隔，以年为单位
N = dat.size # 时间序列的长度
t = np.arange(0, N) * dt + t0 # 构造时间序列数组
# We write the following code to detrend and normalize the input data by its
# standard deviation. Sometimes detrending is not necessary and simply
# removing the mean value is good enough. However, if your dataset has a well
# defined trend, such as the Mauna Loa CO\ :sub:`2` dataset available in the
# above mentioned website, it is strongly advised to perform detrending.
# Here, we fit a one-degree polynomial function and then subtract it from the
# original data.
p = np.polyfit(t - t0, dat, 1) # 线性拟合
dat_notrend = dat - np.polyval(p, t - t0) # 去趋势
std = dat_notrend.std() # 标准差
var = std ** 2 # 方差
dat_norm = dat_notrend / std # 标准化
# plt.figure()
# plt.plot(range(1961, 2021), dat_norm)
# The next step is to define some parameters of our wavelet analysis. We
# select the mother wavelet, in this case the Morlet wavelet with
# :math:`\omega_0=6`.
mother = wavelet.Morlet(6) # Monther Wavelet: Morlet
s0 = -1
dj = 1 / 12 # Twelve sub-octaves per octaves
J = -1
alpha, _, _ = wavelet.ar1(dat) # Lag-1 autocorrelation for red noise
# The following routines perform the wavelet transform and inverse wavelet
# transform using the parameters defined above. Since we have normalized our
# input time-series, we multiply the inverse transform by the standard
# deviation.
wave, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(dat_norm, dt, dj=1 / 12, s0=-1,
J=-1, wavelet='morlet', freqs=None)
iwave = wavelet.icwt(wave, scales, dt, dj, mother) * std
# We calculate the normalized wavelet and Fourier power spectra, as well as
# the Fourier equivalent periods for each wavelet scale.
power = (np.abs(wave)) ** 2
fft_power = np.abs(fft) ** 2
period = 1 / freqs
# We could stop at this point and plot our results. However we are also
# interested in the power spectra significance test. The power is significant
# where the ratio ``power / sig95 > 1``.
signif, fft_theor = wavelet.significance(1.0, dt, scales, 0, alpha,
significance_level=0.95,
wavelet=mother)
sig95 = np.ones([1, N]) * signif[:, None]
sig95 = power / sig95
# Then, we calculate the global wavelet spectrum and determine its
# significance level.
glbl_power = power.mean(axis=1)
dof = N - scales # Correction for padding at edges
glbl_signif, tmp = wavelet.significance(var, dt, scales, 1, alpha,
significance_level=0.95, dof=dof,
wavelet=mother)
# We also calculate the scale average between 2 years and 8 years, and its
# significance level.
sel = find((period >= 2) & (period < 8))
Cdelta = mother.cdelta
scale_avg = (scales * np.ones((N, 1))).transpose()
scale_avg = power / scale_avg # As in Torrence and Compo (1998) equation 24
scale_avg = var * dj * dt / Cdelta * scale_avg[sel, :].sum(axis=0)
scale_avg_signif, tmp = wavelet.significance(var, dt, scales, 2, alpha,
significance_level=0.95,
dof=[scales[sel[0]],
scales[sel[-1]]],
wavelet=mother)
return t, period, power, glbl_signif, glbl_power, sig95, coi, var
def drwa_cwt(t, period, power, glbl_signif, glbl_power, sig95, coi, var, text_i, x, h):
dt = 1
ax = plt.axes([x, h, 0.2, 0.15])
levels = [0.125, 0.25, 0.5, 1, 2, 4, 8]
cs = ax.contourf(t, np.log2(period), np.log2(power), np.log2(levels),
extend='both', cmap=plt.cm.viridis)
plt.colorbar(cs)
extent = [t.min(), t.max(), 0, max(period)]
ax.contour(t, np.log2(period), sig95, [-99, 1], colors='k', linewidths=1, extent=extent)
ax.fill(np.concatenate([t, t[-1:] + dt, t[-1:] + dt, t[:1] - dt, t[:1] - dt]),
np.concatenate([np.log2(coi), [1e-9], np.log2(period[-1:]), np.log2(period[-1:]), [1e-9]]),
'k', alpha=0.4, hatch='x')
ax.set_ylabel('Period (years)')
ax.set_xlabel('Time')
Yticks = 2 ** np.arange(np.ceil(np.log2(period.min())),
np.ceil(np.log2(period.max())))
ax.set_yticks(np.log2(Yticks))
ax.set_yticklabels(Yticks)
ax.set_xlim([t.min(), t.max()])
ax.set_xticks(range(1961, 2021, 10))
plt.text(0.01, 0.9, text_i, transform=ax.transAxes, c='w')
# 方差图
bx = pyplot.axes([x + 0.21, h, 0.1, 0.15], sharey=ax)
bx.plot(glbl_signif, np.log2(period), 'k--')
bx.plot(var * glbl_power, np.log2(period), 'k-', linewidth=1.5)
bx.set_xlabel('Power (mm^2)')
bx.set_ylim(np.log2([period.min(), period.max()]))
bx.set_yticks(np.log2(Yticks))
bx.set_yticklabels(Yticks)
pyplot.setp(bx.get_yticklabels(), visible=False)
def main_code(fn):
'''
:param fn: filename
:return:
'''
# files1 = glob.glob(r'..\data\各区降水指数\逐站\*年降水量.csv')
dataCn1 = np.zeros([60, 1])
text = ['(a)', '(b)', '(c)', '(d)', '(e)', '(f)', '(g)', '(h)']
i = 0
for file in fn:
dataRead1 = pd.read_csv(file)
data = dataRead1.iloc[3:, 2:]
dataArr = np.array(data)
dataCn1 = np.hstack([dataCn1, dataArr])
dat = np.array(data)
dat = np.median(dat, 1)
t, period, power, glbl_signif, glbl_power, sig95, coi, var = cwt_calculate(dat)
# plt.figure()
position = np.array([[0.15, 0.8], [0.6, 0.8],
[0.15, 0.58], [0.6, 0.58],
[0.15, 0.36], [0.6, 0.36],
[0.15, 0.14], [0.6, 0.14]])
p_x, p_y = position[i, 0], position[i, 1]
drwa_cwt(t, period, power, glbl_signif, glbl_power, sig95, coi, var, text[i], p_x, p_y)
i += 1
dataCn1 = np.delete(dataCn1, 0, 1) # delete column 0,全国数据
dataCn1_median = np.median(dataCn1, 1)
t, period, power, glbl_signif, glbl_power, sig95, coi, var = cwt_calculate(dataCn1_median)
p_x, p_y = position[i, 0], position[i, 1]
drwa_cwt(t, period, power, glbl_signif, glbl_power, sig95, coi, var, text[i], p_x, p_y)
if __name__ == '__main__':
plt.figure(1)
text = ['(a)', '(b)', '(c)', '(d)', '(e)', '(f)', '(g)', '(h)']
files1 = glob.glob(r'..\data\各区降水指数\逐站\*年降水量.csv')
main_code(files1)
plt.show()
plt.figure(2)
files2 = glob.glob(r'..\data\各区降水指数\逐站\*5日降水量.csv')
main_code(files2)
plt.show()
plt.figure(3)
files3 = glob.glob(r'..\data\各区降水指数\逐站\*95*日数.csv')
main_code(files3)
plt.show()
结果：
修改y轴区间
# -*- coding: utf-8 -*-
"""
@Features: 计算并画小波和方差图;按审稿意见修改
@Author: L.F. Zhan
@Date：2022-05-06;2024-06-14
"""
import glob
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import pycwt as wavelet
from pycwt.helpers import find
plt.rcParams['font.family'] = 'Times New Roman, SimSun' # 设置字体族，中文为SimSun，英文为Times New Roman
plt.rcParams['mathtext.fontset'] = 'stix' # 设置数学公式字体为stix
def cwt_cal(dat):
'''
Parameters
----------
dat: 一维数据
Returns
-------
'''
t0 = 1961
dt = 1 # In years
# We also create sin_value time array in years.
N = dat.size
t = np.arange(0, N) * dt + t0
# We write the following code to detrend and normalize the input data by its
# standard deviation. Sometimes detrending is not necessary and simply
# removing the mean value is good enough. However, if your dataset has sin_value well
# defined trend, such as the Mauna Loa CO\ :sub:`2` dataset available in the
# above mentioned website, it is strongly advised to perform detrending.
# Here, we fit sin_value one-degree polynomial function and then subtract it from the
# original data.
p = np.polyfit(t - t0, dat, 1)
dat_notrend = dat - np.polyval(p, t - t0)
std = dat_notrend.std() # Standard deviation
var = std ** 2 # Variance
dat_norm = dat_notrend / std # Normalized dataset
# The next step is to define some parameters of our wavelet analysis. We
# select the mother wavelet, in this case the Morlet wavelet with
# :math:`\omega_0=6`.
mother = wavelet.Morlet(6)
s0 = 2 * dt # Starting scale, in this case 2 * 0.25 years = 6 months
dj = 1 / 12 # Twelve sub-octaves per octaves
J = 7 / dj # Seven powers of two with dj sub-octaves
alpha, _, _ = wavelet.ar1(dat) # Lag-1 autocorrelation for red noise
# The following routines perform the wavelet transform and inverse wavelet
# transform using the parameters defined above. Since we have normalized our
# input time-series, we multiply the inverse transform by the standard
# deviation.
wave, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(dat_norm, dt, dj, s0, J, mother)
iwave = wavelet.icwt(wave, scales, dt, dj, mother) * std
# We calculate the normalized wavelet and Fourier power spectra, as well as
# the Fourier equivalent periods for each wavelet scale.
power = (np.abs(wave)) ** 2
fft_power = np.abs(fft) ** 2
period = 1 / freqs
# We could stop at this point and plot our results. However we are also
# interested in the power spectra significance test. The power is significant
# where the ratio ``power / sig95 > 1``.
signif, fft_theor = wavelet.significance(1.0, dt, scales, 0, alpha,
significance_level=0.95,
wavelet=mother)
sig95 = np.ones([1, N]) * signif[:, None]
sig95 = power / sig95
# Then, we calculate the global wavelet spectrum and determine its
# significance level.
glbl_power = power.mean(axis=1)
dof = N - scales # Correction for padding at edges
glbl_signif, tmp = wavelet.significance(var, dt, scales, 1, alpha,
significance_level=0.95, dof=dof,
wavelet=mother)
# We also calculate the scale average between 2 years and 8 years, and its
# significance level.
sel = find((period >= 2) & (period < 8))
Cdelta = mother.cdelta
scale_avg = (scales * np.ones((N, 1))).transpose()
scale_avg = power / scale_avg # As in Torrence and Compo (1998) equation 24
scale_avg = var * dj * dt / Cdelta * scale_avg[sel, :].sum(axis=0)
scale_avg_signif, tmp = wavelet.significance(var, dt, scales, 2, alpha,
significance_level=0.95,
dof=[scales[sel[0]],
scales[sel[-1]]],
wavelet=mother)
y_32 = 45 #设置Y轴32年的位置索引，控制Y轴显示
return t, period[:y_32], power[:y_32, :], glbl_signif[:y_32], glbl_power[:y_32], sig95[:y_32,
:], coi, var
def drwa_cwt(t, period, power, glbl_signif, glbl_power,
sig95, coi, var, text_i, x, h):
dt = 1
ax = plt.axes([x, h, 0.56, 0.15])
levels = [0.125, 0.25, 0.5, 1, 2, 4, 8] #图例等级划分
# levels = [1, 2, 4, 8, 16] # 新的等高线级别
cs = ax.contourf(t, np.log2(period), np.log2(power), np.log2(levels), extend='both', cmap=plt.cm.RdBu_r)
plt.colorbar(cs)
extent = [t.min(), t.max(), 0, max(period)]
ax.contour(t, np.log2(period), sig95, [-99, 1], colors='k', linewidths=1, extent=extent)
ax.fill(np.concatenate([t, t[-1:] + dt, t[-1:] + dt, t[:1] - dt, t[:1] - dt]),
np.concatenate([np.log2(coi), [1e-9], np.log2(period[-1:]), np.log2(period[-1:]), [1e-9]]), 'k',
alpha=0.6, hatch='x')
ax.set_ylabel('周期 (年)')
ax.set_xlabel('年份')
Yticks = 2 ** np.arange(np.ceil(np.log2(period.min())),
np.ceil(np.log2(period.max())))
ax.set_yticks(np.log2(Yticks))
ax.set_yticklabels(Yticks)
ax.set_xlim([t.min(), t.max()])
ax.set_xticks(range(1961, 2021, 10))
plt.text(0, 1.1, text_i, transform=ax.transAxes, c='k', size=10)
# 方差图
bx = plt.axes([x + 0.58, h, 0.15, 0.15], sharey=ax)
bx.plot(glbl_signif, np.log2(period), 'k--')
bx.plot(var * glbl_power, np.log2(period), 'k-', linewidth=1.5)
bx.set_xlabel('能量 (mm^2)')
bx.set_ylim(np.log2([period.min(), period.max()]))
bx.set_yticks(np.log2(Yticks))
bx.set_yticklabels(Yticks)
plt.setp(bx.get_yticklabels(), visible=False)
def main_code(fn):
'''
:param fn: filename
:return:
'''
# files1 = glob.glob(r'..\data\各区降水指数\逐站\*年降水量.csv')
dataCn1 = np.zeros([60, 1])
text = ['(a)干旱区', '(b)半干旱区', '(c)半湿润区', '(d)湿润区']
i = 0
for file in fn:
print(file)
dataRead1 = pd.read_csv(file)
data = dataRead1['1']
dataArr = np.array(data)
t, period, power, glbl_signif, glbl_power, sig95, coi, var = cwt_cal(dataArr)
# plt.figure()图像位置
x_position, y_position, move_position = 0.18, 0.8, 0.23
position = np.array([[x_position, y_position],
[x_position, y_position - move_position],
[x_position, y_position - move_position * 2],
[x_position, y_position - move_position * 3]])
p_x, p_y = position[i, 0], position[i, 1]
drwa_cwt(t, period, power, glbl_signif, glbl_power, sig95, coi, var, text[i], p_x, p_y)
i += 1
if __name__ == '__main__':
plt.figure(1)
plt.figure(figsize=(5, 8))
files1 = glob.glob(r'..\data\4区降水指数\逐站\*_全域_年降水量.csv')
main_code(files1)
plt.savefig(r'F:\cwt_prcptot.svg', dpi=700)
plt.show()
plt.figure(2)
plt.figure(figsize=(5, 8))
files2 = glob.glob(r'..\data\4区降水指数\逐站\*_全域_max5day.csv')
main_code(files2)
plt.savefig(r'F:\cwt_rx5day.svg', dpi=700)
plt.show()
plt.figure(3)
plt.figure(figsize=(5, 8))
files3 = glob.glob(r'..\data\4区降水指数\逐站\*_全域_*95*日数.csv')
main_code(files3)
plt.savefig(r'F:\cwt_r95p.svg', dpi=700)
plt.show()
扫一扫
4302
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
