NCL求教，提取3000米高度上的格点计算平均位势高度，求大神赐教QAQ

hlyahs

非常感谢大家
本人想计算近20年来夏季逐月的，在550hPa上，高原一定范围且高度在3000米以上的格点的位势平均高度，奈何ncl新手菜鸟找不到合适的做法，所以希望各位能指点一二
以下是lz想出的傻白甜又错误百出的脚本QAQ，报错显示if这一块是错的

begin

    f = addfile("E:/ncldata/550mon678/zhishu/678mon.nc", "r") ;位势高度资料
    f2 = addfile("E:/ncldata/550mon678/dixinggaodu/dixinggaodu.nc", "r")；地形高度资料


    wsgd=short2flt(f->z(0:59,:,:))
    dixing=short2flt(f2->z(0,:,:))
    dixing = dixing/9.8

    if(dixing.gt.3000) then
    x  = dim_avg_n_Wrap(wsgd(:,{28:38},{75:100}),(/1,2/)) ;qu yu pingjun
    end if

    x  = x/9.8
    printVarSummary(x)
    print(x)

希望大家能给与指导，非常感谢！英语一般，还没有在官网找到类似的例子，如果有的话谢谢大家分享QAQ

songwei

使用where函数，将不满足条件的设置为缺测，再使用 x  = dim_avg_n_Wrap(wsgd(:,{28:38},{75:100}),(/1,2/))

hlyahs

songwei 发表于 2017-9-18 10:38
使用where函数，将不满足条件的设置为缺测，再使用 x  = dim_avg_n_Wrap(wsgd(:,{28:38},{75:100}),(/1,2/) ...

非常感谢您的回答，但是我在使用where函数时，显示错误：
fatal:where: condition variable (parameter 0) dimension mismatch with parameter 1
fatal:["Execute.c":8575]:Execute: Error occurred at or near line 21 in file /cygdrive/e/ncldata/550mon678/dixinggaodu/3000mzhishu.ncl

脚本修改如下：

    f = addfile("E:/ncldata/550mon678/zhishu/678mon.nc", "r")
    f2 = addfile("E:/ncldata/550mon678/dixinggaodu/dixinggaodu.nc", "r")

    wsgd=short2flt(f->z(0:59,:,:))
    dixing=short2flt(f2->z(0,:,:))
    dixing = dixing/9.8
  
    xwsgd = where(dixing .ge.3000, wsgd, wsgd@_FillValue)（这是添加的where语句）
    xave  = dim_avg_n_Wrap(xwsgd(:,{27:38},{76:97}),(/1,2/)) ;qu yu pingjun

    x  = xave/9.8
    printVarSummary(x)
    print(x)

where参考的脚本是ncl官网提供的例子：
orog = f1->orog
sftlf = f2->sftlf
t = f3->tas
u = f4->uas
v = f5->vas
land_only = orog
land_only = where(sftlf .ge.10, orog, orog@_FillValue)

看错误好像是维数不对，请问这该如何解决呢？
非常非常感谢各位大神帮助！

songwei

hlyahs 发表于 2017-9-18 13:28
非常感谢您的回答，但是我在使用where函数时，显示错误：
fatal:where: condition variable (parameter  ...

wsgd是三维的，dixing是二维的，可以自己定义一个dixing数组也是三维和wsgd对应起来即可。

hlyahs

songwei 发表于 2017-9-18 13:44
wsgd是三维的，dixing是二维的，可以自己定义一个dixing数组也是三维和wsgd对应起来即可。

嗯，正在努力研究维数定义中……担心一直不回复大神您跑掉了，要是待会折腾不出来，可能还要来叨扰您，大神请暂时不要抛弃我，谢谢您

hlyahs

songwei 发表于 2017-9-18 13:44
wsgd是三维的，dixing是二维的，可以自己定义一个dixing数组也是三维和wsgd对应起来即可。

大神大神，我好像算出来了
脚本如下：
    f = addfile("E:/ncldata/550mon678/zhishu/678mon.nc", "r")
    f2 = addfile("E:/ncldata/550mon678/dixinggaodu/dixinggaodu.nc", "r")



    wsgd=short2flt(f->z(0:59,:,:))
    dixing=short2flt(f2->z(0,:,:))

   
    xindixing=new((/60,241,480/),float)
    ndim = (/1,2/)
    tconform = conform(xindixing,dixing,ndim)
    tconform = tconform/9.8

   
    xwsgd = wsgd
    xwsgd = where(tconform .ge.3000, wsgd, wsgd@_FillValue)
    xave  = dim_avg_n_Wrap(xwsgd(:,{27:38},{76:97}),(/1,2/)) ;qu yu pingjun
   

    x  = xave/9.8
    printVarSummary(x)
    print(x)

大神，请问是这样算的吗？这样的数据就是3000米以上的选定区域位势高度平均了咩？

songwei

hlyahs 发表于 2017-9-18 16:07
大神大神，我好像算出来了
脚本如下：
    f = addfile("E:/ncldata/550mon678/zhishu/678mon.nc", "r" ...

你把新地形print一下 如果xindixing(0,:,:)和xindixing(1,:,:)一样就没问题了

hlyahs

songwei 发表于 2017-9-20 14:22
你把新地形print一下 如果xindixing(0,:,:)和xindixing(1,:,:)一样就没问题了

嗯嗯，是对的了，超级感谢大神！！{:eb303:}{:eb303:}

寻路人

谢谢楼主分享

寻路人

谢谢楼主分享