grads里的cdiff函数直接出图和输出格点再出图为什么差距那么大

自言自语

下面三张图第一张是用过cdiff函数后直接出图，第二张是我自己编的，第三张是先用cdiff函数出格点数据，再打开格点数据出图。第一、二图除了地图边缘细微差别，别的几乎都一样。第三张图跟前两张就差得多
三张图的gs文件分别是：
第一张：
'reinit'
'sdfopen e:\ncep\air.2019.nc'

'set t 1424'
'set z 3'

'd cdiff(air,x)/(cos(lat*3.1416/180)*5*3.1416*6371000/180)'
*****************************************************

第二张：
'reinit'
'sdfopen e:\ncep\air.2019.nc'
'set gxout fwrite'
'set fwrite e:\ncep\air.20191222wx.dat'
'set t 1424'
'set z 3'
i=2
while(i<73)
'set y 'i
'define WJ=cos((2.5*'i'-92.5)*3.1416/180)*5*3.1416*6371000/180'
j=2
while(j<144)
'set x 'j
'define WXTD=(air(x+1)-air(x-1))/WJ'
'd WXTD'
j=j+1
endwhile
i=i+1
endwhile
'disable fwrite'


配的描述文件是：
dset e:\ncep\air.20191222wx.dat
undef -9.96921e+36
xdef 142 linear 2.5 2.5
ydef 71 linear -87.5 2.5
zdef 1 levels 850
tdef 1 linear 18Z22dec2019 1mo
vars 1
WXTD 1 99  wxtd
endvars

******************************************************
第三张：
'reinit'
'sdfopen e:\ncep\air.2019.nc'
'set gxout fwrite'
'set fwrite e:\ncep\air.20191222wx.1.dat'
'set t 1424'
'set z 3'
'define WXTD=cdiff(air,x)/(cos(lat*3.1416/180)*5*3.1416*6371000/180)'
'd WXTD'
'disable fwrite'


配的描述文件：
dset e:\ncep\air.20191222wx.1.dat
undef -9.99e+08
xdef 142 linear 2.5 2.5
ydef 71 linear -87.5 2.5
zdef 1 levels 850
tdef 1 linear 18Z22dec2019 1mo
vars 1
WXTD 1 99  wxtd
endvars

lqouc

自言自语 发表于 2020-12-29 11:41
是不一样，大了2kb。我改成了这样：
dset e:\ncep\air.20191222wx.1.dat
undef -9.99e+8

那在输出数据的时候设置上set x 1 144; set y 1 73

自言自语

第三张图用set gxout stat后的提示是“undef value=-9.99e+08”，缺省值应该没问题吧

river

自言自语 发表于 2020-12-28 17:39
第三张图用set gxout stat后的提示是“undef value=-9.99e+08”，缺省值应该没问题吧

你用 q undef 命令，用那个缺测值试一下

自言自语

river 发表于 2020-12-28 21:43
你用 q undef 命令，用那个缺测值试一下

提示“output undef value is set to -999000000.000000”,我把描述文件里面的缺省值改成它画图还是那个样子

自言自语

river 发表于 2020-12-28 21:43
你用 q undef 命令，用那个缺测值试一下

我刚拿era5数据测了一下就没这问题，用了cdiff直接画图和出格点值再画图两个是一模一样的

lqouc

你看下air.20191222wx.1.dat和air.20191222wx.dat这两个文件一样大吗？
我感觉第三张图的ctl应该是144*73啊，只不过边界是缺省值。

自言自语

lqouc 发表于 2020-12-29 10:32
你看下air.20191222wx.1.dat和air.20191222wx.dat这两个文件一样大吗？
我感觉第三张图的ctl应该是144*73 ...

是不一样，大了2kb。我改成了这样：
dset e:\ncep\air.20191222wx.1.dat
undef -9.99e+8
xdef 144 linear 0 2.5
ydef 73 linear -90 2.5
zdef 1 levels 850
tdef 1 linear 18Z22dec2019 1mo
vars 1
WXTD 1 99  wxtd
endvars
出来的图还是跟前面的不一样，但比之前的要正常一些。之前那张图相当于以南极为支点弯了60°，现在差不多弯了45°

自言自语

lqouc 发表于 2020-12-29 13:39
那在输出数据的时候设置上set x 1 144; set y 1 73

我把这句加到gs文件里，再用出来的格点数据画图就正常了
'reinit'
'sdfopen e:\ncep\air.2019.nc'
'set gxout fwrite'
'set fwrite e:\ncep\air.20191222wx.1.dat'
'set x 1 144'
'set y 1 73'
'set t 1424'
'set z 3'
'define WXTD=cdiff(air,x)/(cos(lat*3.1416/180)*5*3.1416*6371000/180)'
'd WXTD'
'disable fwrite'