grads处理nc雷达数据，谢谢

lxe3399

有个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可以将该数据打开，有如下信息：

再用q ctlinfo可以有如下信息


根据以上信息写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'


出图如下：

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


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


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

river

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

出图的脚本：

'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
;

river

你去看一下http://bbs.06climate.com/forum.p ... p;extra=&page=1
先从描述文件入手

lxe3399

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
这部分写法直接影响了数据的读取，
具体该什么改
还请大神明示！
非常感谢！

苏拉苏拉

感觉数据本身的描述文件用的是levels，你改用linear，是不是有影响。另外-1,40,1你是怎么来的，不过这里的东西一般倒也没有影响。

river

lxe3399 发表于 2021-4-1 13:52
谢谢大神！
看了一下，感觉这个数据特殊的地方在于nc里数据类型：data_type = unsigned byte，
因此ctl ...

咱们一步一步来，前面没有问题吗？
unpack scale_factor add_offset 1 32066  你看到这一项了吗？先把这个研究一下看看吧

lxe3399

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）

lxe3399

苏拉苏拉 发表于 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）
这个数据里这个还真是很重要

river

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

lxe3399

river 发表于 2021-4-12 00:11
dset G:/test/1/1606290700.nc
title radar
*dtype netcdf

谢谢尝试！非常感谢
的确很困惑