用fwrite 读nc数据转为dat，读出来再反回去画图发现两张图不一样，是什么原因，请...

心情小雨

我想读取nc的月平均资料，再插值到站点上，读出十进制数据，步骤是这样的
首先，我用gs文件读出年大气可降水量资料，gs是这样的
'reinit'
'set grid off'
'set grads off'
'sdfopen I:\NCEPmonthly\yin\shum.mon.mean.nc'
'sdfopen I:\NCEPmonthly\yin\pres.mon.mean.nc'
'sdfopen I:\NCEPmonthly\yin\uwnd.mon.mean.nc'

*j=1
*while(j<=12)
*'set t 'j''
*'set z 1 8'
'define w1=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=181,t=790,12)'
'define w2=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=182,t=790,12)'
'define w3=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=183,t=790,12)'
'define w4=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=184,t=790,12)'
'define w5=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=185,t=790,12)'
'define w6=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=186,t=790,12)'
'define w7=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=187,t=790,12)'
'define w8=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=188,t=790,12)'
'define w9=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=189,t=790,12)'
'define w10=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=190,t=790,12)'
'define w11=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=191,t=790,12)'
'define w12=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=192,t=790,12)'
*'set z 1'
*'define w6=ave(w,t=192,t=790,12)'
*'define w7=ave(w,t=188,t=790,12)'
*'define w8=ave(w,t=189,t=790,12)'
*'define w9=ave(w,t=190,t=790,12)'
'define sw=(w1+w2+w12)'
*'define sw=(w4+w5+w3)'
*'define sw=(w8+w6+w7)'
*'define sw=(w1+w2+w3+w4+w5+w6+w7+w8+w9+w10+w11+w12)'
*'define sw=(w11+w9+w10)'
*'define ww=aave(sw,lon=19,lon=32,lat=96,lat=108)'
'set lat 18 32'
'set lon 94 110'
'set parea off'
'set csmooth on'
*'set mpdset cnworld'
'set mpdset yunnan'
*'set mpdset hires'
'set map 1 1 6'
'set xlopts 1 100 0.23'
'set ylopts 1 100 0.23'
'set gxout shaded'
*'set gxout contour'
*'set ccolor 1'
'set rbcols auto'
'set clopts 1 5 0.1'
'set cthick 4'
'set cstyle 1'
'set clab forced'
'set clskip 1'
*'set cint 500'
'd sw'
'C:\OpenGrADS\Contents\Resources\Scripts\cbarn.gs'
'cbarn 1 1'
'printim I:\NCEPmonthly\yin\pri-win1.gif gif x1000 y800 white'
*'enable print  I:\NCEPmonthly\yin\pri-year.gmf'
'disable print'
'c'
*j=j+1
*endwhile
*'clear'
*'reinit'

这个时候可以出图，图形如下：
然后，我同样用这个gs文件读出dat数据，程序是这样的：
'reinit'
'set grid off'
'set grads off'
'set gxout fwrite'
'set fwrite I:\NCEPmonthly\yin\tpri-year.dat'
'sdfopen I:\NCEPmonthly\yin\shum.mon.mean.nc'
'sdfopen I:\NCEPmonthly\yin\pres.mon.mean.nc'
'sdfopen I:\NCEPmonthly\yin\uwnd.mon.mean.nc'

*j=1
*while(j<=12)
*'set t 'j''
*'set z 1 8'
'define w1=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=181,t=790,12)'
'define w2=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=182,t=790,12)'
'define w3=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=183,t=790,12)'
'define w4=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=184,t=790,12)'
'define w5=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=185,t=790,12)'
'define w6=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=186,t=790,12)'
'define w7=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=187,t=790,12)'
'define w8=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=188,t=790,12)'
'define w9=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=189,t=790,12)'
'define w10=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=190,t=790,12)'
'define w11=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=191,t=790,12)'
'define w12=ave(vint(pres.2(z=1),shum.1,300)/33.3,t=192,t=790,12)'
*'set z 1'
*'define w6=ave(w,t=192,t=790,12)'
*'define w7=ave(w,t=188,t=790,12)'
*'define w8=ave(w,t=189,t=790,12)'
*'define w9=ave(w,t=190,t=790,12)'
*'define sw=(w1+w2+w12)'
*'define sw=(w4+w5+w3)'
*'define sw=(w8+w6+w7)'
'define sw=(w1+w2+w3+w4+w5+w6+w7+w8+w9+w10+w11+w12)'
*'define sw=(w11+w9+w10)'
*'define ww=aave(sw,lon=19,lon=32,lat=96,lat=108)'
'set lat -90 90'
'set lon 0 360'
*'set parea off'
*'set csmooth on'
*'set mpdset cnworld'
*'set mpdset yunnan'
*'set mpdset hires'
*'set map 1 1 6'
*'set xlopts 1 100 0.23'
*'set ylopts 1 100 0.23'
*'set gxout shaded'
*'set gxout contour'
*'set ccolor 1'
*'set rbcols auto'
*'set clopts 1 5 0.1'
*'set cthick 4'
*'set cstyle 1'
*'set clab forced'
*'set clskip 1'
*'set cint 500'
'd sw'
*'C:\OpenGrADS\Contents\Resources\Scripts\cbarn.gs'
*'cbarn 1 1'
*'printim I:\NCEPmonthly\yin\tpri-year.gif gif x1000 y800 white'
*'enable print  I:\NCEPmonthly\yin\tpri-year.gmf'
*'disable print'
'disable fwrite'
'c'
*j=j+1
*endwhile
'clear'
*'reinit'

为了检查读出的数据的正确性，我写了ctl文件和gs文件如下
：ctl文件
DSET   I:\NCEPmonthly\yin\tpri-year.dat
Title  925 bottom pressure
undef  -32767
xdef  144 linear 0  2.5
ydef  73 linear -90 2.5
zdef  1 levels 1000
tdef 1 linear 00Z01JAN1948 1dy
vars 1
pri   12  -999 pri
endvars

再用gs文件画图比较：
'reinit'
'set grid off'
'set grads off'
*'set gxout fwrite'
*'set fwrite I:\NCEPmonthly\yin\tpri-year.dat'
'open I:\NCEPmonthly\yin\tpri-year.ctl'
'set lat 18 32'
'set lon 94 110'
'set parea off'
'set csmooth on'
*'set mpdset cnworld'
'set mpdset yunnan'
*'set mpdset hires'
'set map 1 1 6'
'set xlopts 1 100 0.23'
'set ylopts 1 100 0.23'
'set gxout shaded'
*'set gxout contour'
*'set ccolor 1'
'set rbcols auto'
'set clopts 1 5 0.1'
'set cthick 4'
'set cstyle 1'
'set clab forced'
'set clskip 1'
*'set cint 500'
'd pri'
'C:\OpenGrADS\Contents\Resources\Scripts\cbarn.gs'
'cbarn 1 1'
'printim I:\NCEPmonthly\yin\tpri-year.gif gif x1000 y800 white'
'enable print  I:\NCEPmonthly\yin\tpri-year.gmf'
'disable print'
*'disable fwrite'
'c'
*j=j+1
*endwhile
'clear'
*'reinit'

结果出来的图是这样的！
检查了好长时间，不知道哪里出问题了，请各位大侠帮忙指导，不胜感谢！

一个不小心

程序太长了，不想看。你可以检查一下是不是读或写的时候把经纬度弄反了，或者计算过程中公式运用错误了

心情小雨

一个不小心 发表于 2014-12-2 07:51
程序太长了，不想看。你可以检查一下是不是读或写的时候把经纬度弄反了，或者计算过程中公式运用错误了

已经解决了，谢谢

caolaomiao

请问是怎样解决的呢？求助！！