ncl 求特定时段标准差

deeli

数据资料：1974.6.1~2013.8.31逐日olr资料，时间格式19740601
我现在想计算每一年5月~10月的olr标准差，求标准差的语句我知道了，用stddev，但是这个时间循环还是没弄明白，怎么才能每一年5~10月的提取出来分别算呢？time维数直接是14337~新手求教~~

longlivehj

; 一维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月数据，则能取多少取多少

longlivehj

另外，生成类似(/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)

deeli

longlivehj 发表于 2014-2-26 16:20

                               
登录/注册后可看大图

另外，生成类似(/19740501, ..., 20131031/)数组可以用下面的方法：

julrange = greg2jul((/1974, 2013/ ...

这个生成的是从19740501~20131013之间的所有月份的日期吧？我只要这些年5月~10月的呢~

longlivehj

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日啊！

longlivehj

deeli 发表于 2014-2-27 16:35

                               
登录/注册后可看大图

这个生成的是从19740501~20131013之间的所有月份的日期吧？我只要这些年5月~10月的呢~

哦，不好意思，刚明白你什么意思。程序里面下面这句是有问题。
olr2d&time = (/19740501, ..., 20131031/) % 10000
应该是
olr2d&time = (/19740501, ..., 19741031/) % 10000

deeli

longlivehj 发表于 2014-2-27 20:39

                               
登录/注册后可看大图

哦，不好意思，刚明白你什么意思。程序里面下面这句是有问题。
olr2d&time = (/19740501, ..., 20131031 ...

查了一些资料才知道%是取整，可是，还是不太明白，如果%10000到底是取的哪一部分的整呢？

longlivehj

deeli 发表于 2014-2-27 20:48

                               
登录/注册后可看大图

查了一些资料才知道%是取整，可是，还是不太明白，如果%10000到底是取的哪一部分的整呢？

%是取余，比如19740501 % 10000，结果应该是501。olr2d&time = (/19740501, ..., 19741031/) % 10000的结果是(/501, 502, ..., 1030, 1031/)，实际上就是想表示月和日。

deeli

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) 结果不太对,我知道数据类型不匹配，可是那些转换语句又不能巧妙的用进来~

deeli

longlivehj 发表于 2014-2-27 21:17

                               
登录/注册后可看大图

%是取余，比如19740501 % 10000，结果应该是501。olr2d&time = (/19740501, ..., 19741031/) % 10000的结 ...

基础太不扎实了，百度了几篇ncl的手册，总感觉不够用，不知道前辈有没有什么ncl的教程，可否分享呢~