- 积分
- 780
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-7-5
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 我喜欢风 于 2015-1-17 01:04 编辑
之前的程序有问题,于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程序如下,完整程序见附件 - ;=====start day==========
- y=1966
- m=10
- d=4
- ;=====end day============
- ye=2015
- me=1
- de=5
- ;========================
- md=(/-1,31,28,31,30,31,30,31,31,30,31,30,31/)
- i=0 ;time index
- nn=0 ;number of weeks contained in a month
- rn=0 ;residual of a week
- rn2=0 ;residual of the last week of a month
- ;=====deal with the first month======
- if(d.gt.1)
- md(2)=where(isleapyear(y),29,28) ;leap year
- nn=(md(m)-d+1)/7
- rn=(md(m)-d+1-nn*7)
- m=m+1
- y=where(m.gt.12,y+1,y)
- m=where(m.gt.12,m-12,m)
- end if
- ;=====deal with the last month======
- md(2)=where(isleapyear(ye),29,28)
- if(de.lt.md(me))
- me=me-1
- ye=where(me.lt.1,ye-1,ye)
- me=where(me.lt.1,me+12,me)
- end if
- ;====================================
- ;=====creat the new var for monthly data======
- ;!!!!!!depended on the dimsizes of weekly var
- months=(ye-y)*12+me-m+1
- MyTime=ispan(m,months-1+m,1)
- MyTime_temp=(y+MyTime/12)*100+MyTime-MyTime/12*12
- MyTime=where((MyTime-MyTime/12*12).eq.0,MyTime_temp-88,MyTime_temp)
- delete(MyTime_temp) ;Change time format to yyyymm, for example, 200402
- m_var=new((/months,d_var(1),d_var(2)/),"float")
- copy_VarAtts(var,m_var)
- m_var@valid_time=y+tostring(m)+"_"+ye+tostring(me)+"_monthly_from_weekly"
- m_var@method="montly values are calulated by weighting the weekly areas "+\
- "according to the number of days of a map week falling in the given month."
- ;=============================================
- ;=====calculate monthly from weekly===========
- wi=nn ;weekly start index
- do i=0,months-1
- md(2)=where(isleapyear(y),29,28)
- nn=(md(m)-(7-rn))/7
- rn2=md(m)-nn*7-(7-rn)
- temp=(7-rn)*1.0/md(m)*var(wi,:,:) ;week start of a month
- temp=temp+dim_sum_n(var(wi+1:wi+nn,:,:),0)*7.0/md(m) ;week middle
- wi=wi+nn+1
- temp=temp+rn2*1.0/md(m)*var(wi,:,:) ;week end
- m_var(i,:,:)=temp
- m=m+1
- y=where(m.gt.12,y+1,y)
- m=where(m.gt.12,m-12,m)
- rn=rn2
- ; print("y="+y+" m="+m+" md2="+md(2)+" nn="+nn+" rn="+rn+" wi="+wi)
- end do
- ;=============================================
- ;=====Coordinations===========================
- m_var!0="MyTime"
- m_var&MyTime=MyTime
- ;=============================================
复制代码
检验下结果: (右图为网站提供的月平均积雪覆盖,左图为由weekly资料计算得到的月平均积雪覆盖)
|
评分
-
查看全部评分
|