Hi everyone,
When I used Example 5 to in http://www.ncl.ucar.edu/Document/Functions/Built-in/specx_anal.shtml to calculate the mean power spectrum, there is still an error message.
In my code, I defined a matrix pw, which is the same as pc2dT:
-------------------------------------------------------------
ncl 217> printVarSummary(pc2dT)
Variable: pc2dT
Type: double
Total Size: 29520 bytes
3690 values
Number of Dimensions: 2
Dimensions and sizes: [year | 30] x [day | 123]
Coordinates:
Number Of Attributes: 1
_FillValue : -999000000
---------------------------------------------------------------
The difference between my code and Example is: I defined r1 and r1zsum before the loop, otherwise the loop cannot continue. Also I defined splt before using splt = specx_ci(df, 0.05, 0.95), because it always said splt is undefined.
Here is my code:
;===============================
; power spectra
;===============================
;x(nseg,ntim), nseg is nyrs here, ntim is nday here;
ntim = nday
nseg = nyrs
pw = pc2dT
d = 0 ; detrending opt: 0=>remove mean 1=>remove mean + detrend
sm = 21 ; smooth: should be at least 3 and odd
pct = 0.10 ; percent taper: (0.0 <= pct <= 1.0) 0.10 common.
;************************************************
; calculate mean spectrum spectrum and lag1 auto cor
;************************************************
; loop over each segment of length ntim
r1=new(1,double)
r1=0.0
r1zsum=new(1,double)
r1zsum=0.0
spcavg = new ( ntim/2, typeof(pw))
spcavg = 0.0
do n=0,nseg-1
dof = specx_anal(pw(n,:),d,sm,pct) ; current segment spc
spcavg = spcavg + dof@spcx ; sum spc of each segment
r1 = dof@xlag1 ; extract segment lag-1
r1zsum = r1zsum + 0.5*(log((1.0+r1)/(1.0-r1))) ; sum the Fischer Z
end do
r1z = r1zsum/nseg ; average r1z
r1 = (exp(2*r1z)-1)/(exp(2*r1z)+1) ; transform back
; this is the mean r1
spcavg = spcavg/nseg ; average spectrum
;************************************************
; Assign mean spectrum to data object
;************************************************
df = 2.0*nseg ; deg of freedom
; all segments
df@spcx = spcavg ; assign the mean spc
df@frq = dof@frq
spcavg@xlag1 = r1 ; assign mean lag-1
;************************************************
; plotting
;************************************************
wks = gsn_open_wks("ps","spec") ; Opens a ps file
res = True
res@tiMainString = "Mean Spectra: "+nseg+" segments, dof="+df ; title
res@tiXAxisString = "Frequency (cycles/month)" ; xaxis
res@tiYAxisString = "Variance" ; yaxis
splt = new((/4,dimsizes(df@spcx)/),float)
splt = specx_ci(df, 0.05, 0.95) ; confidence interval
plot = gsn_csm_xy(wks, df@frq, splt ,res) ; create plot
;===================================
However, The code is still stuck in splt = specx_ci(df, 0.05, 0.95). The error message is
--------------------------
warning:Attempt to reference attribute (xlag1) which is undefined
fatal:Assignment type mismatch, right hand side can't be coerced to type of left hand side
fatal:["Execute.c":8128]:Execute: Error occurred at or near line 1348 in file $NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl
fatal:["Execute.c":8126]:Execute: Error occurred at or near line 216
----------------------------------------------------------
?Thank you very much.
Best regards
Yan
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sun Apr 20 15:38:58 2014
This archive was generated by hypermail 2.1.8 : Tue Apr 29 2014 - 09:04:20 MDT