weekly数据处理成monthly数据

我喜欢风

    之前的程序有问题，于2015-01-17_01:00修改
    前言：lz由于工作需要，要用到积雪资料（snow cover），老板推荐使用Rutgers University Global Snow Lab的资料，于是就愉快地去下载了。 http://climate.rutgers.edu/snowcover/docs.php?target=datareq        资料下载之后发现只有weekly的资料，没有monthly的资料。在网站上纠结了很久之后，lz放弃了，就打算自己处理下。

    处理方法参照网站给的思路（如上），monthly资料由weekly资料进行权重平均得到，权重为某一week包含在该month里的天数除以该month的总天数。例如：某week在2001年6月里有3天，这个月总共有31天，那么其权重就是3/31；若某week在2001年6月里有7天，则权重为7/31.

    NCL程序如下，完整程序见附件

1. ;=====start day==========
   
2. y=1966
   
3. m=10
   
4. d=4
   
5. ;=====end day============
   
6. ye=2015
   
7. me=1
   
8. de=5
   
9. ;========================
   
10. md=(/-1,31,28,31,30,31,30,31,31,30,31,30,31/)
    
11. 
    
12. i=0     ;time index
    
13. nn=0    ;number of weeks contained in a month
    
14. rn=0    ;residual of a week
    
15. rn2=0   ;residual of the last week of a month
    
16. 
    
17. ;=====deal with the first month======
    
18. if(d.gt.1)
    
19.     md(2)=where(isleapyear(y),29,28)  ;leap year
    
20.     nn=(md(m)-d+1)/7
    
21.     rn=(md(m)-d+1-nn*7)
    
22.     m=m+1
    
23.     y=where(m.gt.12,y+1,y)
    
24.     m=where(m.gt.12,m-12,m)
    
25. end if
    
26. 
    
27. ;=====deal with the last month======
    
28. md(2)=where(isleapyear(ye),29,28)
    
29. if(de.lt.md(me))
    
30.     me=me-1
    
31.     ye=where(me.lt.1,ye-1,ye)
    
32.     me=where(me.lt.1,me+12,me)
    
33. end if
    
34. ;====================================
    
35. 
    
36. ;=====creat the new var for monthly data======
    
37. ;!!!!!!depended on the dimsizes of weekly var
    
38. months=(ye-y)*12+me-m+1
    
39. MyTime=ispan(m,months-1+m,1)
    
40. MyTime_temp=(y+MyTime/12)*100+MyTime-MyTime/12*12
    
41. MyTime=where((MyTime-MyTime/12*12).eq.0,MyTime_temp-88,MyTime_temp)
    
42. delete(MyTime_temp)      ;Change time format to yyyymm, for example, 200402
    
43. 
    
44. m_var=new((/months,d_var(1),d_var(2)/),"float")
    
45. copy_VarAtts(var,m_var)
    
46. m_var@valid_time=y+tostring(m)+"_"+ye+tostring(me)+"_monthly_from_weekly"
    
47. m_var@method="montly values are calulated by weighting the weekly areas "+\
    
48. "according to the number of days of a map week falling in the given month."
    
49. ;=============================================
    
50. 
    
51. ;=====calculate monthly from weekly===========
    
52. wi=nn  ;weekly start index
    
53. do i=0,months-1
    
54.     md(2)=where(isleapyear(y),29,28)
    
55.     nn=(md(m)-(7-rn))/7
    
56.     rn2=md(m)-nn*7-(7-rn)
    
57.     temp=(7-rn)*1.0/md(m)*var(wi,:,:) ;week start of a month
    
58.     temp=temp+dim_sum_n(var(wi+1:wi+nn,:,:),0)*7.0/md(m) ;week middle
    
59.     wi=wi+nn+1
    
60.     temp=temp+rn2*1.0/md(m)*var(wi,:,:)  ;week end
    
61.     m_var(i,:,:)=temp
    
62.     m=m+1
    
63.     y=where(m.gt.12,y+1,y)
    
64.     m=where(m.gt.12,m-12,m)
    
65.     rn=rn2
    
66. ;    print("y="+y+" m="+m+" md2="+md(2)+" nn="+nn+" rn="+rn+" wi="+wi)
    
67. end do
    
68. ;=============================================
    
69. 
    
70. ;=====Coordinations===========================
    
71. m_var!0="MyTime"
    
72. m_var&MyTime=MyTime
    
73. ;=============================================

复制代码

    检验下结果： （右图为网站提供的月平均积雪覆盖，左图为由weekly资料计算得到的月平均积雪覆盖）

sun_shine_Xia

nice~给楼主32个赞之后marl学习一下！

沙颖凯

是好方法，我对楼主所说的权重的理解是把weekly的值先转化成日值，作为可以度量月的基本单位，然后转化成月值。

我喜欢风

沙颖凯 发表于 2015-1-16 13:45
是好方法，我对楼主所说的权重的理解是把weekly的值先转化成日值，作为可以度量月的基本单位，然后转化成月 ...

确实如此。处理的时候是将日值用周值表示出来。
从网站上给出的信息来看，weekly数据是从daily数值平均后得到的，所以反过来用weekly数值代表每天的数据也是可取的。

我喜欢风

程序编写的时候让我头疼的主要是第一个月和最后一个月的处理，还有构建YYYYMM这样的时间格式MyTime数组。细节需要调试，可能程序对某些情况还有疏忽，希望对同样想把weekly数据处理成monthly数据的朋友有所帮助

wybnu

留名，感谢分享

gxlbxzx

发邮件问他们要  是可以得到月平均的   没必要自己处理

我喜欢风

gxlbxzx 发表于 2015-1-21 10:36
发邮件问他们要  是可以得到月平均的   没必要自己处理

果然如此么。。当时也想发邮件来着，后来懒得弄，看处理方式不难，就自己编了一个-0-

shaoyanying

批量拆分txt格式的气象数据

六六的奶黄包

感谢分享，太赞了