爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 24119|回复: 29

[其他] ncl 求特定时段标准差

[复制链接]
发表于 2014-2-26 14:52:21 | 显示全部楼层 |阅读模式

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

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

x
数据资料:1974.6.1~2013.8.31逐日olr资料,时间格式19740601
我现在想计算每一年5月~10月的olr标准差,求标准差的语句我知道了,用stddev,但是这个时间循环还是没弄明白,怎么才能每一年5~10月的提取出来分别算呢?time维数直接是14337~新手求教~~
密码修改失败请联系微信:mofangbao
发表于 2014-2-26 15:30:17 | 显示全部楼层
本帖最后由 longlivehj 于 2014-2-26 15:57 编辑

; 一维olr数组
olr1d = (/.../)
; 如果资料是从比如nc文件里面读出,可能已经有时间坐标变量,  下面两句可以不用
olr1d!0 = "time"
olr1d&time = (/19740601, ..., 20130831/)

; 构建二维olr数组,其中184是5-10月总天数,40是1974-2013总年数
olr2d = new((/184, 40/), typeof(olr))
olr2d!0 = "time"
olr2d&time = (/19740501, ..., 20131031/) % 10000
; 采用坐标方式索引olr1d
do y = 1974, 2013
    temp =  olr1d({y * 10000 + 501 : y * 10000 + 1031})
    olr2d({temp&time % 10000}, y - 1974) = temp
end do

; 求逐日标准差
std = dim_stddev(olr2d)

注意两点:
1. olr2d新建的时候是带缺测值(_FillValue)的,而dim_stddev在计算时自动考虑缺测情况
2. 采用坐标方式索引数组,如果坐标越界,例如1974年没有5月数据,则能取多少取多少
密码修改失败请联系微信:mofangbao
发表于 2014-2-26 16:20:14 | 显示全部楼层
另外,生成类似(/19740501, ..., 20131031/)数组可以用下面的方法:

julrange = greg2jul((/1974, 2013/), (/5, 10/), (/1, 31/), (/-1, -1/))
jul = ispan(julrange(0), julrange(1), 1)
greg = jul2greg(jul)

result = greg(:, 0) * 10000 + greg(:, 1) * 100 + greg(:, 2)
密码修改失败请联系微信:mofangbao
 楼主| 发表于 2014-2-27 16:35:02 | 显示全部楼层

这个生成的是从19740501~20131013之间的所有月份的日期吧?我只要这些年5月~10月的呢~
密码修改失败请联系微信:mofangbao
发表于 2014-2-27 20:28:48 | 显示全部楼层
deeli 发表于 2014-2-27 16:35
这个生成的是从19740501~20131013之间的所有月份的日期吧?我只要这些年5月~10月的呢~

do y = 1974, 2013
    temp =  olr1d({y * 10000 + 501 : y * 10000 + 1031})
    olr2d({temp&time % 10000}, y - 1974) = temp
end do

这个循环就是在提取每年5月1日到10月31日啊!
密码修改失败请联系微信:mofangbao
发表于 2014-2-27 20:39:07 | 显示全部楼层
deeli 发表于 2014-2-27 16:35
这个生成的是从19740501~20131013之间的所有月份的日期吧?我只要这些年5月~10月的呢~

哦,不好意思,刚明白你什么意思。程序里面下面这句是有问题。
olr2d&time = (/19740501, ..., 20131031/) % 10000
应该是
olr2d&time = (/19740501, ..., 19741031/) % 10000
密码修改失败请联系微信:mofangbao
 楼主| 发表于 2014-2-27 20:48:53 | 显示全部楼层
longlivehj 发表于 2014-2-27 20:39
哦,不好意思,刚明白你什么意思。程序里面下面这句是有问题。
olr2d&time = (/19740501, ..., 20131031 ...

查了一些资料才知道%是取整,可是,还是不太明白,如果%10000到底是取的哪一部分的整呢?
密码修改失败请联系微信:mofangbao
发表于 2014-2-27 21:17:30 | 显示全部楼层
本帖最后由 longlivehj 于 2014-2-27 21:19 编辑
deeli 发表于 2014-2-27 20:48
查了一些资料才知道%是取整,可是,还是不太明白,如果%10000到底是取的哪一部分的整呢?

%是取余,比如19740501 % 10000,结果应该是501。olr2d&time = (/19740501, ..., 19741031/) % 10000的结果是(/501, 502, ..., 1030, 1031/),实际上就是想表示月和日。
密码修改失败请联系微信:mofangbao
 楼主| 发表于 2014-2-27 21:27:06 | 显示全部楼层
longlivehj 发表于 2014-2-27 21:17
%是取余,比如19740501 % 10000,结果应该是501。olr2d&time = (/19740501, ..., 19741031/) % 10000的结 ...

额,还有一个应该很弱的问题,比如:我数据里面的olr的time是double型,然后是
units :hours since 1800-01-01 00:00:0.0  我可以用cd_calendar得到time的年月日格式,有可能将这种格式重新赋给olr&time吗?我尝试赋了下,直接是用的olr&time:=cd_calendar(olr&time,-2) 结果不太对,我知道数据类型不匹配,可是那些转换语句又不能巧妙的用进来~
密码修改失败请联系微信:mofangbao
 楼主| 发表于 2014-2-27 21:29:26 | 显示全部楼层
longlivehj 发表于 2014-2-27 21:17
%是取余,比如19740501 % 10000,结果应该是501。olr2d&time = (/19740501, ..., 19741031/) % 10000的结 ...

基础太不扎实了,百度了几篇ncl的手册,总感觉不够用,不知道前辈有没有什么ncl的教程,可否分享呢~
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

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

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

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