grads2ascii提取固定经纬度数值不对

边缘人糊涂乐

用grads2ascii提取固定经纬度数值不对：
原本的gs脚本画出的水汽通量散度图如下：附脚本：
'reinit'
'open d:\xkt\FNL\20190814.ctl'

'set mpdset cnworld'
*'set mpdset liaon'
'set map 15 1 4'
'set grads off'
'set grid off'

i=1
while(i<=12)
'set t 'i''

'set lat 0 90 '
'set lon 0 180'
'set lev 850'
'set xlint 5'
'set ylint 5'

'define p1=lev'
'define e=(6.112*exp((17.67*(TMPprs-273.16))/(TMPprs-29.65)))'
'define qs=(0.622*e/(p1-0.378*e))'
'define q3=qs*RHprs/100'
'define qdiv=hdivg(q3*UGRDprs/9.8,q3*VGRDprs/9.8)'

'set gxout shaded'
'set lat 30 55'
'set lon 110 140'
'set xlopts 1 12 0.18'
'set ylopts 1 12 0.18'
'set grads off'
'set grid off'
*'set xlint 2'
*'set ylint 1'
'set cint 10'
'set cmax -10'
'd smth9(qdiv*1e9)'
'cbarn 0.8 1  9.7 4.3'

'set gxout contour'
'set lat 30 55'
'set lon 110 140'
'set xlopts 1 12 0.18'
'set ylopts 1 12 0.18'
'set grads off'
'set grid off'
*'set xlint 2'
*'set ylint 1'
'set cint 10'
'set cmax -10'
'd smth9(qdiv*1e9)'

'define p=lev'
'define e=exp(17.67*(TMPprs-273.15)/(TMPprs-29.65))*(RHprs/100)*6.112'
'define q=622*e/(p-0.378*e)'
'define q1=UGRDprs*q/9.8'
'define q2=VGRDprs*q/9.8'
'set gxout vector'
'set ccolor 3'
*'set cint 0.5'
'set lat 30 55'
'set lon 110 140'
'set ccolor 1'
'set xlopts 1 12 0.18'
'set ylopts 1 12 0.18'
'set grads off'
'set grid off'
*'set xlint 2'
*'set ylint 1'
'd q1;q2'
*'d mag(q1,q2)'

'q time'
rst=substr(result,8,12)
*'set strsiz 0.1 0.1'
*'draw title  'rst''
'printim D:\xkt\shuiqi\'rst'.png  png white'
'c'
i=i+1
endwhile
'print'
'disable print'
*'print'
;

水汽通量散度图

图1 2019年08月14日00时850hPa水汽通量散度图


我想提取经纬度为42°N，124E这一点的水汽通量散度数值。
编程脚本如下：
'reinit'
'open d:\xkt\FNL\20190814.ctl'

'set lat 0 90'
'set lon 0 180'
'set lev 850'
'set t 1 12'
'define p1=lev'
'define e=(6.112*exp((17.67*(TMPprs-273.16))/(TMPprs-29.65)))'
'define qs=(0.622*e/(p1-0.378*e))'
'define q3=qs*RHprs/100'
'define qdiv=hdivg(q3*UGRDprs/9.8,q3*VGRDprs/9.8)'

'define p=lev'
'define e=exp(17.67*(TMPprs-273.15)/(TMPprs-29.65))*(RHprs/100)*6.112'
'define q=622*e/(p-0.378*e)'
'define q1=UGRDprs*q/9.8'
'define q2=VGRDprs*q/9.8'
'define qq=mag(q1,q2)'

'set lat 42'
'set lon 124'
'set lev 850'
'set t 1 12'

'grads2ascii qdiv*1e9 D:\xkt\qdiv.txt a %g 1 1'
;

可以运行得出txt文件，发现txt文本里的数值为-64.79。但是从图中看42°N，124E这一点的数值应该在-20左右，为什么数值对不上呢？

这张图的第一个数值应该对应上一张平面图42°N，124E这一点的数值，但是对不上。

CTL如下：

dset D:/xkt/FNL/fnl_201908%d2_%h2_00.grib2
options template
index D:/xkt/FNL/fnl_20190814_00_00.grib2.idx
undef 9.999E+20
title D:/xkt/FNL/fnl_20190814_00_00.grib2
* produced by g2ctl v0.1.1
* command line options: D:/xkt/FNL/fnl_20190814_00_00.grib2
* griddef=1:0:(360 x 181):grid_template=0:winds(N/S): lat-lon grid:(360 x 181) units 1e-06 input WE:NS output WE:SN res 48 lat 90.000000 to -90.000000 by 1.000000 lon 0.000000 to 359.000000 by 1.000000 #points=65160:winds(N/S)

dtype grib2
ydef 181 linear -90.000000 1
xdef 360 linear 0.000000 1.000000
tdef 12 linear 00Z14aug2019 6hr
* PROFILE hPa
zdef 34 levels 100000 97500 95000 92500 90000 85000 80000 75000 70000 65000 60000 55000 50000 45000 40000 35000 30000 25000 20000 15000 10000 7000 5000 4000 3000 2000 1500 1000 700 500 300 200 100 40
options pascals
vars 147
no4LFTXsfc   0,1,0   0,7,193 ** surface Best (4 layer) Lifted Index [K]
no5WAVH500mb   0,100,50000   0,3,193 ** 500 mb 5-Wave Geopotential Height [gpm]
ABSVprs    34,100  0,2,10 ** (1000 975 950 925 900.. 30 20 15 10 40) Absolute Vorticity [1/s]
CAPEsfc   0,1,0   0,7,6 ** surface Convective Available Potential Energy [J/kg]
CAPE180_0mb  0,108,18000,0   0,7,6 ** 180-0 mb above ground Convective Available Potential Energy [J/kg]
CAPE255_0mb  0,108,25500,0   0,7,6 ** 255-0 mb above ground Convective Available Potential Energy [J/kg]
CFRZRsfc   0,1,0   0,1,193 ** surface Categorical Freezing Rain [-]
CICEPsfc   0,1,0   0,1,194 ** surface Categorical Ice Pellets [-]
CINsfc   0,1,0   0,7,7 ** surface Convective Inhibition [J/kg]
CIN180_0mb  0,108,18000,0   0,7,7 ** 180-0 mb above ground Convective Inhibition [J/kg]
CIN255_0mb  0,108,25500,0   0,7,7 ** 255-0 mb above ground Convective Inhibition [J/kg]
CLWMRprs    34,100  0,1,22 ** (1000 975 950 925 900.. 250 200 150 100 50) Cloud Mixing Ratio [kg/kg]
CLWMRhy1   0,105,1   0,1,22 ** 1 hybrid level Cloud Mixing Ratio [kg/kg]
CPOFPsfc   0,1,0   0,1,39 ** surface Percent frozen precipitation [%]
CRAINsfc   0,1,0   0,1,192 ** surface Categorical Rain [-]
CSNOWsfc   0,1,0   0,1,195 ** surface Categorical Snow [-]
CWATclm   0,200,0   0,6,6 ** entire atmosphere (considered as a single layer) Cloud Water [kg/m^2]
DPT2m   0,103,2   0,0,6 ** 2 m above ground Dew Point Temperature [K]
DZDTprs    21,100  0,2,9 ** (1000 975 950 925 900.. 300 250 200 150 100) Vertical Velocity (Geometric) [m/s]
FLDCPsfc   0,1,0   2,3,203 ** surface Field Capacity [Fraction]
GRLEprs    34,100  0,1,32 ** (1000 975 950 925 900.. 250 200 150 100 50) Graupel [kg/kg]
GRLEhy1   0,105,1   0,1,32 ** 1 hybrid level Graupel [kg/kg]
GUSTsfc   0,1,0   0,2,22 ** surface Wind Speed (Gust) [m/s]
HGTsfc   0,1,0   0,3,5 ** surface Geopotential Height [gpm]
HGTprs    34,100  0,3,5 ** (1000 975 950 925 900.. 5 3 2 1 40) Geopotential Height [gpm]
HGT2pv   0,109,2e-06   0,3,5 ** PV=2e-06 (Km^2/kg/s) surface Geopotential Height [gpm]
HGTneg2pv   0,109,-2e-06   0,3,5 ** PV=-2e-06 (Km^2/kg/s) surface Geopotential Height [gpm]
HGTtop0C   0,204,0   0,3,5 ** highest tropospheric freezing level Geopotential Height [gpm]
HGT0C   0,4,0   0,3,5 ** 0C isotherm Geopotential Height [gpm]
HGTmwl   0,6,0   0,3,5 ** max wind Geopotential Height [gpm]
HGTtrop   0,7,0   0,3,5 ** tropopause Geopotential Height [gpm]
HINDEXsfc   0,1,0   2,4,2 ** surface Haines Index [Numeric]
HLCY3000_0m  0,103,3000,0   0,7,8 ** 3000-0 m above ground Storm Relative Helicity [m^2/s^2]
HPBLsfc   0,1,0   0,3,196 ** surface Planetary Boundary Layer Height [m]
ICAHTmwl   0,6,0   0,3,3 ** max wind ICAO Standard Atmosphere Reference Height [m]
ICAHTtrop   0,7,0   0,3,3 ** tropopause ICAO Standard Atmosphere Reference Height [m]
ICECsfc   0,1,0   10,2,0 ** surface Ice Cover [Proportion]
ICMRprs    34,100  0,1,23 ** (1000 975 950 925 900.. 250 200 150 100 50) Ice Water Mixing Ratio [kg/kg]
ICMRhy1   0,105,1   0,1,23 ** 1 hybrid level Ice Water Mixing Ratio [kg/kg]
LANDsfc   0,1,0   2,0,0 ** surface Land Cover (0=sea, 1=land) [Proportion]
LFTXsfc   0,1,0   0,7,192 ** surface Surface Lifted Index [K]
MSLETmsl   0,101,0   0,3,192 ** mean sea level MSLP (Eta model reduction) [Pa]
O3MRprs    34,100  0,14,192 ** (1000 850 700 500 400.. 5 3 2 1 40) Ozone Mixing Ratio [kg/kg]
PLPL255_0mb  0,108,25500,0   0,3,200 ** 255-0 mb above ground Pressure of level from which parcel was lifted [Pa]
POTsig995   0,104,0.995   0,0,2 ** 0.995 sigma level Potential Temperature [K]
PRATEsfc   0,1,0   0,1,7 ** surface Precipitation Rate [kg/m^2/s]
PRESsfc   0,1,0   0,3,0 ** surface Pressure [Pa]
PRES80m   0,103,80   0,3,0 ** 80 m above ground Pressure [Pa]
PRES2pv   0,109,2e-06   0,3,0 ** PV=2e-06 (Km^2/kg/s) surface Pressure [Pa]
PRESneg2pv   0,109,-2e-06   0,3,0 ** PV=-2e-06 (Km^2/kg/s) surface Pressure [Pa]
PRESmwl   0,6,0   0,3,0 ** max wind Pressure [Pa]
PREStrop   0,7,0   0,3,0 ** tropopause Pressure [Pa]
PRMSLmsl   0,101,0   0,3,1 ** mean sea level Pressure Reduced to MSL [Pa]
PWATclm   0,200,0   0,1,3 ** entire atmosphere (considered as a single layer) Precipitable Water [kg/m^2]
REFCclm   0,10,0   0,16,196 ** entire atmosphere Composite reflectivity [dB]
RHprs    34,100  0,1,1 ** (1000 975 950 925 900.. 7 5 3 2 1) Relative Humidity [%]
RH2m   0,103,2   0,1,1 ** 2 m above ground Relative Humidity [%]
RHsg330_1000  0,104,0.33,1   0,1,1 ** 0.33-1 sigma layer Relative Humidity [%]
RHsg440_1000  0,104,0.44,1   0,1,1 ** 0.44-1 sigma layer Relative Humidity [%]
RHsg720_940  0,104,0.72,0.94   0,1,1 ** 0.72-0.94 sigma layer Relative Humidity [%]
RHsg440_720  0,104,0.44,0.72   0,1,1 ** 0.44-0.72 sigma layer Relative Humidity [%]
RHsig995   0,104,0.995   0,1,1 ** 0.995 sigma level Relative Humidity [%]
RH30_0mb  0,108,3000,0   0,1,1 ** 30-0 mb above ground Relative Humidity [%]
RHclm   0,200,0   0,1,1 ** entire atmosphere (considered as a single layer) Relative Humidity [%]
RHtop0C   0,204,0   0,1,1 ** highest tropospheric freezing level Relative Humidity [%]
RH0C   0,4,0   0,1,1 ** 0C isotherm Relative Humidity [%]
RWMRprs    34,100  0,1,24 ** (1000 975 950 925 900.. 250 200 150 100 50) Rain Mixing Ratio [kg/kg]
RWMRhy1   0,105,1   0,1,24 ** 1 hybrid level Rain Mixing Ratio [kg/kg]
SNMRprs    34,100  0,1,25 ** (1000 975 950 925 900.. 250 200 150 100 50) Snow Mixing Ratio [kg/kg]
SNMRhy1   0,105,1   0,1,25 ** 1 hybrid level Snow Mixing Ratio [kg/kg]
SNODsfc   0,1,0   0,1,11 ** surface Snow Depth [m]
SOILW0_10cm  0,106,0,0.1   2,0,192 ** 0-0.1 m below ground Volumetric Soil Moisture Content [Fraction]
SOILW10_40cm  0,106,0.1,0.4   2,0,192 ** 0.1-0.4 m below ground Volumetric Soil Moisture Content [Fraction]
SOILW40_100cm  0,106,0.4,1   2,0,192 ** 0.4-1 m below ground Volumetric Soil Moisture Content [Fraction]
SOILW100_200cm  0,106,1,2   2,0,192 ** 1-2 m below ground Volumetric Soil Moisture Content [Fraction]
SPFH2m   0,103,2   0,1,0 ** 2 m above ground Specific Humidity [kg/kg]
SPFH80m   0,103,80   0,1,0 ** 80 m above ground Specific Humidity [kg/kg]
SPFH30_0mb  0,108,3000,0   0,1,0 ** 30-0 mb above ground Specific Humidity [kg/kg]
SUNSDsfc   0,1,0   0,6,201 ** surface Sunshine Duration
TCDCprs    34,100  0,6,1 ** (1000 975 950 925 900.. 250 200 150 100 50) Total Cloud Cover [%]
TMPsfc   0,1,0   0,0,0 ** surface Temperature [K]
TMPprs    34,100  0,0,0 ** (1000 975 950 925 900.. 5 3 2 1 40) Temperature [K]
TMP_1829m   0,102,1829   0,0,0 ** 1829 m above mean sea level Temperature [K]
TMP_2743m   0,102,2743   0,0,0 ** 2743 m above mean sea level Temperature [K]
TMP_3658m   0,102,3658   0,0,0 ** 3658 m above mean sea level Temperature [K]
TMP2m   0,103,2   0,0,0 ** 2 m above ground Temperature [K]
TMP80m   0,103,80   0,0,0 ** 80 m above ground Temperature [K]
TMP100m   0,103,100   0,0,0 ** 100 m above ground Temperature [K]
TMPsig995   0,104,0.995   0,0,0 ** 0.995 sigma level Temperature [K]
TMP30_0mb  0,108,3000,0   0,0,0 ** 30-0 mb above ground Temperature [K]
TMP2pv   0,109,2e-06   0,0,0 ** PV=2e-06 (Km^2/kg/s) surface Temperature [K]
TMPneg2pv   0,109,-2e-06   0,0,0 ** PV=-2e-06 (Km^2/kg/s) surface Temperature [K]
TMPmwl   0,6,0   0,0,0 ** max wind Temperature [K]
TMPtrop   0,7,0   0,0,0 ** tropopause Temperature [K]
TOZNEclm   0,200,0   0,14,0 ** entire atmosphere (considered as a single layer) Total Ozone [DU]
TSOIL0_10cm  0,106,0,0.1   2,0,2 ** 0-0.1 m below ground Soil Temperature    Validation to deprecate [K]
TSOIL10_40cm  0,106,0.1,0.4   2,0,2 ** 0.1-0.4 m below ground Soil Temperature    Validation to deprecate [K]
TSOIL40_100cm  0,106,0.4,1   2,0,2 ** 0.4-1 m below ground Soil Temperature    Validation to deprecate [K]
TSOIL100_200cm  0,106,1,2   2,0,2 ** 1-2 m below ground Soil Temperature    Validation to deprecate [K]
UGRDprs    34,100  0,2,2 ** (1000 975 950 925 900.. 7 5 3 2 1) U-Component of Wind [m/s]
UGRD_1829m   0,102,1829   0,2,2 ** 1829 m above mean sea level U-Component of Wind [m/s]
UGRD_2743m   0,102,2743   0,2,2 ** 2743 m above mean sea level U-Component of Wind [m/s]
UGRD_3658m   0,102,3658   0,2,2 ** 3658 m above mean sea level U-Component of Wind [m/s]
UGRD10m   0,103,10   0,2,2 ** 10 m above ground U-Component of Wind [m/s]
UGRD20m   0,103,20   0,2,2 ** 20 m above ground U-Component of Wind [m/s]
UGRD30m   0,103,30   0,2,2 ** 30 m above ground U-Component of Wind [m/s]
UGRD40m   0,103,40   0,2,2 ** 40 m above ground U-Component of Wind [m/s]
UGRD50m   0,103,50   0,2,2 ** 50 m above ground U-Component of Wind [m/s]
UGRD80m   0,103,80   0,2,2 ** 80 m above ground U-Component of Wind [m/s]
UGRD100m   0,103,100   0,2,2 ** 100 m above ground U-Component of Wind [m/s]
UGRDsig995   0,104,0.995   0,2,2 ** 0.995 sigma level U-Component of Wind [m/s]
UGRD30_0mb  0,108,3000,0   0,2,2 ** 30-0 mb above ground U-Component of Wind [m/s]
UGRD2pv   0,109,2e-06   0,2,2 ** PV=2e-06 (Km^2/kg/s) surface U-Component of Wind [m/s]
UGRDneg2pv   0,109,-2e-06   0,2,2 ** PV=-2e-06 (Km^2/kg/s) surface U-Component of Wind [m/s]
UGRDpbl   0,220,0   0,2,2 ** planetary boundary layer U-Component of Wind [m/s]
UGRDmwl   0,6,0   0,2,2 ** max wind U-Component of Wind [m/s]
UGRDtrop   0,7,0   0,2,2 ** tropopause U-Component of Wind [m/s]
USTM6000_0m  0,103,6000,0   0,2,194 ** 6000-0 m above ground U-Component Storm Motion [m/s]
VGRDprs    34,100  0,2,3 ** (1000 975 950 925 900.. 7 5 3 2 1) V-Component of Wind [m/s]
VGRD_1829m   0,102,1829   0,2,3 ** 1829 m above mean sea level V-Component of Wind [m/s]
VGRD_2743m   0,102,2743   0,2,3 ** 2743 m above mean sea level V-Component of Wind [m/s]
VGRD_3658m   0,102,3658   0,2,3 ** 3658 m above mean sea level V-Component of Wind [m/s]
VGRD10m   0,103,10   0,2,3 ** 10 m above ground V-Component of Wind [m/s]
VGRD20m   0,103,20   0,2,3 ** 20 m above ground V-Component of Wind [m/s]
VGRD30m   0,103,30   0,2,3 ** 30 m above ground V-Component of Wind [m/s]
VGRD40m   0,103,40   0,2,3 ** 40 m above ground V-Component of Wind [m/s]
VGRD50m   0,103,50   0,2,3 ** 50 m above ground V-Component of Wind [m/s]
VGRD80m   0,103,80   0,2,3 ** 80 m above ground V-Component of Wind [m/s]
VGRD100m   0,103,100   0,2,3 ** 100 m above ground V-Component of Wind [m/s]
VGRDsig995   0,104,0.995   0,2,3 ** 0.995 sigma level V-Component of Wind [m/s]
VGRD30_0mb  0,108,3000,0   0,2,3 ** 30-0 mb above ground V-Component of Wind [m/s]
VGRD2pv   0,109,2e-06   0,2,3 ** PV=2e-06 (Km^2/kg/s) surface V-Component of Wind [m/s]
VGRDneg2pv   0,109,-2e-06   0,2,3 ** PV=-2e-06 (Km^2/kg/s) surface V-Component of Wind [m/s]
VGRDpbl   0,220,0   0,2,3 ** planetary boundary layer V-Component of Wind [m/s]
VGRDmwl   0,6,0   0,2,3 ** max wind V-Component of Wind [m/s]
VGRDtrop   0,7,0   0,2,3 ** tropopause V-Component of Wind [m/s]
VISsfc   0,1,0   0,19,0 ** surface Visibility [m]
VRATEpbl   0,220,0   0,2,224 ** planetary boundary layer Ventilation Rate [m^2/s]
VSTM6000_0m  0,103,6000,0   0,2,195 ** 6000-0 m above ground V-Component Storm Motion [m/s]
VVELprs    21,100  0,2,8 ** (1000 975 950 925 900.. 300 250 200 150 100) Vertical Velocity (Pressure) [Pa/s]
VVELsig995   0,104,0.995   0,2,8 ** 0.995 sigma level Vertical Velocity (Pressure) [Pa/s]
VWSH2pv   0,109,2e-06   0,2,192 ** PV=2e-06 (Km^2/kg/s) surface Vertical Speed Shear [1/s]
VWSHneg2pv   0,109,-2e-06   0,2,192 ** PV=-2e-06 (Km^2/kg/s) surface Vertical Speed Shear [1/s]
VWSHtrop   0,7,0   0,2,192 ** tropopause Vertical Speed Shear [1/s]
WEASDsfc   0,1,0   0,1,13 ** surface Water Equivalent of Accumulated Snow Depth [kg/m^2]
WILTsfc   0,1,0   2,0,201 ** surface Wilting Point [Fraction]
var00212m   0,103,2   0,0,21 ** 2 m above ground desc [unit]
ENDVARS

river

你画的是哪个时次的图，你一共算了12个时次的，会不会是时次没对上？全部是-64.79吗？

边缘人糊涂乐

river 发表于 2020-11-12 19:55
你画的是哪个时次的图，你一共算了12个时次的，会不会是时次没对上？全部是-64.79吗？

是12个时次的，出来是12个数，我对第一个数据和画出来的第一张图进行了对比，不一样，输出一个时次图和第一个数据对比不一样，我画的是850hPa的图。

学业有成

边缘人糊涂乐 发表于 2020-11-13 09:07
是12个时次的，出来是12个数，我对第一个数据和画出来的第一张图进行了对比，不一样，输出一个时次图和第 ...

你好，我最近也学grads编程，你能帮帮我吗？忙着写论文，很着急，看见加我一下qq：455528540

river

边缘人糊涂乐 发表于 2020-11-13 09:07
是12个时次的，出来是12个数，我对第一个数据和画出来的第一张图进行了对比，不一样，输出一个时次图和第 ...

你要不试一下增加写入，就是循环起来，算一个写出来一个看看，用append

边缘人糊涂乐

river 发表于 2020-11-19 12:47
你要不试一下增加写入，就是循环起来，算一个写出来一个看看，用append

这个问题算是变通的解决了，我求了区域平均导出来的数值是对的，我用区域平均代替点了。但是还是不知道为啥点提取出的值不对。

river

边缘人糊涂乐 发表于 2020-12-7 14:05
这个问题算是变通的解决了，我求了区域平均导出来的数值是对的，我用区域平均代替点了。但是还是不知道为 ...

我觉得你可以用set gxout grid，看一下网格里的值是多少。可能是散度插值以后再赋值给格点过程中有误差？而且你两个脚本插值的经纬度范围也不同，是不是跟这些有关系，可以再研究一下

边缘人糊涂乐

river 发表于 2020-12-8 12:05
我觉得你可以用set gxout grid，看一下网格里的值是多少。可能是散度插值以后再赋值给格点过程中有误差？ ...

嗯嗯，我再研究一下，谢谢。