ncl计算R95P

旺仔女孩

分享一个计算R95P的程序
采用的R95P定义：一年之中大于等于1mm的日降水量RR（雨日）排序，取95百分位数为RR95,把 RR>RR95的降水量累加起来
百分位计算方法：这里采用的是算出来如果是X.x,则采用第x位和第x+1位线性插值出百分位数的方法（代码是白嫖家园一个帖子里的，谢谢TA）
本小白不一定是对的，仅供参考！！！！（错了不要来骂我），如果有大神发现哪里有错误，请一定赐教
CDO计算R95P的方法有没有大神可以分享一下，求求了，我算出来全是缺省值

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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "/wind1/home/tiansht3/shapefile_utils.ncl"
begin
;................计算百分位函数......................
undef("percent_to_value")
function percent_to_value( \
  i_data           : numeric, \
  i_percentiles: numeric  \
)
local None
begin
  retVal = new( dimsizes(i_percentiles), float )

  data1d = ndtooned( i_data )

if (all(ismissing(data1d)))then
  retVal=data1d@_FillValue
else

  notMissing = data1d( ind(.not.ismissing(data1d) ) )

  qsort(notMissing)

  do p = 0, dimsizes(i_percentiles)-1
;   pInd = round( i_percentiles(p) * .01 * dimsizes(notMissing) + 0.5, 3 ) -1
;   pInd = where( pInd.ge.dimsizes(notMissing), dimsizes(notMissing)-1, pInd )
    floatInd = i_percentiles(p) * .01 * dimsizes(notMissing) - 0.5
    floorInd = toint( floor(floatInd) )
    floorInd = where( floorInd.lt.0, 0, floorInd )
    ceilInd = toint( ceil(floatInd) )
    ceilInd = where( ceilInd.ge.dimsizes(notMissing), \
        dimsizes(notMissing)-1, ceilInd )
;   print(pInd + " " + dimsizes(notMissing))
    if( ceilInd.eq.floorInd ) then
      retVal(p) = notMissing(floorInd)
    else
      retVal(p) = notMissing(floorInd) * ( ceilInd - floatInd ) \
          + notMissing(ceilInd) * ( floatInd - floorInd )
    end if
  end do
end if
;........
  return(retVal)
end ; percent_to_value
;.............................................................
f = addfile("/wind1/home/tiansht3/PRE1/pr_cma_05_05_daily.nc","r")  ；这里是多年的日降水数据
time=f->time
TIME=cd_calendar(time,-2)
year=ispan(1961,2020,1)
it_s=new(60,string)
it_e=new(60,string)
rec_s=new(60,integer)
rec_e=new(60,integer)
R95Pyear=new((/60,72,128/),float) ;R95Pyear根据自己的降水数据自己设

do j=0,59
  it_s(j)=year(j)+"0101"
  it_e(j)=year(j)+"1231"
rec_s(j)=ind(it_s(j).eq.TIME)
rec_e(j)=ind(it_e(j).eq.TIME)
pr:=f->pr(rec_s(j):rec_e(j),:,:)      ;这里pr有可能是365，有可能是366，所以加个冒号
pr=where(pr.ge.1,pr,9.96921e+36)      ;先把雨日挑出来

size=dimsizes(pr)
R95P=new((/size(1),size(2)/),float)
per95=new((/size(1),size(2)/),float)
R95P=0
  do m=0,size(1)-1
    do n=0,size(2)-1
per95(m,n)=percent_to_value(pr(:,m,n),95)   ;这里是算95百分位是多少，算99分位就把95改成99
      do i=0,size(0)-1
if(ismissing(per95(m,n)))then
  R95P(m,n)=9.96921e+36
  else if(.not.ismissing(pr(i,m,n)).and.pr(i,m,n).gt.per95(m,n))then
  R95P(m,n)=R95P(m,n)+pr(i,m,n)         ;把大于95百分位数的降水量累加起来
  end if
end if
      end do
    end do
  end do
R95Pyear(j,:,:)=R95P
end do
printVarSummary(R95Pyear)            ；R95Pyear是60*72*128 每个格点上每一年的R95P
end

旺仔女孩

之前程序R95P的定义搞错了，现在的R95P的计算方法为：
是先计算气候态的极端降水阈值，取30年（我取的是1961-1990）的日降水量大于1mm的进行排序，算出第95%分位值，这个数就是气候态的极端降水阈值。然后算每一年的R95P的时候就是把每一年中大于气候态的极端降水阈值的降水量累加起来，就是每一年的R95P了
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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "/wind1/home/tiansht3/shapefile_utils.ncl"
begin
;................计算百分位函数......................
undef("percent_to_value")
function percent_to_value( \
  i_data           : numeric, \
  i_percentiles: numeric  \
)
local None
begin
  retVal = new( dimsizes(i_percentiles), float )

  data1d = ndtooned( i_data )

if (all(ismissing(data1d)))then
  retVal=data1d@_FillValue
else

  notMissing = data1d( ind(.not.ismissing(data1d) ) )

  qsort(notMissing)

  do p = 0, dimsizes(i_percentiles)-1
    floatInd = i_percentiles(p) * .01 * dimsizes(notMissing) - 0.5
    floorInd = toint( floor(floatInd) )
    floorInd = where( floorInd.lt.0, 0, floorInd )
    ceilInd = toint( ceil(floatInd) )
    ceilInd = where( ceilInd.ge.dimsizes(notMissing), \
        dimsizes(notMissing)-1, ceilInd )
;   print(pInd + " " + dimsizes(notMissing))
    if( ceilInd.eq.floorInd ) then
      retVal(p) = notMissing(floorInd)
    else
      retVal(p) = notMissing(floorInd) * ( ceilInd - floatInd ) \
          + notMissing(ceilInd) * ( floatInd - floorInd )
    end if
  end do
end if
;........
  return(retVal)
end ; percent_to_value
;.............................................................
f = addfile("/wind1/home/tiansht3/PRE1/pr_cma_05_05_daily.nc","r")  ;这里是多年的日降水数据
time=f->time
TIME=cd_calendar(time,-2)
year=ispan(1961,2020,1)
it_s=new(60,string)
it_e=new(60,string)
rec_s=new(60,integer)
rec_e=new(60,integer)
R95Pyear=new((/60,72,128/),float) ;R95Pyear根据自己的降水数据自己设

;.........计算1961_1990气候态95%阈值.............................
it_s1=new(1,string)
it_e1=new(1,string)
rec_s1=new(1,integer)
rec_e1=new(1,integer)
it_s1="19610101"
it_e1="19901231"
rec_s1=ind(it_s1.eq.TIME)
rec_e1=ind(it_e1.eq.TIME)
pr1:=f->pr(rec_s1:rec_e1,:,:)
pr1=where(pr1.ge.1,pr1,9.96921e+36)               ;先把雨日挑出来
size1=dimsizes(pr1)
per95=new((/size1(1),size1(2)/),float)
do m=0,size1(1)-1
    do n=0,size1(2)-1
per95(m,n)=percent_to_value(pr1(:,m,n),95)       ;;这个是经度*纬度的气候态的95%阈值
    end do
end do
print(per95(30:40,60:70))            
printVarSummary(per95)
;................................................................
do j=0,5                                      ;j是在循环年
  it_s(j)=year(j)+"0101"
  it_e(j)=year(j)+"1231"
rec_s(j)=ind(it_s(j).eq.TIME)
rec_e(j)=ind(it_e(j).eq.TIME)
pr:=f->pr(rec_s(j):rec_e(j),:,:)              ;这里pr有可能是365，有可能是366，所以加个冒号
; pr=where(pr.ge.1,pr,9.96921e+36)            ;先把雨日挑出来
size=dimsizes(pr)
R95P=new((/size(1),size(2)/),float)
R95P=0
  do m=0,size(1)-1
    do n=0,size(2)-1
      do i=0,size(0)-1                        ;i是在循环天
if(ismissing(per95(m,n)))then
  R95P(m,n)=9.96921e+36
  else if(.not.ismissing(pr(i,m,n)).and.pr(i,m,n).gt.per95(m,n))then     ;per95不参与时间循环     
  R95P(m,n)=R95P(m,n)+pr(i,m,n)               ;把大于95百分位数的降水量累加起来
  end if
end if
      end do
    end do
  end do
  print(R95P(35,60))                        ;输出某个格点的R95P
R95Pyear(j,:,:)=R95P                        
end do
printVarSummary(R95Pyear)                  ;R95Pyear是60*72*128 每个格点上每一年的R95P
print(R95Pyear(:,35,60))                   ;输出某个格点的前6年的R95P
end

旺仔女孩

我还算了其他极端降水的指数，还有八个，但是别的计算程序比较简单，如果有需要可以给我发私信

wangwangxiao

感谢大佬，我来学习一下

wangwangxiao

兄弟资料是哪里下载的啊，我之前需要几十年的找到的都是分开的资料

旺仔女孩

wangwangxiao 发表于 2021-8-18 16:20
兄弟资料是哪里下载的啊，我之前需要几十年的找到的都是分开的资料

降水数据是别人给我的，我也不是太清楚

又又又加零

楼主你好，你程序中设置的R95Pyear=new((/60,72,128/),float)的后两个数字是什么呢？是经纬度吗？

旺仔女孩

又又又加零 发表于 2021-10-11 21:56
楼主你好，你程序中设置的R95Pyear=new((/60,72,128/),float)的后两个数字是什么呢？是经纬度吗？

是经纬度，是格点数据

旺仔女孩

f1= addfile("/wind1/home/tiansht3/CN05.01/CN05.01_PRE_1961_2020_daily.nc","r")  
year=ispan(1980,2020,1)                                                        
time=f1->time
time:=cd_calendar(time,2)
it_s=new(41,string)
it_e=new(41,string)
rec_s=new(41,integer)
rec_e=new(41,integer)
RX1day    =new((/41,81,169/),float)
RX1d      =new((/81,169/),float)
CWD       =new(dimsizes(RX1day),float)
CDD       =new(dimsizes(RX1day),float)
PRECPTOT  =new(dimsizes(RX1day),float)
R25       =new(dimsizes(RX1day),float)
R10       =new(dimsizes(RX1day),float)
SDII      =new(dimsizes(RX1day),float)
cont2     =new(dimsizes(RX1day),float)
RX5d      =new(dimsizes(RX1day),float)
PRECPTOT=0
SDII=0
R25=0
R10=0  
do j=0,40
it_s(j)=year(j)+"0101"
it_e(j)=year(j)+"1231"
rec_s(j)=ind(it_s(j).eq.time)
rec_e(j)=ind(it_e(j).eq.time)
pr:=f1->pre(rec_s(j):rec_e(j),{30:50},{73:115})            ;先挑出一年的降水数据 365*81*169
pr1:=where(pr.ge.0.1,pr,-9999)                             ;雨日
size=dimsizes(pr)
;........CDW&CDD&RX5d&RX1day start..............                                 
x1:=new(dimsizes(pr),integer)
x2:=new(dimsizes(pr),integer)
cont1=new((/size(1),size(2)/),integer)
cont3=new((/size(1),size(2)/),integer)
total:=new(dimsizes(pr),float)
cont1=0             ;CDW的记数
cont3=0             ;CDD的记数
total=0
x1=0
x2=0
do m=0,80
    do n=0,168     
     ; .......RX1d............
        RX1d(m,n)=max(pr(:,m,n))  
     ;.......CDW.............   
      do i=0,size(0)-1
        if(.not.ismissing(pr(i,m,n)).and.pr(i,m,n).ge.1)
            cont1(m,n)=cont1(m,n)+1
        else
            x1(i,m,n)=cont1(m,n)
            cont1(m,n)=0
        end if
      ;.......CDD.............  
        if(.not.ismissing(pr(i,m,n)).and.pr(i,m,n).lt.1)
             cont3(m,n)=cont3(m,n)+1
        else
            x2(i,m,n)=cont3(m,n)
            cont3(m,n)=0
        end if
        end do
        CWD(j,m,n)=max(x1(:,m,n))
        CDD(j,m,n)=max(x2(:,m,n))
    end do
end do
; ........CWD&CDD&RX5d&RX1D end..............

;.........PRECPTOT&SDII&R25&R10 start..............                          
cont2=0                    
do m=0,80
    do n=0,168
        do i=0,size(0)-1
        if(.not.ismissing(pr1(i,m,n)))then ;用的是PR1，已经是雨日了
        cont2(j,m,n)=cont2(j,m,n)+1
        PRECPTOT(j,m,n)=pr1(i,m,n)+PRECPTOT(j,m,n)
        end if
        end do
        if(cont2(j,m,n).eq.0)then
            SDII(j,m,n)=0
        else
         SDII(j,m,n)=PRECPTOT(j,m,n)/cont2(j,m,n)
        end if

        do i=0,size(0)-1
        if(.not.ismissing(pr1(i,m,n)).and.pr1(i,m,n).ge.25)then
        R25(j,m,n)=R25(j,m,n)+1
        end if
        end do

        do i=0,size(0)-1
        if(.not.ismissing(pr1(i,m,n)).and.pr1(i,m,n).ge.10)then
        R10(j,m,n)=R10(j,m,n)+1
        end if  
        end do
    end do
end do
; ; ;............PRECPTOT&SDII&R25&R10 end.....................
RX1day(j,:,:)=RX1d
print(j)
end do
end

旺仔女孩

旺仔女孩 发表于 2021-12-24 10:34
f1= addfile("/wind1/home/tiansht3/CN05.01/CN05.01_PRE_1961_2020_daily.nc","r")  
year=ispan(1980,20 ...

以上是 PRECPTOT,RX1day,CWD,CDD,R25,R10,SDII的计算，如果有错误请指出