NCL求含缺测值数据的90%分位数

837078493

求取12天数据的全部空间分布的90%分位数，思路是，先将12天的数据合并为一组数据，在将这组数据的缺测值删掉，保留剩余值，将剩余值按大小顺序排列，求数据个数*90%序号所对应的数值。但是程序总是提示错误：黄色部分提示fatal:VarVarWrite: Number of dimensions on left hand side does not match right hand side

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
;************************************************
begin
refa1 = new((/12,123,124/),float)
refa2 = new((/12,123,124/),float)
banda1 = new((/12,123,124/),float)
banda2 = new((/12,123,124/),float)
ndvia = new((/12,123,124/),float)

;************************************************
do xuhao = 1,12
a  = addfile(""+xuhao+".hdf","r")
refa1(xuhao-1,:,:) = a->sur_refl_b01(:,:)
refa2(xuhao-1,:,:) = a->sur_refl_b02(:,:)
banda1= refa1*0.0001
banda2= refa2*0.0001
ndvia=(banda2-banda1)/(banda2+banda1)
end do


ndvi = new(183024,"float")
ndvi=ndtooned(ndvia)

ndvi2=new(182845,"float")
ndvi2(i)=ndvi(ind(.not.ismissing(ndvi)))
n=0
do i=0,182844
if (ndvi2(i).gt.0) then
ndvi3(n)=ndvi2(i)
n=n+1
end if
end do
print(n)

x =new(1,"float")
temp = new(1,"float")
temp(0) =0
x(0)=0
do i=1,n-1
x(0)=ndvi3(i)
do j=i+1,n
if (x(0).gt.ndvi3(j))  then
temp(0)=ndvi3(j)
ndvi3(j)=ndvi3(i)
ndvi3(i)=temp(0)
x(0)=ndvi3(i)
end if
end do
end do
print(ndvi3(n*0.05))
print(ndvi3(n*0.9))

另外还有一个问题：在读取12天数据之后输出ndvi及banda1的属性为9.96921e+36，而只读取1天数据时这两者属性为-28672，这与refa1是一致的，想知道这是为什么？

zby

ndvi2(i)=ndvi(ind(.not.ismissing(ndvi)))应该写成ndvi2=ndvi(ind(.not.ismissing(ndvi)))，你那个i是什么意思。。。。

zby

你这个做的有点麻烦了，你将数组转化成一维不含缺测之后，可以直接用ncl中的qsort命令就可以自动排序了

kisscate

学习了