请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 8959|回复: 19

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

[复制链接]

新浪微博达人勋

发表于 2019-7-10 10:57:57 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
如题 我有个dat文件,根据ctl描述用grads写了一个变量的grd主要是为了用ncl画图,用ncl画图的时候没有报错,而且数据应该也是有的,我print了一些看了下,虽然也有缺测,但大部分是有数据的,但是图是空白的,不知道为什么,想了好久也不知道怎么回事,请大神指教。。
10.png
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
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-7-10 11:58:58 | 显示全部楼层
直接拿wrfout.nc来画图
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-7-10 12:48:53 | 显示全部楼层
werewolf 发表于 2019-7-10 11:58
直接拿wrfout.nc来画图

这是别人给的dat数据,我没有wrfout.nc...........
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-7-10 13:46:16 | 显示全部楼层
subtropical 发表于 2019-7-10 12:48
这是别人给的dat数据,我没有wrfout.nc...........

那就用grads处理
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-7-10 16:57:46 | 显示全部楼层
本帖最后由 subtropical 于 2019-7-10 17:31 编辑

觉得grads画的不好看,改用ncl画画看,主要是想弄明白为什么这个没有值显示出来。。。。。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-7-11 08:40:44 | 显示全部楼层
subtropical 发表于 2019-7-10 16:57
觉得grads画的不好看,改用ncl画画看,主要是想弄明白为什么这个没有值显示出来。。。。。

用grads把数据输出为nc。然后用ncl读格点数据
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-7-11 10:34:17 | 显示全部楼层
werewolf 发表于 2019-7-11 08:40
用grads把数据输出为nc。然后用ncl读格点数据

grads能输出nc格式的吗?我只看到写成grd的,请大神指教。谢谢
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-7-11 15:17:25 | 显示全部楼层
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)
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 成长值: 19710
发表于 2019-7-11 15:53:40 | 显示全部楼层
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-7-11 19:25:41 | 显示全部楼层
本帖最后由 subtropical 于 2019-7-11 19:32 编辑

谢谢兰溪大神指点 我看了之后 试着写了 也显示生成了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文件好像是正常的,哎。。。。请大神帮忙看看。。





11-1.png
11.png


密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表