拆分模式的daily数据为每月数据

王磊

（一）问题描述：
       前几天接到老师的一个任务，绘制春夏秋冬各个季节的降水概率密度函数。所给的数据为气候模式daily数据（9年资料，存储在一个文件，时间维长365*9），模式的calendar为noleap（每年都是积分365天）。采用ncl的内置函数pdfx来计算概率密度函数。需要做的工作是从模式中挑出各个季节的数据参与计算。
（二）解决思路
      分割原始数据文件，每年的一个月存储为一个文件，然后再读取各个季节的数据计算概率密度函数。这样拆分也将使得我由各个月的资料计算各个季节的平均值（空间分布）变得简单许多。
（三）程序思路及要求

     （a）确定数据的年份长度，每月的天数（2月都是28天）
     （b）使用ut_calendar来获取所读取的月份的第一天的日期，作为输出文件名的一部分
     （c）最好具备批量处理功能，使得一次可以处理多个数据。

（四）程序代码


load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"


begin


FileIn="../../data/c384lin54._ctrl_PR.91-01.nc"
fi=addfile(FileIn, "r")

Ninterval  =1                            ;when your input data is too large,you can set this one bigger to get some piece of the data.
time   =fi->time(:)
lat    =fi->lat(::Ninterval)
lon    =fi->lon(::Ninterval)
Nleap  =365                                 ;in our model data, calendar is Noleap
Nyear  =dimsizes(time)/Nleap
Nmonth =12

Nlat   =dimsizes(lat)
Nlon   =dimsizes(lon)
Monthday  =(/0,31,28,31,30,31,30,31,31,30,31,30,31/)
Month     =(/"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"/)


casename="c384lin54"
VarInName=(/"PR"/)
VarOutName=(/"prect"/)          ;you can quote different varible to batch processing


;Nyear=1                               ;you can set Nyear equals 1 to test the first year
prec=new((/Nmonth*Nyear,Nlat,Nlon/), float)
OutTime=new(Nmonth*Nyear, double)
OutTime@units= time@units
OutTime@calendar =365       ;this is a must to get the actual time(such as 19991203)


Nvar=dimsizes(VarOutName)
do Ivar=0,Nvar-1                ;the number of the varible you want to split
n=0                                   ;the first month
do Iyear = 1, Nyear
    do Imon=1,Nmonth

        time1=get_cpu_time()        

    prect    = fi->PR(365*(Iyear-1)+sum(Monthday(:Imon-1)) : 365*(Iyear-1)+sum(Monthday(:Imon))-1,::Ninterval,::Ninterval)
    ;printVarSummary(prect)

    YearMonth=ut_calendar(prect&time(0), -1)     ;get the year and month ,like 199109
    OutFileName =casename +"."+ YearMonth +"."+ VarOutName(Ivar) +".daily.nc"
    system("/bin/rm -rf " + OutFileName)
    cdf_file = addfile(OutFileName,"c")
    cdf_file->$VarOutName(Ivar)$          = prect

    prec(n,:,:)=(/dim_avg_n_Wrap(prect, 0)/)     ;calculat  and store Month Mean
    OutTime(n) =prect&time(0)                                           ;store the first day of the month
        n=n+1
    delete([/prect,OutFileName,cdf_file/])
        
        time2=get_cpu_time()
        delt=time2-time1
        print("the time consuming in this loop is : "+ delt +"s")

    end  do
end do


;============================================================================
;store the Monthly meandatasets
prec!0="time"
prec!1="lat"
prec!2="lon"

prec&time=OutTime
prec&lat=lat
prec&lon=lon OutFileName2 =casename +"."+ YearMonth +"."+ VarOutName(Ivar) +".nc"
system("/bin/rm -rf " + OutFileName2)
cdf_file2 = addfile(OutFileName2,"c")
cdf_file2->$VarOutName(Ivar)$          = prec
end do

end




说明：上述程序处理的过程中，每年都是等长的，如果是观测（有闰年）又改如何处理呢？后面将给出代码的调整过程，但是思路十分相似。

王磊

对于观测的处理过程如下：
文件为每年一个。考虑闰年。


load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"


begin

   
FilePath ="../../data/OBS/prect/daily/3B42_daily_0.25X0.25/"
      
Files    = systemfunc ("/bin/ls " + FilePath)
;print(Files)
FilesIn  = addfiles(FilePath + Files, "r")   
ListSetType (FilesIn, "cat")   
Nfiles   =dimsizes(Files)      
;print(Nfiles)



Nyear=Nfiles                            ;the data are store by year
Ninterval  =1                        ;when your input data is too large,you can set this one bigger to get some piece of the data.

time   =FilesIn[0]->time(0)
lat    =FilesIn[0]->lat(::Ninterval)
lon    =FilesIn[0]->lon(::Ninterval)
Nlat   =dimsizes(lat)
Nlon   =dimsizes(lon)
Month     =(/"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"/)

;Nyear=1                                 ;you can set Nyear equals 1 to test the first year
Nmonth =12
year=new(Nyear,integer)


prec=new((/Nmonth*Nyear,Nlat,Nlon/), float)
OutTime=new(Nmonth*Nyear, double)
OutTime@units= time@units
OutTime@calendar ="standard"           ;this is a must to get the actual time(such as 19991203)


casename="Trmm_3B43"
VarInName=(/"r"/)
VarOutName=(/"prect"/)          ;you can quote different varible to batch processing
Nvar=dimsizes(VarInName)



do Ivar=0,Nvar-1
n=0
do Iyear = 0,Nyear-1

        time   =FilesIn[Iyear]->time(0)
    year(Iyear)=ut_calendar(time(0), -1)/100
   
        if (isleapyear(year(Iyear))) then
        Monthday  =(/0,31,29,31,30,31,30,31,31,30,31,30,31/)   ;for leap year ,there are 29 days in February.
        print(year(Iyear) + "is leap year")
    else
        Monthday  =(/0,31,28,31,30,31,30,31,31,30,31,30,31/)
        print(year(Iyear) + "is not leap year")
    end if


    do Imon=1,Nmonth
    time1=get_cpu_time()      
        prect    = FilesIn[Iyear]->r(sum(Monthday(:Imon-1)) : sum(Monthday(:Imon))-1,::Ninterval,::Ninterval)
                ;printVarSummary(prect)
                ;print(dimsizes(prect&time))
        YearMonth=ut_calendar(prect&time(0), -1)     ;get the year and month ,like 199109
        OutFileName =casename +"."+ YearMonth +"."+  VarOutName(Ivar)  +".daily.nc"
        system("/bin/rm -rf " + OutFileName)
        cdf_file = addfile(OutFileName,"c")
        cdf_file->$VarOutName(Ivar)$          = prect

        prec(n,:,:)=(/dim_avg_n_Wrap(prect, 0)/)     ;calculat  and store Month Mean
        OutTime(n) =prect&time(0)                                           ;store the first day of the month
        n=n+1

        delete([/prect,YearMonth,cdf_file/])
    time2=get_cpu_time()
        delt=time2-time1
        ;print("the time consuming in this loop is : "+ delt +"s")       
        end do
end do


;============================================================================
;store the Monthly meandatasets
prec!0="time"
prec!1="lat"
prec!2="lon"

prec&time=OutTime
prec&lat=lat
prec&lon=lon
OutFileName2 =casename +"."+ year(0) +"_to_"+ year(Nyear-1) +"."+ VarOutName(Ivar) +".Monthly_Mean.nc"
system("/bin/rm -rf " + OutFileName2)
cdf_file2 = addfile(OutFileName2,"c")
cdf_file2->$VarOutName(Ivar)$          = prec

end do

end

兰溪之水

直接用现成的函数不好？https://www.ncl.ucar.edu/Documen ... onthly_values.shtml

王磊

兰溪之水 发表于 2016-12-23 00:13
直接用现成的函数不好？https://www.ncl.ucar.edu/Document/Functions/Contributed/calculate_monthly_valu ...

师兄说的是，不过这段代码主要是想把年资料转化为按照月份存储，资料实际上还是daily资料，所以程序主要是为了拆分数据。如果是计算月平均值，肯定是不需要这么麻烦的，用这个帖子中的函数就能快速实现。

王磊

兰溪之水 发表于 2016-12-23 00:13
直接用现成的函数不好？https://www.ncl.ucar.edu/Document/Functions/Contributed/calculate_monthly_valu ...

这个函数链接给的好，以后计算气候态，绘制空间分布就很好用了。

兰溪之水

王磊 发表于 2016-12-23 08:25
师兄说的是，不过这段代码主要是想把年资料转化为按照月份存储，资料实际上还是daily资料，所以程序主要 ...

拆分可以试试cdo或者nco

andrewsoong

非常棒的工作！

王磊

andrewsoong 发表于 2016-12-23 11:40
非常棒的工作！

查阅了cdo手册，使用selmon，selseas 两个参数可以快速实现上述代码功能

张张会长大

请问大家daily函数怎么报错Undefined identifier: (calculate_daily_values) is undefined, can't continue