grads子程序调用

漂流的屋顶

请问一下，我想对eof第一模态的时间序列和nino3指数做一下滞后相关。现在下载了一个滞后相关的gs文件
请问这个文件怎么调用啊

漂流的屋顶

sorry，忘了下载要花钱了
gs文件如下
***************************************************************************************
*        $Id: tlag.gs,v 1.9 2013/04/17 01:36:51 bguan Exp $
*        Copyright (C) 2004-2013 Bin Guan.
*        Distributed under GNU/GPL.
***************************************************************************************
function main(arg)
*
* Lag regression/correlation/composite.
*
rc=gsfallow('on')

switch=subwrd(arg,1)

input1=subwrd(arg,2)
input2=subwrd(arg,3)
output=subwrd(arg,4)
_lags=subwrd(arg,5)
_lage=subwrd(arg,6)
if(output='')
  usage()
  return
endif
if(_lags='')
  _lags=0
endif
if(_lage='')
  _lage=0
endif
lags=_lage-_lags+1

*
* Ensure x-coordinates are integers and there are no redundant grid points.
*
qdims()
_xs_old=_xs
_xe_old=_xe
if(math_int(_xs)!=_xs | math_int(_xe)!=_xe)
  xs_new=math_nint(_xs)
  xe_new=math_nint(_xe)
  'set x 'xs_new' 'xe_new
  qdims()
endif
if(_lone-_lons>=360)
  rddnt_points=(_lone-_lons-360)/_dlon+1
  'set x '_xs' '_xe-rddnt_points
  qdims()
endif

*
* Ensure y-coordinates are integers.
*
qdims()
_ys_old=_ys
_ye_old=_ye
if(math_int(_ys)!=_ys | math_int(_ye)!=_ye)
  ys_new=math_nint(_ys)
  ye_new=math_nint(_ye)
  'set y 'ys_new' 'ye_new
  qdims()
endif

if(_ts=_te)
  say '[tlag ERROR] Time not varying. Use "set t ..." or "set time ..." to set a varying time dimension.'
  return
endif

prompt '[tlag info] Temporary directory is set to '
'printenv /tmp/bGASL-$USER'
mytmpdir=sublin(result,1)
'!mkdir -p 'mytmpdir

prompt '[tlag info] Process ID is '
'printenv $$'
randnum=sublin(result,1)

writectl(mytmpdir,lags,output,switch,randnum)

'set t '_ts' '_te
'lgcrtmp1='input1

'set gxout fwrite'
'set fwrite 'mytmpdir'/tlag.dat.'randnum

lag=_lags
while(lag<=_lage)
  zcnt=_zs
  while(zcnt<=_ze)
    'set z 'zcnt
    'set t '_ts' '_te
    'lgcrtmp1lagged=tloop(lgcrtmp1(t-'lag'))'
    'set t 1'
    if(switch!='comp')
      'lgcrtmp3=t'switch'(lgcrtmp1lagged,'input2',t='_ts',t='_te')'
    else
      'lgcrtmp3=ave(maskout('input2',lgcrtmp1lagged),t='_ts',t='_te')'
    endif
    'display const(lgcrtmp3,'dfile_undef()',-u)'
    zcnt=zcnt+1
  endwhile
  lag=lag+1
endwhile

'disable fwrite'
'undefine lgcrtmp1'
'undefine lgcrtmp1lagged'
'undefine lgcrtmp3'
'set gxout contour'

'open 'mytmpdir'/tlag.ctl.'randnum
file_num=file_number()
'set x '_xs_old' '_xe_old
* The above line is needed to ensure that there will not be a gap near the prime meridian in global maps if unintended.
'set y '_ys_old' '_ye_old
'set z '_zs' '_ze
'set t '_lags' '_lage
*offset=_ts-_lags
offset=1-_lags
output'=tloop('output'.'file_num'(t+'offset'))'
'close 'file_num

*'set x '_xs_old' '_xe_old
*'set y '_ys_old' '_ye_old
*'set z '_zs' '_ze
*These were already set above.
'set t '_ts' '_te

'!unlink 'mytmpdir'/tlag.dat.'randnum

return
***************************************************************************************
function file_number()
*
* Get the number of files opened.
*
'q files'
line1=sublin(result,1)
if(line1='No files open')
  return 0
endif

lines=1
while(sublin(result,lines+1)!='')
  lines=lines+1
endwhile

return lines/3
***************************************************************************************
function dfile_undef()
*
* Get undef value from the default .ctl file. (Not 'q undef', which is for output.)
*
'query ctlinfo'
if(result='No Files Open')
  return 'unknown'
endif

lines=1
while(1)
  lin=sublin(result,lines)
  if(subwrd(lin,1)='undef'|subwrd(lin,1)='UNDEF')
    return subwrd(lin,2)
  endif
  lines=lines+1
endwhile
***************************************************************************************
function default_tims()
*
* Get the beginning time step of the default file.
*
'set t 1'
'query dims'
lt=sublin(result,5)
tims=subwrd(lt,6)

return tims
***************************************************************************************
function writectl(mytmpdir,lags,output,switch,randnum)
*
* Write the .ctl file for the temporary .dat file
*
lines=9
line.1='DSET ^tlag.dat.'randnum
line.2='UNDEF 'dfile_undef()
if(_cal='')
  line.3='*Intentionally left blank.'
else
  line.3='OPTIONS '_cal
endif
line.4=_xdef
line.5=_ydef
line.6=_zdef
*line.7='TDEF 'lags' LINEAR '_tims' '_dtim
line.7='TDEF 'lags' LINEAR 'default_tims()' '_dtim
line.8='VARS 1'
line.9='ENDVARS'
cnt=1
while(cnt<=lines-1)
  status=write(mytmpdir'/tlag.ctl.'randnum,line.cnt)
  cnt=cnt+1
endwhile
cnt=1
while(cnt<=1)
  varline=output' '_nz0' 99 Add description here.'
  status=write(mytmpdir'/tlag.ctl.'randnum,varline)
  cnt=cnt+1
endwhile
status=write(mytmpdir'/tlag.ctl.'randnum,line.lines)
status=close(mytmpdir'/tlag.ctl.'randnum)

return
***************************************************************************************
function usage()
*
* Print usage information.
*
say 'Lag regression/correlation/composite.'
say ''
say 'USAGE: tlag regr|corr|comp <input1> <input2> <output> [<lag_start> [<lag_end>]]'
say '  regr|corr|comp: regr for regression, corr for correlation, comp for composite.'
say '  <input1>: Independent variable, or mask for compositing (mask>=0 where conditions are met). Can be any GrADS expression. Can have a vertical dimension. Cannot have horizontal dimensions.'
say '  <input2>: Dependent variable. Can be any GrADS expression. Can have vertical and horizontal dimensions.'
say '  <output>: Output variable. E.g., "set t 0" and then <output>(t-3) will be lag regression/correlation/composite when <input2> leads <input1> by 3 time steps (<input2> is shifted forward by 3 time steps for this calculation).'
say '  <lag_start>: Beginning lag. E.g., <lag1>=-3 for <input2> leading <input1> by 3 time steps. Default=0.'
say '  <lag_end>: Ending lag. E.g., <lag2>=3 for <input2> lagging <input2> by 3 time steps. Default=0.'
say ''
say 'NOTE 1: <input2> must be on a grid consistent with the default file. If not, use "set dfile" to change the default file.'
say ''
say 'NOTE 2: "set t 0" and use <output>(t+number) to get the value at lag(number), or "set t <lag1> <lag2>" to get values at all lags.'
say ''
say 'EXAMPLE 1: calculate and display auto-correlation function of Nino3.4 index.'
say '  tlag corr nino34 nino34 out -24 24'
say '  set t -24 24'
say '  display out'
say ''
say 'EXAMPLE 2: calculate and display lag correlation between Nino3.4 index and global precipitation.'
say '  tlag regr nino34 precip out -12 12'
say '  set t 0'
say '  subplot 8 1'
say '  display out(t-12)'
say '  subplot 8 2'
say '  display out(t-8)'
say '  subplot 8 3'
say '  display out(t-4)'
say '  ...'
say ''
say 'DEPENDENCIES: qdims.gsf'
say ''
say 'COPYRIGHT (C) 2004-2013 Bin Guan.'
say 'Distributed under GNU/GPL.'
return

river

漂流的屋顶 发表于 2013-5-9 20:37

                               
登录/注册后可看大图

sorry，忘了下载要花钱了
gs文件如下
*************************************************************** ...

gs最后一部分不就是说明用法的么，楼主自己看一下，试试。
say 'USAGE: tlag regr|corr|comp <input1> <input2> <output> [<lag_start> [<lag_end>]]'
say '  regr|corr|comp: regr for regression, corr for correlation, comp for composite.'
say '  <input1>: Independent variable, or mask for compositing (mask>=0 where conditions are met). Can be any GrADS expression. Can have a vertical dimension. Cannot have horizontal dimensions.'
say '  <input2>: Dependent variable. Can be any GrADS expression. Can have vertical and horizontal dimensions.'
say '  <output>: Output variable. E.g., "set t 0" and then <output>(t-3) will be lag regression/correlation/composite when <input2> leads <input1> by 3 time steps (<input2> is shifted forward by 3 time steps for this calculation).'
say '  <lag_start>: Beginning lag. E.g., <lag1>=-3 for <input2> leading <input1> by 3 time steps. Default=0.'
say '  <lag_end>: Ending lag. E.g., <lag2>=3 for <input2> lagging <input2> by 3 time steps. Default=0.'
say ''
say 'NOTE 1: <input2> must be on a grid consistent with the default file. If not, use "set dfile" to change the default file.'
say ''
say 'NOTE 2: "set t 0" and use <output>(t+number) to get the value at lag(number), or "set t <lag1> <lag2>" to get values at all lags.'
say ''
say 'EXAMPLE 1: calculate and display auto-correlation function of Nino3.4 index.'
say '  tlag corr nino34 nino34 out -24 24'
say '  set t -24 24'
say '  display out'
say ''
say 'EXAMPLE 2: calculate and display lag correlation between Nino3.4 index and global precipitation.'
say '  tlag regr nino34 precip out -12 12'
say '  set t 0'
say '  subplot 8 1'
say '  display out(t-12)'
say '  subplot 8 2'
say '  display out(t-8)'
say '  subplot 8 3'
say '  display out(t-4)'

chaijing

学习了~~~~~

somnuszhouying

真是学习了，顶一下

漂流的屋顶

river 发表于 2013-5-9 22:13

                               
登录/注册后可看大图

gs最后一部分不就是说明用法的么，楼主自己看一下，试试。
say 'USAGE: tlag regr|corr|comp    [ []]'
...

嗯，谢谢。我再看看吧。想问一下，滞后相关后得到的应该是个一维数组吗？

river

漂流的屋顶 发表于 2013-5-9 22:31

                               
登录/注册后可看大图

嗯，谢谢。我再看看吧。想问一下，滞后相关后得到的应该是个一维数组吗？

看你做的是什么的吧，有单点对单点的，有单点对场的，有场对场的

漂流的屋顶

river 发表于 2013-5-9 23:06

                               
登录/注册后可看大图

看你做的是什么的吧，有单点对单点的，有单点对场的，有场对场的

是单点对单点的哈，那是不是就是一维数组啊？

river

漂流的屋顶 发表于 2013-5-9 23:08

                               
登录/注册后可看大图

是单点对单点的哈，那是不是就是一维数组啊？

你这个应该是和eof第一模态的时间序列长度一样的一维数组

漂流的屋顶

river 发表于 2013-5-10 07:56

                               
登录/注册后可看大图

你这个应该是和eof第一模态的时间序列长度一样的一维数组

嗯，谢谢大神