画相关系数分布

天道酬勤

最近做了海温跟某区域降水序列的相关分析，用surfer画出的图，效果还算满意，不过用grads画出的图，就非常不满意了，这是出来的图，满屏幕的扑 克 牌！！！     请指点一下，是哪里出了问题，我把ctl也贴出来了。。谢谢大家！！！grd文件就是每个格点的序列跟研究区域降水序列的相关系数。

1. dset e:\sst12.grd
   
2. TITLE sst12
   
3. UNDEF 32767
   
4. XDEF 181 LINEAR 0 2
   
5. YDEF 89 LINEAR -88 2
   
6. ZDEF 1 LEVELS 0 0
   
7. TDEF 1 linear 00z01jan2010 1mo
   
8. vars 1
   
9. sst12 1 -999
   
10. endvars
    
11. 
    

dset e:\sst12.grd<br /> TITLE sst12<br /> UNDEF 32767<br /> XDEF 181 LINEAR 0 2<br /> YDEF 89 LINEAR -88 2<br /> ZDEF 1 LEVELS 0 0<br /> TDEF 1 linear 00z01jan2010 1mo<br /> vars 1<br /> sst12 1 -999 <br /> endvars<br /> <br />

做个霸气的木头

天道酬勤 发表于 2012-6-17 13:27

                               
登录/注册后可看大图

鼓捣了一上午，终于出来了。原来是最后读入相关系数时出错了。
open (3,file='e:\sst12xg.grd')
do iz= ...

grads读取的数据需要是二进制的，比如grd，open选项里必须有form=‘binary’，读写这种数据不需要格式说明如：read(100)，write(100)，这种数据是用记事本看不到
txt的数据为十进制数，可以看到，如果要grads读取需要转换成grd数据才行。

天道酬勤

这是我的程序，实在找不出哪里出了毛病
parameter(nx=180,ny=89,nz=1,n=35,undef=-9.96921e+36)
dimension ty(n),tw(nx,ny,nz,n)
real ya,wa(nx,ny,nz),syw(nx,ny,nz),sy,sw(nx,ny,nz),r(nx,ny,nz),t(nx,ny,nz)
open(1,file='e:\1.txt')
read(1,*) (ty(it),it=1,n)
close(1)

open(2,file='e:\sst12.grd',form='binary')
read(2) ((((tw(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz),it=1,n)
close(2)

do iz=1,nz
do i=1,nx
do j=1,ny
ya=0.0
wa(i,j,iz)=0.0
do it=1,n
ya=ya+ty(it)
wa(i,j,iz)=wa(i,j,iz)+tw(i,j,iz,it)
end do
ya=ya/real(n)
wa(i,j,iz)=wa(i,j,iz)/real(n)

syw(i,j,iz)=0.0
sy=0.0
sw(i,j,iz)=0.0
do it=1,n
syw(i,j,iz)=syw(i,j,iz)+(ty(it)-ya)*(tw(i,j,iz,it)-wa(i,j,iz))
sy=sy+(ty(it)-ya)*(ty(it)-ya)
sw(i,j,iz)=sw(i,j,iz)+(tw(i,j,iz,it)-wa(i,j,iz))*(tw(i,j,iz,it)-wa(i,j,iz))
end do
syw(i,j,iz)=syw(i,j,iz)/real(n)
sy=sy/real(n)
sw(i,j,iz)=sw(i,j,iz)/real(n)
r(i,j,iz)=0.0
r(i,j,iz)=syw(i,j,iz)/sqrt(sy*sw(i,j,iz))
enddo;enddo;enddo
do it=1,n
do iz=1,nz
do i=1,nx
do j=1,ny
if(tw(i,j,iz,it)==undef)then
r(i,j,iz)=0.0
endif
enddo;enddo;enddo;enddo
open (3,file='e:\sst12xg.grd')
do iz=1,nz
do i=1,nx
do j=1,ny
write (3,*) r(i,j,iz)
enddo;enddo;enddo;
end
之后我就编了上面的ctl。结果就是这样。nx到底是181还是180呢，缺测值如果换成32767的话，得出的相关系数全部为零。

mofangbao

可能是缺测值设置不对，要不然就是你的数据有问题了，这种图明显不是因为用GrADS画才造成的，是你自己的数据和文件配置问题

做个霸气的木头

应该是缺测值或者精度设置的问题，你好好检查下，是不是跟数据相对应~

做个霸气的木头

缺测值是不能换的，就是-9.96921e+36，你的数据水平方向是多少个格点，就是多少啊，经度为2的画应该是（360/2+1）-1=180才是；
再者你的程序显示，你是提取的 1文件和sst12的数据计算，最后输出到sst12xg.grd中，那你这里画图怎么画的是sst12呢？

天道酬勤

那是发上来的时候没有改过来。我在想是不是毛病出在。我的1文件是站点资料，sst12是格点资料。因为我用做出的高度场做了一遍，出来的相关系数图，也是那个毛病。。。真心不知道问题出在哪里？求答疑，求解惑。另外给楼上一个链接。
http://blog.sina.com.cn/s/blog_68c2b7570100ycdh.html

做个霸气的木头

天道酬勤 发表于 2012-6-16 21:49

                               
登录/注册后可看大图

那是发上来的时候没有改过来。我在想是不是毛病出在。我的1文件是站点资料，sst12是格点资料。因为我用做出 ...

站点、格点什么没关系，你1号文件的数据只有35个，，sst12的较多。不如你这样，你把读取的ty、tw、r三个数组分别输出到屏幕上，自己看看数据读入和计算有木有问题，没有问题的话，再照着你给的链接改一改试试。ps：你敢不敢点个回复呢，我都不知道你回复了。

天道酬勤

做个霸气的木头 发表于 2012-6-16 22:00

                               
登录/注册后可看大图

站点、格点什么没关系，你1号文件的数据只有35个，，sst12的较多。不如你这样，你把读取的ty、tw、r三个数 ...

鼓捣了一上午，终于出来了。原来是最后读入相关系数时出错了。
open (3,file='e:\sst12xg.grd')
do iz=1,nz
do i=1,nx
do j=1,ny
write (3,*) r(i,j,iz)
enddo;enddo;enddo;
end
写入.txt文件和写入.grd文件完全是不同的，具体原因，求解释啊。我是写入.txt之后，再转化为.grd。就可以正常画图了。不过谢谢楼上的热心帮助啊！！！

天道酬勤

刚才试了试，果然是哇，，，学习了。 非常感谢你的帮助！