求教利用ncl赋值问题

Jaquelinluo

假如存在一个nc文件，其中包含变量var（lat,lon）。现在我希望通过写ncl脚本，创建另一个nc文件，其中包含变量temp（time,lat,lon）,而变量temp所对应的每个值是相应变量var对应的每个值的0.2倍。为了实现这个，我写的一段是这样的：
f->temp = (/fi->$var$/)*0.2
但是运行的时候，总是会提示出错：
fatal:Dimension sizes of left hand side do not match right hand side
fatal:["Execute.c":8128]:Execute: Error occurred at or near line 56 in file ch4.ncl
我知道是两边变量的维数不一致，但是我想不到其他方法实现这一功能，所以想请教各位大神，这种情况的话该怎么处理呢？谢谢大家啦～

longlivehj

看你的描述好像是想把多个var合并成一个temp，增加一个time维。那为什么不把temp数组做好了，一次性写到nc文件里面去呢？

Jaquelinluo

longlivehj 发表于 2014-2-24 15:11

                               
登录/注册后可看大图

看你的描述好像是想把多个var合并成一个temp，增加一个time维。那为什么不把temp数组做好了，一次性写到nc文 ...

其实我就是想把var的每个值乘0.2，再赋给temp，但是var是二维的，temp是三维的。

longlivehj

Jaquelinluo 发表于 2014-2-24 15:17

                               
登录/注册后可看大图

其实我就是想把var的每个值乘0.2，再赋给temp，但是var是二维的，temp是三维的。

do i = ...
    ...
    fi = addfile(...)
    temp(i, :, :) = (/fi->$var$/)*0.2
    ....
end do

f->temp = temp

longlivehj

Jaquelinluo 发表于 2014-2-24 15:17

                               
登录/注册后可看大图

其实我就是想把var的每个值乘0.2，再赋给temp，但是var是二维的，temp是三维的。

或者，用addfiles把所有包含var的nc文件一次性打开，可以直接得到三维的temp数组。比如：

p = ... ; 所有输入nc文件路径的字符串数组
fi = addfiles(p, "r")
ListSetType(fi, "join")
temp = fi[:]->$var$
f->temp = temp

Jaquelinluo

longlivehj 发表于 2014-2-24 16:06

                               
登录/注册后可看大图

或者，用addfiles把所有包含var的nc文件一次性打开，可以直接得到三维的temp数组。比如：

p = ... ; 所 ...

哦哦，这样，好的。我想再请教下，nc里面是不是有变量个数的限制呢？我在里面写了15个变量，到第11个变量的时候就会显示出错：
One or more variable sizes violate format constraints.
是不是最多只能写10个变量呢？

longlivehj

Jaquelinluo 发表于 2014-2-24 16:51

                               
登录/注册后可看大图

哦哦，这样，好的。我想再请教下，nc里面是不是有变量个数的限制呢？我在里面写了15个变量，到第11个变量 ...

你按下面链接里面的提示试试，看行不行。
http://www.ncl.ucar.edu/FAQ/#file_io_004

Jaquelinluo

longlivehj 发表于 2014-2-24 17:00

                               
登录/注册后可看大图

你按下面链接里面的提示试试，看行不行。
http://www.ncl.ucar.edu/FAQ/#file_io_004

哦哦，好的。谢谢～

Einstein01

循环读入文件并增加时间维有个问题，图形在循环内画和在循环外画，差别很大。。

simon7

longlivehj 发表于 2014-2-24 17:00
你按下面链接里面的提示试试，看行不行。
http://www.ncl.ucar.edu/FAQ/#file_io_004

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 "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
    f = addfile("../wrfout_d02_2009-09-16_06:00:00.nc","r")
    times  = wrf_user_list_times(f)  ; get times in the filewww.mnmuc.org)
    ntimes = dimsizes(times)         ; number of times in the file
    tc = wrf_user_getvar(f,"tc",-1)
    p  = wrf_user_getvar(f, "pressure",-1) ;
    rh = wrf_user_getvar(f,"rh",-1)
    z    = wrf_user_getvar(f, "z",-1)
    pressure_levels = (/ 1000.,950.,900.,850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,350.,300.,250.,200.,150.,100./)
    nlevels         = dimsizes(pressure_levels)
    tc_plane = wrf_user_intrp3d(tc,p,"h",pressure_levels,0.,False)
    z_plane  = wrf_user_intrp3d(z,p,"h",pressure_levels,0.,False)
    rh_plane = wrf_user_intrp3d(rh,p,"h",pressure_levels,0.,False)
    p_plane  = wrf_user_intrp3d(p,p,"h",pressure_levels,0.,False)

    ;t         = f->tc
    ;rh        = f->rh
    ;pre_7     = f->lv_ISBL3                   ;(100,1000)   氓<8d><95>盲陆<8d>hPa
;#################################################################################################
;猫庐隆莽庐<97>氓<81><87>莽<9b>赂氓陆<93>盲陆<8d>忙赂?
;==========================
   ;printVarSummary(rh)
   ;printVarSummary(t)
   tk=tc_plane+273.166
  printVarSummary(pressure_levels)
   p = conform(tc_plane, pressure_levels, 1)
   a1 = where(tk .gt. 263.0, 0.622*6.11*exp(17.26*(tk-273.16)/(tk-35.86)), \
                     0.622*6.11*exp(21.87*(tk-273.16)/(tk-7.66)))   ;
   b1= where(tk .gt. 263.0, P-0.278*exp(17.26*(tk-273.16)/(tk-35.86)),\
                 P-0.278*exp(21.87*(tk-273.16)/(tk-7.66)))
   ;prs   = conform(t,pre_7,0)
   ;printVarSummary(p)
   ;es    = 6.112*exp(17.67*(t-273.15)/(t-29.65))
   ;q     = rh*(0.62197*es/(prs-es))/100.
   q = a1/b1
   q1 = q*rh_plane
   e     = p*q1/100./(0.62197+q1/100.)+1e-10
   tlcl  = 55.0+2840.0/(3.5* log(tk)-log(e)-4.805)
   theta = tk*((1000./p)^(0.2854*(1.0-0.28*q1/100.)))
   eqt   = theta*exp(((3376./tlcl)-2.54)*q1/100.0*(1.0+0.81*q/100.0))
   copy_VarCoords(t,eqt)
   printVarSummary(eqt)

   ; system("rm -f eqttest.nc")   ; remove any pre-existing file
   ;    ncdf     = addfile("eqttest.nc" ,"c")  ; open output netCDF file
   ;    ncdf->eqt    = eqt
;##########################################################
;PLOT
;=========================
  wks = gsn_open_wks("x11","h_lat_7-21-00-test")
  do it =0, ntimes-1
  res                      = True
  res@gsnDraw              = False         ; Do not draw plot
  res@gsnFrame             = False         ; Do not advance frome
  res@gsnMaximize      = True
res@cnLevelSelectionMode = "ManualLevels"      ; manual contouring
  res@cnMinLevelValF       = 326               ; set min contour level
  res@cnMaxLevelValF       = 408               ; set max contour level
  res@cnLevelSpacingF      = 2                   ; set contour spacing
  res@cnInfoLabelOn        = False
  res@cnLineThicknessF     = 2
  res@cnLineLabelDensityF  = 0.7

  ;res@cnSmoothingOn        = True
  ;res@cnSmoothingDistanceF = 0.005
  ;res@cnSmoothingTensionF  = -2.5


  plot3  = gsn_csm_pres_hgt(wks,eqt(it,:,{30.5:32.5},{120}),res)
  draw(plot3)
  frame(wks)
  end do
end


请教一下，这个脚本画假相当位温剖面图，总是报错，报错的信息为fatal:Dimension sizes of left hand side and right hand side of assignment do not match
fatal:["Execute.c":8567]:Execute: Error occurred at or near line 31 in file test_theta.ncl，
31行为标红的那行：p = conform(tc_plane, pressure_levels, 1)
请问该如何解决呢，麻烦了啊