用grads的tcorr画相关系数的空间分布

夜未央

我有一组一维数，共140个，对应140个月
想画它和位势高度场的相关系数二维图
我已经把数改成dat并配置了ctl，grads试过了，可以单独打开它
写了个gs，但是不知道为什么总说我beyond....
以下是一维数的ctl和我写的gs
求各位帮助】
dset D:\srdp\xiangguan\2D\pdsijuping.dat
undef -999.0
title pdsi-hdpingjun
xdef  1 linear 1 1
ydef  1 linear 1  1
zdef 1 levels 1000
tdef 140 linear 01MAR1999 1MO
vars 1
pdsi 1 99 timeseries rain data
endvars








'open  d:\srdp\xiangguan\2D\w-pjuping.ctl'
*打开年平均降水的时间序列
'open d:\srdp\xiangguan\2D\hgt.mon.mean.ctl'

*打开.nc场数据
lev1=1000
lev2=850
*设定层次
nt=139
*时间时序列长度
'set lon 0'
'set lat 0'
'set t  1  140'
'set lev ' lev1
'define a=pdsi.1'
'set display color white'
'set lon 0 360'
'set lat -90 90'
'set lev 'lev2
'set t  1'
'd  tcorr(a,hgt.2,t=1,t=140)'
*'set gxout print'
*file='d:/test/test1_xg.txt'
'set gxout fwrite'
'set fwrite d:\srdp\xiangguan\2D\nc_xg_Ravg.grd'
'd  tcorr(a,hgt.2,t=1,t=140)'
*rc=write(file,result,append)
'disable fwrite'
*'set gxout contour'
*'d tcorr(a,hgt.2,t=1,t='nt')'
*lons=lons+2.5
*endwhile
*lats=lats+2.5
*endwhile
;

mofangbao

你想这样做的话就把两个文件打开的顺序换一下，因为第二次你其实是是设置的第一个文件的维度，当然会超出范围，解决的方法还有很多，你可以自己再想想

小蝌蚪

我觉得可能是你时间序列定义不对啊？hgt.mon.mean.ctl中时间序列为多少，要和年平均降水的时间序列
一致。

夜未央

小蝌蚪 发表于 2013-4-7 09:28

                               
登录/注册后可看大图

我觉得可能是你时间序列定义不对啊？hgt.mon.mean.ctl中时间序列为多少，要和年平均降水的时间序列
一致。

d  tcorr(a,hgt.2,t=1,t=140)  这个应该就是设置hgt的时间序列吧？

小蝌蚪

夜未央 发表于 2013-4-7 18:47

                               
登录/注册后可看大图

d  tcorr(a,hgt.2,t=1,t=140)  这个应该就是设置hgt的时间序列吧？

hgt.mon.mean.ctl里面t也是140吗？

夜未央

小蝌蚪 发表于 2013-4-7 22:12

                               
登录/注册后可看大图

hgt.mon.mean.ctl里面t也是140吗？

是的 一共140个数，按列排的
program pdsiTOBinary

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!把一列数写成dat
implicit none
integer,parameter:: tn=140
  real :: val(tn),irec
  integer i
!tn为年数，时间序列长度
!读源数据

  open(1,file ='pdsijuping.txt',status='old')
   

open(3,file='pdsijuping.dat',form='unformatted',access='direct',recl=1)

  
irec=0
do i=1,tn
   irec=irec+1
   read(1,*)  val(i)  
   write(3,rec=irec) val(i)
end do
close(1)
close(3)


stop
end

这个是写成dat 的 fortran程序

一闭一睁

回个帖！

839168298

我也在整这个，整不对啊。。。老提示维数设置不对

夜未央

839168298 发表于 2013-6-9 15:38

                               
登录/注册后可看大图

我也在整这个，整不对啊。。。老提示维数设置不对

我没弄出来  后来转MATLAB了

839168298

夜未央 发表于 2013-6-17 00:12

                               
登录/注册后可看大图

我没弄出来  后来转MATLAB了

后来我整出来了！
下面三我的gs
你参考一下：不过有点迟了，mothod1不对的
************
'reinit'
'open c:\18\swt-summermatch.ctl'
'open c:\18\wy-match.ctl'
'set grid off'
'set grads off'
'set dfile 2'
'set x 1'
'set y 1'
'set z 1'
'set t 1 44'
*****************************************method1
*'set gxout fwrite'
*'set fwrite c:\18\rr.dat'
*y=1
*while(y<=61)
*x=1
*while(x<=72)
*'define aw=ave(wy.2,t=1,t=44)'
*'define as=ave(swt.1(x='x',y='y'),t=1,t=44)'
*'define xy=sumg((wy-aw)*(swt.1(x='x',y='y')-as),t=1,t=44)'
*'define x2=sumg((wy-aw)*(wy-aw),t=1,t=44)'
*'define y2=sumg((swt.1(x='x',y='y')-as)*(swt.1(x='x',y='y')-as),t=1,t=44)'
*'define b=xy/(pow(x2,0.5)*pow(y2,0.5))'
*'d b'
*x=x+1
*endwhile
*y=y+1
*endwhile
*'disable fwrite'

*****************************************method2
'define a=wy.2'
'set dfile 1'
'set lon 0 360'
'set lat -60 60'
'set t 1'
'define b=tcorr(a,swt.1,t=1,t=44)'
*******************************************draw it
'set gxout shaded'
'set t 1'
***********自由度n-m-1=41
***********|r|>0.30079通过95%可信检验，|r|>0.38868通过99%可信检验
*'set black -0.38868 0.38868'
'set black -0.29 0.29'
*'set black -0.25079 0.25079'
'set lon 50 250'
'set lat -20 45'
'd b'
'set font 1'
'draw title W-Y INDEX'
'set gxout contour'
'set cthick 5'
'set csmooth on'
'set clab on'
'set cint 0.3'
'd b'
'run cbarn.gs'
'draw line 0.5 3.5 10.5 3.5'
'printim c:\18\1.jpg white x2000 y1000'