grads利用格点降水和温度数据计算干旱指数

幽居寒舍

大家好！
我想利用格点降水和温度数据计算上年10月到当年5月的平均干旱指数,时间从哪1950年10月开始（t=598），温度和降水数据的经纬度和空间分辨率是完全一致的。
公式为mt=p/（tm+10）；p为降水，tm为月均温度。
代码如下（'d mt'下面的代码可以忽略不计）
'reinit'
'set grads off'
'set grid off'
'sdfopen e:\grads\pre.nc'
'sdfopen e:\grads\temp.nc'
'set lat 21 30'
'set lon 97 112'
'define pre=(ave(pre,t=598,t=1380,12)+ave(pre,t=599,t=1380,12))+ave(pre,t=600,t=1380,12)+ave(pre,t=601,t=1386,12)+ave(pre,t=602,t=1386,12)+ave(pre,t=603,t=1386,12)+ave(pre,t=604,t=1386,12)+ave(pre,t=605,t=1386,12))'
'define tmp=(ave(tmp,t=598,t=1380,12)+ave(tmp,t=599,t=1381,12))+ave(tmp,t=600,t=1382,12)+ave(tmp,t=601,t=1383,12)+
ave(tmp,t=602,t=1384,12)+ave(tmp,t=603,t=1385,12)+ave(tmp,t=604,t=1386,12)+ave(tmp,t=605,t=1387,12))/8'
'define mt=pre/(tmp+10)'
'set clevs -2 -1.5 -1 -0.5 0 0.5 1 1.5 2'
'e:\grads\opengrads\output2.gs'
'set gxout shaded'
'set gxout line'
'set mpdset country1'
'set xlpos -4'
'set ylpos -4'
'set parea 0 11 1 7.8'
'd mt'
'run axis.gs -type L  -lsize 0.15 -lfont 5 -position o -interval 3 -start 23 -end 29'
'run axis.gs -type T  -lsize 0.15 -lfont 5 -position o -interval 5 -start 95 -end 115'
'cbarn 1'
'gxprint e:\grads\picture\mt.png png white'
;
但是显示出错，如下图
求大神帮忙看看问题出在哪儿，感激不尽！下图分别是温度数据和降水数据的数据信息




经过自己摸索已经弄明白了，问题出在没有set dfile
下面是修改后的代码，仅供参考：
'reinit'
'sdfopen e:\grads\asianpre.nc'
'sdfopen e:\grads\asiantempmonmean.nc'
'set dfile 1'
'set lat 21.1 30'
'set lon 97 112'
'define p=(ave(pre,t=598,t=1379,12)+ave(pre,t=599,t=1380,12))+ave(pre,t=600,t=1381,12)+ave(pre,t=601,t=1382,12)+ave(pre,t=602,t=1383,12)+ave(pre,t=603,t=1384,12)+ave(pre,t=604,t=1385,12)+ave(pre,t=605,t=1386,12))/8'
'set dfile 2'
'set lat 21 30'
'set lon 97 112'
'define t=(ave(tmp,t=598,t=1380,12)+ave(tmp,t=599,t=1381,12))+ave(tmp,t=600,t=1382,12)+ave(tmp,t=601,t=1383,12)+ave(tmp,t=602,t=1384,12)+ave(tmp,t=603,t=1385,12)+ave(tmp,t=604,t=1386,12)+ave(tmp,t=605,t=1387,12))/8'
'set mpdset country1 yunn guiz guangx'
'set gxout shaded'
'set xlpos -4'
'set ylpos -4'
'set parea 0 11 1 7.8'
'set clevs 1 1.5 2 2.5 3 3.5 4 4.5 5'
'e:\grads\opengrads\output2.gs'
'd p/(t+10)'
'set gxout contour'
'set clab off'
'set cthick 3'
'd p/(t+10)'
'run axis.gs -type L -lsize 0.15 -lfont 5 -position o -interval 3 -start 23 -end 29'
'run axis.gs -type T -lsize 0.15 -lfont 5 -position o -interval 5 -start 95 -end 115'
'cbarn 1'
'gxprint e:\grads\picture\drypt.png png white'
;

gwk

它不是说了吗，tmp不是变量名，检查一下nc文件里面变量名是啥吧

幽居寒舍

gwk 发表于 2018-9-13 21:08
它不是说了吗，tmp不是变量名，检查一下nc文件里面变量名是啥吧

我检查过了，变量名确实是tmp，所以有点懵，不懂为什么出错

gwk

幽居寒舍 发表于 2018-9-13 21:11
我检查过了，变量名确实是tmp，所以有点懵，不懂为什么出错

define pre=(ave........
'define tmp=(ave........
把定义的变量名换一下再试一试吧（我也没看出什么毛病），定义的变量名和原变量名一样，看着头大。

幽居寒舍

gwk 发表于 2018-9-13 21:20
define pre=(ave........
'define tmp=(ave........
把定义的变量名换一下再试一试吧（我也没看出什么 ...

改了，还是一样的问题

小其其格

'sdfopen e:\grads\pre.nc'
'sdfopen e:\grads\temp.nc'

你同时打开了2个nc文件
是不是tmp要写成tmp.2才行？

幽居寒舍

小其其格 发表于 2018-9-20 18:17
'sdfopen e:\grads\pre.nc'
'sdfopen e:\grads\temp.nc'

'reinit'
'sdfopen e:\grads\pre.nc'
'sdfopen e:\grads\temp.nc'
'set dfile 1'
'set lat 21.1 30'
'set lon 97 112'
'define p=(ave(pre,t=598,t=1379,12)+ave(pre,t=599,t=1380,12))+ave(pre,t=600,t=1381,12)+ave(pre,t=601,t=1382,12)+ave(pre,t=602,t=1383,12)+ave(pre,t=603,t=1384,12)+ave(pre,t=604,t=1385,12)+ave(pre,t=605,t=1386,12))/8'
'set dfile 2'
'set lat 21 30'
'set lon 97 112'
'define t=(ave(tmp,t=598,t=1380,12)+ave(tmp,t=599,t=1381,12))+ave(tmp,t=600,t=1382,12)+ave(tmp,t=601,t=1383,12)+ave(tmp,t=602,t=1384,12)+ave(tmp,t=603,t=1385,12)+ave(tmp,t=604,t=1386,12)+ave(tmp,t=605,t=1387,12))/8'
'set mpdset country1 yunn guiz guangx'
'set gxout shaded'
'set xlpos -4'
'set ylpos -4'
'set parea 0 11 1 7.8'
'set clevs 1 1.5 2 2.5 3 3.5 4 4.5 5'
'e:\grads\opengrads\output2.gs'
'd p/(t+10)'
'set gxout contour'
'set clab off'
'set cthick 3'
'd p/(t+10)'
'run axis.gs -type L -lsize 0.15 -lfont 5 -position o -interval 3 -start 23 -end 29'
'run axis.gs -type T -lsize 0.15 -lfont 5 -position o -interval 5 -start 95 -end 115'
'cbarn 1'
'e:\grads\samplingsites.gs'
'gxprint e:\grads\picture\drypt.png png white'

幽居寒舍

小其其格 发表于 2018-9-20 18:17
'sdfopen e:\grads\pre.nc'
'sdfopen e:\grads\temp.nc'

'reinit'
'sdfopen e:\grads\pre.nc'
'sdfopen e:\grads\temp.nc'
'set dfile 1'
'set lat 21.1 30'
'set lon 97 112'
'define p=(ave(pre,t=598,t=1379,12)+ave(pre,t=599,t=1380,12))+ave(pre,t=600,t=1381,12)+ave(pre,t=601,t=1382,12)+ave(pre,t=602,t=1383,12)+ave(pre,t=603,t=1384,12)+ave(pre,t=604,t=1385,12)+ave(pre,t=605,t=1386,12))/8'
'set dfile 2'
'set lat 21 30'
'set lon 97 112'
'define t=(ave(tmp,t=598,t=1380,12)+ave(tmp,t=599,t=1381,12))+ave(tmp,t=600,t=1382,12)+ave(tmp,t=601,t=1383,12)+ave(tmp,t=602,t=1384,12)+ave(tmp,t=603,t=1385,12)+ave(tmp,t=604,t=1386,12)+ave(tmp,t=605,t=1387,12))/8'
'set mpdset country1 yunn guiz guangx'
'set gxout shaded'
'set xlpos -4'
'set ylpos -4'
'set parea 0 11 1 7.8'
'set clevs 1 1.5 2 2.5 3 3.5 4 4.5 5'
'e:\grads\opengrads\output2.gs'
'd p/(t+10)'
'set gxout contour'
'set clab off'
'set cthick 3'
'd p/(t+10)'
'run axis.gs -type L -lsize 0.15 -lfont 5 -position o -interval 3 -start 23 -end 29'
'run axis.gs -type T -lsize 0.15 -lfont 5 -position o -interval 5 -start 95 -end 115'
'cbarn 1'
'e:\grads\samplingsites.gs'
'gxprint e:\grads\picture\drypt.png png white'

幽居寒舍

小其其格 发表于 2018-9-20 18:17
'sdfopen e:\grads\pre.nc'
'sdfopen e:\grads\temp.nc'

问题出在图层设置上，需要设置set dfile就好了