请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 2788|回复: 3

脚本:离散功率谱估计,双X轴如何实现,请老师帮忙看下。

[复制链接]

新浪微博达人勋

发表于 2018-4-22 19:53:06 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
黄嘉佑老师编写教课书,离散功率谱估计,脚本如下:
#打开1月上旬温度数据文件
table = readtable('f:/temp.csv', delimiter=',', format='%s%f')
xt = table['temp']
#print xt
n = len(xt)
m = n/2
ind = arange(m+1)
ak = zeros(m+1)
bk = zeros(m+1)
sksq = zeros(m+1)
ak[0] = sum(xt)/n
for k in arange(1,m+1):
    for t in arange(0,n):   
        ak[k] = ak[k] + xt[t]*cos((2*pi*k*t)/n)
        bk[k] = bk[k] + xt[t]*sin((2*pi*k*t)/n)
    ak[k] = ak[k]*2/n
    bk[k] = bk[k]*2/n
    sksq[k] = (ak[k]*ak[k]+bk[k]*bk[k])/2
#print sksq
width = 0.05
ax1=axes()
y = arange(0,2.0,0.5)
x1 = arange(1,10,2.0)
x2 = 20/x1
bar(ind,sksq,width,color='black',encoding='gbk')
yaxis(tickin=False,tickfontsize=17)
xaxis(tickin=False,tickfontsize=17)
yaxis(location='bottom',tickin=True,tickfontsize=18)
yticks(y)
xticks(x1)
xlabel('k',fontsize=16)
ylabel(r'$\rm(S^2)$')',fontsize=16)

#ax2 = twinx(ax1)
#xaxis(location='top',tickin=False,tickfontsize=18)
#xlabel('T(year)',fontsize=18,bold=False)
#antialias(True)
图形的顶部如何加刻度,呈现双X轴。



s.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-4-22 22:10:48 | 显示全部楼层
我把数据贴过来

temp.csv

236 Bytes, 下载次数: 2, 下载积分: 金钱 -5

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-4-23 17:45:41 | 显示全部楼层

增加了twiny函数,可以在MeteoInfo QQ群下载最新版本的MeteoInfo,参考此脚本:

  1. table = readtable('D:/Temp/ascii/temp.csv', delimiter=',', format='%s%f')
  2. xt = table['temp']
  3. #print xt
  4. n = len(xt)
  5. m = n/2
  6. ind = arange(m+1)
  7. ak = zeros(m+1)
  8. bk = zeros(m+1)
  9. sksq = zeros(m+1)
  10. ak[0] = sum(xt)/n
  11. for k in arange(1,m+1):
  12.     for t in arange(0,n):   
  13.         ak[k] = ak[k] + xt[t]*cos((2*pi*k*t)/n)
  14.         bk[k] = bk[k] + xt[t]*sin((2*pi*k*t)/n)
  15.     ak[k] = ak[k]*2/n
  16.     bk[k] = bk[k]*2/n
  17.     sksq[k] = (ak[k]*ak[k]+bk[k]*bk[k])/2
  18. #print sksq
  19. width = 0.05
  20. ax1=axes()
  21. y = arange(0,2.0,0.5)
  22. x1 = arange(1,12,2.0)
  23. x2 = 20/x1
  24. bar(ind,sksq,width,color='black')
  25. yaxis(tickin=False)
  26. xaxis(tickin=False)
  27. yticks(y)
  28. xticks(x1)
  29. xlabel('k')
  30. ylabel(r'$\rm(S^2)$', fontname='Arial')

  31. ax2 = twiny(ax1)
  32. xaxis(location='top',tickin=False)
  33. xlabel('T(year)',bold=False)
  34. title('xx test')
  35. #antialias(True)


密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-4-23 21:59:04 | 显示全部楼层
王老师,我看了下miplot的函数。改成脚本如下,发现两个问题1、只能改写底部,顶部不知怎么写;2、xlabel('k')语句报错,错误信息附后:
table = readtable('f:/temp.csv', delimiter=',', format='%s%f')
xt = table['temp']
#print xt
n = len(xt)
m = n/2
ind = arange(m+1)
ak = zeros(m+1)
bk = zeros(m+1)
sksq = zeros(m+1)
ak[0] = sum(xt)/n
for k in arange(1,m+1):
    for t in arange(0,n):   
        ak[k] = ak[k] + xt[t]*cos((2*pi*k*t)/n)
        bk[k] = bk[k] + xt[t]*sin((2*pi*k*t)/n)
    ak[k] = ak[k]*2/n
    bk[k] = bk[k]*2/n
    sksq[k] = (ak[k]*ak[k]+bk[k]*bk[k])/2
#print sksq
width = 0.05
ax1=axes()
y = zeros(5)
x1 = zeros(5)
x2 = zeros(5)
y = arange(0,2.0,0.5)
x1 = arange(1,11,2.0)
x2 = 20/x1
xlbl = []
j = 0
for i in x1:
    xlbl.append(str(format(x2[j],'0.1f')))
    j = j+1
print xlbl  
bar(ind,sksq,width,color='black')
yaxis(tickin=False)
xaxis(tickin=False)
yticks(y)
xticks(x1,xlbl)
#xlabel('k')
ylabel(r'$\rm(S^2)$', fontname='Arial')
ax2 = twiny(ax1)
#xaxis(ax2,shift = 10,color = 'r')
xaxis(ax2)
#print ax2
xaxis(location='top',tickin=False)

#xlabel(u'T(year)',bold=False)
#xlabel(u'T(年)',fontname=u'黑体',bold=False)
title('Discrete power spectrum estimation')
#antialias(True)

报错信息
xlabel('k')
TypeError: 'list' object is not callable


密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表