ncl处理grd数据格式 但是图是空白的

subtropical

如题 我有个dat文件，根据ctl描述用grads写了一个变量的grd主要是为了用ncl画图，用ncl画图的时候没有报错，而且数据应该也是有的，我print了一些看了下，虽然也有缺测，但大部分是有数据的，但是图是空白的，不知道为什么，想了好久也不知道怎么回事，请大神指教。。

dat文件的ctl如下：
dset D:\work2019\wrf_dat_read\wrf02d2018072612.dat
options  byteswapped
undef 1.e30
title  OUTPUT FROM WRF V3.8 MODEL
pdef  288 288 lcc  32.482  117.397  144.500  144.500  60.00000  30.00000  117.30000   3000.000   3000.000
xdef  747 linear  112.35242   0.01351351
ydef  606 linear   28.35665   0.01351351
zdef   17 levels  
1000.00000
975.00000
950.00000
925.00000
900.00000
850.00000
800.00000
750.00000
700.00000
650.00000
600.00000
550.00000
500.00000
400.00000
300.00000
200.00000
100.00000
tdef   25 linear 12Z26JUL2018      60MN      
VARS   19
XLAT           1  0  LATITUDE, SOUTH IS NEGATIVE (degree_north)
XLONG          1  0  LONGITUDE, WEST IS NEGATIVE (degree_east)
U             17  0  x-wind component (m s-1)
V             17  0  y-wind component (m s-1)
W             17  0  z-wind component (m s-1)
T2             1  0  TEMP at 2 M (K)
HGT            1  0  Terrain Height (m)
RAINC          1  0  ACCUMULATED TOTAL CUMULUS PRECIPITATION (mm)
RAINNC         1  0  ACCUMULATED TOTAL GRID SCALE PRECIPITATION (mm)
tc            17  0  Temperature (C)
td            17  0  Dewpoint Temperature (C)
td2            1  0  Dewpoint Temperature at 2m (C)
rh            17  0  Relative Humidity (%)
rh2            1  0  Relative Humidity at 2m (%)
u10m           1  0  Rotated wind component (m s-1)
v10m           1  0  Rotated wind component (m s-1)
slp            1  0  Sea Levelp Pressure (hPa)
dbz           17  0  Reflectivity (-)
cape          17  0  CAPE (J/kg)
ENDVARS
脚本如下：
begin
min_lat =  28.5
max_lat =35.0
min_lon = 113.0
max_lon = 121.0

diri   = "D:/work2019/wrf_dat_read/"
fili1   = "rh2.grd"
    f1=fbindirread(diri+fili1,0,(/10,747,606/),"float")
      
   mlon =747                              ; XDEF
   nlat =606                              ; YDEF
   ntim = 10                              ; time steps
  time  = new (ntim, float) ; generate a "time" variable,"No_FillValue"
  time = ispan(0,ntim-1,1)*1.
  time!0         = "time"
  time@long_name = "time"
  time@units     = "day"           ; "yyyy.fraction_of_year"
  time&time      = time
                                         ; generate lat/lon
  lon = ispan(0,mlon-1,1)*0.0135 +112.35
  lon!0 = "lon"
  lon@long_name = "longitude"
  lon@units     = "degrees_east"
  lat = ispan(0,nlat-1,1)*0.0135+28.35;-90.
  lat!0 = "lat"
  lat@long_name = "latitude"
  lat@units     = "degrees_north"
     rh2=new((/10,747,606/),"float")
      rh2=f1(:,:,:)     
         rh2!0="time"  
          rh2!1="lon"
          rh2!2="lat"
         rh2&time=time
         rh2&lon=lon
         rh2&lon@units= "degrees_east"
          rh2&lat=lat  
           rh2&lat@units= "degrees_north"
          rh2@_FillValue=-9.99e+08        
   printVarSummary(rh2)  
   wks = gsn_open_wks("ps","rh2")  
gsn_define_colormap(wks,"precip_11lev")
  
res=True
res@gsnDraw                        = False
res@gsnFrame                       = False
;res@gsnMaximize                    = True
res@tmXBLabelFont                               = 21
res@tmYLLabelFont                               = 21
res@tmXBLabelFontHeightF                        = 0.02
res@tmYLLabelFontHeightF                        = 0.02
res@tmYROn                                      = False
res@tmXTOn                                      = False
res@tmYRLabelsOn                                = False
res@tmXTLabelsOn                                = False
res@tmYRMode                                    = "Automatic"
res@tmXBMinorOn                                 = False   
res@tmBorderThicknessF                          = 4
;res1@tmXBMajorThicknessF                         = 1.8
;res1@tmYLMajorThicknessF                         = 1.8
res@tiYAxisOn                                   = True
res@tiYAxisString                               = ""
res@tiYAxisFont                                 = 21
res@tiYAxisFontHeightF                          = 0.022
res@tiXAxisString                               = ""  ;Latitude
res@tiXAxisFont                                 = 21
res@tiXAxisFontHeightF                          = 0.022
res@gsnLeftString=""
res@gsnRightString=""
res@trGridType="TriangularMesh"
mpres                             = res
mpres@mpDataSetName               = "Earth..4"
mpres@mpDataBaseVersion           = "Ncarg4_1"
mpres@mpOutlineOn                 = False
mpres@mpMaskAreaSpecifiers       = (/"China"/)
mpres@mpGeophysicalLineThicknessF = 1
mpres@mpNationalLineThicknessF    = 1
mpres@mpCenterLonF         = 180.0
mpres@mpMinLatF         = min_lat
mpres@mpMaxLatF         = max_lat
mpres@mpMinLonF         = min_lon
mpres@mpMaxLonF         = max_lon
map= gsn_csm_map(wks,mpres)
resp= res
resp@gsnDraw  = False
resp@gsnFrame= False         
resp@cnFillPalette = "precip_11lev"
resp@cnFillOn = True     
;resp@cnLevelSelectionMode = "ManualLevels"      ; manually set cn levels
;resp@cnMinLevelValF       = 50              ; min level
;resp@cnMaxLevelValF       =  100               ; max level
;resp@cnLevelSpacingF=10
resp@cnLinesOn =True  
resp@cnLineLabelsOn = False     ; do not draw contour lines
resp@gsnLeftString=""
res@gsnLeftString=""
resp@gsnRightString=""
resp@lbLabelBarOn =True     ; Turn off labelbar
resp@cnInfoLabelOn = True
resp@lbOrientation="vertical"
; Turn off informational label
plot1   = gsn_csm_contour(wks,rh2(2,:,:),resp)
print(rh2(1,{114:120},{30:35}))
cnres           = True
cnres@china     = True       ;draw china map or not
cnres@river     = False       ;draw changjiang&huanghe or not
cnres@province  = True      ;draw province boundary or not
cnres@nanhai    = False       ;draw nanhai or not
cnres@diqu      = False       ; draw diqujie or not
chinamap = add_china_map(wks,plot1,cnres)
overlay(map,plot1)
draw(map)
frame(wks)
end

werewolf

直接拿wrfout.nc来画图

subtropical

werewolf 发表于 2019-7-10 11:58
直接拿wrfout.nc来画图

这是别人给的dat数据，我没有wrfout.nc...........

werewolf

subtropical 发表于 2019-7-10 12:48
这是别人给的dat数据，我没有wrfout.nc...........

那就用grads处理

subtropical

werewolf 发表于 2019-7-10 13:46
那就用grads处理

觉得grads画的不好看，改用ncl画画看，主要是想弄明白为什么这个没有值显示出来。。。。。

werewolf

subtropical 发表于 2019-7-10 16:57
觉得grads画的不好看，改用ncl画画看，主要是想弄明白为什么这个没有值显示出来。。。。。

用grads把数据输出为nc。然后用ncl读格点数据

subtropical

werewolf 发表于 2019-7-11 08:40
用grads把数据输出为nc。然后用ncl读格点数据

grads能输出nc格式的吗？我只看到写成grd的，请大神指教。谢谢

subtropical

werewolf 发表于 2019-7-11 08:40
用grads把数据输出为nc。然后用ncl读格点数据

之前是这样写grd的，
'reinit'
'open D:\work2019\wrf_dat_read\wrf02d2018072612.ctl'
'set fwrite d:\work2019\wrf_dat_read\rh2.grd
'set gxout fwrite'
'set z 1'
'set t 1 10'
'd rh2'
'disable fwrite'
;
我试着把第三句改成rh2.nc,但是ncl报错
fatal:NetCDF: Unknown file format
fatal:Could not open (d:/work2019/wrf_dat_read/rh2.nc)

兰溪之水

参考这个http://cola.gmu.edu/grads/gadoc/gradcomdsdfwrite.html

subtropical

兰溪之水 发表于 2019-7-11 15:53
参考这个http://cola.gmu.edu/grads/gadoc/gradcomdsdfwrite.html

谢谢兰溪大神指点 我看了之后 试着写了 也显示生成了nc文件，但是我用ncl画图的时候还是不对
这是写一个时次一个变量nc的gs：
'reinit'
'open D:\work2019\wrf_dat_read\wrf02d2018072612.ctl'
'set gxout fwrite'
'set x 1 747'
'set y 1 606'
'set lev 850'
'set t 1'
'define dbz1 = dbz'
'set sdfwrite d:\work2019\wrf_dat_read\dbz1.nc'
'sdfwrite dbz1'
'clear sdfwrite'
;
这是ncl读取脚本：
begin
min_lat =  28.5
max_lat =35.0
min_lon = 113.0
max_lon = 121.0
f1=addfile("d:/work2019/wrf_dat_read/dbz1.nc","r")
print(f1)
  lat=f1->latitude
lon=f1->longitude
;time   = f1->time
dbz = f1->dbz1(:,:)
dbz@_FillValue= -999000000
print(dbz({30:35},{114:120}))
print的结果数据是错的感觉，f1文件好像是正常的，哎。。。。请大神帮忙看看。。