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

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 22688|回复: 15

grads处理nc雷达数据,谢谢

[复制链接]

新浪微博达人勋

发表于 2021-4-1 09:26:44 | 显示全部楼层 |阅读模式
10金钱
本帖最后由 lxe3399 于 2021-4-1 09:31 编辑

有个nc格式的雷达数据,折腾了半天  用grads没画出正确的图来,恳请版上高人帮解决,非常感谢!下面是该数据的头文件:
netcdf 1606290700 {
dimensions:
longtitude = 1019 ;
latitude = 881 ;
layer = 20 ;

variables:
float longtitude(longtitude) ;
  longtitude:units = "degrees_east" ;
float latitude(latitude) ;
  latitude:units = "degrees_north" ;
long layer(layer) ;
  layer:units = "meter" ;
byte reflectivity(layer, latitude, longtitude) ;
  reflectivity:name = "Reflectivity" ;
  reflectivity:units = "dBZ" ;
  reflectivity:scale_factor = 0.5f ;
  reflectivity:add_offset = -33s ;
  reflectivity:data_type = "unsigned byte" ;
  reflectivity:formula = "dBZ=-33+0.5*N" ;
  reflectivity:valid_min = 0s ;
  reflectivity:valid_max = 255s ;
  reflectivity:no_data = 0s ;
  reflectivity:bad_data = 1s ;
  reflectivity:no_observation = 2s ;
  reflectivity:censored_data = 3s ;

// global attributes:
  :Mosaic_version = 1.f ;
  :QC_version = 1.f ;
  :Netcdf_mosaic_format_version = 1.f ;
  :Name = "B08_Mosaic" ;
  :Height_datum = "ASL" ;
  :NumOfRadars = 5s ;
  :Date = "20160629" ;
  :Time = "070000" ;
  :Year = 2016s ;
  :Month = 6s ;
  :Day = 29s ;
  :Hour = 7s ;
  :Minute = 0s ;
  :Second = 0s ;
  :FirstLon = 111.69f ;
  :FirstLat = 35.348999f ;
  :LastLon = 121.88f ;
  :LastLat = 44.159f ;
  :LonRes = 0.0099999998f ;
  :LatRes = 0.0099999998f ;
  :SiteName-1 = "BeiJing" ;
  :SiteID-1 = " " ;
  :SiteType-1 = "SA" ;
  :SiteLon-1 = 116.47f ;
  :SiteLat-1 = 39.810001f ;
  :SiteHgt-1 = 31.f ;
  :ObsRange-1 = 300s ;
  :SiteName-2 = "TianJin" ;
  :SiteID-2 = " " ;
  :SiteType-2 = "SA" ;
  :SiteLon-2 = 117.72f ;
  :SiteLat-2 = 39.040001f ;
  :SiteHgt-2 = 70.f ;
  :ObsRange-2 = 300s ;
  :SiteName-3 = "ShiJiaZhuang" ;
  :SiteID-3 = " " ;
  :SiteType-3 = "SA" ;
  :SiteLon-3 = 114.71f ;
  :SiteLat-3 = 38.349998f ;
  :SiteHgt-3 = 134.f ;
  :ObsRange-3 = 300s ;
  :SiteName-4 = "ZhangBei" ;
  :SiteID-4 = " " ;
  :SiteType-4 = "CB" ;
  :SiteLon-4 = 114.69f ;
  :SiteLat-4 = 41.150002f ;
  :SiteHgt-4 = 1425.f ;
  :ObsRange-4 = 300s ;
  :SiteName-5 = "QinHuangDao" ;
  :SiteID-5 = " " ;
  :SiteType-5 = "SA" ;
  :SiteLon-5 = 118.88f ;
  :SiteLat-5 = 39.880001f



用sdfopen可以将该数据打开,有如下信息:
sdfopen.png
再用q ctlinfo可以有如下信息
qctlinfo1.png qctlinfo2.png qctlinfo3.png

根据以上信息写ctl如下 aa.ctl:
dset 1606290700.nc
title aa
undef -9.99e+33
xdef 1019 linear  111.69   0.01
ydef 881 linear   35.349   0.01
zdef 20 linear    500     500
tdef   1 linear 11:10Z01JUL2019      10MN      
VARS   1
reflectivity 20 -1,40,1   Reflectivity
ENDVARS


画图gs如下:
'reinit'
'open aa.ctl'
'set z 8'
'set mpdset cnworld'
'set gxout shaded'
'd reflectivity/2-33'
'cbarn'
'reinit'


出图如下: tu.png

此图看似没错  实际上回波大概向左偏移了3个经度,
此时正确的单部s波段雷达回波如下:
s.png

为了检验   将所读数据画z=1 的图  用等值线表示  出图如下:
z=1.png

明显在下面有错误数据出现,因此说明数据读取的有问题,整了半天也没整明白到底是哪里的问题
参照了版上用grads画卫星中心TBB数据的那些帖子 也没找出自己的这个问题在哪,请大家帮忙指正,非常感谢!
将数据文件也上传上来,供参考  谢谢!




sdfopen.png

1606290700.nc

17.13 MB, 下载次数: 17

最佳答案

查看完整内容

我修改了你的描述文件,增加了fileheader项,基本能画出你要的效果,和原图在纬度上稍微差一点点。fileheader项的数值是我试出来的,目前我也说不清楚道理,你再研究研究吧。 dset G:/test/1/1606290700.nc *dtype netcdf title radar undef 9.96921e+36 unpack scale_factor add_offset 0.5f -33s fileheader 720 *options yrev xdef 1019 linear 111.69 0.0099999998 ydef 881 linear 35.348999 0.0099999998 z ...
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2021-4-1 09:26:45 | 显示全部楼层
lxe3399 发表于 2021-4-14 14:33
再次感谢您的尝试!
是的编码值是0-255,读出来这个编码值之后,再除以2减去33即为正确的雷达回波dbz值 ...

我修改了你的描述文件,增加了fileheader项,基本能画出你要的效果,和原图在纬度上稍微差一点点。fileheader项的数值是我试出来的,目前我也说不清楚道理,你再研究研究吧。

dset G:/test/1/1606290700.nc
*dtype netcdf
title radar
undef  9.96921e+36
unpack scale_factor add_offset 0.5f -33s
fileheader 720
*options yrev
xdef 1019 linear 111.69 0.0099999998
ydef 881 linear 35.348999 0.0099999998
zdef 20 levels 500 1000 1500 2000 2500 3000 3500 4000
4500 5000 5500 6000 7000 8000 9000 10000 12000 14000
17000 19000
tdef 1 linear 00Z01JAN0001 1mn
vars 1
reflectivity 20 -1,40,1 Reflectivity
endvars


出来的图:

radar3

radar3


出图的脚本:

'reinit'
'open g:/test/1/1606290700.ctl'
i=1
while(i<=20)
'set grads off'
'set mpdset cnhimap'
'set z 'i''
'define dBZ= -33+0.5*reflectivity'
'set gxout shade2b'
'g:/test/1/15colors.gs'
'set clevs   0  5 10 15 20 25 30 35 40 45 50 55 60 65 70'
'set ccols  -1 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34'
'd dBZ'
'cbar'
'gxprint g:/test/1/radar/radar'i'.png white'
'c'
i=i+1
endwhile
;



密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2021-4-1 12:19:48 | 显示全部楼层
你去看一下http://bbs.06climate.com/forum.p ... p;extra=&page=1
先从描述文件入手
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-4-1 13:52:51 | 显示全部楼层
river 发表于 2021-4-1 12:19
你去看一下http://bbs.06climate.com/forum.php?mod=viewthread&tid=18358&extra=&page=1
先从描述文件入 ...

谢谢大神!
看了一下,感觉这个数据特殊的地方在于nc里数据类型:data_type = unsigned byte,
因此ctl里
VARS   1
reflectivity 20 -1,40,1   Reflectivity
ENDVARS
这部分写法直接影响了数据的读取,
具体该什么改
还请大神明示!
非常感谢!
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2021-4-1 14:46:50 | 显示全部楼层
感觉数据本身的描述文件用的是levels,你改用linear,是不是有影响。另外-1,40,1你是怎么来的,不过这里的东西一般倒也没有影响。
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2021-4-2 18:23:55 | 显示全部楼层
lxe3399 发表于 2021-4-1 13:52
谢谢大神!
看了一下,感觉这个数据特殊的地方在于nc里数据类型:data_type = unsigned byte,
因此ctl ...

咱们一步一步来,前面没有问题吗?
unpack scale_factor add_offset 1 32066  你看到这一项了吗?先把这个研究一下看看吧
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-4-7 08:44:45 | 显示全部楼层
river 发表于 2021-4-2 18:23
咱们一步一步来,前面没有问题吗?
unpack scale_factor add_offset 1 32066  你看到这一项了吗?先把这 ...

谢谢大神回复,我看到了这个信息了,也仿照写了如下的ctl:
dset 1606290700.nc
title aa
dtype netcdf
undef -9.99e+33
unpack scale_factor add_offset 0.5 -33
xdef 1019 linear  111.69   0.01
ydef 881 linear   35.349   0.01
zdef 20  linear    500     500
tdef   1 linear 11:10Z01JUL2019      10MN      
VARS   1
reflectivity 20 z,y,x  Reflectivity
ENDVARS
能画出图,但值依然是原始的值(非回波编码值,正确的编码值应为0-255)
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-4-7 08:54:04 | 显示全部楼层
苏拉苏拉 发表于 2021-4-1 14:46
感觉数据本身的描述文件用的是levels,你改用linear,是不是有影响。另外-1,40,1你是怎么来的,不过这里的 ...

谢谢回复!
我把linear改成了levels,依然是同样的问题,因此应该不是这个问题
-1,40,1是依据头文件“reflectivity:data_type = "unsigned byte" ;”这个来的,表示数据是整形,1-byte unsigned chars(0-255)
这个数据里这个还真是很重要
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2021-4-12 00:11:42 | 显示全部楼层
本帖最后由 river 于 2021-4-12 00:15 编辑
lxe3399 发表于 2021-4-7 08:54
谢谢回复!
我把linear改成了levels,依然是同样的问题,因此应该不是这个问题
-1,40,1是依据头文件 ...

dset G:/test/1/1606290700.nc
title radar
*dtype netcdf
undef  9.96921e+36
unpack scale_factor add_offset 0.5 -33
xdef 1019 linear 111.69 0.0099999998
ydef 881 linear 35.348999 0.0099999998
zdef 20 levels 500 1000 1500 2000 2500 3000 3500 4000
4500 5000 5500 6000 7000 8000 9000 10000 12000 14000
17000 19000
tdef 1 linear 00Z01JAN0001 1mn
vars 1
reflectivity 20 -1,40,1 Reflectivity
endvars

这样写数值是正常了,但是经度上的偏移还没有解决。它也不是说偏移,把图上最右边的补到左边就一样了,不知道为啥莫名其妙从那个地方被分成两半了

radar

radar


密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-4-12 16:48:21 | 显示全部楼层
river 发表于 2021-4-12 00:11
dset G:/test/1/1606290700.nc
title radar
*dtype netcdf

谢谢尝试!非常感谢
的确很困惑
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

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

本版积分规则

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

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

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