关于维数定义问题的请教

MacTaylor

我是NCL初学者，写脚本的过程中碰到一个问题，请教大家，其中有问题的一段如下：

f1 = addfile("hgt.mon.mean.nc","r")
  ztemp = short2flt(f1->z)
  hgt = ztemp/9.80665
  hgt!0 = "time"
  hgt!1 = "lev"
  hgt!2 = "lat"
  hgt!3 = "lon"
  printVarSummary(hgt)
  delete(ztemp)
  epi = month_to_season(hgt,"JJA")
  printVarSummary(epi)
  delete(hgt)

然后我在后续计算的时候都没有出问题，但是生产的nc文件用grads打不开，提示信息：no discernable X corrdinate，我看了兰溪关于这个的帖子，知道是写成的nc数据不规范问题，我想请教大家的是：想输出标准定义的nc文件，那么这部分数组维数在写ncl的脚本中应该怎么定义？

mofangbao

帮顶了，好多人放假回家了，估计问题解决的不会很及时

hzwjy

从你这里的代码看，你只是对维度做了命名（name），但是并没有赋值。
虽然NCL在读取nc文件时会保留元数据，但是你中途处理的时候有元数据丢失：

  f1 = addfile("hgt.mon.mean.nc","r")
  ztemp = short2flt(f1->z)  ; 这里还有ztemp元数据
  hgt = ztemp/9.80665  ; 这里hgt就丢失了元数据
  hgt!0 = "time"
  hgt!1 = "lev"
  hgt!2 = "lat"
  hgt!3 = "lon"
  printVarSummary(hgt)
  delete(ztemp)
  epi = month_to_season(hgt,"JJA")
  printVarSummary(epi)
  delete(hgt)

其实可以用NCL自带函数将元数据复制给新变量
  copy_VarCoords(ztemp, hgt)
这样结果应该就OK了。

或者你可以自己设置time, lev, lat, lon数组，然后给维度赋值：
  lon = lonGlobeF(360, "lon", "lon", "degrees_east")
  hgt!3 = "lon"
  hgt&lon = lon

MacTaylor

hzwjy 发表于 2012-7-13 11:07

                               
登录/注册后可看大图

从你这里的代码看，你只是对维度做了命名（name），但是并没有赋值。
虽然NCL在读取nc文件时会保留元数据， ...

十分感谢！不好意思，还想多请教你一些，这段我自己改了下，结果倒是没问题了，我是这样写的：

  f1   = addfile("hgt.mon.mean.nc","r")
  hgt = short2flt(f1->z)
  hgt = floattoshort(hgt/9.80665)
  printVarSummary(hgt)
  epi = month_to_season(hgt,"JJA")
  printVarSummary(epi)

根据你说的，copy_VarCoords(ztemp, hgt)和我hgt = short2flt(f->z)实质上是一个意思，就是先将元数据数组的数值和属性纬度等都赋值给新变量hgt，这样新变量具有了原变量的所有属性，我的理解对吗？
另外你说的两种方法与我自己改的从脚本的书写标准角度上来说哪种更常用，更好一些呢？
多谢你的解答。

hzwjy

MacTaylor 发表于 2012-7-13 14:10

                               
登录/注册后可看大图

十分感谢！不好意思，还想多请教你一些，这段我自己改了下，结果倒是没问题了，我是这样写的：

  f1   ...

  hgt = short2flt(f1->z)
  hgt = floattoshort(hgt/9.80665)

你这里hgt先从文件中读取出来，然后short2flt把数据转换成float类型，这样是没问题的。
但你第二步的floattoshort这就不对了。float类型转换为short型会损失数据精度，所以你这样得到的结果是错的。

  hgt = short2flt(f1->z)
  hgt = hgt/9.80665
这样就可以了，如果hgt这里元数据没有丢失，就可以用来输出了。

顺带说一下，short2flt的作用是读取short的数据，并且读取scale和offset属性，
然后根据这两个属性计算后得到float的数据，换算公式是：
Float = Short*scale + offset
(如果没有数据中没scale和offset属性，那么上式中scale=1.  offset=0.)
这个函数的说明可以参考
http://www.ncl.ucar.edu/Document/Functions/Contributed/short2flt.shtml
以及 $NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl

MacTaylor

hzwjy 发表于 2012-7-13 16:30

                               
登录/注册后可看大图

hgt = short2flt(f1->z)
  hgt = floattoshort(hgt/9.80665)

十分感谢，高手啊，多多指教！

wsdeqx

受教了，谢神仙们

felin

hzwjy 发表于 2012-7-13 11:07
从你这里的代码看，你只是对维度做了命名（name），但是并没有赋值。
虽然NCL在读取nc文件时会保留元数据 ...

您好！请问您一个问题，我新建了一个格点数组，并按照您给的方法给维度赋值
                avetasi        = new((/nlati,nloni/),"float")       
                lat                  = latGlobeF(nlati, "lat", "latitude", "degrees_north")
                print(lat)
                avetasi!0        = "lat"
                avetasi&lat = lat
                lon                 = lonGlobeF(nloni, "lon", "longitude", "degrees_east")
                avetasi!1        = "lon"
                avetasi&lon = lon
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;               
                if (i .eq. 4)then        
                        avetasi                 = dim_avg_n(tasi(:1592,:,:),0)-273.15
                else
                        avetasi                 = dim_avg_n(tasi,0)-273.15
                end if
                printVarSummary(avetasi)
;                print(avetasi)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;                调整平均后的数据经度位置适合地图       
                newtasi                = new((/nlati,nloni/),"float")               
                lni                                = nloni - 1               ;less 1 nloni
                hni                                = nloni / 2               ;half nloni
                hlni                        = hni - 1                 ;half less 1 nloni
;                print(hlni)
;                print(hni)
;                print(lni)
                newtasi(:,0:hlni)                = avetasi(:,hni:lni)
                newtasi(:,hni:lni)        = avetasi(:,0:hlni)
                printVarSummary(newtasi)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;                取东亚面积               
                EA_newtasi                   = newtasi({20:50},{75:150})
                printVarSummary(EA_newtasi)               
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

其中之所以要调整数据经度位置，是因为我的数据和ncl默认世界地图在经度上不匹配，相差180度
而我把经度调整后，执行、报错如下
Variable: newtasi
Type: float
Total Size: 111360 bytes
            27840 values
Number of Dimensions: 2
Dimensions and sizes:        [lat | 145] x [lon | 192]
Coordinates:
            lat: [-90..90]
            lon: [180..178.125]
Number Of Attributes: 1
  _FillValue :        9.96921e+36
fatal:NclOneDValGetRangeIndex: Non-monotonic coordinate value being used, can't complete coordinate subscript
fatal:Could not obtain coordinate indexes, unable to perform subscript
fatal:["Execute.c":8640]:Execute: Error occurred at or near line 143 in file draw_EA_tas.ncl

错误位置也就是选取东亚经纬度的这一句
                EA_newtasi                   = newtasi({20:50},{75:150})
之所以是选取经纬度而不是格点，是应为各种模拟结果的格点不一
而在之前没有给数组赋值维度，报错总是出现newtasi的第一维没有命名以致变量坐标有误。
所以想请您帮忙，在脚本已经写成这样的情况下，如何解决上述错误？
同时也想请您指点，如何解决一直困扰我的数据和地图相差180°的问题？

木人

牛厉害了我的哥

QQtangxian

学习了，谢谢分享