［已解决］画的出梯度图画不出1/梯度的图？

echo要读书

RT，此为我计算中的一部分，因为1/GA画不出来导致后面的步骤也画不出。
求助。。求解决方法。。

截取计算梯度的部分：


'reinit'

'open z:\prog\grd\fro\tempmean.ctl'
'open z:\prog\grd\fro\rhummean.ctl'
'open z:\prog\grd\fro\uwndmean.ctl'
'open z:\prog\grd\fro\vwndmean.ctl'
'open z:\prog\grd\fro\850uwndmean.ctl'
'open z:\prog\grd\fro\850vwndmean.ctl'

'set display color white'
'c'

'set xlopts 1 8 0.25'
'set ylopts 1 8 0.25'
'set clopts 1 5 0.15'

'set mpdset cnworld'
'set xlevs 90 95 100 105 110 115 120 125 130'
'set ylevs 15 20 25 30 35 40 45'

'set lon 80 130'
'set lat 10 50'

'define m=0.0'
'define n=0.0'
'define uu=0.0'
'define vv=0.0'

i=1
while(i<=46)

'set t 'i''
'define tem=air.1'
'define rh=rhum.2'
'define wu=uwnd.3'
'define wv=vwnd.4'
'define ku=uwnd.5'
'define kv=vwnd.6'
****mixing ratio*****

'define es=(6.112*exp(17.67*(tem-273.15)/(tem-36)))'
'define pot=rh*(0.62197*es/(1000-0.378*es))/100.'
****gradient***

'define r=6.37e6'
'define dtx = cdiff(pot,x)'
'define dty = cdiff(pot,y)'
'define dx = cdiff(lon,x)*r*cos(lat*3.1416/180)*3.1416/180'
'define dy = cdiff(lat,y)*r*3.1416/180'
'define ga=mag(dtx/dx,dty/dy)''m=m+ga'
i＝i+1
endwhile

'set gxout shaded'
'set grid off'
'set grads off'
*'d m'
'd 1/m'
'run cbarn'

；

echo要读书

但是我把1/GA展开写成CA=(dx*dy)/(mag(dtx*dy,dty*dx))有的时候能画出来部分。。WHY。。。

虫儿飞

乘以一个系数，将值控制在1左右试试

echo要读书

虫儿飞 发表于 2014-6-26 21:46
乘以一个系数，将值控制在1左右试试

我在最后画图的那步乘了，还是会缺值的。

虫儿飞

echo要读书 发表于 2014-6-26 22:34
我在最后画图的那步乘了，还是会缺值的。

你在计算梯度那里就乘以一个系数试试，这类问题我没有遇到过，我的猜测是倒数值和系统默认的缺测值一致，但是得要多试试才行

瞌睡宝宝

楼主，您好。我把你的程序下载了调试了下，数据用的你提供网站的测试数据。循环测试可以看出是第4天和第46天存在这一问题，取了第四天发现发现你图中空白地区的原始计算值mag(dtx*dy,dty*dx)值为0，取倒数后grads估计是自身问题，没去研究，导致了空白的产生，我做了小处理 mag(dtx*dy,dty*dx)+0.0000001，可以正常出图，而且对你原始图没有偏差，带入46天循环后也基本正常。Ps:我加的小数位数已经算是最大了，再添加精度会继续出现空白缺值区。附上结果


调试程序代码： a.gs (2.38 KB, 下载次数: 29)

2014-6-27 21:00 上传

点击文件名下载附件

这样应该可以解决你的问题了吧

冬日的阳光

想basemap？

斌行乾坤

.............